JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
TMCFMUtils.cc
Go to the documentation of this file.
1 #include "MELAStreamHelpers.hh"
2 #include "TMCFMUtils.hh"
3 #include "TUtilHelpers.hh"
4 #include "TMath.h"
5 
6 
7 using MELAStreamHelpers::MELAout;
8 using MELAStreamHelpers::MELAerr;
9 using namespace std;
10 using namespace PDGHelpers;
11 using namespace TNumericUtil;
12 
13 namespace TMCFMUtils{
14  const std::vector<intQuad_t> MCFMHash_QQVVQQAny = Hash_QQVVQQAny();
15 }
16 
17 
18 void TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rSel, int sSel, int order[2]){
19  const std::vector<intQuad_t>& hash = MCFMHash_QQVVQQAny;
20  bool outFound=false;
21  for (intQuad_t const& hashel:hash){
22  if (
23  !(
24  (PDGHelpers::isAnUnknownJet(iSel) || iSel==hashel[0])
25  &&
26  (PDGHelpers::isAnUnknownJet(jSel) || jSel==hashel[1])
27  )
28  ) continue;
29  // Final particles are q
30  if (PDGHelpers::isAJet(rSel) && PDGHelpers::isAJet(sSel)){
31  if (
32  (PDGHelpers::isAnUnknownJet(rSel) || rSel==hashel[2])
33  &&
34  (PDGHelpers::isAnUnknownJet(sSel) || sSel==hashel[3])
35  ){
36  order[0]=0;
37  order[1]=1;
38  outFound=true;
39  //MELAout << "Hash requested outgoing "<< hashel[2] << " " << hashel[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
40  }
41  else if (
42  (PDGHelpers::isAnUnknownJet(rSel) || rSel==hashel[3])
43  &&
44  (PDGHelpers::isAnUnknownJet(sSel) || sSel==hashel[2])
45  ){
46  order[0]=1;
47  order[1]=0;
48  outFound=true;
49  //MELAout << "Hash requested outgoing "<< hashel[2] << " " << hashel[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
50  }
51  }
52  // Final particles l/nu
53  else if ((PDGHelpers::isALepton(rSel) || PDGHelpers::isANeutrino(rSel)) && (PDGHelpers::isALepton(sSel) || PDGHelpers::isANeutrino(sSel))){
54  if (abs(hashel[0])==abs(hashel[1]) && abs(hashel[0])==abs(hashel[2]) && abs(hashel[0])==abs(hashel[3])) continue; // Do not consider the ordering in uquq_uquq or dqdq_dqdq
55 
56  if (
57  (
58  TMath::Sign(1, rSel)==TMath::Sign(1, hashel[2]) &&
59  ((PDGHelpers::isALepton(rSel) && PDGHelpers::isDownTypeQuark(hashel[2])) || (PDGHelpers::isANeutrino(rSel) && PDGHelpers::isUpTypeQuark(hashel[2])))
60  )
61  &&
62  (
63  TMath::Sign(1, sSel)==TMath::Sign(1, hashel[3]) &&
64  ((PDGHelpers::isALepton(sSel) && PDGHelpers::isDownTypeQuark(hashel[3])) || (PDGHelpers::isANeutrino(sSel) && PDGHelpers::isUpTypeQuark(hashel[3])))
65  )
66  ){
67  order[0]=0;
68  order[1]=1;
69  outFound=true;
70  //MELAout << "Hash requested outgoing "<< hashel[2] << " " << hashel[3] << ", unswapped r, s = " << rSel << " " << sSel << endl;
71  }
72  else if (
73  (
74  TMath::Sign(1, rSel)==TMath::Sign(1, hashel[3]) &&
75  ((PDGHelpers::isALepton(rSel) && PDGHelpers::isDownTypeQuark(hashel[3])) || (PDGHelpers::isANeutrino(rSel) && PDGHelpers::isUpTypeQuark(hashel[3])))
76  )
77  &&
78  (
79  TMath::Sign(1, sSel)==TMath::Sign(1, hashel[2]) &&
80  ((PDGHelpers::isALepton(sSel) && PDGHelpers::isDownTypeQuark(hashel[2])) || (PDGHelpers::isANeutrino(sSel) && PDGHelpers::isUpTypeQuark(hashel[2])))
81  )
82  ){
83  order[0]=1;
84  order[1]=0;
85  outFound=true;
86  //MELAout << "Hash requested outgoing "<< hashel[2] << " " << hashel[3] << ", swapped r, s = " << rSel << " " << sSel << endl;
87  }
88  if (PDGHelpers::getCoupledVertex(rSel, sSel)!=PDGHelpers::getCoupledVertex(hashel[2], hashel[3])) outFound=false;
89  }
90  if (outFound) break;
91  }
92  if (!outFound){ for (unsigned int ip=0; ip<2; ip++) order[ip]=-1; }
93 }
94 std::vector<intQuad_t> TMCFMUtils::Hash_QQVVQQAny(){
95  std::vector<intQuad_t> pcfg;
96  std::vector<intQuad_t> hash_qqvvqq = TMCFMUtils::Hash_QQVVQQ();
97  std::vector<intQuad_t> hash_qqvvqqstrong = TMCFMUtils::Hash_QQVVQQStrong();
98  TUtilHelpers::copyVector(hash_qqvvqq, pcfg);
99  TUtilHelpers::copyVector(hash_qqvvqqstrong, pcfg);
100  return pcfg;
101 }
102 std::vector<intQuad_t> TMCFMUtils::Hash_QQVVQQ(){
103  /*
104  Based on the following cases in MCFM:
105  parameter(
106  & uqcq_uqcq=1,uquq_uquq=2,dqsq_dqsq=3,
107  & dqdq_dqdq=4,uqbq_uqbq=5,dqcq_dqcq=6,
108  & dquq_dquq=7,dqcq_uqsq=8,uqsq_dqcq=9
109  & jmax=12)
110  integer,parameter:: j1(jmax)=(/1,2,8,8, 7,2,7,1, 1,7,2,7/)
111  integer,parameter:: j2(jmax)=(/2,1,7,7, 2,7,1,7, 7,1,7,2/)
112  integer,parameter:: j7(jmax)=(/7,7,2,1, 1,8,2,8, 2,8,1,8/)
113  integer,parameter:: j8(jmax)=(/8,8,1,2, 8,1,8,2, 8,2,8,1/)
114  */
115  std::vector<intQuad_t> base_cfg; base_cfg.reserve(17);
116  // uc_uc
117  base_cfg.emplace_back(2, 4, 2, 4);
118  // ds_ds
119  base_cfg.emplace_back(1, 3, 1, 3);
120  base_cfg.emplace_back(1, 5, 1, 5);
121  base_cfg.emplace_back(3, 5, 3, 5);
122  // ub_ub
123  base_cfg.emplace_back(2, 3, 2, 3);
124  base_cfg.emplace_back(2, 5, 2, 5);
125  base_cfg.emplace_back(4, 5, 4, 5);
126  // dc_dc
127  base_cfg.emplace_back(1, 4, 1, 4);
128  // du_du
129  base_cfg.emplace_back(1, 2, 1, 2);
130  base_cfg.emplace_back(3, 4, 3, 4);
131  // dc_us
132  base_cfg.emplace_back(1, 4, 2, 3);
133  // us_dc
134  base_cfg.emplace_back(2, 3, 1, 4);
135  // uu_uu
136  base_cfg.emplace_back(2);
137  base_cfg.emplace_back(4);
138  // dd_dd
139  base_cfg.emplace_back(1);
140  base_cfg.emplace_back(3);
141  base_cfg.emplace_back(5);
142 
143  std::vector<intQuad_t> jcfg; jcfg.reserve(12);
144  jcfg.emplace_back(0, 1, 2, 3);
145  jcfg.emplace_back(1, 0, 2, 3);
146  jcfg.emplace_back(3, 2, 1, 0);
147  jcfg.emplace_back(3, 2, 0, 1);
148  jcfg.emplace_back(2, 1, 0, 3);
149  jcfg.emplace_back(1, 2, 3, 0);
150  jcfg.emplace_back(2, 0, 1, 3);
151  jcfg.emplace_back(0, 2, 3, 1);
152  jcfg.emplace_back(0, 2, 1, 3);
153  jcfg.emplace_back(2, 0, 3, 1);
154  jcfg.emplace_back(1, 2, 0, 3);
155  jcfg.emplace_back(2, 1, 3, 0);
156 
157  std::vector<intQuad_t> pcfg;
158  for (unsigned int j=0; j<jcfg.size(); j++){
159  for (unsigned int p=0; p<base_cfg.size(); p++){
160  intQuad_t cfg;
161  for (unsigned int ipos=0; ipos<4; ipos++){
162  int idpos = jcfg.at(j)[ipos];
163  int idAssigned = base_cfg.at(p)[ipos];
164  if ((idpos<2 && ipos>=2) || (idpos>=2 && ipos<2)) idAssigned = -idAssigned;
165  cfg[jcfg.at(j)[ipos]] = idAssigned;
166  }
167  if (!TUtilHelpers::checkElementExists(cfg, pcfg)) pcfg.push_back(cfg);
168  }
169  }
170  // Uncommenting the lines below prints out the hash when the library is loaded.
171  /*
172  for (unsigned int ic=0; ic<pcfg.size(); ic++) MELAout
173  << "TMCFMUtils::Hash_QQVVQQAny: Hash configuration " << ic << " requests ids=( "
174  << pcfg.at(ic)[0] << ", " << pcfg.at(ic)[1] << ", " << pcfg.at(ic)[2] << ", " << pcfg.at(ic)[3]
175  << ")" << endl;
176  */
177  return pcfg;
178 }
179 std::vector<intQuad_t> TMCFMUtils::Hash_QQVVQQStrong(){
180  /*
181  Based on the following cases in MCFM:
182  call qq4lggampf(1,2,3,4,5,6,7,8,3,4,za,zb,msqgg)
183  msq(1,-1)=msq(1,-1)+stat*aveqq*msqgg(1)
184  call qq4lggampf(2,1,3,4,5,6,7,8,3,4,za,zb,msqgg)
185  msq(-1,1)=msq(-1,1)+stat*aveqq*msqgg(1)
186  call qq4lggampf(7,8,3,4,5,6,1,2,3,4,za,zb,msqgg)
187  msq(0,0)=msq(0,0)+avegg*(3d0*msqgg(1)+2d0*msqgg(2))
188  call qq4lggampf(7,2,3,4,5,6,1,8,3,4,za,zb,msqgg)
189  msq(0,-1)=msq(0,-1)+aveqg*msqgg(1)
190  call qq4lggampf(7,1,3,4,5,6,2,8,3,4,za,zb,msqgg)
191  msq(-1,0)=msq(-1,0)+aveqg*msqgg(1)
192  call qq4lggampf(2,8,3,4,5,6,1,7,3,4,za,zb,msqgg)
193  msq(0,1)=msq(0,1)+aveqg*msqgg(1)
194  call qq4lggampf(1,8,3,4,5,6,2,7,3,4,za,zb,msqgg)
195  msq(1,0)=msq(1,0)+aveqg*msqgg(1)
196  */
197  std::vector<intQuad_t> base_cfg; base_cfg.reserve(5);
198  // Start with qqb_gg
199  for (int iq=1; iq<=5; iq++) base_cfg.emplace_back(iq, -iq, 21, 21);
200 
201  std::vector<intQuad_t> jcfg; jcfg.reserve(7);
202  jcfg.emplace_back(0, 1, 2, 3); // qqb->gg
203  jcfg.emplace_back(1, 0, 2, 3); // qbq->gg
204  jcfg.emplace_back(2, 3, 0, 1); // gg->qbq
205  jcfg.emplace_back(2, 1, 0, 3); // gqb->qbg
206  jcfg.emplace_back(2, 0, 1, 3); // qbg->qbg
207  jcfg.emplace_back(1, 3, 2, 0); // gq->gq
208  jcfg.emplace_back(0, 3, 2, 1); // qg->gq
209 
210  std::vector<intQuad_t> pcfg;
211  for (unsigned int j=0; j<jcfg.size(); j++){
212  for (unsigned int p=0; p<base_cfg.size(); p++){
213  intQuad_t cfg;
214  for (unsigned int ipos=0; ipos<4; ipos++){
215  int idpos = jcfg.at(j)[ipos];
216  int idAssigned = base_cfg.at(p)[ipos];
217  if (((idpos<2 && ipos>=2) || (idpos>=2 && ipos<2)) && !PDGHelpers::isAGluon(idAssigned)) idAssigned = -idAssigned;
218  cfg[jcfg.at(j)[ipos]] = idAssigned;
219  }
220  if (!TUtilHelpers::checkElementExists(cfg, pcfg)) pcfg.push_back(cfg);
221  }
222  }
223  // Uncommenting the lines below prints out the hash when the library is loaded.
224  /*
225  for (unsigned int ic=0; ic<pcfg.size(); ic++) MELAout
226  << "TMCFMUtils::Hash_QQVVQQAny: Hash configuration " << ic << " requests ids=( "
227  << pcfg.at(ic)[0] << ", " << pcfg.at(ic)[1] << ", " << pcfg.at(ic)[2] << ", " << pcfg.at(ic)[3]
228  << ")" << endl;
229  */
230  return pcfg;
231 }
PDGHelpers::isALepton
bool isALepton(const int id)
Definition: PDGHelpers.cc:62
TUtilHelpers::checkElementExists
bool checkElementExists(T const &element, std::vector< U > const &elementlist)
Definition: TUtilHelpers.hh:42
PDGHelpers
Definition: PDGHelpers.h:12
TMCFMUtils::Hash_QQVVQQStrong
std::vector< TNumericUtil::intQuad_t > Hash_QQVVQQStrong()
Definition: TMCFMUtils.cc:179
TMCFMUtils::Hash_QQVVQQAny
std::vector< TNumericUtil::intQuad_t > Hash_QQVVQQAny()
Definition: TMCFMUtils.cc:94
PDGHelpers::isAnUnknownJet
bool isAnUnknownJet(const int id)
Definition: PDGHelpers.cc:22
PDGHelpers::isANeutrino
bool isANeutrino(const int id)
Definition: PDGHelpers.cc:67
TMCFMUtils
Definition: TMCFMUtils.hh:13
MELAStreamHelpers::MELAout
MELAOutputStreamer MELAout
PDGHelpers::isDownTypeQuark
bool isDownTypeQuark(const int id)
Definition: PDGHelpers.cc:45
TNumericUtil
Definition: MELAAccumulators.h:6
PDGHelpers::getCoupledVertex
int getCoupledVertex(const int idfirst, const int idsecond, int *hel=0, int *useAHcoupl=0)
Definition: PDGHelpers.cc:235
TUtilHelpers::copyVector
void copyVector(std::vector< T > const &input, std::vector< T > &target)
Definition: TUtilHelpers.hh:24
TMCFMUtils::MCFMHash_QQVVQQAny
const std::vector< intQuad_t > MCFMHash_QQVVQQAny
Definition: TMCFMUtils.cc:14
MELAStreamHelpers::MELAerr
MELAOutputStreamer MELAerr
TUtilHelpers.hh
PDGHelpers::isAGluon
bool isAGluon(const int id)
Definition: PDGHelpers.cc:58
TNumericUtil::quadruplet
Definition: TNumericUtil.hh:27
MELAStreamHelpers.hh
TMCFMUtils.hh
TMCFMUtils::Hash_QQVVQQ
std::vector< TNumericUtil::intQuad_t > Hash_QQVVQQ()
Definition: TMCFMUtils.cc:102
TMCFMUtils::AssociatedParticleOrdering_QQVVQQAny
void AssociatedParticleOrdering_QQVVQQAny(int iSel, int jSel, int rSel, int sSel, int order[2])
Definition: TMCFMUtils.cc:18
PDGHelpers::isUpTypeQuark
bool isUpTypeQuark(const int id)
Definition: PDGHelpers.cc:40
PDGHelpers::isAJet
bool isAJet(const int id)
Definition: PDGHelpers.cc:18