#include "PHASIC++/Process/Process_Info.H" #include "PHASIC++/Process/Virtual_ME2_Base.H" #include "MODEL/Main/Running_AlphaS.H" #include "AddOns/MCFM/MCFM_Wrapper.H" namespace MCFM { // README: // For Higgs production, choose model: MODEL = SMEHC // It is important for the Higgs production to have all five flavours // in the initial state, but the Yukawa coupling of the b must be // switched off: YUKAWA[5] = 0. // Also, MCFM acts in the limit of mt->infinity, // thus a further correction term has been introduced // // We could actually also extend this to BSM models. class MCFM_gg_hg: public PHASIC::Virtual_ME2_Base { private: int m_pID; double * p_p, *p_msqv; MODEL::Running_AlphaS * p_as; double m_mh2,m_Gh2,m_mZ2,m_GZ2,m_mW2,m_GW2; double m_ehcscale2,m_cplcorr,m_normcorr; double CallMCFM(const int & i,const int & j); public: MCFM_gg_hg(const int & pID, const PHASIC::Process_Info& pi, const ATOOLS::Flavour_Vector& flavs); ~MCFM_gg_hg(); void Calc(const ATOOLS::Vec4D_Vector& momenta); double Eps_Scheme_Factor(const ATOOLS::Vec4D_Vector& mom); }; }// end of namespace MCFM extern "C" { void gg_hg_v_(double *p,double *msqv); void gg_hwwg_v_(double *p,double *msqv); void gg_hzzg_v_(double *p,double *msqv); void gg_hgagag_v_(double *p,double *msqv); } #include "MODEL/Main/Model_Base.H" #include "ATOOLS/Org/Run_Parameter.H" using namespace MCFM; using namespace PHASIC; using namespace ATOOLS; MCFM_gg_hg::MCFM_gg_hg(const int & pID,const Process_Info& pi, const Flavour_Vector& flavs) : Virtual_ME2_Base(pi,flavs), m_pID(pID), p_as((MODEL::Running_AlphaS *) MODEL::s_model->GetScalarFunction(std::string("alpha_S"))), m_mh2(sqr(Flavour(kf_h0).Mass())), m_Gh2(sqr(Flavour(kf_h0).Width())), m_mZ2(sqr(Flavour(kf_Z).Mass())), m_GZ2(sqr(Flavour(kf_Z).Width())), m_mW2(sqr(Flavour(kf_Wplus).Mass())), m_GW2(sqr(Flavour(kf_Wplus).Width())), m_ehcscale2(MODEL::s_model->ScalarConstant(std::string("EHC_SCALE2"))), m_cplcorr(ewcouple_.vevsq/ sqr(MODEL::s_model->ScalarConstant(std::string("vev")))* sqr((*p_as)(m_ehcscale2)/qcdcouple_.as)* sqr(MODEL::s_model->ScalarConstant(std::string("h0_gg_fac"))/ (2./3.))), m_normcorr(4.*9./qcdcouple_.ason2pi) { rpa->gen.AddCitation (1,"The NLO matrix elements have been taken from MCFM \\cite{}."); switch (m_pID) { case 204: m_cplcorr *= sqr(Flavour(kf_tau).Yuk())/masses_.mbsq/3. * 4.*sqr(masses_.wmass)/ewcouple_.gwsq/ sqr(MODEL::s_model->ScalarConstant(std::string("vev"))); break; case 208: case 209: m_cplcorr *= pow(4.*M_PI*MODEL::s_model->ScalarConstant("alpha_QED")/ std::abs(MODEL::s_model->ComplexConstant("csin2_thetaW"))/ ewcouple_.gwsq,3.); break; case 210: m_cplcorr *= sqr(MODEL::s_model->ScalarFunction(std::string("alpha_QED"),m_ehcscale2)/ (ewcouple_.esq/(4.*M_PI))) * 1./(sqrt(2.)*ewcouple_.Gf* sqr(MODEL::s_model->ScalarConstant(std::string("vev")))); break; } p_p = new double[4*MCFM_NMX]; p_msqv = new double[sqr(2*MCFM_NF+1)]; m_mode=1; } MCFM_gg_hg::~MCFM_gg_hg() { delete [] p_p; delete [] p_msqv; } double MCFM_gg_hg::CallMCFM(const int & i,const int & j) { //msg_Out()<<METHOD<<"("<<i<<", "<<j<<") for "<<m_pID<<".\n"; switch (m_pID) { case 204: gg_hg_v_(p_p,p_msqv); break; case 208: gg_hwwg_v_(p_p,p_msqv); break; case 209: gg_hzzg_v_(p_p,p_msqv); break; case 210: gg_hgagag_v_(p_p,p_msqv); break; } return p_msqv[mr(i,j)]; } void MCFM_gg_hg::Calc(const Vec4D_Vector &p) { double corrfactor(m_cplcorr*m_normcorr*(*p_as)(m_mur2)/qcdcouple_.as); double sh((m_pID==204)? (p[2]+p[3]).Abs2() : (p[2]+p[3]+p[4]+p[5]).Abs2()); corrfactor *= (sqr(sh-sqr(masses_.hmass))+ sqr(masses_.hmass*masses_.hwidth))/ (sqr(sh-m_mh2)+m_mh2*m_Gh2); if (m_pID==208 || m_pID==209) { double s1((p[2]+p[3]).Abs2()),s2((p[4]+p[5]).Abs2()); if (m_pID==208) { corrfactor *= (sqr(s1-sqr(masses_.wmass))+ sqr(masses_.wmass*masses_.wwidth))/ (sqr(s1-m_mW2)+m_mW2*m_GW2); corrfactor *= (sqr(s2-sqr(masses_.wmass))+ sqr(masses_.wmass*masses_.wwidth))/ (sqr(s2-m_mW2)+m_mW2*m_GW2); } else { corrfactor *= (sqr(s1-sqr(masses_.zmass))+ sqr(masses_.zmass*masses_.zwidth))/ (sqr(s1-m_mZ2)+m_mZ2*m_GZ2); corrfactor *= (sqr(s2-sqr(masses_.zmass))+ sqr(masses_.zmass*masses_.zwidth))/ (sqr(s2-m_mZ2)+m_mZ2*m_GZ2); } } for (int n(0);n<2;++n) GetMom(p_p,n,-p[n]); for (size_t n(2);n<p.size();++n) GetMom(p_p,n,p[n]); long int i(m_flavs[0]), j(m_flavs[1]); if (i==21) { i=0; corrfactor *= 8./3.; } if (j==21) { j=0; corrfactor *= 8./3.; } scale_.musq=m_mur2; scale_.scale=sqrt(scale_.musq); epinv_.epinv=epinv2_.epinv2=0.0; double res(CallMCFM(i,j) * corrfactor); epinv_.epinv=1.0; double res1(CallMCFM(i,j) * corrfactor); epinv2_.epinv2=1.0; double res2(CallMCFM(i,j) * corrfactor); m_res.Finite() = res; m_res.IR() = (res1-res); m_res.IR2() = (res2-res1); //msg_Out()<<METHOD<<" yields "<<m_res.Finite() // <<" + 1/eps * "<<m_res.IR() // <<" + 1/eps^2 * "<<m_res.IR2() // <<" for mb = "<<masses_.mb // <<" and "<<nflav_.nflav<<" active flavours" // <<" in "<<scheme_.scheme<<" ... .\n"; } double MCFM_gg_hg::Eps_Scheme_Factor(const Vec4D_Vector& mom) { return 4.*M_PI; } extern "C" { void chooser_(); } DECLARE_VIRTUALME2_GETTER(MCFM_gg_hg,"MCFM_gg_hg") Virtual_ME2_Base *ATOOLS::Getter <Virtual_ME2_Base,Process_Info,MCFM_gg_hg>:: operator()(const Process_Info &pi) const { return NULL; if (pi.m_loopgenerator!="MCFM") return NULL; else msg_Out()<<".\n"; if (MODEL::s_model->Name()!=std::string("SMEHC") || Flavour(kf_b).Yuk()>0. || !Flavour(kf_h0).IsOn()) return NULL; if (!(pi.m_fi.m_nlotype&nlo_type::loop)) return NULL; if (pi.m_fi.m_nlocpl[1]!=0.) return NULL; Flavour_Vector fl(pi.ExtractFlavours()); if (!(fl[0].Strong() && fl[1].Strong())) return NULL; if (pi.m_fi.m_ps.size()!=2) return NULL; Flavour flh(pi.m_fi.m_ps[0].m_fl[0]); if (flh != Flavour(kf_h0)) return NULL; msg_Out()<<"Check numbers: " <<fl.size()<<" external particles, " <<pi.m_fi.m_ps.size()<<" props"; if (pi.m_fi.m_ps.size()==0) msg_Out()<<".\n"; else { msg_Out()<<":\n"; for (size_t i=0;i<pi.m_fi.m_ps.size();i++) { msg_Out()<<" "<<pi.m_fi.m_ps[i].m_ps.size() <<" final state particles for propagator "<<i<<" " <<"with flavour "<<pi.m_fi.m_ps[i].m_fl[0]<<".\n"; for (size_t j=0;j<pi.m_fi.m_ps[i].m_ps.size();j++) { msg_Out()<<" "<<pi.m_fi.m_ps[i].m_ps[j].m_ps.size() <<" final state particles for propagator "<<i<<"["<<j<<"] " <<"with flavour "<<pi.m_fi.m_ps[i].m_ps[j].m_fl[0]<<".\n"; } } } int pID(0); ////////////////////////////////////////////////////////////////// // tau tau final states ////////////////////////////////////////////////////////////////// if (fl.size()==5 && pi.m_fi.m_ps.size()==2 && fl[4].Strong() && fl[2]==fl[3].Bar() && fl[2].Kfcode()==15) { if (Flavour(kf_tau).Yuk()<=0.) { msg_Error()<<"Error in "<<METHOD<<":"<<std::endl <<" Setup for gg->[h->tau tau] (+jet), but tau Yukawa = 0." <<std::endl; THROW(fatal_error,"Inconsistent setup."); } pID = 204; } ////////////////////////////////////////////////////////////////// // photon photon final states ////////////////////////////////////////////////////////////////// else if (fl.size()==5 && pi.m_fi.m_ps.size()==2 && fl[4].Strong() && fl[2]==fl[3] && fl[2].Kfcode()==22) { msg_Error()<<"Error in "<<METHOD<<":"<<std::endl <<" gg->[h->gamma gamma] not yet debugged.\n"; THROW(fatal_error,"Not yet working."); if (fl.size()==4 && pi.m_fi.m_ps.size()==1) pID = 210; } ////////////////////////////////////////////////////////////////// // Z Z final states ////////////////////////////////////////////////////////////////// else if (fl.size()==7 && pi.m_fi.m_ps.size()==2 && fl[6].Strong() && pi.m_fi.m_ps[0].m_ps[0].m_fl[0].Kfcode()==23 && pi.m_fi.m_ps[0].m_ps[0].m_ps.size()==2 && pi.m_fi.m_ps[0].m_ps[1].m_fl[0].Kfcode()==23 && pi.m_fi.m_ps[0].m_ps[1].m_ps.size()==2) { Flavour Z11,Z12,Z21,Z22; Z11 = pi.m_fi.m_ps[0].m_ps[0].m_ps[0].m_fl[0]; Z12 = pi.m_fi.m_ps[0].m_ps[0].m_ps[1].m_fl[0]; Z21 = pi.m_fi.m_ps[0].m_ps[1].m_ps[0].m_fl[0]; Z22 = pi.m_fi.m_ps[0].m_ps[1].m_ps[1].m_fl[0]; msg_Out()<<" test for ZZ: \n" <<" {"<<Z11<<", "<<Z12<<"} {"<<Z21<<", "<<Z22<<"}.\n"; if (!(Z11==Z12.Bar() && Z21==Z22.Bar())) { msg_Error()<<"Error in "<<METHOD<<":"<<std::endl <<" wrong flavour pairs: " <<Z11<<" & "<<Z12<<", "<<Z21<<" & "<<Z22<<".\n"; THROW(fatal_error,"wrong process."); return NULL; } if (Z11.Kfcode()==15 || Z21.Kfcode()==15) { msg_Error()<<"Error in "<<METHOD<<":"<<std::endl <<" process should not involve taus: " <<Z11<<" & "<<Z12<<", "<<Z21<<" & "<<Z22<<".\n"; THROW(fatal_error,"wrong process."); return NULL; } if (Z11.Charge()!=0 && Z21.Charge()!=0) pID = 209; else if (Z11.Charge()!=0 && Z21.Charge()==0) { msg_Error()<<"Error in "<<METHOD<<":"<<std::endl <<" cannot deal with Z->nu nu yet.\n" <<" Must rescale couplings, will be done soon.\n"; THROW(fatal_error,"wrong process."); } } ////////////////////////////////////////////////////////////////// // W W final states ////////////////////////////////////////////////////////////////// else if (fl.size()==7 && pi.m_fi.m_ps.size()==2 && fl[6].Strong() && pi.m_fi.m_ps[0].m_ps[0].m_fl[0].Kfcode()==24 && pi.m_fi.m_ps[0].m_ps[0].m_ps.size()==2 && pi.m_fi.m_ps[0].m_ps[1].m_fl[0].Kfcode()==24 && pi.m_fi.m_ps[0].m_ps[1].m_ps.size()==2) { Flavour W11,W12,W21,W22; W11 = pi.m_fi.m_ps[0].m_ps[0].m_ps[0].m_fl[0]; W12 = pi.m_fi.m_ps[0].m_ps[0].m_ps[1].m_fl[0]; W21 = pi.m_fi.m_ps[0].m_ps[1].m_ps[0].m_fl[0]; W22 = pi.m_fi.m_ps[0].m_ps[1].m_ps[1].m_fl[0]; msg_Out()<<" test for WW: \n" <<" {"<<W11<<", "<<W12<<"} {"<<W21<<", "<<W22<<"}.\n"; if (!W11.IsAnti() && W12.IsAnti() && !W21.IsAnti() && W22.IsAnti() && ((W11.Kfcode()==11 && W12.Kfcode()==12) || (W11.Kfcode()==13 && W12.Kfcode()==14)) && ((W21.Kfcode()==12 && W22.Kfcode()==11) || (W21.Kfcode()==14 && W22.Kfcode()==13))) pID = 208; } ////////////////////////////////////////////////////////////////// // check for 5 light flavours, this is hard-coded in MCFM ////////////////////////////////////////////////////////////////// Flavour lq(kf_quark); double nlf(0.); for (size_t i(0); i<lq.Size(); ++i) if (!lq[i].IsMassive()) nlf++; nlf/=2.; if (nlf!=5.) { msg_Error()<<METHOD<<"(): MCFM only available for five light quark " <<"flavours for this process"<<std::endl; THROW(fatal_error, "Five-flavour calculations with MCFM only.") } if (pID>0) { zerowidth_.zerowidth=true; if (nproc_.nproc>=0) { if (nproc_.nproc!=pID) THROW(not_implemented, "Only one process class allowed when using MCFM"); } nproc_.nproc=pID; chooser_(); msg_Info()<<"Initialise MCFM with nproc = "<<nproc_.nproc<<"\n"; return new MCFM_gg_hg(pID,pi,fl); } return NULL; }