JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
MELAPConstant.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <cmath>
4 #include "MELAStreamHelpers.hh"
5 #include "MelaPConstant.h"
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TString.h"
9 
10 
11 using namespace std;
12 using TVar::simple_event_record;
13 using MELAStreamHelpers::MELAout;
14 using MELAStreamHelpers::MELAerr;
15 
16 
17 MelaPConstant::MelaPConstant(
18  TVar::MatrixElement me_,
19  TVar::Production prod_,
20  TVar::Process proc_,
21  const char* path,
22  const char* spname
23  ) :
24  processME(me_),
25  processProd(prod_),
26  processProc(proc_),
27  fcnLow(0),
28  fcnHigh(0),
29  fcnMid(0),
30  fname(path)
31 {
32  GetFcnFromFile(path, spname);
33 }
34 
35 MelaPConstant::~MelaPConstant(){}
36 
37 void MelaPConstant::GetFcnFromFile(const char* path, const char* spname){
38  TString spname_core = spname;
39  spname_core.Append("_Smooth");
40  spname_core.Prepend("sp_tg_");
41 
42  TFile* fin = TFile::Open(path, "read");
43  gROOT->cd();
44 
45  if (fin!=0 && !fin->IsZombie() && fin->IsOpen()){
46  TString spname_low = spname_core; spname_low.Prepend("lowFcn_");
47  fcnLow = (TF1*)fin->Get(spname_low);
48 
49  TString spname_high = spname_core; spname_high.Prepend("highFcn_");
50  fcnHigh = (TF1*)fin->Get(spname_high);
51 
52  TString spname_mid = spname_core;
53  fcnMid = (TSpline3*)fin->Get(spname_mid);
54  }
55  else MELAerr << "MelaPConstant::GetFcnFromFile: Failed to open file in path " << path << endl;
56  if (fin!=0 && fin->IsOpen()) fin->Close();
57 }
58 
59 double MelaPConstant::Eval(const MelaIO* RcdME, TVar::VerbosityLevel verbosity)const{
60  if (verbosity>=TVar::DEBUG) MELAout << "Begin MelaPConstant::Eval" << endl;
61 
62  double result=1;
63  if (RcdME->melaCand->genStatus==-1) return result;
64 
65  double candMass = RcdME->melaCand->m();
66  if (verbosity>=TVar::DEBUG) MELAout << "MelaPConstant::Eval: Candidate mass is " << candMass << endl;
67  if (candMass<=0.) return 0.;
68  else if (fcnLow!=0 && candMass<fcnLow->GetXmax()) result = fcnLow->Eval(candMass);
69  else if (fcnHigh!=0 && candMass>fcnHigh->GetXmin()) result = fcnHigh->Eval(candMass);
70  else if (fcnMid!=0){
71  double var = candMass;
72  if (var<fcnMid->GetXmin()) var = fcnMid->GetXmin();
73  else if (var>fcnMid->GetXmax()) var = fcnMid->GetXmax();
74  result = fcnMid->Eval(candMass);
75  }
76  else return result;
77 
78  if (verbosity>=TVar::DEBUG) MELAout << "MelaPConstant::Eval: Spline evaluated to " << result << endl;
79 
80  bool multiplyALR=false;
81  bool multiplyHprop=false;
82  bool multiplyAVjjprop=false;
83  int sigmaZ4lprop=0;
84  unsigned int powAlphaSatMZ=0;
85  result = pow(10., result);
86  if (
87  processME==TVar::JHUGen
88  &&
89  processProd==TVar::ZZINDEPENDENT
90  &&
91  processProc==TVar::HSMHiggs
92  ){
93  multiplyALR=true;
94  }
95  else if (
96  (processME==TVar::JHUGen || processME==TVar::MCFM)
97  &&
98  processProd==TVar::ZZGG
99  &&
100  processProc==TVar::HSMHiggs
101  ){
102  multiplyALR=true;
103  powAlphaSatMZ=2;
104  multiplyHprop=true;
105  }
106  else if (
107  processME==TVar::MCFM
108  &&
109  processProd==TVar::ZZGG
110  &&
111  processProc==TVar::bkgZZ
112  ){
113  multiplyALR=true;
114  powAlphaSatMZ=2;
115  }
116  else if (
117  processME==TVar::MCFM
118  &&
119  processProd==TVar::ZZQQB
120  &&
121  processProc==TVar::bkgZZ
122  ){
123  sigmaZ4lprop=4;
124  multiplyALR=true;
125  }
126  else if (
127  processME==TVar::JHUGen
128  &&
129  processProd==TVar::JQCD
130  &&
131  processProc==TVar::HSMHiggs
132  ){
133  powAlphaSatMZ=3;
134  }
135  else if (
136  processME==TVar::JHUGen
137  &&
138  processProd==TVar::JJQCD
139  &&
140  processProc==TVar::HSMHiggs
141  ){
142  powAlphaSatMZ=4;
143  }
144  else if (
145  processME==TVar::JHUGen
146  &&
147  (processProd==TVar::JJVBF || processProd==TVar::Had_WH || processProd==TVar::Had_ZH)
148  &&
149  processProc==TVar::HSMHiggs
150  ){
151  multiplyHprop=(processProd==TVar::Had_WH || processProd==TVar::Had_ZH);
152  //multiplyAVjjprop=(processProd==TVar::Had_WH || processProd==TVar::Had_ZH); // Disabled to preserve discrimination power
153  }
154  else if (
155  processME==TVar::MCFM
156  &&
157  (
158  (processProc==TVar::HSMHiggs && (processProd==TVar::JJVBF_S || processProd==TVar::Had_WH_S || processProd==TVar::Had_ZH_S))
159  ||
160  (processProc==TVar::bkgZZ && (processProd==TVar::JJVBF || processProd==TVar::Had_WH || processProd==TVar::Had_ZH))
161  )
162  ){
163  multiplyALR=true;
164  multiplyHprop=(processProc==TVar::HSMHiggs);
165  sigmaZ4lprop=4*(processProc==TVar::bkgZZ);
166  //multiplyAVjjprop=(processProd==TVar::Had_WH || processProd==TVar::Had_ZH || processProd==TVar::Had_WH_S || processProd==TVar::Had_ZH_S); // Disabled to preserve discrimination power
167  }
168  else if (
169  processME==TVar::MCFM
170  &&
171  processProd==TVar::JJQCD
172  &&
173  processProc==TVar::bkgZZ
174  ){
175  multiplyALR=true;
176  sigmaZ4lprop=4;
177  powAlphaSatMZ=2;
178  }
179  else if (
180  processME==TVar::MCFM
181  &&
182  processProd==TVar::JJQCD
183  &&
184  processProc==TVar::bkgZJets
185  ){
186  multiplyALR=true;
187  }
188  else{
189  if (verbosity>=TVar::DEBUG) MELAout << "MelaPConstant::Eval: Multiplying " << result << " by nothing." << endl;
190  }
191 
192  if (multiplyALR) result *= this->GetVDaughterCouplings(RcdME, verbosity);
193  if (multiplyHprop) result *= this->GetHPropagator(RcdME, verbosity);
194  if (sigmaZ4lprop!=0) result *= this->GetZPropagator(RcdME, sigmaZ4lprop, verbosity);
195  if (powAlphaSatMZ!=0) result *= this->GetAlphaSatMZ(RcdME, powAlphaSatMZ, verbosity);
196  if (multiplyAVjjprop!=0) result *= this->GetAssociatedVjjPropagator(RcdME, verbosity);
197 
198  if (verbosity>=TVar::DEBUG) MELAout << "End MelaPConstant::Eval with " << result << endl;
199  return result;
200 }
201 
202 double MelaPConstant::GetHPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{
203  double propagator=1.;
204  double sh = pow(RcdME->melaCand->m(), 2);
205  double mh, gah;
206  RcdME->getHiggsMassWidth(mh, gah, 0);
207  propagator = 1./(pow(sh-pow(mh, 2), 2) + pow(mh*gah, 2));
208  if (verbosity>=TVar::DEBUG) MELAout
209  << "MelaPConstant::GetHPropagator: "
210  << "propagatorH=" << propagator << " "
211  << endl;
212  return propagator;
213 }
214 double MelaPConstant::GetZPropagator(const MelaIO* RcdME, const int& sigmaZ4lprop, const TVar::VerbosityLevel& verbosity)const{
215  double propagator=1.;
216  double candMass = RcdME->melaCand->m();
217  double sh = pow(candMass, 2);
218  double mz = TUtil::GetMass(23);
219  double gaz = TUtil::GetMass(23);
220  double prop_sh = 1./(pow(sh-pow(mz, 2), 2) + pow(mz*gaz, 2));
221  if (sigmaZ4lprop>0.){
222  const double dsigmaZ4lprop=sigmaZ4lprop;
223  if (fabs(candMass-mz)<=dsigmaZ4lprop*gaz){
224  double shdn = pow(mz-dsigmaZ4lprop*gaz, 2);
225  double shup = pow(mz+dsigmaZ4lprop*gaz, 2);
226  double prop_shdn = 1./(pow(shdn-pow(mz, 2), 2) + pow(mz*gaz, 2));
227  double prop_shup = 1./(pow(shup-pow(mz, 2), 2) + pow(mz*gaz, 2));
228  double fsh = (sh-shdn)/(shup-shdn);
229  propagator = prop_sh / (prop_shdn*(1.-fsh) + prop_shup*fsh);
230  }
231  }
232  else if (sigmaZ4lprop<0.) propagator = prop_sh;
233  if (verbosity>=TVar::DEBUG) MELAout
234  << "MelaPConstant::GetZPropagator: "
235  << "propagatorZ=" << propagator << " "
236  << endl;
237  return propagator;
238 }
239 double MelaPConstant::GetAssociatedVjjPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{
240  double propagator=1.;
241  if (
242  processProd==TVar::Had_WH || processProd==TVar::Had_WH_S
243  ||
244  processProd==TVar::Had_ZH || processProd==TVar::Had_ZH_S
245  ){
246  // Get a simple event record, safest way to handle jet mass corrections
247  int nRequested_AssociatedJets=2;
248  int AssociationVCompatibility=0;
249  int partIncCode=TVar::kUseAssociated_Jets;
250  double mv=0, gav=0;
251  if (processProd==TVar::Had_WH || processProd==TVar::Had_WH_S){
252  AssociationVCompatibility=24;
253  mv = TUtil::GetMass(24);
254  gav = TUtil::GetMass(24);
255  }
256  else{
257  AssociationVCompatibility=23;
258  mv = TUtil::GetMass(23);
259  gav = TUtil::GetMass(23);
260  }
261  simple_event_record mela_event;
262  mela_event.AssociationCode=partIncCode;
263  mela_event.AssociationVCompatibility=AssociationVCompatibility;
264  mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets;
265  TUtil::GetBoostedParticleVectors(
266  RcdME->melaCand,
267  mela_event,
268  verbosity
269  );
270 
271  float mJJval=-1;
272  vector<TLorentzVector> pJets;
273  const SimpleParticleCollection_t& pAssociated = mela_event.pAssociated;
274  for (auto& part : pAssociated){
275  if (PDGHelpers::isAJet(part.first)) pJets.push_back(part.second);
276  if (pJets.size()==2) break;
277  }
278  if (pJets.size()==2) mJJval = (pJets[0] + pJets[1]).M();
279  if (mJJval>=0.) propagator = 1./(pow(pow(mJJval, 2)-pow(mv, 2), 2) + pow(mv*gav, 2));
280  }
281  if (verbosity>=TVar::DEBUG) MELAout
282  << "MelaPConstant::GetAssociatedVjjPropagator: "
283  << "propagatorV=" << propagator << " "
284  << endl;
285  return propagator;
286 }
287 double MelaPConstant::GetVDaughterCouplings(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{
288  double result=1;
289  double aL1, aR1, aL2, aR2;
290  RcdME->getVDaughterCouplings(aL1, aR1, 0);
291  RcdME->getVDaughterCouplings(aL2, aR2, 1);
292  if (verbosity>=TVar::DEBUG) MELAout
293  << "MelaPConstant::GetVDaughterCouplings: "
294  << "L**2+R**2 couplings=" << pow(aL1, 2)+pow(aR1, 2) << " " << pow(aL2, 2)+pow(aR2, 2) << " "
295  << endl;
296  if (fabs(aL1)>0. || fabs(aR1)>0.) result *= pow(aL1, 2)+pow(aR1, 2);
297  if (fabs(aL2)>0. || fabs(aR2)>0.) result *= pow(aL2, 2)+pow(aR2, 2);
298  return result;
299 }
300 double MelaPConstant::GetAlphaSatMZ(const MelaIO* RcdME, const unsigned int& powAlphaSatMZ, const TVar::VerbosityLevel& verbosity)const{
301  double result;
302  double alphasVal = RcdME->getAlphaSatMZ();
303  result = pow(alphasVal, powAlphaSatMZ);
304  if (verbosity>=TVar::DEBUG) MELAout
305  << "MelaPConstant::GetAlphaSatMZ: "
306  << "alphas(MZ)=" << alphasVal << " (**" << powAlphaSatMZ << ") "
307  << endl;
308  return result;
309 }
310 
311 
MelaPConstant::fcnMid
TSpline3 * fcnMid
Definition: MelaPConstant.h:18
TVar::JJVBF
@ JJVBF
Definition: TVar.hh:72
TVar::Had_ZH_S
@ Had_ZH_S
Definition: TVar.hh:86
MelaIO
Definition: MelaIO.h:8
TVar::VerbosityLevel
VerbosityLevel
Definition: TVar.hh:47
MelaPConstant::fcnHigh
TF1 * fcnHigh
Definition: MelaPConstant.h:17
MelaPConstant.h
hto_masses::mz
real *8, parameter mz
Definition: CALLING_cpHTO.f:76
TVar::simple_event_record::nRequested_AssociatedJets
int nRequested_AssociatedJets
Definition: TVar.hh:229
MelaPConstant::MelaPConstant
MelaPConstant(TVar::MatrixElement me_, TVar::Production prod_, TVar::Process proc_, const char *path, const char *spname)
Definition: MELAPConstant.cc:17
MelaPConstant::processProc
TVar::Process processProc
Definition: MelaPConstant.h:14
MelaIO::getVDaughterCouplings
void getVDaughterCouplings(double &left, double &right, int iv) const
Definition: MelaIO.h:151
TVar::bkgZJets
@ bkgZJets
Definition: TVar.hh:164
MelaPConstant::~MelaPConstant
virtual ~MelaPConstant()
Definition: MELAPConstant.cc:35
MelaPConstant::processProd
TVar::Production processProd
Definition: MelaPConstant.h:13
TVar::ZZINDEPENDENT
@ ZZINDEPENDENT
Definition: TVar.hh:65
TVar::Process
Process
Definition: TVar.hh:124
TVar::JJQCD
@ JJQCD
Definition: TVar.hh:71
MelaPConstant::GetAssociatedVjjPropagator
double GetAssociatedVjjPropagator(const MelaIO *RcdME, const TVar::VerbosityLevel &verbosity) const
Definition: MELAPConstant.cc:239
TVar::ZZQQB
@ ZZQQB
Definition: TVar.hh:63
TVar::kUseAssociated_Jets
@ kUseAssociated_Jets
Definition: TVar.hh:33
TVar::simple_event_record::AssociationVCompatibility
int AssociationVCompatibility
Definition: TVar.hh:228
MelaIO::getAlphaSatMZ
double getAlphaSatMZ() const
Definition: MelaIO.h:140
TVar::ZZGG
@ ZZGG
Definition: TVar.hh:62
MELAStreamHelpers::MELAout
MELAOutputStreamer MELAout
MELAParticle::m
double m() const
Definition: MELAParticle.h:66
TVar::simple_event_record::AssociationCode
int AssociationCode
Definition: TVar.hh:227
MelaPConstant::Eval
double Eval(const MelaIO *RcdME, TVar::VerbosityLevel verbosity) const
Definition: MELAPConstant.cc:59
TVar::HSMHiggs
@ HSMHiggs
Definition: TVar.hh:125
MELAParticle::genStatus
int genStatus
Definition: MELAParticle.h:20
MelaPConstant::GetFcnFromFile
void GetFcnFromFile(const char *path, const char *spname)
Definition: MELAPConstant.cc:37
MelaPConstant::GetHPropagator
double GetHPropagator(const MelaIO *RcdME, const TVar::VerbosityLevel &verbosity) const
Definition: MELAPConstant.cc:202
TVar::MCFM
@ MCFM
Definition: TVar.hh:56
MelaIO::getHiggsMassWidth
void getHiggsMassWidth(double &mass_, double &width_, int jh) const
Definition: MelaIO.h:144
MelaPConstant::GetZPropagator
double GetZPropagator(const MelaIO *RcdME, const int &sigmaZ4lprop, const TVar::VerbosityLevel &verbosity) const
Definition: MELAPConstant.cc:214
TVar::MatrixElement
MatrixElement
Definition: TVar.hh:55
TVar::JHUGen
@ JHUGen
Definition: TVar.hh:57
TVar::Had_ZH
@ Had_ZH
Definition: TVar.hh:75
TUtil::GetMass
double GetMass(int ipart)
Definition: TUtil.cc:1582
MELAStreamHelpers::MELAerr
MELAOutputStreamer MELAerr
TVar::Had_WH_S
@ Had_WH_S
Definition: TVar.hh:87
TVar::bkgZZ
@ bkgZZ
Definition: TVar.hh:165
TVar::simple_event_record::pAssociated
SimpleParticleCollection_t pAssociated
Definition: TVar.hh:238
TVar::JQCD
@ JQCD
Definition: TVar.hh:69
MelaPConstant::fcnLow
TF1 * fcnLow
Definition: MelaPConstant.h:16
MelaPConstant::GetAlphaSatMZ
double GetAlphaSatMZ(const MelaIO *RcdME, const unsigned int &powAlphaSatMZ, const TVar::VerbosityLevel &verbosity) const
Definition: MELAPConstant.cc:300
TVar::DEBUG
@ DEBUG
Definition: TVar.hh:51
MELAStreamHelpers.hh
SimpleParticleCollection_t
std::vector< SimpleParticle_t > SimpleParticleCollection_t
Definition: TVar.hh:25
modvhiggs::propagator
complex(8) function propagator(inv_mass, mass, width)
Definition: mod_VHiggs.F90:1357
TVar::Had_WH
@ Had_WH
Definition: TVar.hh:76
MelaPConstant::GetVDaughterCouplings
double GetVDaughterCouplings(const MelaIO *RcdME, const TVar::VerbosityLevel &verbosity) const
Definition: MELAPConstant.cc:287
MelaIO::melaCand
MELACandidate * melaCand
Definition: MelaIO.h:29
TVar::Production
Production
Definition: TVar.hh:61
MelaPConstant::processME
TVar::MatrixElement processME
Definition: MelaPConstant.h:12
TVar::JJVBF_S
@ JJVBF_S
Definition: TVar.hh:83
TVar::simple_event_record
Definition: TVar.hh:226
TUtil::GetBoostedParticleVectors
void GetBoostedParticleVectors(MELACandidate *melaCand, TVar::simple_event_record &mela_event, TVar::VerbosityLevel verbosity=TVar::DEBUG)
Definition: TUtil.cc:8212
PDGHelpers::isAJet
bool isAJet(const int id)
Definition: PDGHelpers.cc:18