diff --git a/current_generator/jV_jqqbar.frm b/current_generator/jV_jqqbar.frm new file mode 100644 index 0000000..6adf26b --- /dev/null +++ b/current_generator/jV_jqqbar.frm @@ -0,0 +1,71 @@ +*/** +* \brief Contraction of vector boson current with extremal qqbar current +* +* \authors The HEJ collaboration (see AUTHORS for details) +* \date 2020 +* \copyright GPLv2 or later +*/ +#include- include/helspin.frm +#include- include/write.frm + +* Obtained from jV_juno.frm using crossing symmetry (see eq:crossing in developer manual) + +v pqbar,pq,pa,p1,pb,pl,plbar,q2,q3,p,pr; +i mu,nu,sigma; +s h; + +#do h2={+,-} + l [U1X_jV `h2'] = ( + + i_*Current(`h2'1, pq, nu, pb)*i_*i_*Current(`h2'1, pb, mu, pqbar) + + 2*pq(nu)*i_*Current(`h2'1, pq, mu, pqbar) + )/m2(pq - pb); + l [U2X_jV `h2'] = ( + - 2*i_*Current(`h2'1, pq, mu, pqbar)*pqbar(nu) + - i_*Current(`h2'1, pq, mu, pb)*i_*i_*Current(`h2'1, pb, nu, pqbar) + )/m2(-pqbar + pb); + l [LX_jV `h2'] = ( + - (q2(nu) + q3(nu))*i_*Current(`h2'1, pq, mu, pqbar) + + 2*pb(mu)*i_*Current(`h2'1, pq, nu, pqbar) + - 2*i_*Current(`h2'1, pq, pb, pqbar)*d_(mu, nu) + + m2(q2)*i_*Current(`h2'1, pq, mu, pqbar)*(p1(nu)/m2(p1-pb) + pa(nu)/m2(pa-pb)) + )/m2(q3); +#enddo +.sort +drop; + +* multiply with polarisation vector and other currents +#do h1={+,-} + #do hl={+,-} + #do h2={+,-} + #do hg={+,-} + #do TENSOR={U1X_jV,U2X_jV,LX_jV} + l [`TENSOR' `h1'`hl'`h2'`hg'] = ( + [`TENSOR' `h2'] + *JV(`h1'1, `hl'1, mu, pa, p1, plbar, pl) + *Eps(`hg'1, nu, pr) + ); + #enddo + #enddo + #enddo + #enddo +#enddo + +* choice of auxiliary vector +id Eps(h?, mu?, pr)*Current(h?, ?a) = Eps(h, mu, pb, pq)*Current(h, ?a); +also Eps(h?, mu?, pr) = Eps(h, mu, pb, pqbar); + +multiply replace_( + q2,pq-pb+pqbar, + q3,pq+pqbar +); +.sort +#call ContractCurrents +.sort +format O4; +format c; +#call WriteHeader(`OUTPUT') +#call WriteOptimised(`OUTPUT',U1X_jV,4,pa,pb,p1,pq,pqbar,pl,plbar) +#call WriteOptimised(`OUTPUT',U2X_jV,4,pa,pb,p1,pq,pqbar,pl,plbar) +#call WriteOptimised(`OUTPUT',LX_jV,4,pa,pb,p1,pq,pqbar,pl,plbar) +#call WriteFooter(`OUTPUT') +.end diff --git a/current_generator/jW_jqqbar.frm b/current_generator/jW_jqqbar.frm deleted file mode 100644 index 7f6a9a3..0000000 --- a/current_generator/jW_jqqbar.frm +++ /dev/null @@ -1,48 +0,0 @@ -*/** -* \brief Contraction of W current with extremal qqbar emission current -* -* TODO: transcribe formulas to developer manual -* -* \authors The HEJ collaboration (see AUTHORS for details) -* \date 2020 -* \copyright GPLv2 or later -*/ -#include- include/helspin.frm -#include- include/write.frm - -v pb,p2,p3,pa,p1,pl,plbar,pr,p; -i mu,nu,sigma,alpha; -s h; - -* Equation 3.23 in James Cockburn's Thesis. -l qggm1 = Current(-1, p2, mu, p3-pb, nu, p3); -l qggm2 = -Current(-1, p2, nu, p2-pb, mu, p3); -l qggm3 = -2*i_/m2(p2+p3)*( - (p2(nu) + p3(nu))*d_(mu, sigma) - - pb(mu)*d_(nu, sigma) - + pb(sigma)*d_(mu, nu) -)*Current(-1, p2, sigma, p3); -.sort -drop; - -#do i=1,3 - #do h1={+,-} - #do hg={+,-} - l [jW_qggm`i' `h1'`hg'] = qggm`i'*Eps(`hg'1, nu, pa, pb)*JV( - `h1'1, -1, mu, pa, p1, plbar, pl - ); - #enddo - #enddo -#enddo -id Current(h?, p1?, mu?, p?, nu?, p2?) = Current(h, p1, mu, p, nu, p2)/m2(p); - -#call ContractCurrents -.sort -format O4; -format c; -#call WriteHeader(`OUTPUT') -#call WriteOptimised(`OUTPUT',jW_qggm1,2,pb,p2,p3,pa,p1,pl,plbar) -#call WriteOptimised(`OUTPUT',jW_qggm2,2,pb,p2,p3,pa,p1,pl,plbar) -#call WriteOptimised(`OUTPUT',jW_qggm3,2,pb,p2,p3,pa,p1,pl,plbar) -#call WriteFooter(`OUTPUT') -.end diff --git a/include/HEJ/Wjets.hh b/include/HEJ/Wjets.hh index 41e7218..cab1cc8 100644 --- a/include/HEJ/Wjets.hh +++ b/include/HEJ/Wjets.hh @@ -1,182 +1,182 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ /** \file * \brief Functions computing the square of current contractions in W+Jets. * * This file contains all the W+Jet specific components to compute * the current contractions for valid HEJ processes. */ #pragma once #include #include "CLHEP/Vector/LorentzVector.h" namespace HEJ { struct ParticleProperties; namespace currents { using HLV = CLHEP::HepLorentzVector; //! Square of qQ->qenuQ W+Jets Scattering Current /** * @param p1out Momentum of final state quark * @param plbar Momentum of final state anti-lepton * @param pl Momentum of final state lepton * @param p1in Momentum of initial state quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param wprop Mass and width of the W boson * @returns Square of the current contractions for qQ->qenuQ * * This returns the square of the current contractions in qQ->qenuQ scattering * with an emission of a W Boson off the quark q. */ double ME_W_qQ(HLV const & p1out, HLV const & plbar, HLV const & pl, HLV const & p1in, HLV const & p2out, HLV const & p2in, ParticleProperties const & wprop); - //! qQg Wjets Unordered backwards opposite leg to W + //! W+uno same leg /** * @param p1out Momentum of final state quark a * @param p1in Momentum of initial state quark a * @param p2out Momentum of final state quark b * @param p2in Momentum of intial state quark b * @param pg Momentum of final state unordered gluon * @param plbar Momentum of final state anti-lepton * @param pl Momentum of final state lepton * @param wprop Mass and width of the W boson - * @returns Square of the current contractions for qQ->qQg Scattering + * @returns Square of the current contractions for qQ->qQg * - * This returns the square of the current contractions in qQg->qQg scattering + * This returns the square of the current contractions in gqQ->gqQ scattering * with an emission of a W Boson. */ - double ME_W_unob_qQ(HLV const & p1out, HLV const & p1in, - HLV const & p2out, HLV const & p2in, HLV const & pg, - HLV const & plbar, HLV const & pl, - ParticleProperties const & wprop); + double ME_Wuno_qQ(HLV const & p1out, HLV const & p1in, + HLV const & p2out, HLV const & p2in, HLV const & pg, + HLV const & plbar, HLV const & pl, + ParticleProperties const & wprop); - //! W+uno same leg + //! qQg Wjets Unordered backwards opposite leg to W /** * @param p1out Momentum of final state quark a * @param p1in Momentum of initial state quark a * @param p2out Momentum of final state quark b * @param p2in Momentum of intial state quark b * @param pg Momentum of final state unordered gluon * @param plbar Momentum of final state anti-lepton * @param pl Momentum of final state lepton * @param wprop Mass and width of the W boson - * @returns Square of the current contractions for qQ->qQg + * @returns Square of the current contractions for qQ->qQg Scattering * - * This returns the square of the current contractions in gqQ->gqQ scattering + * This returns the square of the current contractions in qQg->qQg scattering * with an emission of a W Boson. */ - double ME_Wuno_qQ(HLV const & p1out, HLV const & p1in, - HLV const & p2out, HLV const & p2in, HLV const & pg, - HLV const & plbar, HLV const & pl, - ParticleProperties const & wprop); + double ME_W_unob_qQ(HLV const & p1out, HLV const & p1in, + HLV const & p2out, HLV const & p2in, HLV const & pg, + HLV const & plbar, HLV const & pl, + ParticleProperties const & wprop); //! W+Extremal qqbar. qqbar+Q /** * @param pgin Momentum of initial state gluon * @param pqbarout Momentum of final state anti-quark a * @param plbar Momentum of final state anti-lepton * @param pl Momentum of final state lepton * @param pqout Momentum of final state quark a * @param p2out Momentum of initial state anti-quark b * @param p2in Momentum of final state gluon b * @param wprop Mass and width of the W boson * @returns Square of the current contractions for ->qqbarQ * * Calculates the square of the current contractions with extremal qqbar pair * production. This is calculated through the use of crossing symmetry. */ double ME_WExqqbar_qqbarQ(HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl, HLV const & pqout, HLV const & p2out, HLV const & p2in, ParticleProperties const & wprop); //! W+Extremal qqbar. gg->qqbarg. qqbar on forwards leg, W emission backwards leg. /** * @param pa Momentum of initial state (anti-)quark * @param pb Momentum of initial state gluon * @param p1 Momentum of final state (anti-)quark (after W emission) - * @param p2 Momentum of final state anti-quark - * @param p3 Momentum of final state quark + * @param pq Momentum of final state quark + * @param pqbar Momentum of final state anti-quark * @param plbar Momentum of final state anti-lepton * @param pl Momentum of final state lepton * @param aqlinepa Is opposite extremal leg to qqbar a quark or antiquark line * @param wprop Mass and width of the W boson * @returns Square of the current contractions for gq->qqbarqW * * Calculates the square of the current contractions with extremal qqbar pair * production. This is calculated via current contraction of existing * currents. Assumes qqbar split from forwards leg, W emission from backwards * leg. Switch input (pa<->pb, p1<->pn) if calculating forwards qqbar. */ double ME_W_Exqqbar_QQq(HLV const & pa, HLV const & pb, HLV const & p1, HLV const & p2, HLV const & p3, HLV const & plbar, HLV const & pl, bool aqlinepa, ParticleProperties const & wprop); //! W+Jets qqbarCentral. qqbar W emission. /** * @param pa Momentum of initial state particle a * @param pb Momentum of initial state particle b * @param pl Momentum of final state lepton * @param plbar Momentum of final state anti-lepton * @param partons Vector of outgoing parton momenta * @param aqlinepa True= pa is anti-quark * @param aqlinepb True= pb is anti-quark * @param qqbar_marker Ordering of the qqbar pair produced (q-qbar vs qbar-q) * @param nabove Number of lipatov vertices "above" qqbar pair * @param wprop Mass and width of the W boson * @returns Square of the current contractions for qq>qQQbarWq * * Calculates the square of the current contractions with extremal qqbar pair * production. This is calculated through the use of crossing symmetry. */ double ME_WCenqqbar_qq(HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar, std::vector const & partons, bool aqlinepa, bool aqlinepb, bool qqbar_marker, int nabove, ParticleProperties const & wprop); //! W+Jets qqbarCentral. W emission from backwards leg. /** * @param pa Momentum of initial state particle a * @param pb Momentum of initial state particle b * @param pl Momentum of final state lepton * @param plbar Momentum of final state anti-lepton * @param partons outgoing parton momenta * @param aqlinepa True= pa is anti-quark * @param aqlinepb True= pb is anti-quark * @param qqbar_marker Ordering of the qqbar pair produced (q-qbar vs qbar-q) * @param nabove Number of lipatov vertices "above" qqbar pair * @param nbelow Number of lipatov vertices "below" qqbar pair * @param forwards Swap to emission off front leg * TODO: remove so args can be const * @param wprop Mass and width of the W boson * @returns Square of the current contractions for qq>qQQbarWq * * Calculates the square of the current contractions with extremal qqbar pair * production. This is calculated through the use of crossing symmetry. */ double ME_W_Cenqqbar_qq(HLV pa, HLV pb, HLV pl, HLV plbar, std::vector partons, bool aqlinepa, bool aqlinepb, bool qqbar_marker, int nabove, int nbelow, bool forwards, ParticleProperties const & wprop); } // namespace currents } // namespace HEJ diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index 6d7fcdb..c9d0f0f 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,2475 +1,2475 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2022 * \copyright GPLv2 or later */ #include "HEJ/MatrixElement.hh" #include #include #include #include #include #include #include #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Event.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/Hjets.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Particle.hh" #include "HEJ/WWjets.hh" #include "HEJ/Wjets.hh" #include "HEJ/Zjets.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/jets.hh" #include "HEJ/utility.hh" namespace HEJ { namespace { // Colour acceleration multiplier for gluons // see eq:K_g in developer manual 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()); } // Colour acceleration multiplier for quarks // see eq:K_q in developer manual constexpr double K_q = C_F; // Colour acceleration multipliers double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ){ if(type == pid::gluon) return K_g(pout, pin); return K_q; } } // namespace double MatrixElement::omega0( double alpha_s, double mur, fastjet::PseudoJet const & q_j ) const { const double lambda = param_.regulator_lambda; const double result = - alpha_s*N_C/M_PI*std::log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; return ( 1. + alpha_s/(4.*M_PI)*BETA0*std::log(mur*mur/(q_j.perp()*lambda)) )*result; } Weights MatrixElement::operator()(Event const & event) const { std::vector tree_kin_part=tree_kin(event); std::vector virtual_part=virtual_corrections(event); if(tree_kin_part.size() != virtual_part.size()) { throw std::logic_error("tree and virtuals have different sizes"); } Weights sum = Weights{0., std::vector(event.variations().size(), 0.)}; for(size_t i=0; i tree_kin_part=tree_kin(event); double sum = 0.; for(double i : tree_kin_part) { sum += i; } return tree_param(event)*sum; } Weights MatrixElement::tree_param(Event const & event) const { if(! is_resummable(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; } std::vector MatrixElement::virtual_corrections(Event const & event) const { if(!event.valid_hej_state(param_.soft_pt_regulator)) { return {Weights{0., std::vector(event.variations().size(), 0.)}}; } // only compute once for each renormalisation scale std::unordered_map > known_vec; std::vector central_vec=virtual_corrections(event, event.central().mur); known_vec.emplace(event.central().mur, central_vec); for(auto const & var: event.variations()) { const auto ME_it = known_vec.find(var.mur); if(ME_it == end(known_vec)) { known_vec.emplace(var.mur, virtual_corrections(event, var.mur)); } } // At this stage known_vec contains one vector of virtual corrections for each mur value // Now put this into a vector of Weights std::vector result_vec; for(size_t i=0; isecond.at(i)); } result_vec.emplace_back(result); } return result_vec; } template std::vector MatrixElement::virtual_corrections_interference( InputIterator begin_parton, InputIterator end_parton, fastjet::PseudoJet const & q0_t, fastjet::PseudoJet const & q0_b, const double mur ) const{ const double alpha_s = alpha_s_(mur); auto qi_t = q0_t; auto qi_b = q0_b; double sum_top = 0.; double sum_bot = 0.; double sum_mix = 0.; for(auto parton_it = begin_parton; parton_it != end_parton; ++parton_it){ Particle parton = *parton_it; Particle parton_next = *(parton_it+1); const double dy = parton_next.rapidity() - parton.rapidity(); const double tmp_top = omega0(alpha_s, mur, qi_t)*dy; const double tmp_bot = omega0(alpha_s, mur, qi_b)*dy; sum_top += tmp_top; sum_bot += tmp_bot; sum_mix += (tmp_top + tmp_bot) / 2.; qi_t -= parton_next.p; qi_b -= parton_next.p; } if (param_.nlo.enabled){ return {(sum_top), (sum_bot), (sum_mix)}; } return {exp(sum_top), exp(sum_bot), exp(sum_mix)}; } double MatrixElement::virtual_corrections_W( Event const & event, const double mur, Particle const & WBoson ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(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(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - partons[0].p; std::size_t first_idx = 0; std::size_t last_idx = partons.size() - 1; #ifndef NDEBUG bool wc = true; #endif bool wqq = false; // With extremal qqbar or unordered gluon outside the extremal // partons then it is not part of the FKL ladder and does not // contribute to the virtual corrections. W emitted from the // most backward leg must be taken into account in t-channel if (event.type() == event_type::unob) { q -= partons[1].p; ++first_idx; if (in[0].type != partons[1].type ){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } else if (event.type() == event_type::qqbar_exb) { q -= partons[1].p; ++first_idx; if (std::abs(partons[0].type) != std::abs(partons[1].type)){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } else { if(event.type() == event_type::unof || event.type() == event_type::qqbar_exf){ --last_idx; } if (in[0].type != partons[0].type ){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } std::size_t first_idx_qqbar = last_idx; std::size_t last_idx_qqbar = last_idx; //if qqbarMid event, virtual correction do not occur between qqbar pair. if(event.type() == event_type::qqbar_mid){ const auto backquark = std::find_if( begin(partons) + 1, end(partons) - 1 , [](Particle const & s){ return (s.type != pid::gluon); } ); if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0; if(std::abs(backquark->type) != std::abs((backquark+1)->type)) { wqq=true; #ifndef NDEBUG wc=false; #endif } last_idx = std::distance(begin(partons), backquark); first_idx_qqbar = last_idx+1; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(std::size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q)*( partons[j+1].rapidity() - partons[j].rapidity() ); q -=partons[j+1].p; } // End Loop one if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p; if (wqq) q -= WBoson.p; for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){ exponent += omega0(alpha_s, mur, q)*( partons[j+1].rapidity() - partons[j].rapidity() ); q -= partons[j+1].p; } #ifndef NDEBUG if (wc) q -= WBoson.p; assert( nearby(q, -1*pb, norm) || is_AWZH_boson(partons.back().type) || event.type() == event_type::unof || event.type() == event_type::qqbar_exf ); #endif if (param_.nlo.enabled){ double nlo_virtual =1.; //Only apply virtual corrections to a nlo order event. if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent; return nlo_virtual; } return std::exp(exponent); } std::vector MatrixElement::virtual_corrections_WW( Event const & event, const double mur ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); assert(event.decays().size() == 2); std::vector plbar; std::vector pl; for(auto const & decay_pair : event.decays()) { auto const decay = decay_pair.second; if(decay.at(0).type < 0) { plbar.emplace_back(decay.at(0).p); pl .emplace_back(decay.at(1).p); } else { pl .emplace_back(decay.at(0).p); plbar.emplace_back(decay.at(1).p); } } fastjet::PseudoJet q_t = pa - partons[0].p - pl[0] - plbar[0]; fastjet::PseudoJet q_b = pa - partons[0].p - pl[1] - plbar[1]; auto const begin_parton = cbegin(partons); auto const end_parton = cend(partons) - 1; if (param_.nlo.enabled){ std::vector virt_corrections_nlo {1.0,1.0,1.0}; //set virtual corrections to 1. // Only apply virtual corrections to a nlo order event. if (partons.size() == param_.nlo.nj) { std::vector virt_corrections_nlo_interference = virtual_corrections_interference( begin_parton, end_parton, q_t, q_b, mur ); assert( virt_corrections_nlo_interference.size() == virt_corrections_nlo.size() ); for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){ virt_corrections_nlo[i] += virt_corrections_nlo_interference[i]; } } return virt_corrections_nlo; } return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur); } std::vector MatrixElement::virtual_corrections_Z_qq( Event const & event, const double mur, Particle const & ZBoson ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p; fastjet::PseudoJet q_b = pa - partons[0].p; auto begin_parton = cbegin(partons); auto end_parton = cend(partons) - 1; // Unordered gluon does not contribute to the virtual corrections if (event.type() == event_type::unob) { // Gluon is partons[0] and is already subtracted // partons[1] is the backward quark q_t -= partons[1].p; q_b -= partons[1].p; ++begin_parton; } else if (event.type() == event_type::unof) { // End sum at forward quark --end_parton; } if (param_.nlo.enabled){ //set virtual corrections to 1. std::vector virt_corrections_nlo {1.0,1.0,1.0}; // Only apply virtual corrections to a nlo order event. if (partons.size() == param_.nlo.nj) { std::vector virt_corrections_nlo_interference = virtual_corrections_interference( begin_parton, end_parton, q_t, q_b, mur ); assert( virt_corrections_nlo.size() == virt_corrections_nlo_interference.size() ); for(std::size_t i = 0; i < virt_corrections_nlo.size(); ++i){ virt_corrections_nlo[i] += virt_corrections_nlo_interference[i]; } } return virt_corrections_nlo; } return virtual_corrections_interference(begin_parton, end_parton, q_t, q_b, mur); } double MatrixElement::virtual_corrections_Z_qg( Event const & event, const double mur, Particle const & ZBoson, const bool is_gq_event ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); // If this is a gq event, don't subtract the Z momentum from first q fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p); size_t first_idx = 0; size_t last_idx = partons.size() - 1; // Unordered gluon does not contribute to the virtual corrections if (event.type() == event_type::unob) { // Gluon is partons[0] and is already subtracted // partons[1] is the backward quark q -= partons[1].p; ++first_idx; } else if (event.type() == event_type::unof) { // End sum at forward quark --last_idx; } double sum=0.; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity() - partons[j].rapidity()); q -= partons[j+1].p; } if (param_.nlo.enabled){ double nlo_virtual =1.; //Only apply virtual corrections to a nlo order event. if (partons.size() == param_.nlo.nj) nlo_virtual+=sum; return nlo_virtual; } return exp(sum); } std::vector MatrixElement::virtual_corrections( Event const & event, const 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 std::vector bosons = filter_AWZH_bosons(out); if(event.jets().size() != param_.nlo.nj && param_.nlo.enabled) { throw std::logic_error{ "Input event has number of jets different to stated NLO " "input in config file: " + std::to_string(event.jets().size()) + " vs " +std::to_string(param_.nlo.nj) + "\n" }; } if(bosons.size() > 2) { throw not_implemented("Emission of >2 bosons is unsupported"); } if(bosons.size() == 2) { if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) { return virtual_corrections_WW(event, mur); } else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) { return virtual_corrections_WW(event, mur); } throw not_implemented("Emission of bosons of unsupported type"); } if(bosons.size() == 1) { const auto AWZH_boson = bosons[0]; if(std::abs(AWZH_boson.type) == pid::Wp){ return {virtual_corrections_W(event, mur, AWZH_boson)}; } if(AWZH_boson.type == pid::Z_photon_mix){ if(is_gluon(in.back().type)){ // This is a qg event return {virtual_corrections_Z_qg(event, mur, AWZH_boson, false)}; } if(is_gluon(in.front().type)){ // This is a gq event return {virtual_corrections_Z_qg(event, mur, AWZH_boson, true)}; } // This is a qq event return virtual_corrections_Z_qq(event, mur, AWZH_boson); } } 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; std::size_t first_idx = 0; std::size_t last_idx = out.size() - 1; // if there is a Higgs boson _not_ emitted off an incoming gluon, // extremal qqbar 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 && in.front().type != pid::gluon) || event.type() == event_type::unob || event.type() == event_type::qqbar_exb){ q -= out[1].p; ++first_idx; } if((out.back().type == pid::Higgs && in.back().type != pid::gluon) || event.type() == event_type::unof || event.type() == event_type::qqbar_exf){ --last_idx; } std::size_t first_idx_qqbar = last_idx; std::size_t last_idx_qqbar = last_idx; //if central qqbar event, virtual correction do not occur between q-qbar. if(event.type() == event_type::qqbar_mid){ const auto backquark = std::find_if( begin(out) + 1, end(out) - 1 , [](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); } ); if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.}; last_idx = std::distance(begin(out), backquark); first_idx_qqbar = last_idx+1; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(std::size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p; for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){ exponent += omega0(alpha_s, mur, q)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } assert( nearby(q, -1*pb, norm) || (out.back().type == pid::Higgs && in.back().type != pid::gluon) || event.type() == event_type::unof || event.type() == event_type::qqbar_exf ); if (param_.nlo.enabled){ const auto partons = filter_partons(event.outgoing()); double nlo_virtual =1.; //Only apply virtual corrections to a nlo order event. if (partons.size() == param_.nlo.nj) nlo_virtual+=exponent; return {nlo_virtual}; } return {std::exp(exponent)}; } namespace { //! Lipatov vertex for partons emitted into extremal jets CLHEP::HepLorentzVector CLipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ) { const CLHEP::HepLorentzVector p5 = qav-qbv; const CLHEP::HepLorentzVector CL = -(qav+qbv) + 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; } double C2Lipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ){ const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets double C2Lipatovots( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2())); const double kperp=(qav-qbv).perp(); if (kperp>lambda) return Cls; return Cls-4./(kperp*kperp); } double C2Lipatov_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ) { const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2); const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2); return -CL_t.dot(CL_b); } double C2Lipatovots_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2) / sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2()); const double kperp = (qav_t - qbv_t).perp(); if (kperp > lambda){ return Cls; } return Cls - 4.0 / (kperp * kperp); } CLHEP::HepLorentzVector CLipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ){ const CLHEP::HepLorentzVector p5 = qav-qbv; const CLHEP::HepLorentzVector CL = -(qav+qbv) + 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; } //! Lipatov vertex double C2Lipatov( // B CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ){ const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction double C2Lipatovots( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2())); const double kperp=(qav-qbv).perp(); if (kperp>lambda) return Cls; return Cls-4./(kperp*kperp); } double C2Lipatov_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ) { const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop); const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop); return -CL_t.dot(CL_b); } double C2Lipatovots_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2) / sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2()); const double kperp = (qav_t - qbv_t).perp(); if (kperp > lambda) { return Cls; } return Cls - 4.0 / (kperp * kperp); } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pg Unordered gluon momentum * @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 * * @note The unof contribution can be calculated by reversing the argument ordering. */ double ME_uno_current( [[maybe_unused]] ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ assert(aptype!=pid::gluon); // aptype cannot be gluon const double t1 = (pa - p1 - pg).m2(); const double t2 = (pb - pn).m2(); return K_q * K(bptype, pn, pb)*currents::ME_unob_qq(pg, p1, pa, pn, pb) / (t1 * t2); } /** Matrix element squared for tree-level current-current scattering * @param bptype Particle b PDG ID * @param pgin Incoming gluon momentum * @param p1 More backward (anti-)quark from splitting Momentum * @param p2 Less backward (anti-)quark from splitting Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @returns ME Squared for Tree-Level Current-Current Scattering * * @note The forward qqbar contribution can be calculated by reversing the * argument ordering. */ double ME_qqbar_current( ParticleID bptype, CLHEP::HepLorentzVector const & pgin, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb ){ const double t1 = (pgin - p1 - p2).m2(); const double t2 = (pn - pb).m2(); return K(bptype, pn, pb)*currents::ME_qqbar_qg(pb, pgin, pn, p2, p1) / (t1 * t2); } /* \brief Matrix element squared for central qqbar tree-level current-current * scattering * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param nabove Number of gluons emitted before central qqbarpair * @param nbelow Number of gluons emitted after central qqbarpair * @param pa Initial state a Momentum * @param pb Initial state b Momentum * @param pq Final state qbar Momentum * @param pqbar Final state q Momentum * @param partons Vector of all outgoing partons * @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering */ double ME_qqbar_mid_current( ParticleID aptype, ParticleID bptype, int nabove, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, std::vector const & partons ){ using namespace currents; // CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards) const bool swap_qqbar=pqbar.rapidity() < pq.rapidity(); CLHEP::HepLorentzVector const & p1 = partons.front(); CLHEP::HepLorentzVector const & pn = partons.back(); const double t1 = (p1 - pa).m2(); const double t2 = (pb - pn).m2(); return K(aptype, p1, pa) *K(bptype, pn, pb) *ME_Cenqqbar_qq( pa, pb, partons, swap_qqbar, nabove ) / (t1 * t2); } /** 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( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ const double t1 = (p1 - pa).m2(); const double t2 = (pb - pn).m2(); return K(aptype, p1, pa) * K(bptype, pn, pb) * currents::ME_qq(p1, pa, pn, pb)/(t1 * t2); } /** Matrix element squared for tree-level current-current scattering With W+Jets * @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 wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_W_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; assert(!(aptype==pid::gluon && bptype==pid::gluon)); if(aptype == pid::gluon || wc) { // swap currents to ensure that the W is emitted off the first one return ME_W_current(bptype, aptype, p1, pa, pn, pb, plbar, pl, false, Wprop); } // we assume that the W is emitted off a quark line // if this is not the case, we have to apply CP conjugation, // which is equivalent to swapping lepton and antilepton momenta const double current_contr = is_quark(aptype)? ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop): ME_W_qQ(p1,pl,plbar,pa,pn,pb,Wprop); const double t1 = (pa - p1 - pl - plbar).m2(); const double tn = (pn - pb).m2(); return K(aptype, p1, pa) * K(bptype, pn, pb) * current_contr/(4.*(N_C*N_C - 1) * t1 * tn); } /** Matrix element squared for backwards uno tree-level current-current * scattering With W+Jets * * @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 pg Unordered gluon momentum * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for unob Tree-Level Current-Current Scattering * * @note The unof contribution can be calculated by reversing the argument ordering. */ double ME_W_uno_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; assert(bptype != pid::gluon || aptype != pid::gluon); if(aptype == pid::gluon || wc) { // emission off pb -- pn line // we assume that the W is emitted off a quark line // if this is not the case, we have to apply CP conjugation, // which is equivalent to swapping lepton and antilepton momenta if(is_antiquark(bptype)) std::swap(plbar, pl); const double t1 = (pa - p1 - pg).m2(); const double tn = (pb - pn - plbar - pl).m2(); return K_q*K_q*ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(4.*(N_C*N_C - 1) * t1 * tn); } // emission off pa -- p1 line if(is_antiquark(aptype)) std::swap(plbar, pl); const double t1 = (pa - p1 - pg - plbar - pl).m2(); const double tn = (pb - pn).m2(); return K(bptype, pn, pb)/C_F*ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop)/(t1 * tn); } /** \brief Matrix element squared for backward qqbar tree-level current-current * scattering With W+Jets * * @param bptype Particle b PDG ID * @param pa Initial state a Momentum * @param pb Initial state b Momentum * @param pq Final state q Momentum * @param pqbar Final state qbar Momentum * @param pn Final state n Momentum * @param plbar Final state anti-lepton momentum * @param pl Final state lepton momentum * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for qqbarb Tree-Level Current-Current Scattering * * @note calculate forwards qqbar contribution by reversing argument ordering. */ double ME_W_qqbar_current( ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, - CLHEP::HepLorentzVector pq, - CLHEP::HepLorentzVector pqbar, + CLHEP::HepLorentzVector const & pq, + CLHEP::HepLorentzVector const & pqbar, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; if(is_anyquark(bptype) && wc) { // W Must be emitted from forwards leg. const double t1 = (pa - pq - pqbar).m2(); const double tn = (pb - pn - pl - plbar).m2(); - return ME_W_Exqqbar_QQq( + return K_q * K_q * ME_W_Exqqbar_QQq( pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop - )/(t1 * tn); + ) / (4.*(N_C*N_C - 1) * t1 * tn); } const double t1 = (pa - pl - plbar - pq - pqbar).m2(); const double tn = (pb - pn).m2(); return K(bptype, pn, pb)/C_F * ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop) / (t1 * tn); } /* \brief Matrix element squared for central qqbar tree-level current-current * scattering With W+Jets * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param nabove Number of gluons emitted before central qqbarpair * @param nbelow Number of gluons emitted after central qqbarpair * @param pa Initial state a Momentum * @param pb Initial state b Momentum\ * @param pq Final state qbar Momentum * @param pqbar Final state q Momentum * @param partons Vector of all outgoing partons * @param plbar Final state anti-lepton momentum * @param pl Final state lepton momentum * @param wqq Boolean. True siginfies W boson is emitted from Central qqbar * @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false. * @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering */ double ME_W_qqbar_mid_current( ParticleID aptype, ParticleID bptype, int nabove, int nbelow, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wqq, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; // CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards) const bool swap_qqbar=pqbar.rapidity() < pq.rapidity(); double wt=1.; if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F; if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F; if(wqq) return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype),is_antiquark(bptype), swap_qqbar, nabove, Wprop); return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype), is_antiquark(bptype), swap_qqbar, nabove, nbelow, wc, Wprop); } /** Matrix element squared for tree-level current-current scattering With Z+Jets * @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 plbar Final state positron momentum * @param pl Final state electron momentum * @param Zprop Z properties * @param stw2 Value of sin(theta_w)^2 * @param ctw Value of cos(theta_w) * @returns ME Squared for Tree-Level Current-Current Scattering */ std::vector ME_Z_current( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, ParticleProperties const & Zprop, const double stw2, const double ctw ){ using namespace currents; const auto pZ = pl + plbar; const double pref = K(aptype, p1, pa) * K(bptype, pn, pb)/(4.*(N_C*N_C-1)); // we know they are not both gluons assert(!is_gluon(aptype) || !is_gluon(bptype)); if(is_anyquark(aptype) && is_gluon(bptype)){ // This is a qg event const double t1 = (pa-p1-pZ).m2(); const double tn = (pb-pn ).m2(); return { pref*ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw)/(t1 * tn) }; } if(is_gluon(aptype) && is_anyquark(bptype)){ // This is a gq event const double t1 = (pa-p1 ).m2(); const double tn = (pb-pn-pZ).m2(); return { pref*ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw)/(t1 * tn) }; } // This is a qq event assert(is_anyquark(aptype) && is_anyquark(bptype)); const double t1_top = (pa-p1-pZ).m2(); const double t2_top = (pb-pn ).m2(); const double t1_bot = (pa-p1 ).m2(); const double t2_bot = (pb-pn-pZ).m2(); std::vector res = ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw); assert(res.size() == 3); res[0] *= pref/(t1_top * t2_top); res[1] *= pref/(t1_bot * t2_bot); res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot); return res; } /** Matrix element squared for backwards uno tree-level current-current * scattering With Z+Jets * * @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 pg Unordered gluon momentum * @param plbar Final state positron momentum * @param pl Final state electron momentum * @param Zprop Z properties * @param stw2 Value of sin(theta_w)^2 * @param ctw Value of cos(theta_w) * @returns ME Squared for unob Tree-Level Current-Current Scattering * * @note The unof contribution can be calculated by reversing the argument ordering. */ std::vector ME_Z_uno_current( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, ParticleProperties const & Zprop, const double stw2, const double ctw ){ using namespace currents; const auto pZ = pl + plbar; const double pref = K(aptype, p1, pa)/C_F * K(bptype, pn, pb)/C_F; // we know they are not both gluons assert(!is_gluon(aptype) || !is_gluon(bptype)); if (is_anyquark(aptype) && is_gluon(bptype)) { // This is a qg event const double t1 = (pa-p1-pg-pZ).m2(); const double tn = (pb-pn ).m2(); return { pref*ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw)/(t1 * tn) }; } if (is_gluon(aptype) && is_anyquark(bptype)) { // This is a gq event const double t1 = (pa-p1 ).m2(); const double tn = (pb-pn-pg-pZ).m2(); return { pref*ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw)/(t1 * tn) }; } // This is a qq event assert(is_anyquark(aptype) && is_anyquark(bptype)); const double t1_top = (pa-pg-p1-pZ).m2(); const double t2_top = (pb-pn ).m2(); const double t1_bot = (pa-pg-p1).m2(); const double t2_bot = (pb-pn-pZ).m2(); std::vector res = ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw); assert(res.size() == 3); res[0] *= pref/(t1_top * t2_top); res[1] *= pref/(t1_bot * t2_bot); res[2] *= pref/sqrt(t1_top * t2_top * t1_bot * t2_bot); return res; } /** \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( ParticleID aptype, ParticleID 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, double vev ){ const double t1 = (pa - p1).m2(); const double tn = (pb - pn).m2(); return K(aptype, p1, pa) *K(bptype, pn, pb) *currents::ME_H_qQ( pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev ) / (t1 * tn * qH.m2() * qHp1.m2()); } /** \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 pg 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 * * @note This function assumes unordered gluon backwards from pa-p1 current. * For unof, reverse call order */ double ME_Higgs_current_uno( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pg, 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, double vev ){ const double t1 = (pa - p1 - pg).m2(); const double tn = (pn - pb).m2(); return K(aptype, p1, pa) *K(bptype, pn, pb) *currents::ME_H_unob_qQ( pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev ) / (t1 * qH.m2() * qHp1.m2() * tn); } /** Matrix element squared for tree-level scattering with WW+Jets * @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 pl1bar Particle l1bar Momentum * @param pl1 Particle l1 Momentum * @param pl2bar Particle l2bar Momentum * @param pl2 Particle l2 Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ std::vector ME_WW_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pl1bar, CLHEP::HepLorentzVector const & pl1, CLHEP::HepLorentzVector const & pl2bar, CLHEP::HepLorentzVector const & pl2, ParticleProperties const & Wprop ){ using namespace currents; if (aptype > 0 && bptype > 0) return ME_WW_qQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop); if (aptype < 0 && bptype > 0) return ME_WW_qbarQ(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop); if (aptype > 0 && bptype < 0) return ME_WW_qQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop); if (aptype < 0 && bptype < 0) return ME_WW_qbarQbar(p1, pl1bar, pl1, pa, pn, pl2bar, pl2, pb, Wprop); throw std::logic_error("unreachable"); } void validate(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 MatrixElement::MatrixElement( std::function alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } std::vector MatrixElement::tree_kin( Event const & ev ) const { if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.}; if(!ev.valid_incoming()){ throw std::invalid_argument{ "Invalid momentum for one or more incoming particles. " "Incoming momenta must have vanishing mass and transverse momentum." }; } std::vector bosons = filter_AWZH_bosons(ev.outgoing()); if(bosons.empty()) { return {tree_kin_jets(ev)}; } if(bosons.size() == 1) { switch(bosons[0].type){ case pid::Higgs: return {tree_kin_Higgs(ev)}; case pid::Wp: case pid::Wm: return {tree_kin_W(ev)}; case pid::Z_photon_mix: return tree_kin_Z(ev); // TODO case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } if(bosons.size() == 2) { if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp){ return tree_kin_WW(ev); } else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm){ return tree_kin_WW(ev); } throw not_implemented("Emission of bosons of unsupported type"); } throw not_implemented("Emission of >2 bosons is unsupported"); } 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 lambda ){ 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, lambda)*C_A; } else{ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A; } qi = qip1; } return wt; } template std::vector FKL_ladder_weight_mix( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn, const double lambda ){ double wt_top = 1; double wt_bot = 1; double wt_mix = 1; auto qi_t = q0_t; auto qi_b = q0_b; 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_t = qi_t - g; const auto qip1_b = qi_b - g; if(treat_as_extremal(*gluon_it)){ wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A; wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A; wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A; } else{ wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A; wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A; wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A; } qi_t = qip1_t; qi_b = qip1_b; } return {wt_top, wt_bot, wt_mix}; } std::vector tag_extremal_jet_partons( Event const & ev ){ 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; } auto const & jets = ev.jets(); std::vector extremal_jets; if(! is_backward_g_to_h(ev)) { auto most_backward = begin(jets); // skip jets caused by unordered emission or qqbar if(ev.type() == event_type::unob || ev.type() == event_type::qqbar_exb){ assert(jets.size() >= 2); ++most_backward; } extremal_jets.emplace_back(*most_backward); } if(! is_forward_g_to_h(ev)) { auto most_forward = end(jets) - 1; if(ev.type() == event_type::unof || ev.type() == event_type::qqbar_exf){ assert(jets.size() >= 2); --most_forward; } extremal_jets.emplace_back(*most_forward); } const auto extremal_jet_indices = ev.particle_jet_indices( extremal_jets ); assert(extremal_jet_indices.size() == out_partons.size()); for(std::size_t i = 0; i < out_partons.size(); ++i){ assert(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 tree_kin_jets_qqbar_mid( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, double lambda ){ CLHEP::HepLorentzVector pq; CLHEP::HepLorentzVector pqbar; const auto backmidquark = std::find_if( begin(partons)+1, end(partons)-1, [](Particle const & s){ return s.type != pid::gluon; } ); assert(backmidquark!=end(partons)-1); if (is_quark(backmidquark->type)){ pq = to_HepLorentzVector(*backmidquark); pqbar = to_HepLorentzVector(*(backmidquark+1)); } else { pqbar = to_HepLorentzVector(*backmidquark); pq = to_HepLorentzVector(*(backmidquark+1)); } auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto q0 = pa - p1; // t-channel momentum after qqbar auto qqbart = q0; const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder_1 = (backmidquark); const auto begin_ladder_2 = (backmidquark+2); const auto end_ladder = cend(partons) - 1; for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){ qqbart -= to_HepLorentzVector(*parton_it); } const int nabove = std::distance(begin_ladder, backmidquark); std::vector partonsHLV; partonsHLV.reserve(partons.size()); for (std::size_t i = 0; i != partons.size(); ++i) { partonsHLV.push_back(to_HepLorentzVector(partons[i])); } const double current_factor = ME_qqbar_mid_current( aptype, bptype, nabove, pa, pb, pq, pqbar, partonsHLV ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder_1, q0, pa, pb, p1, pn, lambda )*FKL_ladder_weight( begin_ladder_2, end_ladder, qqbart, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, double lambda){ const auto pgin = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto p1 = to_HepLorentzVector(*BeginPart); const auto p2 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon. const double current_factor = ME_qqbar_current( (EndIn-1)->type, pgin, p1, p2, pn, pb )/(4.*(N_C*N_C - 1.)); const double ladder_factor = FKL_ladder_weight( (BeginPart+2), (EndPart-1), pgin-p1-p2, pgin, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, double lambda ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const double current_factor = ME_uno_current( (BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa )/(4.*(N_C*N_C - 1.)); const double ladder_factor = FKL_ladder_weight( (BeginPart+2), (EndPart-1), pa-p1-pg, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } } // namespace double MatrixElement::tree_kin_jets(Event const & ev) const { auto const & incoming = ev.incoming(); const auto partons = tag_extremal_jet_partons(ev); if (ev.type()==event_type::FKL){ 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, param_.regulator_lambda ); } if (ev.type()==event_type::unordered_backward){ return tree_kin_jets_uno(incoming.begin(), incoming.end(), partons.begin(), partons.end(), param_.regulator_lambda); } if (ev.type()==event_type::unordered_forward){ return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(), partons.rbegin(), partons.rend(), param_.regulator_lambda); } if (ev.type()==event_type::extremal_qqbar_backward){ return tree_kin_jets_qqbar(incoming.begin(), incoming.end(), partons.begin(), partons.end(), param_.regulator_lambda); } if (ev.type()==event_type::extremal_qqbar_forward){ return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(), partons.rbegin(), partons.rend(), param_.regulator_lambda); } if (ev.type()==event_type::central_qqbar){ return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type, to_HepLorentzVector(incoming[0]), to_HepLorentzVector(incoming[1]), partons, param_.regulator_lambda); } throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets"); } namespace { double tree_kin_W_FKL( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, double lambda, ParticleProperties const & Wprop ){ auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder = cend(partons) - 1; bool wc = aptype==partons[0].type; //leg b emits w auto q0 = pa - p1; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_current( aptype, bptype, pn, pb, p1, pa, plbar, pl, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_W_uno(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, double lambda, ParticleProperties const & Wprop ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w auto q0 = pa - p1 - pg; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_uno_current( (BeginIn)->type, (BeginIn+1)->type, pn, pb, p1, pa, pg, plbar, pl, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, double lambda, ParticleProperties const & Wprop ){ const bool swap_qqbar=is_quark(*BeginPart); const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1))); const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0))); const auto p1 = to_HepLorentzVector(*(BeginPart)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w auto q0 = pa - pq - pqbar; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_qqbar_current( (BeginIn+1)->type, pa, pb, pq, pqbar, pn, plbar, pl, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } double tree_kin_W_qqbar_mid( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, double lambda, ParticleProperties const & Wprop ){ CLHEP::HepLorentzVector pq; CLHEP::HepLorentzVector pqbar; const auto backmidquark = std::find_if( begin(partons)+1, end(partons)-1, [](Particle const & s){ return s.type != pid::gluon; } ); assert(backmidquark!=end(partons)-1); if (is_quark(backmidquark->type)){ pq = to_HepLorentzVector(*backmidquark); pqbar = to_HepLorentzVector(*(backmidquark+1)); } else { pqbar = to_HepLorentzVector(*backmidquark); pq = to_HepLorentzVector(*(backmidquark+1)); } auto p1 = to_HepLorentzVector(partons.front()); auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; // t-channel momentum after qqbar auto qqbart = q0; bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W bool wc = !wqq && (aptype==partons.front().type); //leg b emits w assert(!wqq || !wc); if(wqq){ // emission from qqbar qqbart -= pl + plbar; } else if(!wc) { // emission from leg a q0 -= pl + plbar; qqbart -= pl + plbar; } const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder_1 = (backmidquark); const auto begin_ladder_2 = (backmidquark+2); const auto end_ladder = cend(partons) - 1; for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){ qqbart -= to_HepLorentzVector(*parton_it); } const int nabove = std::distance(begin_ladder, backmidquark); const int nbelow = std::distance(begin_ladder_2, end_ladder); std::vector partonsHLV; partonsHLV.reserve(partons.size()); for (std::size_t i = 0; i != partons.size(); ++i) { partonsHLV.push_back(to_HepLorentzVector(partons[i])); } const double current_factor = ME_W_qqbar_mid_current( aptype, bptype, nabove, nbelow, pa, pb, pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder_1, q0, pa, pb, p1, pn, lambda )*FKL_ladder_weight( begin_ladder_2, end_ladder, qqbart, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } } // namespace double MatrixElement::tree_kin_W(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); #ifndef NDEBUG // assert that there is exactly one decay corresponding to the W assert(ev.decays().size() == 1); auto const & w_boson{ std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(), [] (Particle const & p) -> bool { return std::abs(p.type) == ParticleID::Wp; }) }; assert(w_boson != ev.outgoing().cend()); assert( static_cast(ev.decays().cbegin()->first) == std::distance(ev.outgoing().cbegin(), w_boson) ); #endif // find decay products of W auto const & decay{ ev.decays().cbegin()->second }; assert(decay.size() == 2); assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) ) || ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) ); // get lepton & neutrino CLHEP::HepLorentzVector plbar; CLHEP::HepLorentzVector pl; if (decay.at(0).type < 0){ plbar = to_HepLorentzVector(decay.at(0)); pl = to_HepLorentzVector(decay.at(1)); } else{ pl = to_HepLorentzVector(decay.at(0)); plbar = to_HepLorentzVector(decay.at(1)); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto partons = tag_extremal_jet_partons(ev); if(ev.type() == FKL){ return tree_kin_W_FKL(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == unordered_backward){ return tree_kin_W_uno(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == unordered_forward){ return tree_kin_W_uno(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == extremal_qqbar_backward){ return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == extremal_qqbar_forward){ return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } assert(ev.type() == central_qqbar); return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } namespace /* WW */ { std::vector tree_kin_WW_FKL( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & pl1bar, CLHEP::HepLorentzVector const & pl1, CLHEP::HepLorentzVector const & pl2bar, CLHEP::HepLorentzVector const & pl2, double lambda, ParticleProperties const & Wprop ){ assert(is_anyquark(aptype)); assert(is_anyquark(bptype)); auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); const std::vector current_factor = ME_WW_current( aptype, bptype, pn, pb, p1, pa, pl1bar, pl1, pl2bar, pl2, Wprop ); auto const begin_ladder = cbegin(partons) + 1; auto const end_ladder = cend(partons) - 1; // pa -> W1 p1, pb -> W2 + pn const auto q0 = pa - p1 - (pl1 + pl1bar); // pa -> W2 p1, pb -> W1 + pn const auto q1 = pa - p1 - (pl2 + pl2bar); const std::vector ladder_factor = FKL_ladder_weight_mix( begin_ladder, end_ladder, q0, q1, pa, pb, p1, pn, lambda ); assert(current_factor.size() == 3); assert(current_factor.size() == ladder_factor.size()); std::vector result(current_factor.size()); for(size_t i=0; i MatrixElement::tree_kin_WW(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); auto const pa = to_HepLorentzVector(incoming[0]); auto const pb = to_HepLorentzVector(incoming[1]); auto const partons = tag_extremal_jet_partons(ev); // W1 & W2 assert(ev.decays().size() == 2); std::vector plbar; std::vector pl; for(auto const & decay_pair : ev.decays()) { auto const decay = decay_pair.second; // TODO: how to label W1, W2 if(decay.at(0).type < 0) { plbar.emplace_back(to_HepLorentzVector(decay.at(0))); pl .emplace_back(to_HepLorentzVector(decay.at(1))); } else { pl .emplace_back(to_HepLorentzVector(decay.at(0))); plbar.emplace_back(to_HepLorentzVector(decay.at(1))); } } if(ev.type() == FKL) { return tree_kin_WW_FKL( incoming[0].type, incoming[1].type, pa, pb, partons, plbar[0], pl[0], plbar[1], pl[1], param_.regulator_lambda, param_.ew_parameters.Wprop() ); } throw std::logic_error("Can only reweight FKL events in WW"); } namespace{ std::vector tree_kin_Z_FKL( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, const double lambda, ParticleProperties const & Zprop, const double stw2, const double ctw ){ const auto p1 = to_HepLorentzVector(partons[0]); const auto pn = to_HepLorentzVector(partons[partons.size() - 1]); const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder = cend(partons) - 1; const std::vector current_factor = ME_Z_current( aptype, bptype, pn, pb, p1, pa, plbar, pl, Zprop, stw2, ctw ); std::vector ladder_factor; if(is_gluon(bptype)){ // This is a qg event const auto q0 = pa-p1-plbar-pl; ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda)); } else if(is_gluon(aptype)){ // This is a gq event const auto q0 = pa-p1; ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda)); } else { // This is a qq event const auto q0 = pa-p1-plbar-pl; const auto q1 = pa-p1; ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder, q0, q1, pa, pb, p1, pn, lambda); } std::vector result; for(size_t i=0; i std::vector tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, const double lambda, ParticleProperties const & Zprop, const double stw2, const double ctw){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const ParticleID aptype = (BeginIn)->type; const ParticleID bptype = (BeginIn+1)->type; const std::vector current_factor = ME_Z_uno_current( aptype, bptype, pn, pb, p1, pa, pg, plbar, pl, Zprop, stw2, ctw ); std::vector ladder_factor; if(is_gluon(bptype)){ // This is a qg event const auto q0 = pa-pg-p1-plbar-pl; ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda)); }else if(is_gluon(aptype)){ // This is a gq event const auto q0 = pa-pg-p1; ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda)); }else{ // This is a qq event const auto q0 = pa-pg-p1-plbar-pl; const auto q1 = pa-pg-p1; ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1, q0, q1, pa, pb, p1, pn, lambda); } std::vector result; for(size_t i=0; i MatrixElement::tree_kin_Z(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); // find decay products of Z auto const & decay{ ev.decays().cbegin()->second }; assert(decay.size() == 2); assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0)) && decay.at(0).type==-decay.at(1).type); // get leptons CLHEP::HepLorentzVector plbar; CLHEP::HepLorentzVector pl; if (decay.at(0).type < 0){ plbar = to_HepLorentzVector(decay.at(0)); pl = to_HepLorentzVector(decay.at(1)); } else{ pl = to_HepLorentzVector(decay.at(0)); plbar = to_HepLorentzVector(decay.at(1)); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto partons = tag_extremal_jet_partons(ev); const double stw2 = param_.ew_parameters.sin2_tw(); const double ctw = param_.ew_parameters.cos_tw(); if(ev.type() == FKL){ return tree_kin_Z_FKL(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } if(ev.type() == unordered_backward){ return tree_kin_Z_uno(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } if(ev.type() == unordered_forward){ return tree_kin_Z_uno(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets"); } double MatrixElement::tree_kin_Higgs(Event const & ev) const { 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); } // kinetic matrix element square for backward Higgs emission // cf. eq:ME_h_jets_peripheral in developer manual, // but without factors \alpha_s and the FKL ladder double MatrixElement::MH2_backwardH( const ParticleID type_forward, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pH, CLHEP::HepLorentzVector const & pn ) const { using namespace currents; const double vev = param_.ew_parameters.vev(); return K(type_forward, pn, pb)*ME_H_gq( pH, pa, pn, pb, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, vev )/(4.*(N_C*N_C - 1)*(pa-pH).m2()*(pb-pn).m2()); } // kinetic matrix element square for backward unordered emission // and forward g -> Higgs double MatrixElement::MH2_unob_forwardH( CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pH ) const { using namespace currents; const double vev = param_.ew_parameters.vev(); constexpr double K_f1 = K_q; constexpr double nhel = 4.; return K_f1*ME_juno_jgH( pg, p1, pa, pH, pb, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, vev )/(nhel*(N_C*N_C - 1)*(pa - p1 - pg).m2()*(pb - pH).m2()); } 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(is_anyquark(incoming.front())) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(ev); } const auto partons = tag_extremal_jet_partons(ev); const auto pa = to_HepLorentzVector(incoming.front()); const auto pb = to_HepLorentzVector(incoming.back()); const auto pH = to_HepLorentzVector(outgoing.front()); const auto end_ladder = end(partons) - ((ev.type() == event_type::unof)?2:1); const auto pn = to_HepLorentzVector(*end_ladder); const double ladder = FKL_ladder_weight( begin(partons), end_ladder, pa - pH, pa, pb, pa, pn, param_.regulator_lambda ); if(ev.type() == event_type::unof) { const auto pg = to_HepLorentzVector(outgoing.back()); return MH2_unob_forwardH( pb, pa, pg, pn, pH )*ladder; } return MH2_backwardH( incoming.back().type, pa, pb, pH, pn )*ladder; } 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(is_anyquark(incoming.back())) { assert(incoming.back().type == outgoing[outgoing.size()-2].type); return tree_kin_Higgs_between(ev); } const auto partons = tag_extremal_jet_partons(ev); const auto pa = to_HepLorentzVector(incoming.front()); const auto pb = to_HepLorentzVector(incoming.back()); const auto pH = to_HepLorentzVector(outgoing.back()); auto begin_ladder = begin(partons) + 1; auto q0 = pa - to_HepLorentzVector(partons.front()); if(ev.type() == event_type::unob) { q0 -= to_HepLorentzVector(*begin_ladder); ++begin_ladder; } const auto p1 = to_HepLorentzVector(*(begin_ladder - 1)); const double ladder = FKL_ladder_weight( begin_ladder, end(partons), q0, pa, pb, p1, pb, param_.regulator_lambda ); if(ev.type() == event_type::unob) { const auto pg = to_HepLorentzVector(outgoing.front()); return MH2_unob_forwardH( pa, pb, pg, p1, pH )*ladder; } return MH2_backwardH( incoming.front().type, pb, pa, pH, p1 )*ladder; } namespace { template double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, CLHEP::HepLorentzVector const & qH, CLHEP::HepLorentzVector const & qHp1, double mt, bool inc_bot, double mb, double vev ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); return ME_Higgs_current_uno( (BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa, qH, qHp1, mt, inc_bot, mb, vev ); } } // namespace 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 = NAN; if(ev.type() == FKL){ 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, param_.ew_parameters.vev() ); } else if(ev.type() == unob){ current_factor = tree_kin_Higgs_uno( begin(incoming), end(incoming), begin(partons), end(partons), qH, qH-pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, param_.ew_parameters.vev() ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(ev.type() == unof){ current_factor = tree_kin_Higgs_uno( rbegin(incoming), rend(incoming), rbegin(partons), rend(partons), qH-pH, qH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, param_.ew_parameters.vev() ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ throw std::logic_error("Can only reweight FKL or uno processes in H+Jets"); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, q0, pa, pb, p1, pn, param_.regulator_lambda )*FKL_ladder_weight( first_after_Higgs, end_ladder, qH - pH, pa, pb, p1, pn, param_.regulator_lambda ); return current_factor/(4.*(N_C*N_C-1.))*ladder_factor; } namespace { double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) { std::vector bosons = filter_AWZH_bosons(ev.outgoing()); if(bosons.empty()) { return 1.; } if(bosons.size() == 1) { switch(bosons[0].type){ case pid::Higgs: return alpha_s*alpha_s; case pid::Wp: case pid::Wm: return alpha_w*alpha_w; case pid::Z_photon_mix: return alpha_w*alpha_w; // TODO case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } if(bosons.size() == 2) { if(bosons[0].type == pid::Wp && bosons[1].type == pid::Wp) { return alpha_w*alpha_w*alpha_w*alpha_w; } else if(bosons[0].type == pid::Wm && bosons[1].type == pid::Wm) { return alpha_w*alpha_w*alpha_w*alpha_w; } throw not_implemented("Emission of bosons of unsupported type"); } throw not_implemented("Emission of >2 bosons is unsupported"); } } // namespace double MatrixElement::tree_param(Event const & ev, double mur) const { assert(is_resummable(ev.type())); const auto begin_partons = ev.begin_partons(); const auto end_partons = ev.end_partons(); const auto num_partons = std::distance(begin_partons, end_partons); const double alpha_s = alpha_s_(mur); const double gs2 = 4.*M_PI*alpha_s; double res = std::pow(gs2, num_partons); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(num_partons >= 2); const auto first_emission = std::next(begin_partons); const auto last_emission = std::prev(end_partons); for(auto parton = first_emission; parton != last_emission; ++parton){ res *= 1. + alpha_s/(2.*M_PI)*BETA0*std::log(mur/parton->perp()); } } return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res; } } // namespace HEJ diff --git a/src/Wjets.cc b/src/Wjets.cc index dadae22..007b3cb 100644 --- a/src/Wjets.cc +++ b/src/Wjets.cc @@ -1,520 +1,524 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2020-2022 * \copyright GPLv2 or later */ #include "HEJ/Wjets.hh" #include #include #include #include #include #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/exceptions.hh" // generated headers #include "HEJ/currents/jV_j.hh" #include "HEJ/currents/jV_juno.hh" #include "HEJ/currents/jVuno_j.hh" -#include "HEJ/currents/jW_jqqbar.hh" +#include "HEJ/currents/jV_jqqbar.hh" #include "HEJ/currents/jW_qqbar_j.hh" #include "HEJ/currents/jWqqbar_j.hh" #include "HEJ/currents/j_Wqqbar_j.hh" namespace HEJ { namespace currents { namespace { using COM = std::complex; // --- Helper Functions --- // FKL W Helper Functions double WProp(const HLV & plbar, const HLV & pl, ParticleProperties const & wprop ){ COM propW = COM(0.,-1.)/( (pl+plbar).m2() - wprop.mass*wprop.mass + COM(0.,1.)*wprop.mass*wprop.width); double PropFactor=(propW*conj(propW)).real(); return PropFactor; } + } // Anonymous Namespace + + //! W+Jets FKL Contributions + double ME_W_qQ( + HLV const & p1out, + HLV const & plbar, HLV const & pl, + HLV const & p1in, + HLV const & p2out, HLV const & p2in, + ParticleProperties const & wprop + ){ + using helicity::minus; + using helicity::plus; + const COM ampm = jV_j(p1in,p1out,p2in,p2out,pl,plbar); + const COM ampp = jV_j (p1in,p1out,p2in,p2out,pl,plbar); + + const double Msqr = std::norm(ampm) + std::norm(ampp); + + return WProp(plbar, pl, wprop) * Msqr; + } + + namespace { /** * @brief Contraction of W + unordered jet current with FKL current * * See eq:wunocurrent in the developer manual for the definition * of the W + unordered jet current */ template - double jM2Wuno( + double amp_jWuno_j( HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl, HLV const & pa, HLV const & p2, HLV const & pb, ParticleProperties const & wprop ){ using helicity::minus; const COM u1 = U1(p1, p2, pa, pb, pg, pl, plbar); const COM u2 = U2(p1, p2, pa, pb, pg, pl, plbar); - const COM l = L(p1, p2, pa, pb, pg, pl, plbar); + const COM l = L (p1, p2, pa, pb, pg, pl, plbar); const COM x = u1 - l; const COM y = u2 + l; // eq:VunoSumAveS in developer manual // TODO: use same form as for other uno currents, // extracting at least a factor of K_q = C_F const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real(); return WProp(plbar, pl, wprop) * amp; } - /** - * @brief Contraction of W + extremal qqbar current with FKL current - * - * See eq:crossed in the developer manual for the definition - * of the W + extremal qqbar current. - * - */ - template - double jM2_W_extremal_qqbar( - HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl, - HLV const & pqbar, HLV const & p3, HLV const & pb, - ParticleProperties const & wprop + // helicity sum helper for jWuno_j(...) + template + double jWuno_j_helsum( + HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl, + HLV const & pa, HLV const & p2, HLV const & pb, + ParticleProperties const & wprop ){ + using helicity::minus; + using helicity::plus; - const COM u1Xcontr = U1X(pg, pq, plbar, pl, pqbar, p3, pb); - const COM u2Xcontr = U2X(pg, pq, plbar, pl, pqbar, p3, pb); - const COM lXcontr = LX(pg, pq, plbar, pl, pqbar, p3, pb); + const double ME2h1pp = amp_jWuno_j( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); - const COM x = u1Xcontr - lXcontr; - const COM y = u2Xcontr + lXcontr; + const double ME2h1pm = amp_jWuno_j( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); - const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real(); - return WProp(plbar, pl, wprop) * amp; + const double ME2h1mp = amp_jWuno_j( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); + + const double ME2h1mm = amp_jWuno_j( + pg, p1, plbar, pl, pa, p2, pb, wprop + ); + + return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm; } } // Anonymous Namespace - //! W+Jets FKL Contributions - double ME_W_qQ( - HLV const & p1out, - HLV const & plbar, HLV const & pl, - HLV const & p1in, - HLV const & p2out, HLV const & p2in, - ParticleProperties const & wprop + double ME_Wuno_qQ( + HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, + HLV const & pg, HLV const & plbar, HLV const & pl, + ParticleProperties const & wprop ){ - using helicity::minus; - using helicity::plus; - - const COM ampm = jV_j(p1in,p1out,p2in,p2out,pl,plbar); - const COM ampp = jV_j(p1in,p1out,p2in,p2out,pl,plbar); - - const double Msqr = std::norm(ampm) + std::norm(ampp); - - return WProp(plbar, pl, wprop) * Msqr; + return jWuno_j_helsum( + pg, p1out, plbar, pl, p1in, p2out, p2in, wprop + )/(4.*C_A*C_A); } namespace { // helicity amplitude squares for contraction of W current with unordered // current template - double ampsq_jW_juno( + double amp_jW_juno( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p2, HLV const & pg, HLV const & pl, HLV const & plbar ){ using helicity::minus; const COM u1 = U1_jV(pa,p1,pb,p2,pg,pl,plbar); const COM u2 = U2_jV(pa,p1,pb,p2,pg,pl,plbar); - const COM l = L_jV(pa,p1,pb,p2,pg,pl,plbar); + const COM l = L_jV (pa,p1,pb,p2,pg,pl,plbar); const COM x = u1 - l; const COM y = u2 + l; return C_F*norm(x + y) - C_A*(x*conj(y)).real(); } } // Anonymous Namespace double ME_W_unob_qQ( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & pg, HLV const & plbar, HLV const & pl, ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; // helicity assignments assume quarks // in the end, this is irrelevant since we sum over all helicities const double ampsq = - + ampsq_jW_juno(p2in,p2out,p1in,p1out,pg,pl,plbar) - + ampsq_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) - + ampsq_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) - + ampsq_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) + + amp_jW_juno(p2in,p2out,p1in,p1out,pg,pl,plbar) + + amp_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) + + amp_jW_juno(p2in,p2out,p1in,p1out,pg,pl,plbar) + + amp_jW_juno (p2in,p2out,p1in,p1out,pg,pl,plbar) ; return WProp(plbar, pl, wprop)*ampsq; } namespace { + /** + * @brief Contraction of W + extremal qqbar current with FKL current + * + * See eq:crossed in the developer manual for the definition + * of the W + extremal qqbar current. + * + */ + template + double amp_jWqqbar_j( + HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl, + HLV const & pqbar, HLV const & p3, HLV const & pb, + ParticleProperties const & wprop + ){ - // helicity sum helper for jWuno_j(...) + const COM u1Xcontr = U1X(pg, pq, plbar, pl, pqbar, p3, pb); + const COM u2Xcontr = U2X(pg, pq, plbar, pl, pqbar, p3, pb); + const COM lXcontr = LX (pg, pq, plbar, pl, pqbar, p3, pb); + + const COM x = u1Xcontr - lXcontr; + const COM y = u2Xcontr + lXcontr; + + const double amp = C_A*C_F*C_F/2.*(norm(x)+norm(y)) - C_F/2.*(x*conj(y)).real(); + return WProp(plbar, pl, wprop) * amp; + } + + // helicity sum helper for jWqqbar_j(...) template - double jWuno_j_helsum( - HLV const & pg, HLV const & p1, HLV const & plbar, HLV const & pl, - HLV const & pa, HLV const & p2, HLV const & pb, - ParticleProperties const & wprop + double jWqqbar_j_helsum( + HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl, + HLV const & pqbar, HLV const & p3, HLV const & pb, + ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; - const double ME2h1pp = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1pp = amp_jWqqbar_j( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1pm = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1pm = amp_jWqqbar_j( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1mp = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1mp = amp_jWqqbar_j( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - - const double ME2h1mm = jM2Wuno( - pg, p1, plbar, pl, pa, p2, pb, wprop + const double amp_h1mm = amp_jWqqbar_j( + pg, pq, plbar, pl, pqbar, p3, pb, wprop ); - return ME2h1pp + ME2h1pm + ME2h1mp + ME2h1mm; + return amp_h1pp + amp_h1pm + amp_h1mp + amp_h1mm; } - } // Anonymous Namespace - double ME_Wuno_qQ( - HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, - HLV const & pg, HLV const & plbar, HLV const & pl, - ParticleProperties const & wprop - ){ - return jWuno_j_helsum( - pg, p1out, plbar, pl, p1in, p2out, p2in, wprop - )/(4.*C_A*C_A); - } - - // helicity sum helper for jWqqbar_j(...) - template - double jWqqbar_j_helsum( - HLV const & pg, HLV const & pq, HLV const & plbar, HLV const & pl, - HLV const & pqbar, HLV const & p3, HLV const & pb, - ParticleProperties const & wprop - ){ - using helicity::minus; - using helicity::plus; - - const double amp_h1pp = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1pm = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1mp = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - const double amp_h1mm = jM2_W_extremal_qqbar( - pg, pq, plbar, pl, pqbar, p3, pb, wprop - ); - - return amp_h1pp + amp_h1pm + amp_h1mp + amp_h1mm; - } - double ME_WExqqbar_qqbarQ( HLV const & pgin, HLV const & pqbarout, HLV const & plbar, HLV const & pl, HLV const & pqout, HLV const & p2out, HLV const & p2in, ParticleProperties const & wprop ){ //Helicity sum and average over initial states. double ME2 = jWqqbar_j_helsum( pgin, pqbarout, plbar, pl, pqout, p2out, p2in, wprop )/(4.*C_A*C_A); //Correct colour averaging after crossing: ME2*=(3.0/8.0); return ME2; } namespace { - template - double jW_jqqbar( - HLV const & pb, - HLV const & p2, - HLV const & p3, - HLV const & pa, - HLV const & p1, - HLV const & pl, - HLV const & plbar + // helicity amplitude squares for contraction of W current with qqbar + // current + template + double amp_jW_jqqbar( + HLV const & pa, HLV const & pb, + HLV const & p1, HLV const & pq, HLV const & pqbar, + HLV const & pl, HLV const & plbar ){ - using std::norm; + using helicity::minus; - static constexpr double cm1m1 = 8./3.; - static constexpr double cm2m2 = 8./3.; - static constexpr double cm3m3 = 6.; - static constexpr double cm1m2 =-1./3.; - static constexpr COM cm1m3 = COM{0.,-3.}; - static constexpr COM cm2m3 = COM{0.,3.}; + const COM u1 = U1X_jV(pa,pb,p1,pq,pqbar,pl,plbar); + const COM u2 = U2X_jV(pa,pb,p1,pq,pqbar,pl,plbar); + const COM l = LX_jV (pa,pb,p1,pq,pqbar,pl,plbar); - const COM m1 = jW_qggm1(pb,p2,p3,pa,p1,pl,plbar); - const COM m2 = jW_qggm2(pb,p2,p3,pa,p1,pl,plbar); - const COM m3 = jW_qggm3(pb,p2,p3,pa,p1,pl,plbar); + const COM x = u1 - l; + const COM y = u2 + l; - return - + cm1m1*norm(m1) - + cm2m2*norm(m2) - + cm3m3*norm(m3) - + 2.*real(cm1m2*m1*conj(m2)) - + 2.*real(cm1m3*m1*conj(m3)) - + 2.*real(cm2m3*m2*conj(m3)) ; + return C_F*norm(x + y) - C_A*(x*conj(y)).real(); + } + + // helicity sum helper for jW_jqqbar(...) + template + double jW_jqqbar_helsum( + HLV const & pa, HLV const & pb, + HLV const & p1, HLV const & pq, HLV const & pqbar, + HLV const & pl, HLV const & plbar + ){ + using helicity::minus; + using helicity::plus; + + return amp_jW_jqqbar(pa,pb,p1,pq,pqbar,pl,plbar) + + amp_jW_jqqbar (pa,pb,p1,pq,pqbar,pl,plbar) + + amp_jW_jqqbar(pa,pb,p1,pq,pqbar,pl,plbar) + + amp_jW_jqqbar (pa,pb,p1,pq,pqbar,pl,plbar); } } // Anonymous Namespace // contraction of W current with extremal qqbar on the other side double ME_W_Exqqbar_QQq( HLV const & pa, HLV const & pb, - HLV const & p1, HLV const & p2, - HLV const & p3, HLV const & plbar, HLV const & pl, bool aqlinepa, + HLV const & p1, HLV const & pq, HLV const & pqbar, + HLV const & plbar, HLV const & pl, bool aqlinepa, ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; const double wPropfact = WProp(plbar, pl, wprop); - const double prefactor = 2.*wPropfact/24./4.; - if(aqlinepa) { - return prefactor*( - jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) - + jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) - ); - } - return prefactor*( - jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) - + jW_jqqbar(pb,p2,p3,pa,p1,pl,plbar) - ); + + const double ME2 = aqlinepa? + jW_jqqbar_helsum(pa,pb,p1,pq,pqbar,pl,plbar): + jW_jqqbar_helsum(pa,pb,p1,pq,pqbar,pl,plbar); + + // correct colour averaging after crossing + const double avg_fac = N_C / (N_C*N_C - 1.); + + return avg_fac * wPropfact * ME2; } namespace { // helper function for matrix element for W + Jets with central qqbar template double amp_WCenqqbar_qq( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p4, HLV const & pq, HLV const & pqbar, HLV const & pl, HLV const & plbar, HLV const & q11, HLV const & q12 ){ using std::norm; const COM sym = M_sym_W( pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12 ); const COM cross = M_cross_W( pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12 ); const COM uncross = M_uncross_W( pa, p1, pb, p4, pq, pqbar, pl, plbar, q11, q12 ); // Colour factors static constexpr double cmsms = 3.; static constexpr double cmumu = 4./3.; static constexpr double cmcmc = 4./3.; static constexpr COM cmsmu = COM{0., 3./2.}; static constexpr COM cmsmc = COM{0.,-3./2.}; static constexpr double cmumc = -1./6.; return +cmsms*norm(sym) +cmumu*norm(uncross) +cmcmc*norm(cross) +2.*real(cmsmu*sym*conj(uncross)) +2.*real(cmsmc*sym*conj(cross)) +2.*real(cmumc*uncross*conj(cross)) ; } } // Anonymous Namespace // matrix element for W + Jets with W emission off central qqbar double ME_WCenqqbar_qq( HLV const & pa, HLV const & pb, HLV const & pl, HLV const & plbar, std::vector const & partons, bool /* aqlinepa */, bool /* aqlinepb */, bool qqbar_marker, int nabove, ParticleProperties const & wprop ){ using helicity::plus; using helicity::minus; CLHEP::HepLorentzVector q1 = pa; for(int i = 0; i <= nabove; ++i){ q1 -= partons[i]; } const auto qq = split_into_lightlike(q1); // since q1.m2() < 0 the following assertion is always true // see documentation for split_into_lightlike assert(qq.second.e() < 0); HLV pq; HLV pqbar; if (qqbar_marker){ pqbar = partons[nabove+1]; pq = partons[nabove+2];} else{ pq = partons[nabove+1]; pqbar = partons[nabove+2]; } const HLV p1 = partons.front(); const HLV p4 = partons.back(); // 4 Different Helicity Choices (Differs from Pure Jet Case, where there is // also the choice in qqbar helicity. // the first helicity label is for aqlinepa == true, // the second one for aqlinepb == true // In principle the corresponding helicity should be flipped // if either of them is false, but we sum up all helicities, // so this has no effect in the end. const double amp_mm = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_mp = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pm = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pp = amp_WCenqqbar_qq( pa, p1, pb, p4, pq, pqbar, pl, plbar, qq.first, -qq.second ); // Divide by t channels, extremal + adjacent central vertex const double ta = (pa-p1).m2(); const double t1 = q1.m2(); const double t3 = (q1-pl-plbar-pq-pqbar).m2(); const double tb = (p4-pb).m2(); const double amp = WProp(plbar, pl, wprop)*( amp_mm+amp_mp+amp_pm+amp_pp )/(9.*4.*ta*t1*t3*tb); return amp; } - // helper function for matrix element for W + Jets with central qqbar - // W emitted off extremal parton - // @TODO: code duplication with amp_WCenqqbar_qq - template - double amp_W_Cenqqbar_qq( - HLV const & pa, HLV const & p1, - HLV const & pb, HLV const & pn, - HLV const & pq, HLV const & pqbar, - HLV const & pl, HLV const & plbar, - HLV const & q11, HLV const & q12 - ){ - using std::norm; + namespace { + // helper function for matrix element for W + Jets with central qqbar + // W emitted off extremal parton + // @TODO: code duplication with amp_WCenqqbar_qq + template + double amp_W_Cenqqbar_qq( + HLV const & pa, HLV const & p1, + HLV const & pb, HLV const & pn, + HLV const & pq, HLV const & pqbar, + HLV const & pl, HLV const & plbar, + HLV const & q11, HLV const & q12 + ){ + using std::norm; - const COM crossed = M_cross( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); - const COM uncrossed = M_qbar( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); - const COM sym = M_sym( - pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 - ); + const COM crossed = M_cross( + pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 + ); + const COM uncrossed = M_qbar( + pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 + ); + const COM sym = M_sym( + pa, p1, pb, pn, pq, pqbar, pl, plbar, q11, q12 + ); - //Colour factors: - static constexpr double cmsms = 3.; - static constexpr double cmumu = 4./3.; - static constexpr double cmcmc = 4./3.; - static constexpr COM cmsmu = COM{0.,3./2.}; - static constexpr COM cmsmc = COM{0.,-3./2.}; - static constexpr double cmumc = -1./6.; - - return - +cmsms*norm(sym) - +cmumu*norm(uncrossed) - +cmcmc*norm(crossed) - +2.*real(cmsmu*sym*conj(uncrossed)) - +2.*real(cmsmc*sym*conj(crossed)) - +2.*real(cmumc*uncrossed*conj(crossed)) - ; - } + //Colour factors: + static constexpr double cmsms = 3.; + static constexpr double cmumu = 4./3.; + static constexpr double cmcmc = 4./3.; + static constexpr COM cmsmu = COM{0.,3./2.}; + static constexpr COM cmsmc = COM{0.,-3./2.}; + static constexpr double cmumc = -1./6.; + + return + +cmsms*norm(sym) + +cmumu*norm(uncrossed) + +cmcmc*norm(crossed) + +2.*real(cmsmu*sym*conj(uncrossed)) + +2.*real(cmsmc*sym*conj(crossed)) + +2.*real(cmumc*uncrossed*conj(crossed)) + ; + } + } // Anonymous Namespace // matrix element for W + Jets with W emission *not* off central qqbar double ME_W_Cenqqbar_qq( HLV pa, HLV pb, HLV pl,HLV plbar, std::vector partons, bool aqlinepa, bool aqlinepb, bool qqbar_marker, int nabove, int nbelow, bool forwards, ParticleProperties const & wprop ){ using helicity::minus; using helicity::plus; if (!forwards){ //If Emission from Leg a instead, flip process. std::swap(pa, pb); std::reverse(partons.begin(),partons.end()); std::swap(aqlinepa, aqlinepb); qqbar_marker = !qqbar_marker; std::swap(nabove,nbelow); } HLV q1=pa; for(int i=0;i( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_mp = amp_W_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pm = amp_W_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); const double amp_pp = amp_W_Cenqqbar_qq( pa, p1, pb, pn, pq, pqbar, pl, plbar, qq.first, -qq.second ); // Divide by t channels, extremal + adjacent central vertex const double ta = (pa-p1).m2(); const double t1 = q1.m2(); const double t3 = (q1 - pq - pqbar).m2(); const double tb = (pn+pl+plbar-pb).m2(); const double amp= WProp(plbar, pl, wprop)*( amp_mm+amp_mp+amp_pm+amp_pp )/(9.*4.*ta*t1*t3*tb); return amp; } } // namespace currents } // namespace HEJ diff --git a/t/ME_data/ME_Wm.dat b/t/ME_data/ME_Wm.dat index 43575ed..63b5b45 100644 --- a/t/ME_data/ME_Wm.dat +++ b/t/ME_data/ME_Wm.dat @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2dbb36e04293c9a49f21a4ef015715bda5e32872e4db928edb46d8c11344a9f2 +oid sha256:328d808a5aa481f2e4ace288398ce83b973a991d46d476751b13f38c7d04cdf9 size 80000 diff --git a/t/ME_data/ME_Wp.dat b/t/ME_data/ME_Wp.dat index 1b16be0..25150b0 100644 --- a/t/ME_data/ME_Wp.dat +++ b/t/ME_data/ME_Wp.dat @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:a257770ad3701f172d7e12549f334571724ff3ccf01188961d4b77a015d158b3 +oid sha256:e0e4ab5364715d7ff32d828bf662590ea0b9b1970a018f71fd1ff76d8b5ac67f size 80000