diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index 019f8e1..4c0b62f 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,853 +1,855 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/MatrixElement.hh" #include #include #include #include #include #include #include #include "CLHEP/Vector/LorentzVector.h" #include "fastjet/ClusterSequence.hh" #include "HEJ/Constants.hh" #include "HEJ/currents.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" #include "HEJ/Particle.hh" #include "HEJ/utility.hh" namespace HEJ{ //cf. last line of eq. (22) in \ref Andersen:2011hs double MatrixElement::omega0( double alpha_s, double mur, fastjet::PseudoJet const & q_j, double lambda ) const { const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; // use alpha_s(sqrt(q_j*lambda)), evolved to mur return ( 1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda)) )*result; } Weights MatrixElement::operator()( Event const & event ) const { return tree(event)*virtual_corrections(event); } Weights MatrixElement::tree( Event const & event ) const { return tree_param(event)*tree_kin(event); } Weights MatrixElement::tree_param( Event const & event ) const { if(! is_HEJ(event.type())) { return Weights{0., std::vector(event.variations().size(), 0.)}; } Weights result; // only compute once for each renormalisation scale std::unordered_map known; result.central = tree_param(event, event.central().mur); known.emplace(event.central().mur, result.central); for(auto const & var: event.variations()) { const auto ME_it = known.find(var.mur); if(ME_it == end(known)) { const double wt = tree_param(event, var.mur); result.variations.emplace_back(wt); known.emplace(var.mur, wt); } else { result.variations.emplace_back(ME_it->second); } } return result; } Weights MatrixElement::virtual_corrections( Event const & event ) const { if(! is_HEJ(event.type())) { return Weights{0., std::vector(event.variations().size(), 0.)}; } Weights result; // only compute once for each renormalisation scale std::unordered_map known; result.central = virtual_corrections(event, event.central().mur); known.emplace(event.central().mur, result.central); for(auto const & var: event.variations()) { const auto ME_it = known.find(var.mur); if(ME_it == end(known)) { const double wt = virtual_corrections(event, var.mur); result.variations.emplace_back(wt); known.emplace(var.mur, wt); } else { result.variations.emplace_back(ME_it->second); } } return result; } double MatrixElement::virtual_corrections( Event const & event, double mur ) const{ auto const & in = event.incoming(); auto const & out = event.outgoing(); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif assert(std::is_sorted(out.begin(), out.end(), rapidity_less{})); assert(out.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - out[0].p; size_t first_idx = 0; size_t last_idx = out.size() - 1; // if there is a Higgs or unordered gluon outside the extremal partons // then it is not part of the FKL ladder and does not contribute // to the virtual corrections if(out.front().type == pid::Higgs || event.type() == event_type::unob){ q -= out[1].p; ++first_idx; } if(out.back().type == pid::Higgs || event.type() == event_type::unof){ --last_idx; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q, CLAMBDA)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } assert( nearby(q, -1*pb, norm) || out.back().type == pid::Higgs || event.type() == event_type::unof ); return exp(exponent); } } // namespace HEJ namespace { //! Lipatov vertex for partons emitted into extremal jets double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2)) - p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2)); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>HEJ::CLAMBDA) return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2())); return Cls-4./(kperp*kperp); } } //! Lipatov vertex double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip, CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2. - qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2. + ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom)) + pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom)) - pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.; return -CL.dot(CL); } //! Lipatov vertex with soft subtraction double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>HEJ::CLAMBDA) return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2())); double temp=Cls-4./(kperp*kperp); return temp; } } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ if (aptype==21&&bptype==21) { return jM2gg(pn,pb,p1,pa); } else if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2qg(pn,pb,p1,pa); else return jM2qbarg(pn,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jM2qg(p1,pa,pn,pb); else return jM2qbarg(p1,pa,pn,pb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2qQ(pn,pb,p1,pa); else return jM2qQbar(pn,pb,p1,pa); } else { if (aptype>0) return jM2qQbar(p1,pa,pn,pb); else return jM2qbarQbar(pn,pb,p1,pa); } } throw std::logic_error("unknown particle types"); } /** \brief Matrix element squared for tree-level current-current scattering with Higgs * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared for Tree-Level Current-Current Scattering with Higgs */ double ME_Higgs_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype==21) // gg initial state return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else if (aptype==21&&bptype!=21) { if (bptype > 0) return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; else return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; } else if (bptype==21&&aptype!=21) { if (aptype > 0) return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; else return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; } else { // they are both quark if (bptype>0) { if (aptype>0) return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } else { if (aptype>0) return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered forward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param punof Unordered Particle Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered forward emission */ double ME_Higgs_current_unof( int aptype, int bptype, CLHEP::HepLorentzVector const & punof, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (aptype>0) return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered backward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param punob Unordered back Particle Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered backward emission */ double ME_Higgs_current_unob( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & punob, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (bptype==21&&aptype!=21) { if (aptype > 0) return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (aptype>0) { if (bptype>0) return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (bptype>0) return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){ return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()}; } void validate(HEJ::MatrixElementConfig const & config) { #ifndef HEJ_BUILD_WITH_QCDLOOP if(!config.Higgs_coupling.use_impact_factors) { throw std::invalid_argument{ "Invalid Higgs coupling settings.\n" "HEJ without QCDloop support can only use impact factors.\n" "Set use_impact_factors to true or recompile HEJ.\n" }; } #endif if(config.Higgs_coupling.use_impact_factors && config.Higgs_coupling.mt != std::numeric_limits::infinity()) { throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } } // namespace anonymous namespace HEJ{ MatrixElement::MatrixElement( std::function alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } double MatrixElement::tree_kin( Event const & ev ) const { if(! is_HEJ(ev.type())) return 0.; auto AWZH_boson = std::find_if( begin(ev.outgoing()), end(ev.outgoing()), [](Particle const & p){return is_AWZH_boson(p);} ); if(AWZH_boson == end(ev.outgoing())){ return tree_kin_jets(ev); } switch(AWZH_boson->type){ case pid::Higgs: { return tree_kin_Higgs(ev); } // TODO case pid::Wp: case pid::Wm: case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } namespace{ constexpr int extremal_jet_idx = 1; constexpr int no_extremal_jet_idx = 0; bool treat_as_extremal(Particle const & parton){ return parton.p.user_index() == extremal_jet_idx; } template double FKL_ladder_weight( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn ){ double wt = 1; auto qi = q0; for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){ assert(gluon_it->type == pid::gluon); const auto g = to_HepLorentzVector(*gluon_it); const auto qip1 = qi - g; if(treat_as_extremal(*gluon_it)){ wt *= C2Lipatovots(qip1, qi, pa, pb)*C_A; } else{ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A; } qi = qip1; } return wt; } } // namespace anonymous std::vector MatrixElement::tag_extremal_jet_partons( Event const & ev ) const{ auto out_partons = filter_partons(ev.outgoing()); if(out_partons.size() == ev.jets().size()){ // no additional emissions in extremal jets, don't need to tag anything for(auto & parton: out_partons){ parton.p.set_user_index(no_extremal_jet_idx); } return out_partons; } // TODO: avoid reclustering fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def()); const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt())); assert(jets.size() >= 2); auto most_backward = begin(jets); auto most_forward = end(jets) - 1; // skip jets caused by unordered emission if(ev.type() == event_type::unob){ assert(jets.size() >= 3); ++most_backward; } else if(ev.type() == event_type::unof){ assert(jets.size() >= 3); --most_forward; } const auto extremal_jet_indices = cs.particle_jet_indices( {*most_backward, *most_forward} ); assert(extremal_jet_indices.size() == out_partons.size()); for(size_t i = 0; i < out_partons.size(); ++i){ assert(HEJ::is_parton(out_partons[i])); const int idx = (extremal_jet_indices[i]>=0)? extremal_jet_idx: no_extremal_jet_idx; out_partons[i].p.set_user_index(idx); } return out_partons; } double MatrixElement::tree_kin_jets( Event const & ev ) const { auto const & incoming = ev.incoming(); const auto partons = tag_extremal_jet_partons(ev); if(is_uno(ev.type())){ throw not_implemented("unordered emission not implemented for pure jets"); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); return ME_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa )/(4*(N_C*N_C - 1))*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, pa - p1, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs( Event const & ev ) const { if(is_uno(ev.type())){ return tree_kin_Higgs_between(ev); } if(ev.outgoing().front().type == pid::Higgs){ return tree_kin_Higgs_first(ev); } if(ev.outgoing().back().type == pid::Higgs){ return tree_kin_Higgs_last(ev); } return tree_kin_Higgs_between(ev); } namespace { // Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113 +#ifdef HEJ_BUILD_WITH_QCDLOOP // TODO: code duplication with currents.cc double K_g(double p1minus, double paminus) { return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; } double K_g( CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(type == ParticleID::gluon) return K_g(pout, pin); return C_F; } +#endif // Colour factor in strict MRK limit double K_MRK(ParticleID type) { return (type == ParticleID::gluon)?C_A:C_F; } } double MatrixElement::MH2_forwardH( CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, ParticleID type2, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double t1, double t2 ) const{ ignore(p2out, p2in); const double shat = p1in.invariantMass2(p2in); // gluon case #ifdef HEJ_BUILD_WITH_QCDLOOP if(!param_.Higgs_coupling.use_impact_factors){ return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH( p1out, p1in, p2out, p2in, pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb )/(4*(N_C*N_C - 1)); } #endif return K_MRK(type2)/C_A*9./2.*shat*shat*( C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH) )/(t1*t2); } double MatrixElement::tree_kin_Higgs_first( Event const & ev ) const { auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); assert(outgoing.front().type == pid::Higgs); if(outgoing[1].type != pid::gluon) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(ev); } const auto pH = to_HepLorentzVector(outgoing.front()); const auto partons = tag_extremal_jet_partons( ev ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); const auto q0 = pa - p1 - pH; const double t1 = q0.m2(); const double t2 = (pn - pb).m2(); return MH2_forwardH( p1, pa, incoming[1].type, pn, pb, pH, t1, t2 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_last( Event const & ev ) const { auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); assert(outgoing.back().type == pid::Higgs); if(outgoing[outgoing.size()-2].type != pid::gluon) { assert(incoming.back().type == outgoing[outgoing.size()-2].type); return tree_kin_Higgs_between(ev); } const auto pH = to_HepLorentzVector(outgoing.back()); const auto partons = tag_extremal_jet_partons( ev ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; const double t1 = q0.m2(); const double t2 = (pn + pH - pb).m2(); return MH2_forwardH( pn, pb, incoming[0].type, p1, pa, pH, t2, t1 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_between( Event const & ev ) const { using namespace event_type; auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); const auto the_Higgs = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return s.type == pid::Higgs; } ); assert(the_Higgs != end(outgoing)); const auto pH = to_HepLorentzVector(*the_Higgs); const auto partons = tag_extremal_jet_partons(ev); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector( partons[(ev.type() == unob)?1:0] ); auto pn = to_HepLorentzVector( partons[partons.size() - ((ev.type() == unof)?2:1)] ); auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing)); assert( (first_after_Higgs == end(partons) && ( (ev.type() == unob) || partons.back().type != pid::gluon )) || first_after_Higgs->rapidity() >= the_Higgs->rapidity() ); assert( (first_after_Higgs == begin(partons) && ( (ev.type() == unof) || partons.front().type != pid::gluon )) || (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity() ); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_Higgs == begin(partons)) ++first_after_Higgs; else if(first_after_Higgs == end(partons)) --first_after_Higgs; // t-channel momentum before Higgs auto qH = pa; for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){ qH -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if(ev.type() == unob){ current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(ev.type() == unof){ current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ current_factor = ME_Higgs_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, q0, pa, pb, p1, pn )*FKL_ladder_weight( first_after_Higgs, end_ladder, qH - pH, pa, pb, p1, pn ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double MatrixElement::tree_param_partons( double alpha_s, double mur, std::vector const & partons ) const{ const double gs2 = 4.*M_PI*alpha_s; double wt = std::pow(gs2, partons.size()); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(partons.size() >= 2); for(size_t i = 1; i < partons.size()-1; ++i){ wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp()); } } return wt; } double MatrixElement::tree_param( Event const & ev, double mur ) const{ assert(is_HEJ(ev.type())); const auto & out = ev.outgoing(); const double alpha_s = alpha_s_(mur); auto AWZH_boson = std::find_if( begin(out), end(out), [](auto const & p){return is_AWZH_boson(p);} ); double AWZH_coupling = 1.; if(AWZH_boson != end(out)){ switch(AWZH_boson->type){ case pid::Higgs: { AWZH_coupling = alpha_s*alpha_s; break; } // TODO case pid::Wp: case pid::Wm: case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } if(ev.type() == event_type::unob){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(out) + 1, end(out)}) ); } if(ev.type() == event_type::unof){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(out), end(out) - 1}) ); } return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(out)); } } // namespace HEJ