JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
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
testME_all.p
p
Definition: testME_all.py:11
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
Generated by   doxygen 1.8.17