diff --git a/MELA/interface/Mela.h b/MELA/interface/Mela.h index 4992fcdc..b295f3bb 100755 --- a/MELA/interface/Mela.h +++ b/MELA/interface/Mela.h @@ -371,14 +371,17 @@ class Mela{ TVar::MatrixElement me_, TVar::Production prod_, TVar::Process proc_, - const char* relpath, - const char* spname + TString relpath, + TString spname, + const bool useSqrts=false ); + void deletePConstantHandle(MelaPConstant*& handle); void computeConstant(float& prob); void setConstant(); float getConstant_JHUGenUndecayed(); float getConstant_4l(); float getConstant_2l2q(); + float getConstant_FourFermionDecay(int decid); }; diff --git a/MELA/interface/MelaPConstant.h b/MELA/interface/MelaPConstant.h index c108b6c5..3d3954b0 100755 --- a/MELA/interface/MelaPConstant.h +++ b/MELA/interface/MelaPConstant.h @@ -7,6 +7,25 @@ class MelaPConstant{ +protected: + + TVar::MatrixElement processME; + TVar::Production processProd; + TVar::Process processProc; + + TF1* fcnLow; + TF1* fcnHigh; + TSpline3* fcnMid; + TString fname; + + void GetFcnFromFile(const char* path, const char* spname); + + double GetHPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const; + double GetZPropagator(const MelaIO* RcdME, const bool restricted, const TVar::VerbosityLevel& verbosity)const; + double GetAssociatedVjjPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const; + double GetVDaughterCouplings(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const; + double GetAlphaSatMZ(const MelaIO* RcdME, const int p, const TVar::VerbosityLevel& verbosity)const; + public: MelaPConstant( @@ -19,19 +38,12 @@ class MelaPConstant{ virtual ~MelaPConstant(); - double Eval(MelaIO* RcdME, TVar::VerbosityLevel verbosity)const; - -private: - - TVar::MatrixElement processME; - TVar::Production processProd; - TVar::Process processProc; + double Eval(const MelaIO* RcdME, TVar::VerbosityLevel verbosity)const; - TF1* fcnLow; - TF1* fcnHigh; - TSpline3* fcnMid; + bool IsValid(){ return bool(fcnMid!=0); } - void GetFcnFromFile(const char* path, const char* spname); + TString GetFileName(){ return fname; } + TString GetSplineName(){ TString sname=""; if (IsValid()) sname=fcnMid->GetName(); return sname; } }; diff --git a/MELA/src/MELAPConstant.cc b/MELA/src/MELAPConstant.cc index 5f4456db..33d789ae 100755 --- a/MELA/src/MELAPConstant.cc +++ b/MELA/src/MELAPConstant.cc @@ -8,6 +8,7 @@ using namespace std; +using TVar::simple_event_record; MelaPConstant::MelaPConstant( @@ -22,7 +23,8 @@ MelaPConstant::MelaPConstant( processProc(proc_), fcnLow(0), fcnHigh(0), - fcnMid(0) + fcnMid(0), + fname(path) { GetFcnFromFile(path, spname); } @@ -51,9 +53,7 @@ void MelaPConstant::GetFcnFromFile(const char* path, const char* spname){ if (fin!=0 && fin->IsOpen()) fin->Close(); } -double MelaPConstant::Eval(MelaIO* RcdME, TVar::VerbosityLevel verbosity)const{ - using TVar::simple_event_record; - +double MelaPConstant::Eval(const MelaIO* RcdME, TVar::VerbosityLevel verbosity)const{ if (verbosity>=TVar::DEBUG) cout << "Begin MelaPConstant::Eval" << endl; double result=1; @@ -432,5 +432,108 @@ double MelaPConstant::Eval(MelaIO* RcdME, TVar::VerbosityLevel verbosity)const{ return result; } +double MelaPConstant::GetHPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{ + double propagator=1.; + double sh = pow(RcdME->melaCand->m(), 2); + double mh, gah; + RcdME->getHiggsMassWidth(mh, gah, 0); + propagator = 1./(pow(sh-pow(mh, 2), 2) + pow(mh*gah, 2)); + if (verbosity>=TVar::DEBUG) cout + << "MelaPConstant::GetHPropagator: " + << "propagatorH=" << propagator << " " + << endl; + return propagator; +} +double MelaPConstant::GetZPropagator(const MelaIO* RcdME, const bool restricted, const TVar::VerbosityLevel& verbosity)const{ + double propagator=1.; + double candMass = RcdME->melaCand->m(); + double sh = pow(candMass, 2); + double mz = TUtil::GetMass(23); + double gaz = TUtil::GetMass(23); + double prop_sh = 1./(pow(sh-pow(mz, 2), 2) + pow(mz*gaz, 2)); + if (restricted){ + if (fabs(candMass-mz)<=4.*gaz){ + double shdn = pow(mz-4.*gaz, 2); + double shup = pow(mz+4.*gaz, 2); + double prop_shdn = 1./(pow(shdn-pow(mz, 2), 2) + pow(mz*gaz, 2)); + double prop_shup = 1./(pow(shup-pow(mz, 2), 2) + pow(mz*gaz, 2)); + double fsh = (sh-shdn)/(shup-shdn); + propagator = prop_sh / (prop_shdn*(1.-fsh) + prop_shup*fsh); + } + } + else propagator = prop_sh; + if (verbosity>=TVar::DEBUG) cout + << "MelaPConstant::GetZPropagator: " + << "propagatorZ=" << propagator << " " + << endl; + return propagator; +} +double MelaPConstant::GetAssociatedVjjPropagator(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{ + double propagator=1.; + if (processProd==TVar::Had_WH || processProd==TVar::Had_ZH){ + // Get a simple event record, safest way to handle jet mass corrections + int nRequested_AssociatedJets=2; + int AssociationVCompatibility=0; + int partIncCode=TVar::kUseAssociated_Jets; + double mv=0, gav=0; + if (processProd==TVar::Had_WH){ + AssociationVCompatibility=24; + mv = TUtil::GetMass(24); + gav = TUtil::GetMass(24); + } + else{ + AssociationVCompatibility=23; + mv = TUtil::GetMass(23); + gav = TUtil::GetMass(23); + } + simple_event_record mela_event; + mela_event.AssociationCode=partIncCode; + mela_event.AssociationVCompatibility=AssociationVCompatibility; + mela_event.nRequested_AssociatedJets=nRequested_AssociatedJets; + TUtil::GetBoostedParticleVectors( + RcdME->melaCand, + mela_event, + verbosity + ); + + float mJJval=-1; + vector pJets; + const SimpleParticleCollection_t& pAssociated = mela_event.pAssociated; + for (auto& part : pAssociated){ + if (PDGHelpers::isAJet(part.first)) pJets.push_back(part.second); + if (pJets.size()==2) break; + } + if (pJets.size()==2) mJJval = (pJets[0] + pJets[1]).M(); + if (mJJval>=0.) propagator = 1./(pow(pow(mJJval, 2)-pow(mv, 2), 2) + pow(mv*gav, 2)); + } + if (verbosity>=TVar::DEBUG) cout + << "MelaPConstant::GetAssociatedVjjPropagator: " + << "propagatorV=" << propagator << " " + << endl; + return propagator; +} +double MelaPConstant::GetVDaughterCouplings(const MelaIO* RcdME, const TVar::VerbosityLevel& verbosity)const{ + double result=1; + double aL1, aR1, aL2, aR2; + RcdME->getVDaughterCouplings(aL1, aR1, 0); + RcdME->getVDaughterCouplings(aL2, aR2, 1); + if (verbosity>=TVar::DEBUG) cout + << "MelaPConstant::GetVDaughterCouplings: " + << "L**2+R**2 couplings=" << pow(aL1, 2)+pow(aR1, 2) << " " << pow(aL2, 2)+pow(aR2, 2) << " " + << endl; + if (fabs(aL1)>0. || fabs(aR1)>0.) result *= pow(aL1, 2)+pow(aR1, 2); + if (fabs(aL2)>0. || fabs(aR2)>0.) result *= pow(aL2, 2)+pow(aR2, 2); + return result; +} +double MelaPConstant::GetAlphaSatMZ(const MelaIO* RcdME, const int p, const TVar::VerbosityLevel& verbosity)const{ + double result; + double alphasVal = RcdME->getAlphaSatMZ(); + result = pow(alphasVal, p); + if (verbosity>=TVar::DEBUG) cout + << "MelaPConstant::GetAlphaSatMZ: " + << "alphas(MZ)=" << alphasVal << " (**" << p << ") " + << endl; + return result; +} diff --git a/MELA/src/Mela.cc b/MELA/src/Mela.cc index dac2c8b2..f3817c80 100755 --- a/MELA/src/Mela.cc +++ b/MELA/src/Mela.cc @@ -1667,26 +1667,26 @@ float Mela::getConstant_JHUGenUndecayed(){ MelaPConstant* pchandle=0; unsigned int iarray=0; - double correction=1; - + if (TUtil::JetMassScheme == TVar::ConserveDifermionMass) iarray=0; // First element points to the case when the difermion invariant mass is conserved in mass removal scheme + else if (TUtil::JetMassScheme == TVar::MomentumToEnergy) iarray=1; // Second element points to the case when the 3-momentum vector magnitude is scaled to energy in mass removal scheme if (myProduction_ == TVar::JQCD){ - if (TUtil::JetMassScheme == TVar::ConserveDifermionMass) iarray=0; // First element points to the case when the difermion invariant mass is conserved in mass removal scheme - else if (TUtil::JetMassScheme == TVar::MomentumToEnergy) iarray=1; // Second element points to the case when the 3-momentum vector magnitude is scaled to energy in mass removal scheme pchandle = pAvgSmooth_JHUGen_JQCD_HSMHiggs[iarray]; } else if (myProduction_ == TVar::JJQCD){ - if (TUtil::JetMassScheme == TVar::ConserveDifermionMass) iarray=0; // First element points to the case when the difermion invariant mass is conserved in mass removal scheme - else if (TUtil::JetMassScheme == TVar::MomentumToEnergy) iarray=1; // Second element points to the case when the 3-momentum vector magnitude is scaled to energy in mass removal scheme pchandle = pAvgSmooth_JHUGen_JJQCD_HSMHiggs[iarray]; } else if (myProduction_ == TVar::JJVBF){ - if (TUtil::JetMassScheme == TVar::ConserveDifermionMass) iarray=0; // First element points to the case when the difermion invariant mass is conserved in mass removal scheme - else if (TUtil::JetMassScheme == TVar::MomentumToEnergy) iarray=1; // Second element points to the case when the 3-momentum vector magnitude is scaled to energy in mass removal scheme pchandle = pAvgSmooth_JHUGen_JJVBF_HSMHiggs[iarray]; } + else if (myProduction_ == TVar::Had_ZH){ + pchandle = pAvgSmooth_JHUGen_Had_ZH_HSMHiggs[iarray]; + } + else if (myProduction_ == TVar::Had_WH){ + pchandle = pAvgSmooth_JHUGen_Had_WH_HSMHiggs[iarray]; + } /* - else if (myProduction_ == TVar::Lep_ZH || myProduction_ == TVar::Had_ZH) - else if (myProduction_ == TVar::Lep_WH || myProduction_ == TVar::Had_WH) + else if (myProduction_ == TVar::Lep_ZH) + else if (myProduction_ == TVar::Lep_WH) else if (myProduction_ == TVar::GammaH) else if (myProduction_ == TVar::ttH) else if (myProduction_ == TVar::bbH) @@ -1694,102 +1694,42 @@ float Mela::getConstant_JHUGenUndecayed(){ else return constant; constant = pchandle->Eval(getIORecord(), myVerbosity_); - if (myProduction_==TVar::JJVBF && LHCsqrts==7.){ - // Fitting is not good enough for region starting at ~105 GeV due to poor statistics. - const double a0=0.67; - const double a1=22.; - const double a2=73.; - double var = melaCand->m(); - if (var>a2) correction = 1.+a0*exp(-pow((var-a2)/a1, 2)); - else correction = 1.+a0; // Smooth by virtue of the correction function itself. - } - else if (myProduction_==TVar::JJVBF && LHCsqrts==8.){ - // Fitting is not good enough for region starting at ~105 GeV due to poor statistics. - const double a0=0.53; - const double a1=21.; - const double a2=73.; - double var = melaCand->m(); - if (var>a2) correction = 1.+a0*exp(-pow((var-a2)/a1, 2)); - else correction = 1.+a0; // Smooth by virtue of the correction function itself. - } - else if (myProduction_==TVar::JJVBF && LHCsqrts==13.){ - // Fitting is not good enough for region starting at ~105 GeV due to poor statistics. - const double a0=0.2; - const double a1=22.; - const double a2=73.; - double var = melaCand->m(); - if (var>a2) correction = 1.+a0*exp(-pow((var-a2)/a1, 2)); - else correction = 1.+a0; // Smooth by virtue of the correction function itself. - } - // - else if (myProduction_==TVar::JJQCD && LHCsqrts==8.){ - // Fitting is not good enough for region starting at ~105 GeV due to poor statistics. - const double a0=-0.24; - const double a1=80.; - const double a2=9.; - const double a3=0.08; - const double a4=100.; - const double a5=20.; - double var = melaCand->m(); - if (var>a1) correction = 1+a0*exp(-pow((var-a1)/a2, 2))+a3*exp(-pow((var-a4)/a5, 2)); - else correction = 1.+a0+a3*exp(-pow((var-a4)/a5, 2)); // Smooth by virtue of the correction function itself. - } - // - else if (myProduction_==TVar::JQCD && LHCsqrts==7.){ - const double a0=-0.5; - const double a1=80.; - const double a2=9.; - const double a3=-0.35352; - const double a4=1500.; - const double a5=268.; - double var = melaCand->m(); - if (var>a1 && vara1) correction = 1+a0*exp(-pow((var-a1)/a2, 2))+a3; - else correction = 1.+a0+a3*exp(-pow((var-a4)/a5, 2)); - } - else if (myProduction_==TVar::JQCD && LHCsqrts==8.){ - const double a0=-0.2; - const double a1=80.; - const double a2=9.; - const double a3=-0.0792; - const double a4=1500.; - const double a5=615.; - double var = melaCand->m(); - if (var>a1 && vara1) correction = 1+a0*exp(-pow((var-a1)/a2, 2))+a3; - else correction = 1.+a0+a3*exp(-pow((var-a4)/a5, 2)); - } - else if (myProduction_==TVar::JQCD && LHCsqrts==13.){ - // 1+[0]*exp(-pow((x-[1])/[2],2))+[3]*exp(-pow((x-[4])/[5],2)) - const double a0=0.15; - const double a1=320.; - const double a2=300.; - const double a3=0.179; - const double a4=1530.; - const double a5=212.; - const double offset=0.8444; - double var = melaCand->m(); - if (vargetSortedV(0)->getDaughter(0)->id)* abs(melaCand->getSortedV(0)->getDaughter(1)->id)* abs(melaCand->getSortedV(1)->getDaughter(0)->id)* abs(melaCand->getSortedV(1)->getDaughter(1)->id); - const bool is4mu = (idprod==28561); - const bool is4e = (idprod==14641 || idprod==50625); // Use 4e for 4tau as well (I don't know why you would do this, but anyway - const bool is2mu2e = (idprod==20449 || idprod==27225 || idprod==38025); // Use 2e2mu for 2e2tau and 2mu2tau as well + return getConstant_FourFermionDecay(decid); +} +float Mela::getConstant_2l2q(){ + float constant = 1; + if (melaCand==0) return constant; + int decid1 = + abs(melaCand->getSortedV(0)->getDaughter(0)->id)* + abs(melaCand->getSortedV(0)->getDaughter(1)->id); + int decid2 = + abs(melaCand->getSortedV(1)->getDaughter(0)->id)* + abs(melaCand->getSortedV(1)->getDaughter(1)->id); + int decid = 1; + if (decid1!=0) decid*=decid1; + if (decid2!=0) decid*=decid2; + return getConstant_FourFermionDecay(decid); +} + +float Mela::getConstant_FourFermionDecay(int decid){ + float constant = 1; + + const bool is4mu = (decid==28561); + const bool is4e = (decid==14641 || decid==50625); // Use 4e for 4tau as well (I don't know why you would do this, but anyway) + const bool is2mu2e = (decid==20449 || decid==27225 || decid==38025 || decid==169 || decid==121); // Use 2e2mu for 2e2tau and 2mu2tau as well + + const unsigned int nPossibleHandles=6; + MelaPConstant* pchandle[nPossibleHandles]={ 0 }; float constant_tmp=0; if (myME_ == TVar::JHUGen){ @@ -1822,6 +1762,7 @@ float Mela::getConstant_4l(){ } } else if (myME_ == TVar::MCFM){ + // ZZQQB if (myProduction_ == TVar::ZZQQB){ if (myModel_ == TVar::bkgZZ){ if (is2mu2e) pchandle[0] = pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e; @@ -1829,6 +1770,7 @@ float Mela::getConstant_4l(){ else if (is4e) pchandle[0] = pAvgSmooth_MCFM_ZZQQB_bkgZZ_4e; } } + // ZZGG else if (myProduction_ == TVar::ZZGG){ if (myModel_ == TVar::bkgZZ){ if (is2mu2e) pchandle[0] = pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e; @@ -1840,6 +1782,7 @@ float Mela::getConstant_4l(){ else if (is4mu) pchandle[0] = pAvgSmooth_MCFM_ZZGG_HSMHiggs_4mu; else if (is4e) pchandle[0] = pAvgSmooth_MCFM_ZZGG_HSMHiggs_4e; } + // Sum signal and bkg else if (myModel_ == TVar::bkgZZ_SMHiggs){ if (is2mu2e){ pchandle[0] = pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e; @@ -1855,57 +1798,64 @@ float Mela::getConstant_4l(){ } } } - else if (myProduction_ == TVar::JJQCD){ - if (myModel_ == TVar::bkgZJets){ - pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q; // Only option at the moment + // JJEW and components + else if (myProduction_ == TVar::JJVBF || myProduction_ == TVar::Had_WH || myProduction_ == TVar::Had_ZH || myProduction_ == TVar::JJEW){ + MelaPConstant* hvbf=0; + MelaPConstant* hwh=0; + MelaPConstant* hzh=0; + MelaPConstant* hvbs=0; + MelaPConstant* hwzz=0; + MelaPConstant* hzzz=0; + if (is2mu2e){ + hvbf = pAvgSmooth_MCFM_JJVBF_HSMHiggs_2mu2e; + hwh = pAvgSmooth_MCFM_Had_WH_HSMHiggs_2mu2e; + hzh = pAvgSmooth_MCFM_Had_ZH_HSMHiggs_2mu2e; + hvbs = pAvgSmooth_MCFM_JJVBF_bkgZZ_2mu2e; + hwzz = pAvgSmooth_MCFM_Had_WH_bkgZZ_2mu2e; + hzzz = pAvgSmooth_MCFM_Had_ZH_bkgZZ_2mu2e; } - } - } - - bool hasNullHandle=true; - for (unsigned int ihandle=0; ihandleEval(getIORecord(), myVerbosity_); hasNullHandle=false; } } - if (hasNullHandle) return constant; - - constant = constant_tmp; - return constant; -} -float Mela::getConstant_2l2q(){ - float constant = 1; - if (melaCand==0) return constant; - - const unsigned int nPossibleHandles=2; - MelaPConstant* pchandle[nPossibleHandles]={ 0 }; - - float constant_tmp=0; - // Most constants use the 2e2mu constant. MelaPConstant scales for the left/right couplings itself based on what is recorded into the MelaIO object. - if (myME_ == TVar::JHUGen){ - if (myProduction_ == TVar::ZZGG){ - if (myModel_ == TVar::HSMHiggs){ - pchandle[0] = pAvgSmooth_JHUGen_ZZGG_HSMHiggs_2mu2e; + else if (is4mu){ + hvbf = pAvgSmooth_MCFM_JJVBF_HSMHiggs_4mu; + hwh = pAvgSmooth_MCFM_Had_WH_HSMHiggs_4mu; + hzh = pAvgSmooth_MCFM_Had_ZH_HSMHiggs_4mu; + hvbs = pAvgSmooth_MCFM_JJVBF_bkgZZ_4mu; + hwzz = pAvgSmooth_MCFM_Had_WH_bkgZZ_4mu; + hzzz = pAvgSmooth_MCFM_Had_ZH_bkgZZ_4mu; } - } - } - else if (myME_ == TVar::MCFM){ - if (myProduction_ == TVar::ZZQQB){ - if (myModel_ == TVar::bkgZZ){ - pchandle[0] = pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e; + else if (is4e){ + hvbf = pAvgSmooth_MCFM_JJVBF_HSMHiggs_4e; + hwh = pAvgSmooth_MCFM_Had_WH_HSMHiggs_4e; + hzh = pAvgSmooth_MCFM_Had_ZH_HSMHiggs_4e; + hvbs = pAvgSmooth_MCFM_JJVBF_bkgZZ_4e; + hwzz = pAvgSmooth_MCFM_Had_WH_bkgZZ_4e; + hzzz = pAvgSmooth_MCFM_Had_ZH_bkgZZ_4e; } - } - else if (myProduction_ == TVar::ZZGG){ - if (myModel_ == TVar::bkgZZ){ - pchandle[0] = pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e; - } - else if (myModel_ == TVar::HSMHiggs){ - pchandle[0] = pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e; + + if (myModel_ == TVar::HSMHiggs || myModel_ == TVar::bkgZZ_SMHiggs){ + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::JJVBF) pchandle[0]=hvbf; + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::Had_ZH) pchandle[1]=hzh; + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::Had_WH) pchandle[2]=hwh; } - else if (myModel_ == TVar::bkgZZ_SMHiggs){ - pchandle[0] = pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e; - pchandle[1] = pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e; + if (myModel_ == TVar::bkgZZ || myModel_ == TVar::bkgZZ_SMHiggs){ + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::JJVBF) pchandle[3]=hvbs; + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::Had_ZH) pchandle[4]=hzzz; + if (myProduction_ == TVar::JJEW || myProduction_ == TVar::Had_WH) pchandle[5]=hwzz; } } else if (myProduction_ == TVar::JJQCD){ if (myModel_ == TVar::bkgZJets){ - pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q; + pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q; // Only option at the moment + } + else if (myModel_ == TVar::bkgZZ){ + if (is2mu2e){ + pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZZ_2mu2e; + } + else if (is4mu){ + pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZZ_4mu; + } + else if (is4e){ + pchandle[0] = pAvgSmooth_MCFM_JJQCD_bkgZZ_4e; + } } } } @@ -1917,33 +1867,22 @@ float Mela::getConstant_2l2q(){ constant = constant_tmp; return constant; } + + void Mela::getPConstantHandles(){ if (myVerbosity_>=TVar::DEBUG) cout << "Begin Mela::getPConstantHandles" << endl; - // Find closest sqrts allowable - const unsigned int npossiblesqrts=3; - const double possible_sqrts[npossiblesqrts]={ 7, 8, 13 }; - unsigned int sqrts_index=0; - double sqrtsdiff = 99.; // Some large number - for (unsigned isq=0; isq=TVar::DEBUG) cout << "End Mela::getPConstantHandles" << endl; @@ -2026,11 +2068,16 @@ MelaPConstant* Mela::getPConstantHandle( TVar::MatrixElement me_, TVar::Production prod_, TVar::Process proc_, - const char* relpath, - const char* spname + TString relpath, + TString spname, + const bool useSqrts ){ if (myVerbosity_>=TVar::DEBUG) cout << "Begin Mela::getPConstantHandle" << endl; + MelaPConstant* pchandle=0; + string cfile_fullpath; + + // Get data/ path if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: relpath and spline name: " << relpath << ", " << spname << endl; #ifdef _melapkgpathstr_ const string MELAPKGPATH = _melapkgpathstr_; @@ -2039,45 +2086,115 @@ MelaPConstant* Mela::getPConstantHandle( assert(0); #endif const string path = MELAPKGPATH + "data/"; - string cfile_fullpath = path; - cfile_fullpath.append(relpath); - cfile_fullpath.append(".root"); if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: path and spline name: " << path << ", " << spname << endl; - if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: Full path and spline name: " << cfile_fullpath << ", " << spname << endl; - MelaPConstant* pchandle = new MelaPConstant(me_, prod_, proc_, cfile_fullpath.c_str(), spname); + if (useSqrts){ // Loop over possible sqrts values to get the closest one + const unsigned int npossiblesqrts=3; + const double possible_sqrts[npossiblesqrts]={ 7, 8, 13 }; + vector trysqrts; + for (unsigned isq=0; isq::iterator it = trysqrts.begin(); itdiff){ + inserted=true; + trysqrts.insert(it, val); + break; + } + } + if (!inserted) trysqrts.push_back(val); + } + for (auto& dsqrts : trysqrts){ + TString strsqrts = Form("%s_%.0f%s", relpath.Data(), dsqrts, "TeV"); + cfile_fullpath = path; + cfile_fullpath.append(strsqrts.Data()); + cfile_fullpath.append(".root"); + pchandle = new MelaPConstant(me_, prod_, proc_, cfile_fullpath.c_str(), spname.Data()); + if (pchandle->IsValid()){ + if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: Full path and spline name: " << cfile_fullpath << ", " << spname << " is valid." << endl; + break; + } + else{ + if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: Full path and spline name: " << cfile_fullpath << ", " << spname << " is invalid." << endl; + deletePConstantHandle(pchandle); + } + } + } + else{ + cfile_fullpath = path; + cfile_fullpath.append(relpath.Data()); + cfile_fullpath.append(".root"); + pchandle = new MelaPConstant(me_, prod_, proc_, cfile_fullpath.c_str(), spname.Data()); + if (!pchandle->IsValid()) deletePConstantHandle(pchandle); + } + if (myVerbosity_>=TVar::DEBUG) cout << "Mela::getPConstantHandle: Full path and spline name: " << cfile_fullpath << ", " << spname << endl; + if (myVerbosity_>=TVar::DEBUG && pchandle==0) cerr << "Mela::getPConstantHandle: Handle of " << spname << " from " << cfile_fullpath << " is invalid!" << endl; if (myVerbosity_>=TVar::DEBUG) cout << "End Mela::getPConstantHandle" << endl; return pchandle; } void Mela::deletePConstantHandles(){ - for (unsigned int isch=0; isch<(unsigned int)(TVar::nFermionMassRemovalSchemes-1); isch++){ - if (pAvgSmooth_JHUGen_JJQCD_HSMHiggs[isch]!=0) delete pAvgSmooth_JHUGen_JJQCD_HSMHiggs[isch]; - // - if (pAvgSmooth_JHUGen_JJVBF_HSMHiggs[isch]!=0) delete pAvgSmooth_JHUGen_JJVBF_HSMHiggs[isch]; - // - if (pAvgSmooth_JHUGen_JQCD_HSMHiggs[isch]!=0) delete pAvgSmooth_JHUGen_JQCD_HSMHiggs[isch]; - // + for (int isch=(unsigned int)(TVar::nFermionMassRemovalSchemes-2); isch>=0; isch--){ + if (isch==0 || pAvgSmooth_JHUGen_JJQCD_HSMHiggs[isch]!=pAvgSmooth_JHUGen_JJQCD_HSMHiggs[0]) deletePConstantHandle(pAvgSmooth_JHUGen_JJQCD_HSMHiggs[isch]); + if (isch==0 || pAvgSmooth_JHUGen_JQCD_HSMHiggs[isch]!=pAvgSmooth_JHUGen_JQCD_HSMHiggs[0]) deletePConstantHandle(pAvgSmooth_JHUGen_JQCD_HSMHiggs[isch]); + if (isch==0 || pAvgSmooth_JHUGen_JJVBF_HSMHiggs[isch]!=pAvgSmooth_JHUGen_JJVBF_HSMHiggs[0]) deletePConstantHandle(pAvgSmooth_JHUGen_JJVBF_HSMHiggs[isch]); + if (isch==0 || pAvgSmooth_JHUGen_Had_ZH_HSMHiggs[isch]!=pAvgSmooth_JHUGen_Had_ZH_HSMHiggs[0]) deletePConstantHandle(pAvgSmooth_JHUGen_Had_ZH_HSMHiggs[isch]); + if (isch==0 || pAvgSmooth_JHUGen_Had_WH_HSMHiggs[isch]!=pAvgSmooth_JHUGen_Had_WH_HSMHiggs[0]) deletePConstantHandle(pAvgSmooth_JHUGen_Had_WH_HSMHiggs[isch]); } // - if (pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q!=0) delete pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q; + deletePConstantHandle(pAvgSmooth_MCFM_JJQCD_bkgZJets_2l2q); + // + deletePConstantHandle(pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4mu); + deletePConstantHandle(pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4e); + deletePConstantHandle(pAvgSmooth_JHUGen_ZZGG_HSMHiggs_2mu2e); // - if (pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4mu!=0) delete pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4mu; - if (pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4e!=0) delete pAvgSmooth_JHUGen_ZZGG_HSMHiggs_4e; - if (pAvgSmooth_JHUGen_ZZGG_HSMHiggs_2mu2e!=0) delete pAvgSmooth_JHUGen_ZZGG_HSMHiggs_2mu2e; + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_HSMHiggs_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_HSMHiggs_4e); + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e); // - if (pAvgSmooth_MCFM_ZZGG_HSMHiggs_4mu!=0) delete pAvgSmooth_MCFM_ZZGG_HSMHiggs_4mu; - if (pAvgSmooth_MCFM_ZZGG_HSMHiggs_4e!=0) delete pAvgSmooth_MCFM_ZZGG_HSMHiggs_4e; - if (pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e!=0) delete pAvgSmooth_MCFM_ZZGG_HSMHiggs_2mu2e; + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_HSMHiggs_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_HSMHiggs_4e); + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_HSMHiggs_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_HSMHiggs_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_HSMHiggs_4e); + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_HSMHiggs_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_HSMHiggs_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_HSMHiggs_4e); + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_HSMHiggs_2mu2e); // - if (pAvgSmooth_MCFM_ZZGG_bkgZZ_4mu!=0) delete pAvgSmooth_MCFM_ZZGG_bkgZZ_4mu; - if (pAvgSmooth_MCFM_ZZGG_bkgZZ_4e!=0) delete pAvgSmooth_MCFM_ZZGG_bkgZZ_4e; - if (pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e!=0) delete pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e; + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_ZZGG_bkgZZ_2mu2e); // - if (pAvgSmooth_MCFM_ZZQQB_bkgZZ_4mu!=0) delete pAvgSmooth_MCFM_ZZQQB_bkgZZ_4mu; - if (pAvgSmooth_MCFM_ZZQQB_bkgZZ_4e!=0) delete pAvgSmooth_MCFM_ZZQQB_bkgZZ_4e; - if (pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e!=0) delete pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e; + deletePConstantHandle(pAvgSmooth_MCFM_ZZQQB_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_ZZQQB_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_ZZQQB_bkgZZ_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_JJVBF_bkgZZ_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_Had_ZH_bkgZZ_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_Had_WH_bkgZZ_2mu2e); + // + deletePConstantHandle(pAvgSmooth_MCFM_JJQCD_bkgZZ_4mu); + deletePConstantHandle(pAvgSmooth_MCFM_JJQCD_bkgZZ_4e); + deletePConstantHandle(pAvgSmooth_MCFM_JJQCD_bkgZZ_2mu2e); // } +void Mela::deletePConstantHandle(MelaPConstant*& handle){ + if (myVerbosity_>=TVar::DEBUG) cout << "Mela::deletePConstantHandle: Deleting PConstant handle " << handle->GetSplineName() << " at " << handle->GetFileName() << endl; + delete handle; handle=0; + if (myVerbosity_>=TVar::DEBUG) cout << "End Mela::deletePConstantHandle." << endl; +} + void Mela::computeDijetConvBW(float& prob){ melaCand = getCurrentCandidate(); diff --git a/MELA/test/calcC_lintolog.c b/MELA/test/calcC_lintolog.c index f0488d39..0bb39cf4 100755 --- a/MELA/test/calcC_lintolog.c +++ b/MELA/test/calcC_lintolog.c @@ -3229,7 +3229,7 @@ void get_PAvgProfile_MCFM_JJPROD_HSMHiggs_13TeV(TString strprod, int sqrts=13, b /* SPECIFIC COMMENT: OUTPUT ME DIVIDED BY - (aL1**2+aR1**2)*(aL2**2+aR2**2) TO REMAIN INDEPENDENT OF CHANNEL -- PROP(Z)/LINE_PROP(Z) AROUND M4L~=MZ (FIXME: NOT IMPLEMENTED YET) +- PROP(Z)/LINE_PROP(Z) AROUND M4L~=MZ - mJJ propagator is taken out in WH or ZH */ void get_PAvgProfile_MCFM_JJPROD_bkgZZ_13TeV(TString strprod, int sqrts=13, bool recalculate = true){ @@ -3602,6 +3602,23 @@ void get_PAvgProfile_MCFM_JJPROD_bkgZZ_13TeV(TString strprod, int sqrts=13, bool propagatorV = 1./(pow(pow(mjj, 2)-pow(mv, 2), 2) + pow(mv*gav, 2)); mesq_conserveDifermMass /= propagatorV; } + + double mz, gaz, propagator; + mz = mela.getPrimaryMass(23); + gaz = mela.getPrimaryWidth(23); + if (fabs(mzz-mz)<=4.*gaz){ + double sh = pow(mzz, 2); + double shdn = pow(mz-4.*gaz, 2); + double shup = pow(mz+4.*gaz, 2); + double prop_sh = 1./(pow(sh-pow(mz, 2), 2) + pow(mz*gaz, 2)); + double prop_shdn = 1./(pow(shdn-pow(mz, 2), 2) + pow(mz*gaz, 2)); + double prop_shup = 1./(pow(shup-pow(mz, 2), 2) + pow(mz*gaz, 2)); + double fsh = (sh-shdn)/(shup-shdn); + propagator = prop_sh / (prop_shdn*(1.-fsh) + prop_shup*fsh); + } + else propagator=1.; + mesq_conserveDifermMass /= propagator; + mela.resetInputEvent(); if (doRewgt){ // LEFT HERE: Implement reweighting diff --git a/MELA/test/testME_v2.c b/MELA/test/testME_v2.c index 2edb0213..f06e1c9f 100755 --- a/MELA/test/testME_v2.c +++ b/MELA/test/testME_v2.c @@ -1069,7 +1069,7 @@ void testME_Dec_MCFM_Ping(int flavor=2, int useMothers=0, bool useConstants=fals tout.close(); } -void testME_VH_JHUGen_Ping(int erg_tev=13){ +void testME_VH_JHUGen_Ping(int erg_tev=13, bool useConstants=false){ TString strtout = Form("testME_VH_JHUGen_%iTeV_Ping.out", erg_tev); ofstream tout(strtout.Data()); streambuf* coutbuf = cout.rdbuf(); @@ -1198,69 +1198,69 @@ void testME_VH_JHUGen_Ping(int erg_tev=13){ float p0mplus=0; mela.setProcess(TVar::HSMHiggs, TVar::JHUGen, prod); - mela.computeProdP_VH(p0mplus, false, false); + mela.computeProdP_VH(p0mplus, false, useConstants); cout << "p0mplus: " << p0mplus << '\n' << endl; float p0g1prime2=0; mela.setProcess(TVar::H0_g1prime2, TVar::JHUGen, prod); - mela.computeProdP_VH(p0g1prime2, false, false); + mela.computeProdP_VH(p0g1prime2, false, useConstants); cout << "p0g1prime2: " << p0g1prime2 << '\n' << endl; float p0hplus=0; mela.setProcess(TVar::H0hplus, TVar::JHUGen, prod); - mela.computeProdP_VH(p0hplus, false, false); + mela.computeProdP_VH(p0hplus, false, useConstants); cout << "p0hplus: " << p0hplus << '\n' << endl; float p0minus=0; mela.setProcess(TVar::H0minus, TVar::JHUGen, prod); - mela.computeProdP_VH(p0minus, false, false); + mela.computeProdP_VH(p0minus, false, useConstants); cout << "p0minus: " << p0minus << '\n' << endl; float p0gzgs1prime2=0; mela.setProcess(TVar::H0_Zgsg1prime2, TVar::JHUGen, prod); - mela.computeProdP_VH(p0gzgs1prime2, false, false); + mela.computeProdP_VH(p0gzgs1prime2, false, useConstants); cout << "p0gzgs1prime2: " << p0gzgs1prime2 << '\n' << endl; float p0hpluszgs=0; mela.setProcess(TVar::H0_Zgs, TVar::JHUGen, prod); - mela.computeProdP_VH(p0hpluszgs, false, false); + mela.computeProdP_VH(p0hpluszgs, false, useConstants); cout << "p0hpluszgs: " << p0hpluszgs << '\n' << endl; // SelfD MEs float p0mplus_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_VV_1][0]=1; - mela.computeProdP_VH(p0mplus_selfD, false, false); + mela.computeProdP_VH(p0mplus_selfD, false, useConstants); cout << "p0mplus_selfD: " << p0mplus_selfD << '\n' << endl; float p0g1prime2_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_VV_1_PRIME2][0]=1; - mela.computeProdP_VH(p0g1prime2_selfD, false, false); + mela.computeProdP_VH(p0g1prime2_selfD, false, useConstants); cout << "p0g1prime2_selfD: " << p0g1prime2_selfD << '\n' << endl; float p0hplus_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_VV_2][0]=1; - mela.computeProdP_VH(p0hplus_selfD, false, false); + mela.computeProdP_VH(p0hplus_selfD, false, useConstants); cout << "p0hplus_selfD: " << p0hplus_selfD << '\n' << endl; float p0minus_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_VV_4][0]=1; - mela.computeProdP_VH(p0minus_selfD, false, false); + mela.computeProdP_VH(p0minus_selfD, false, useConstants); cout << "p0minus_selfD: " << p0minus_selfD << '\n' << endl; float p0gzgs1prime2_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_ZA_1_PRIME2][0]=1; - mela.computeProdP_VH(p0gzgs1prime2_selfD, false, false); + mela.computeProdP_VH(p0gzgs1prime2_selfD, false, useConstants); cout << "p0gzgs1prime2_selfD: " << p0gzgs1prime2_selfD << '\n' << endl; float p0hpluszgs_selfD=0; mela.setProcess(TVar::SelfDefine_spin0, TVar::JHUGen, prod); mela.selfDHzzcoupl[0][gHIGGS_ZA_2][0]=1; - mela.computeProdP_VH(p0hpluszgs_selfD, false, false); + mela.computeProdP_VH(p0hpluszgs_selfD, false, useConstants); cout << "p0hpluszgs_selfD: " << p0hpluszgs_selfD << '\n' << endl; if (prod==TVar::Had_ZH || prod==TVar::Had_WH){