diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index b1507e4..b6eff06 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,996 +1,996 @@ #include "RHEJ/MatrixElement.hh" #include #include #include "RHEJ/Constants.hh" #include "RHEJ/currents.hh" #include "RHEJ/PDG_codes.hh" #include "RHEJ/uno.hh" #include "RHEJ/utility.hh" namespace RHEJ{ //cf. last line of eq. (22) in \ref Andersen:2011hs double MatrixElement::omega0( double alpha_s, double mur, fastjet::PseudoJet const & q_j, double lambda ) const { const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; // use alpha_s(sqrt(q_j*lambda)), evolved to mur return ( 1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda)) )*result; } double MatrixElement::virtual_corrections( double mur, std::array const & in, std::vector const & out ) const{ fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif assert(std::is_sorted(out.begin(), out.end(), rapidity_less{})); assert(out.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - out[0].p; size_t first_idx = 0; size_t last_idx = out.size() - 1; // if there is a Higgs or unordered gluon outside the extremal partons // then it is not part of the FKL ladder and does not contribute // to the virtual corrections if(out.front().type == pid::Higgs || has_unob_gluon(in, out)){ q -= out[1].p; ++first_idx; } if(out.back().type == pid::Higgs || has_unof_gluon(in, out)){ --last_idx; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q, CLAMBDA)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } assert( nearby(q, -1*pb, norm) || out.back().type == pid::Higgs || has_unof_gluon(in, out) ); return exp(exponent); } } // namespace RHEJ namespace { //! Lipatov vertex for partons emitted into extremal jets double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2)) - p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2)); // cout << "#Fadin qa : "<fabs(CL.dot(CL))) // not sufficient! // return 0.; // else return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>RHEJ::CLAMBDA) return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2())); return Cls-4./(kperp*kperp); } } //! Lipatov vertex double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip, CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2. - qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2. + ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom)) + pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom)) - pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.; return -CL.dot(CL); } //! Lipatov vertex with soft subtraction double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>RHEJ::CLAMBDA) return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2())); double temp=Cls-4./(kperp*kperp); return temp; } } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ if (aptype==21&&bptype==21) { return jM2gg(pn,pb,p1,pa); } else if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2qg(pn,pb,p1,pa); else return jM2qbarg(pn,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jM2qg(p1,pa,pn,pb); else return jM2qbarg(p1,pa,pn,pb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2qQ(pn,pb,p1,pa); else return jM2qQbar(pn,pb,p1,pa); } else { if (aptype>0) return jM2qQbar(p1,pa,pn,pb); else return jM2qbarQbar(pn,pb,p1,pa); } } throw std::logic_error("unknown particle types"); } /** 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 * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_W_current( int aptype, int 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 ){ if (aptype==21&&bptype==21) { throw std::logic_error("gg incoming in W+jets, qqx not yet implemented"); } else if (aptype==21&&bptype!=21) { if (bptype > 0) return jMWqg(pn,pl,plbar,pb,p1,pa); else return jMWqbarg(pn,pl,plbar,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jMWqg(p1,pl,plbar,pa,pn,pb); else return jMWqbarg(p1,pl,plbar,pa,pn,pb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jMWqQ(pn,pl,plbar,pb,p1,pa); else return jMWqQbar(pn,pl,plbar,pb,p1,pa); } else { if (aptype>0) return jMWqQbar(p1,pl,plbar,pa,pn,pb); else return jMWqbarQbar(pn,pl,plbar,pb,p1,pa); } } throw std::logic_error("unknown particle types"); } /** \brief Matrix element squared for tree-level current-current scattering with Higgs * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared for Tree-Level Current-Current Scattering with Higgs */ double ME_Higgs_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype==21) // gg initial state return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else if (aptype==21&&bptype!=21) { if (bptype > 0) return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; else return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; } else if (bptype==21&&aptype!=21) { if (aptype > 0) return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; else return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; } else { // they are both quark if (bptype>0) { if (aptype>0) return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } else { if (aptype>0) return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered forward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param punof Unordered Particle Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered forward emission */ double ME_Higgs_current_unof( int aptype, int bptype, CLHEP::HepLorentzVector const & punof, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (aptype>0) return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered backward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param punob Unordered back Particle Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered backward emission */ double ME_Higgs_current_unob( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & punob, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (bptype==21&&aptype!=21) { if (aptype > 0) return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (aptype>0) { if (bptype>0) return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (bptype>0) return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Particle const & particle){ return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()}; } void validate(RHEJ::MatrixElementConfig const & config) { #ifndef RHEJ_BUILD_WITH_QCDLOOP if(!config.Higgs_coupling.use_impact_factors) { throw std::invalid_argument{ "Invalid Higgs coupling settings.\n" "HEJ without QCDloop support can only use impact factors.\n" "Set use_impact_factors to true or recompile HEJ.\n" }; } #endif if(config.Higgs_coupling.use_impact_factors && config.Higgs_coupling.mt != std::numeric_limits::infinity()) { throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } } // namespace anonymous namespace RHEJ{ MatrixElement::MatrixElement( std::function alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } double MatrixElement::operator()( double mur, std::array const & incoming, std::vector const & outgoing, std::unordered_map> const & decays, bool check_momenta ) const { return tree( mur, incoming, outgoing, decays, check_momenta )*virtual_corrections( mur, incoming, outgoing ); } double MatrixElement::tree_kin( std::array const & incoming, std::vector const & outgoing, std::unordered_map> const & decays, bool check_momenta ) const { assert( std::is_sorted( incoming.begin(), incoming.end(), [](Particle o1, Particle o2){return o1.p.pz()type){ case pid::Higgs: { return tree_kin_Higgs(incoming, outgoing, check_momenta); } // TODO case pid::Wp: { return tree_kin_W(incoming, outgoing, decays, true, check_momenta); } case pid::Wm: { return tree_kin_W(incoming, outgoing, decays, false, check_momenta); } case pid::photon: case pid::Z: default: throw std::logic_error("Emission of boson of unsupported type."); } } namespace{ constexpr int extremal_jet_idx = 1; constexpr int no_extremal_jet_idx = 0; bool treat_as_extremal(Particle const & parton){ return parton.p.user_index() == extremal_jet_idx; } template double FKL_ladder_weight( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn ){ double wt = 1; auto qi = q0; for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){ assert(gluon_it->type == pid::gluon); const auto g = to_HepLorentzVector(*gluon_it); const auto qip1 = qi - g; if(treat_as_extremal(*gluon_it)){ wt *= C2Lipatovots(qip1, qi, pa, pb)*C_A; } else{ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A; } qi = qip1; } return wt; } } // namespace anonymous std::vector MatrixElement::tag_extremal_jet_partons( std::array const & incoming, std::vector out_partons, bool check_momenta ) const{ if(!check_momenta){ for(auto & parton: out_partons){ parton.p.set_user_index(no_extremal_jet_idx); } return out_partons; } fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def); const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt)); assert(jets.size() >= 2); auto most_backward = begin(jets); auto most_forward = end(jets) - 1; // skip jets caused by unordered emission if(has_unob_gluon(incoming, out_partons)){ assert(jets.size() >= 3); ++most_backward; } else if(has_unof_gluon(incoming, out_partons)){ assert(jets.size() >= 3); --most_forward; } const auto extremal_jet_indices = cs.particle_jet_indices( {*most_backward, *most_forward} ); assert(extremal_jet_indices.size() == out_partons.size()); for(size_t i = 0; i < out_partons.size(); ++i){ assert(RHEJ::is_parton(out_partons[i])); const int idx = (extremal_jet_indices[i]>=0)? extremal_jet_idx: no_extremal_jet_idx; out_partons[i].p.set_user_index(idx); } return out_partons; } double MatrixElement::tree_kin_jets( std::array const & incoming, std::vector partons, bool check_momenta ) const { partons = tag_extremal_jet_partons(incoming, partons, check_momenta); if(has_unob_gluon(incoming, partons) || has_unof_gluon(incoming, partons)){ throw std::logic_error("unordered emission not implemented for pure jets"); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); return ME_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa )/(4*(N_C*N_C - 1))*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, pa - p1, pa, pb, p1, pn ); } double MatrixElement::tree_kin_W( std::array const & incoming, std::vector const & outgoing, std::unordered_map> const & decays, bool WPlus, bool check_momenta ) const { bool qqxExtremal=false; bool qqxCentral=false; if(has_unob_gluon(incoming, outgoing)){ throw std::logic_error("unordered emission not yet implemented for W+jets"); //return tree_kin_W_unob(incoming, outgoing, check_momenta); } else if(has_unof_gluon(incoming, outgoing)){ throw std::logic_error("unordered emission not yet implemented for W+jets"); // return tree_kin_W_unof(incoming, outgoing, check_momenta); } else if(qqxExtremal){ throw std::logic_error("Extremal qqx not yet implemented for W+jets"); // return tree_kin_W_Exqqx(incoming, outgoing, check_momenta); } else if(qqxCentral){ throw std::logic_error("Central qqx not yet implemented for W+jets"); // return tree_kin_W_qqxCentral(incoming, outgoing, check_momenta); } else{ return tree_kin_W_FKL(incoming, outgoing, decays, WPlus, check_momenta); } } double MatrixElement::tree_kin_W_FKL( std::array const & incoming, std::vector const & outgoing, std::unordered_map> const & decays, bool WPlus, bool check_momenta ) const { const auto the_W = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return abs(s.type) == pid::Wp; } ); HLV plbar, pl; for (auto& x: decays) { if (x.second.at(0).type < 0){ plbar = to_HepLorentzVector(x.second.at(0)); pl = to_HepLorentzVector(x.second.at(1)); } else{ pl = to_HepLorentzVector(x.second.at(0)); plbar = to_HepLorentzVector(x.second.at(1)); } } const auto pW = to_HepLorentzVector(*the_W); std::vector partons(begin(outgoing), the_W); partons.insert(end(partons), the_W + 1, end(outgoing)); partons = tag_extremal_jet_partons(incoming, partons, check_momenta); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto first_after_W = begin(partons) + (the_W-begin(outgoing)); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_W == begin(partons)) ++first_after_W; else if(first_after_W == end(partons)) --first_after_W; // t-channel momentum before W auto qW = pa; for(auto parton_it = begin(partons); parton_it != first_after_W; ++parton_it){ qW -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if (WPlus){ current_factor = ME_W_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, pl, plbar ); } else{ current_factor = ME_W_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, plbar, pl ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_W, q0, pa, pb, p1, pn )*FKL_ladder_weight( begin_ladder, end_ladder, pa - p1, pa, pb, p1, pn ); return current_factor*9./8.*ladder_factor; } double MatrixElement::tree_kin_Higgs( std::array const & incoming, std::vector const & outgoing, bool check_momenta ) const { if(has_uno_gluon(incoming, outgoing)){ return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } if(outgoing.front().type == pid::Higgs){ return tree_kin_Higgs_first(incoming, outgoing, check_momenta); } if(outgoing.back().type == pid::Higgs){ return tree_kin_Higgs_last(incoming, outgoing, check_momenta); } return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } namespace { // Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113 // TODO: code duplication with currents.cc double K_g(double p1minus, double paminus) { return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1/C_A) + 1/C_A; } double K_g( CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(type == ParticleID::gluon) return K_g(pout, pin); return C_F; } // Colour factor in strict MRK limit double K_MRK(ParticleID type) { return (type == ParticleID::gluon)?C_A:C_F; } } double MatrixElement::MH2_forwardH( CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, ParticleID type2, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double t1, double t2 ) const{ ignore(p2out, p2in); const double shat = p1in.invariantMass2(p2in); // gluon case #ifdef RHEJ_BUILD_WITH_QCDLOOP if(!param_.Higgs_coupling.use_impact_factors){ - return K(type2, p2out, p2in)/C_F*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH( + return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH( p1out, p1in, p2out, p2in, pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb - ); + )/(4*(N_C*N_C - 1)); } #endif return K_MRK(type2)/C_A*9./2.*shat*shat*( C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH) )/(t1*t2); } double MatrixElement::tree_kin_Higgs_first( std::array const & incoming, std::vector const & outgoing, bool check_momenta ) const { assert(outgoing.front().type == pid::Higgs); if(outgoing[1].type != pid::gluon) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } const auto pH = to_HepLorentzVector(outgoing.front()); const auto partons = tag_extremal_jet_partons( incoming, std::vector(begin(outgoing) + 1, end(outgoing)), check_momenta ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); const auto q0 = pa - p1 - pH; const double t1 = q0.m2(); const double t2 = (pn - pb).m2(); return MH2_forwardH( p1, pa, incoming[1].type, pn, pb, pH, t1, t2 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_last( std::array const & incoming, std::vector const & outgoing, bool check_momenta ) const { assert(outgoing.back().type == pid::Higgs); if(outgoing[outgoing.size()-2].type != pid::gluon) { assert(incoming.back().type == outgoing[outgoing.size()-2].type); return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } const auto pH = to_HepLorentzVector(outgoing.back()); const auto partons = tag_extremal_jet_partons( incoming, std::vector(begin(outgoing), end(outgoing) - 1), check_momenta ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; const double t1 = q0.m2(); const double t2 = (pn + pH - pb).m2(); return MH2_forwardH( pn, pb, incoming[0].type, p1, pa, pH, t2, t1 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_between( std::array const & incoming, std::vector const & outgoing, bool check_momenta ) const { 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); std::vector partons(begin(outgoing), the_Higgs); partons.insert(end(partons), the_Higgs + 1, end(outgoing)); partons = tag_extremal_jet_partons(incoming, partons, check_momenta); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector( partons[has_unob_gluon(incoming, outgoing)?1:0] ); auto pn = to_HepLorentzVector( partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)] ); auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing)); assert( (first_after_Higgs == end(partons) && ( has_unob_gluon(incoming, outgoing) || partons.back().type != pid::gluon )) || first_after_Higgs->rapidity() >= the_Higgs->rapidity() ); assert( (first_after_Higgs == begin(partons) && ( has_unof_gluon(incoming, outgoing) || partons.front().type != pid::gluon )) || (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity() ); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_Higgs == begin(partons)) ++first_after_Higgs; else if(first_after_Higgs == end(partons)) --first_after_Higgs; // t-channel momentum before Higgs auto qH = pa; for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){ qH -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if(has_unob_gluon(incoming, outgoing)){ current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(has_unof_gluon(incoming, outgoing)){ current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ current_factor = ME_Higgs_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, q0, pa, pb, p1, pn )*FKL_ladder_weight( first_after_Higgs, end_ladder, qH - pH, pa, pb, p1, pn ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double MatrixElement::tree_param_partons( double alpha_s, double mur, std::vector const & partons ) const{ const double gs2 = 4.*M_PI*alpha_s; double wt = std::pow(gs2, partons.size()); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(partons.size() >= 2); for(size_t i = 1; i < partons.size()-1; ++i){ wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp()); } } return wt; } double MatrixElement::tree_param( double mur, std::array const & incoming, std::vector const & outgoing ) const{ const double alpha_s = alpha_s_(mur); auto AWZH_boson = std::find_if( begin(outgoing), end(outgoing), [](auto const & p){return is_AWZH_boson(p);} ); double AWZH_coupling = 1.; if(AWZH_boson != end(outgoing)){ switch(AWZH_boson->type){ case pid::Higgs: { AWZH_coupling = alpha_s*alpha_s; break; } // TODO case pid::Wp:{ AWZH_coupling = alpha_w*alpha_w/2; break; } case pid::Wm:{ AWZH_coupling = alpha_w*alpha_w/2; break; } case pid::photon: case pid::Z: default: throw std::logic_error("Emission of boson of unsupported type"); } } if(has_unob_gluon(incoming, outgoing)){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)}) ); } if(has_unof_gluon(incoming, outgoing)){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1}) ); } return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(outgoing)); } double MatrixElement::tree( double mur, std::array const & incoming, std::vector const & outgoing, std::unordered_map> const & decays, bool check_momenta ) const { return tree_param(mur, incoming, outgoing)*tree_kin( incoming, outgoing, decays, check_momenta ); } } // namespace RHEJ diff --git a/src/currents.cc b/src/currents.cc index 01b94b7..66be799 100644 --- a/src/currents.cc +++ b/src/currents.cc @@ -1,3346 +1,3345 @@ ////////////////////////////////////////////////// ////////////////////////////////////////////////// // This source code is Copyright (2012) of // // Jeppe R. Andersen and Jennifer M. Smillie // // and is distributed under the // // Gnu Public License version 2 // // http://www.gnu.org/licenses/gpl-2.0.html // // You are allowed to distribute and alter the // // source under the conditions of the GPLv2 // // as long as this copyright notice // // is unaltered and distributed with the source // // Any use should comply with the // // MCNET GUIDELINES // // for Event Generator Authors and Users // // as distributed with this source code // ////////////////////////////////////////////////// ////////////////////////////////////////////////// #include "RHEJ/currents.hh" //#include "ZJets/Flags.h" #include "RHEJ/Constants.hh" #include "RHEJ/utility.hh" #include "RHEJ/PDG_codes.hh" const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4); //const double HVE = 246.21845810181637; #ifdef RHEJ_BUILD_WITH_QCDLOOP #include "qcdloop/qcdloop.h" #endif #include namespace { // Loop integrals #ifdef RHEJ_BUILD_WITH_QCDLOOP COM B0DD(CLHEP::HepLorentzVector q, double mq) { static std::vector> result(3); static auto ql_B0 = [](){ ql::Bubble,double,double> ql_B0; ql_B0.setCacheSize(100); return ql_B0; }(); static std::vector masses(2); static std::vector momenta(1); for(auto & m: masses) m = mq*mq; momenta.front() = q.m2(); ql_B0.integral(result, 1, masses, momenta); return result[0]; } COM C0DD(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { static std::vector> result(3); static auto ql_C0 = [](){ ql::Triangle,double,double> ql_C0; ql_C0.setCacheSize(100); return ql_C0; }(); static std::vector masses(3); static std::vector momenta(3); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = (q1+q2).m2(); ql_C0.integral(result, 1, masses, momenta); return result[0]; } COM D0DD(CLHEP::HepLorentzVector q1,CLHEP::HepLorentzVector q2, CLHEP::HepLorentzVector q3, double mq) { static std::vector> result(3); static auto ql_D0 = [](){ ql::Box,double,double> ql_D0; ql_D0.setCacheSize(100); return ql_D0; }(); static std::vector masses(4); static std::vector momenta(6); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = q3.m2(); momenta[3] = (q1+q2+q3).m2(); momenta[4] = (q1+q2).m2(); momenta[5] = (q2+q3).m2(); ql_D0.integral(result, 1, masses, momenta); return result[0]; } COM A1(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt) // As given in Eq. (B.2) of VDD { double q12,q22,Q2; CLHEP::HepLorentzVector Q; double Delta3,mt2; COM ans(COM(0.,0.)); q12=q1.m2(); q22=q2.m2(); Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD Q2=Q.m2(); // std::cout<<"Higgs mass? : "<pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } } CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin) { COM cur[4]; cur[0]=0.; cur[1]=0.; cur[2]=0.; cur[3]=0.; double sqpop=sqrt(pout.plus()); double sqpom=sqrt(pout.minus()); COM poperp=pout.x()+COM(0,1)*pout.y(); if (helout!=helin) { std::cerr<< "void j : Non-matching helicities\n"; } else if (helout==false) { // negative helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } CCurrent temp(cur[0],cur[1],cur[2],cur[3]); return temp; } CCurrent jio (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout) { COM cur[4]; cur[0]=0.; cur[1]=0.; cur[2]=0.; cur[3]=0.; double sqpop=sqrt(pout.plus()); double sqpom=sqrt(pout.minus()); COM poperp=pout.x()+COM(0,1)*pout.y(); if (helout!=helin) { std::cerr<< "void j : Non-matching helicities\n"; } else if (helout==false) { // negative helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } CCurrent temp(cur[0],cur[1],cur[2],cur[3]); return temp; } // Current for void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) { cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helin!=helout){ std::cout<<__LINE__<<" "<<__FILE__< pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp); cur[2] = COM(0,1) * cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = -COM(0,1) * cur[1]; cur[3] = -cur[0]; } } else { if (pin.plus() > pin.minus()) { // if forward double sqpip = sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip*poperp/abs(poperp); cur[2] = -COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * poperp/abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = COM(0,1)*cur[1]; cur[3] = -cur[0]; } } } // Current for void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) { // Zero our current cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helj){ std::cout<<__LINE__<<" "<<__FILE__< void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) { cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helin){ std::cout<<__LINE__<<" "<<__FILE__< pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip * poperp/abs(poperp); cur[2] = -COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom*sqpim*poperp/abs(poperp); cur[1] = -sqpim*sqpop; cur[2] = COM(0,1)*cur[1]; cur[3] = -cur[0]; } } else { if (pin.plus() > pin.minus()) { // if forward double sqpip = sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip*conj(poperp)/abs(poperp); cur[2] = COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * conj(poperp)/abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = -COM(0,1)*cur[1]; cur[3] = -cur[0]; } } } double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) { // std::cerr<<"Current: "< 16 pi mt^2/v alphas, // and we divide by a factor 4 at the amp sqaured level later // which I absorb here (i.e. I divide by 2) /// @TODO move factor 1/2 from S to |ME|^2 => consistent with general notation return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt)); else return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt)) + 8.*M_PI*mb*mb/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb)); } } } // namespace anonymous double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m, q1v, q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); jio (p2in,true,p2out,true,j2p); jio (p2in,false,p2out,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); jio (p2in,true,p2out,true,j2p); jio (p2in,false,p2out,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon) // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); //cout << "Bits in MH2qg: " << Mpp << " " << Mpm << " " << Mmp << " " << Mmm << endl; const double K = K_g(p2out, p2in); double sst=K/C_A*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm)); // Cf*Ca=4 // return 4.*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon) // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM amp,amm,apm,app; app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp); const double K = K_g(p2out, p2in); MH2sum*=K/C_A; // Cf*Ca=4 // return 4.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // g~p1 g~p2 // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM amp,amm,apm,app; app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp); const double K_g1 = K_g(p1out, p1in); const double K_g2 = K_g(p2out, p2in); MH2sum*=K_g1/C_A*K_g2/C_A; // Ca*Ca=9 // return 9.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } // // Z's stuff // void jZ(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) { // // Init current to zero // cur[0] = 0.0; // cur[1] = 0.0; // cur[2] = 0.0; // cur[3] = 0.0; // // Temporary variables // COM temp; // current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2; // // Momentum of virtual gluons aroun weak boson emission site // HLV qa = pout + pep + pem; // HLV qb = pin - pep - pem; // double ta = qa.m2(); // double tb = qb.m2(); // // Out-Out currents: // current Em_Ep, Out_Em, Out_Ep; // // Other currents: // current Out_In, Em_In, Ep_In; // joi(pout, HelPartons, pin, HelPartons, Out_In); // joi(pem, HelLeptons, pin, HelPartons, Em_In); // joi(pep, HelLeptons, pin, HelPartons, Ep_In); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // joo(pout, HelPartons, pem, HelLeptons, Out_Em); // joo(pout, HelPartons, pep, HelLeptons, Out_Ep); // if (HelLeptons == HelPartons) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // temp = cdot(Out_Em, Em_Ep); // cmult(temp / ta , Em_In, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // temp = -cdot(Ep_In, Em_Ep); // cmult(temp / tb, Out_Ep, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // if (HelPartons == true) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // joo(pout, true, pep, true, TempCur1); // joi(pep, true, pin, true, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // joo(pout, true, pem, true, TempCur1); // joi(pem, true, pin, true, TempCur2); // temp = -cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // joo(pout, false, pep, false, TempCur1); // joi(pep, false, pin, false, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta, TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // joo(pout, false, pem, false, TempCur1); // joi(pem, false, pin, false, TempCur2); // temp = -cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // } // } // void jZbar(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) { // // Init current to zero // cur[0] = 0.0; // cur[1] = 0.0; // cur[2] = 0.0; // cur[3] = 0.0; // // Temporary variables // COM temp; // current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2; // // Transfered 4-momenta // HLV qa = pout + pep + pem; // HLV qb = pin - pep - pem; // // The square of the transfered 4-momenta // double ta = qa.m2(); // double tb = qb.m2(); // // Out-Out currents: // current Em_Ep, Em_Out, Ep_Out; // // In-Out currents: // current In_Out, In_Em, In_Ep; // // Safe to use the currents since helicity structure is ok // if (HelPartons == HelLeptons) { // jio(pin, HelPartons, pout, HelPartons, In_Out); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // jio(pin, HelPartons, pem, HelLeptons, In_Em); // jio(pin, HelPartons, pep, HelLeptons, In_Ep); // joo(pem, HelLeptons, pout, HelPartons, Em_Out); // joo(pep, HelLeptons, pout, HelPartons, Ep_Out); // } // else { // jio(pin, HelPartons, pout, HelPartons, In_Out); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // In_Em[0] = 0.0; // In_Em[1] = 0.0; // In_Em[2] = 0.0; // In_Em[3] = 0.0; // In_Ep[0] = 0.0; // In_Ep[1] = 0.0; // In_Ep[2] = 0.0; // In_Ep[3] = 0.0; // Em_Out[0] = 0.0; // Em_Out[1] = 0.0; // Em_Out[2] = 0.0; // Em_Out[3] = 0.0; // Ep_Out[0] = 0.0; // Ep_Out[1] = 0.0; // Ep_Out[2] = 0.0; // Ep_Out[3] = 0.0; // } // if (HelLeptons == HelPartons) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // temp = cdot(Ep_Out, Em_Ep); // cmult(temp / ta, In_Ep, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // temp = - cdot(In_Em, Em_Ep); // cmult(temp / tb, Em_Out, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // if (HelPartons == true) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // joo(pem, true, pout, true, TempCur1); // jio(pin, true, pem, true, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // joo(pep, true, pout, true, TempCur1); // jio(pin, true, pep, true, TempCur2); // temp = - cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // joo(pem, false, pout, false, TempCur1); // jio(pin, false, pem, false, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // joo(pep, false, pout, false, TempCur1); // jio(pin, false, pep, false, TempCur2); // temp = - cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // } // } // // Progagators // COM PZ(double s) { // double MZ, GammaZ; // MZ = 9.118800e+01; // Mass of the mediating gauge boson // GammaZ = 2.441404e+00; // Z peak width // // Return Z Prop value // return 1.0 / (s - MZ * MZ + COM(0.0, 1.0) * GammaZ * MZ); // } // COM PG(double s) { // return 1.0 / s; // } // // Non-gluonic with pa emitting // std::vector jMZqQ (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector VProducts, std::vector < std::vector > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) { // std::vector ScaledWeights; // double Sum; // // Propagator factors // COM PZs = PZ((pep + pem).m2()); // COM PGs = PG((pep + pem).m2()); // // Emitting current initialisation // current j1pptop, j1pmtop; // Emission from top line // current j1ppbot, j1pmbot; // Emission from bottom line // // Non-emitting current initialisation // current j2ptop, j2mtop; // Emission from top line // current j2pbot, j2mbot; // Emission from bottom line // // Currents for top emission // // Upper current calculations // // if a is a quark // if (aptype > 0) { // jZ(pa, p1, pem, pep, true, true, j1pptop); // jZ(pa, p1, pem, pep, true, false, j1pmtop); // } // // if a is an antiquark // else { // jZbar(pa, p1, pem, pep, true, true, j1pptop); // jZbar(pa, p1, pem, pep, true, false, j1pmtop); // } // // Lower current calculations // // if b is a quark // if (bptype > 0) { // joi(p2, true, pb, true, j2ptop); // joi(p2, false, pb, false, j2mtop); // } // // if b is an antiquark // else { // jio(pb, true, p2, true, j2ptop); // jio(pb, false, p2, false, j2mtop); // } // // Currents for bottom emission // // Lower current calculations // if (bptype > 0) { // jZ(pb, p2, pem, pep, true, true, j1ppbot); // jZ(pb, p2, pem, pep, true, false, j1pmbot); // } // else { // jZbar(pb, p2, pem, pep, true, true, j1ppbot); // jZbar(pb, p2, pem, pep, true, false, j1pmbot); // } // // Upper current calculations // if (aptype > 0) { // joi(p1, true, pa, true, j2pbot); // joi(p1, false, pa, false, j2mbot); // } // else { // jio(pa, true, p1, true, j2pbot); // jio(pa, false, p1, false, j2mbot); // } // COM Coeff[2][8]; // if (!Interference) { // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // double ZCharge_b_P = Zq(bptype, true); // double ZCharge_b_M = Zq(bptype, false); // if (BottomLineEmit) { // // Emission from top-line quark (pa/p1 line) // Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // } // else { // // Emission from bottom-line quark (pb/p2 line) // Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot); // Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot); // Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot); // Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot); // Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot)); // Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot)); // Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot)); // Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot)); // } // } // // Else calculate all the possiblities // else { // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // double ZCharge_b_P = Zq(bptype, true); // double ZCharge_b_M = Zq(bptype, false); // // Emission from top-line quark (pa/p1 line) // Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // // Emission from bottom-line quark (pb/p2 line) // Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot); // Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot); // Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot); // Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot); // Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot)); // Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot)); // Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot)); // Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot)); // } // // Find the numbers of scales // int ScaleCount; // #if calcscaleunc // ScaleCount = 20; // #else // ScaleCount = 1; // #endif // // For each scale... // for (int j = 0; j < ScaleCount; j++) { // Sum = 0.0; // // If we want to compare back to the W's code only emit from one quark and only couple to left handed particles // // virtuals arent here since they are calculated and included in weight() call. // if (!Interference) { // if (BottomLineEmit) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[1][i]) * VProducts.at(1); // else for (int i = 0; i < 8; i++) Sum += abs2(Coeff[0][i]) * VProducts.at(0); // } // // Else work out the full interference // else { // // For the full calculation... // if (UseVirtuals) { // for (int i = 0; i < 8; i++) { // Sum += abs2(Coeff[0][i]) * VProducts.at(0) * Virtuals.at(j).at(0) // + abs2(Coeff[1][i]) * VProducts.at(1) * Virtuals.at(j).at(1) // + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2) * Virtuals.at(j).at(2); // } // } // // For the tree level calculation... // else { // for (int i = 0; i < 8; i++) { // Sum += abs2(Coeff[0][i]) * VProducts.at(0) // + abs2(Coeff[1][i]) * VProducts.at(1) // + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2); // } // } // } // // Add this to the vector to be returned with the other factors of C_A and the helicity sum/average factors. // ScaledWeights.push_back(Sum / 18.0); // } // // Return all the scale values // return ScaledWeights; // } // // Semi-gluonic with pa emitting // std::vector jMZqg (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector VProducts, std::vector < std::vector > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) { // COM Coeff[8]; // double Sum; // std::vector ScaledWeights; // COM PZs = PZ((pep + pem).m2()); // COM PGs = PG((pep + pem).m2()); // // Emitting current initialisation - Emission from top line // current j1pptop, j1pmtop; // // Non-emitting current initialisation - Emission from top line // current j2ptop, j2mtop; // // Currents for top emission // // Upper current calculations // if (aptype > 0) { // jZ (pa, p1, pem, pep, true, true, j1pptop); // jZ (pa, p1, pem, pep, true, false, j1pmtop); // } // else { // jZbar(pa, p1, pem, pep, true, true, j1pptop); // jZbar(pa, p1, pem, pep, true, false, j1pmtop); // } // // Lower current calculations // joi(p2, true, pb, true, j2ptop); // joi(p2, false, pb, false, j2mtop); // // Calculate all the possiblities // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // // Emission from top-line quark (pa/p1 line) // Coeff[0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // // Calculate gluon colour accelerated factor // double CAMFactor, z; // // If b is a forward moving gluon define z (C.F. multiple jets papers) // if (pb.pz() > 0) z = p2.plus() / pb.plus(); // else z = p2.minus() / pb.minus(); // CAMFactor = (1.0 - 1.0 / 9.0) / 2.0 * (z + 1.0 / z) + 1.0 / 9.0; // // Find the numbers of scales // int ScaleCount; // #if calcscaleunc // ScaleCount = 20; // #else // ScaleCount = 1; // #endif // // For each scale... // for (int j = 0; j < ScaleCount; j++) { // Sum = 0.0; // // If we dont want the interference // if (!Interference) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0); // // Else work out the full interference // else { // if (UseVirtuals) { // for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0) * Virtuals.at(j).at(0); // } // else { // for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0); // } // } // // Add this to the vector to be returned with the other factors of C_A, the colour accelerated factor and the helicity sum/average factors.: (4/3)*3/32 // ScaledWeights.push_back(CAMFactor * Sum / 8.0); // } // return ScaledWeights; // } // // Electroweak Charge Functions // double Zq (int PID, bool Helcitiy) { // double temp; // // Positive Spin // if (Helcitiy == true) { // if (PID == 1 || PID == 3 || PID == 5) temp = (+ 1.0 * stw2 / 3.0) / ctw; // if (PID == 2 || PID == 4) temp = (- 2.0 * stw2 / 3.0) / ctw; // if (PID == -1 || PID == -3 || PID == -5) temp = (- 1.0 * stw2 / 3.0) / ctw; // if (PID == -2 || PID == -4) temp = (+ 2.0 * stw2 / 3.0) / ctw; // // If electron or positron // if (PID == 7 || PID == -7) temp = Zep; // } // // Negative Spin // else { // if (PID == 1 || PID == 3 || PID == 5) temp = (-0.5 + 1.0 * stw2 / 3.0) / ctw; // if (PID == 2 || PID == 4) temp = ( 0.5 - 2.0 * stw2 / 3.0) / ctw; // if (PID == -1 || PID == -3 || PID == -5) temp = ( 0.5 - 1.0 * stw2 / 3.0) / ctw; // if (PID == -2 || PID == -4) temp = (-0.5 + 2.0 * stw2 / 3.0) / ctw; // // If electron or positron // if (PID == 7 || PID == -7) temp = Zem; // } // return temp; // } // double Gq (int PID) { // if (!VirtualPhoton) return 0.0; // if (PID == -1) return 1.0 * ee / 3.0; // if (PID == -2) return -2.0 * ee / 3.0; // if (PID == -3) return 1.0 * ee / 3.0; // if (PID == -4) return -2.0 * ee / 3.0; // if (PID == -5) return 1.0 * ee / 3.0; // if (PID == 1) return -1.0 * ee / 3.0; // if (PID == 2) return 2.0 * ee / 3.0; // if (PID == 3) return -1.0 * ee / 3.0; // if (PID == 4) return 2.0 * ee / 3.0; // if (PID == 5) return -1.0 * ee / 3.0; // std::cout << "ERROR! No Electroweak Charge Found at line " << __LINE__ << "..." << std::endl; // return 0.0; // } namespace { //@{ /// @brief Higgs vertex contracted with one current CCurrent jH (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j2 = j(pout,helout,pin,helin); CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz()); if(mt == infinity) return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); } } CCurrent jioH (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j2 = jio(pin,helin,pout,helout); CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz()); if(mt == infinity) return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); } } CCurrent jHtop (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j1 = j(pout,helout,pin,helin); CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz()); if(mt == infinity) return ((q1.dot(q2))*j1 - j1.dot(q2)*jq1)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); } } CCurrent jioHtop (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j1 = jio(pin,helin,pout,helout); CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz()); if(mt == infinity) return ((q1.dot(q2))*j1 - j1.dot(q2)*jq1)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); } } //@} } // namespace anonymous double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=RHEJ::C_F*RHEJ::C_F/RHEJ::C_A/RHEJ::C_A; // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; //Higgs coupling is included in Hjets.C ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. return ampsq; } double jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // here we need 2 to match with the normalization // gq is 9./4. times the qQ //Higgs coupling is included in Hjets.C const double K = K_g(p2out, p2in); return ampsq*K/C_A*9./4.; //ca/cf = 9/4 } double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // here we need 2 to match with the normalization // gq is 9./4. times the qQ //Higgs coupling is included in Hjets.C const double K = K_g(p2out, p2in); return ampsq*K/C_F; } double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // std::cout << "####################\n"; // std::cout << "# p1in : "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs const double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=RHEJ::C_F*RHEJ::C_F/(RHEJ::C_A*RHEJ::C_A); // Factor of (Cf/Ca) for each quark to match MH2qQ. return ampsq; } double jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // std::cout << "####################\n"; // std::cout << "# p1in : "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // need twice to match the normalization //Higgs coupling is included in Hjets.C const double K = K_g(p1out, p1in); return ampsq*K/C_F; } double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "< 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C const double K = K_g(p1out, p1in); return ampsq*K/C_F; //ca/cf = 9/4 } // Begin finite mass stuff #ifdef RHEJ_BUILD_WITH_QCDLOOP namespace { // All the stuff needed for the box functions in qg->qgH now... //COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2=-(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) + S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) + 2.*(s34 + S1)*(s34 + S1)/Delta + S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + 2.*(s34 + S1)/Delta + S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) - C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1 + k2, q2, mq)*2.*s34/ S1 - (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 + S1)/(S1*Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(2.*s34*(s34 + S1)*(S1 - S2)/(Delta*Sigma) + 2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S1)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } //COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-S2*D0DD(k1, k2, q2, mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) - s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 + S2)*(S2 - S1)/(Delta*Sigma) + (B0DD( k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S2)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } //COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor*(S2*D0DD(k1, q2, k2, mq)*(Delta/s12/s12 - 4.*mq*mq/s12) - S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) - S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) - C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) - C0DD(k1, q2, mq))*4.*mq*mq/ s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./ s12 + (B0DD(k1 + q2, mq) - B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1))); } //COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (-s12*D0DD(k2, k1, q2, mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) - s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S1),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta + 2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) + s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma))); } //COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (-s12*D0DD(k1, k2, q2, mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) + s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S2),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/Delta + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma))); } //COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor* (-D0DD(k1, q2, k2, mq)*(Delta/s12 + (s12 + S1)/2. - 4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD( k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD( k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*2./(s12 + S2)); } //COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta + 12.*mq*mq*S1*(s34 + S1)/Delta/Delta - 4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) + C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) - 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) + 8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s34*(4.*s12 + 3.*S1 + S2)/(Delta*Sigma) + 8.*s12*s34*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*S1*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } //COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (s12*D0DD(k1, k2, q2, mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta + 12.*mq*mq*s34*(s12 + S1)/Delta/Delta - 4.*s12*pow((s34 + S2),2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*(s34 + S2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) - C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) + 4.*s12*s34*(S2 - S1)/(Delta*Sigma) + 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) - (B0DD( k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) + 8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(-12*s34*(2*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s12*s34*s34/(S2*Delta*Delta) + 4.*s34*S1/(Delta*Sigma) - 4.*s34*(s12*s34*(2.*s12 + S2) - S1*S1*(2.*s12 + S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*(2.*s34 + S1)/(Delta*Sigma) - 8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } //COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. + 4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*4.*(s34 + S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) - B0DD(k2 + q2, mq))*4.*s34/(Delta*S2)); } //COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E1(k1,k2,k3,k4,mq)+F1(k1,k2,k3,k4,mq)+G1(k1,k2,k3,k4,mq); return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq); } //COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E4(k1,k2,k3,k4,mq)+F4(k1,k2,k3,k4,mq)+G4(k1,k2,k3,k4,mq); return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq); } //COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E10(k1,k2,k3,k4,mq)+F10(k1,k2,k3,k4,mq)+G10(k1,k2,k3,k4,mq); return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq); } //COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H1(k2,k1,k3,k4,mq); return -1.*H1(k2,k1,kh,mq); } //COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H4(k2,k1,k3,k4,mq); return -1.*H4(k2,k1,kh,mq); } //COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H10(k2,k1,k3,k4,mq); return -1.*H10(k2,k1,kh,mq); } // FL and FT functions COM FL(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { CLHEP::HepLorentzVector Q = q1 + q2; double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2); return -1./(2.*detQ2)*((2.- 3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) - B0DD(Q, mq)) + (2. - 3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) - B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() + Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD( q1, q2, mq) - 2.); } COM FT(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { CLHEP::HepLorentzVector Q = q1 + q2; double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2); return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) - 2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() - q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) - q1.dot(q2)*FL(q1, q2, mq); } CLHEP::HepLorentzVector ParityFlip(CLHEP::HepLorentzVector p) { CLHEP::HepLorentzVector flippedVector; flippedVector.setE(p.e()); flippedVector.setX(-p.x()); flippedVector.setY(-p.y()); flippedVector.setZ(-p.z()); return flippedVector; } /// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H) void g_gH_HC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns) { current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1, epsHapart2,conjepsH1part1,conjepsH1part2; COM ang1a,sqa1; const double F = 4.*mq*mq/v; // Easier to have the whole thing as current object so I can use cdot functionality. // Means I need to write pa,p1 as current objects to_current(pa, pacur); to_current(p1,p1cur); to_current(pH,pHcur); bool gluonforward = true; if(pa.z() < 0) gluonforward = false; //HEJ gauge jio(pa,false,p1,false,cura1); if(gluonforward){ // sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp) ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp(); // sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp) sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp(); } else { ang1a = sqrt(pa.minus()*p1.plus()); sqa1 = sqrt(pa.minus()*p1.plus()); } const double prop = (pa-p1-pH).m2(); cmult(-1./sqrt(2)/ang1a,cura1,conjeps1); cmult(1./sqrt(2)/sqa1,cura1,epsa); const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2(); const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2(); const COM h4 = H4(p1,-pa,pH,mq); const COM h5 = H5(p1,-pa,pH,mq); const COM h10 = H10(p1,-pa,pH,mq); const COM h12 = H12(p1,-pa,pH,mq); cmult(Fta*pa.dot(pH), epsa, epsHapart1); cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2); cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1); cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2); cadd(epsHapart1, epsHapart2, epsHa); cadd(conjepsH1part1, conjepsH1part2, conjepsH1); const COM aH1 = cdot(pHcur, cura1); current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10; if(gluonforward){ cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1); cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2); } else{ cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus()) *((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1); cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus()) *((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2); } cmult(sqrt(2.)/ang1a*aH1, epsHa, T3); cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4); cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5); cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6); cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7); cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8); cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9); cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10); current ans; for(int i=0;i<4;i++) { ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i]; } retAns[0] = F/prop*ans[0]; retAns[1] = F/prop*ans[1]; retAns[2] = F/prop*ans[2]; retAns[3] = F/prop*ans[3]; } /// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H) void g_gH_HNC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns) { const double F = 4.*mq*mq/v; COM ang1a,sqa1; current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur, p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1, conjepsH1part2; // Find here if pa, meaning the gluon, is forward or backward bool gluonforward = true; if(pa.z() < 0) gluonforward = false; jio(pa,true,p1,true,cura1); j(p1,true,pa,true,cur1a); to_current(pa,pacur); to_current(p1,p1cur); to_current(pH,pHcur); to_current(pa+p1,paplusp1cur); to_current(p1-pa,p1minuspacur); const COM aH1 = cdot(pHcur,cura1); const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a) if(gluonforward){ // sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp) ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp(); // sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp) sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp(); } else { ang1a = sqrt(pa.minus()*p1.plus()); sqa1 = sqrt(pa.minus()*p1.plus()); } const double prop = (pa-p1-pH).m2(); cmult(1./sqrt(2)/sqa1, cur1a, epsa); cmult(-1./sqrt(2)/sqa1, cura1, conjeps1); const COM phase = cdot(conjeps1, epsa); const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2(); const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2(); const COM Falpha = FT(p1-pa,pa-p1-pH,mq); const COM Fbeta = FL(p1-pa,pa-p1-pH,mq); const COM h1 = H1(p1,-pa, pH, mq); const COM h2 = H2(p1,-pa, pH, mq); const COM h4 = H4(p1,-pa, pH, mq); const COM h5 = H5(p1,-pa, pH, mq); const COM h10 = H10(p1,-pa, pH, mq); const COM h12 = H12(p1,-pa, pH, mq); cmult(Fta*pa.dot(pH), epsa, epsHapart1); cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2); cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1); cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2); cadd(epsHapart1, epsHapart2, epsHa); cadd(conjepsH1part1, conjepsH1part2, conjepsH1); current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a, T11b,T12a,T12b,T13; if(gluonforward){ cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1); cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2); } else{ cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp()) *prop/sqa1, conjepsH1, T1); cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp()) *prop/sqa1, epsHa, T2); } const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI; cmult(aH1*sqrt(2.)/sqa1, epsHa, T3); cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4); cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a); cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b); cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6); cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7); cmult(-boxdiagFact*phase*h2, pacur, T8a); cmult(boxdiagFact*phase*h1, p1cur, T8b); cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9); cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10); cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a); cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b); cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a); cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b); cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13); current ans; for(int i=0;i<4;i++) { ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i]; } retAns[0] = F/prop*ans[0]; retAns[1] = F/prop*ans[1]; retAns[2] = F/prop*ans[2]; retAns[3] = F/prop*ans[3]; } } // namespace anonymous // JDC - new amplitude with Higgs emitted close to gluon with full mt effects. Keep usual HEJ-style function call double MH2gq_outsideH(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double mq, bool includeBottom, double mq2) { current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip; current retAns,retAnsb; j(p2out,true,p2in,true,cur2bplus); j(p2out,false,p2in,false,cur2bminus); j(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip); j(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip); COM app1,app2,apm1,apm2; COM app3, app4, apm3, apm4; if(!includeBottom) { g_gH_HC(p1in,p1out,pH,mq,retAns); app1=cdot(retAns,cur2bplus); app2=cdot(retAns,cur2bminus); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); app3=cdot(retAns,cur2bplusFlip); app4=cdot(retAns,cur2bminusFlip); // And non-conserving bits g_gH_HNC(p1in,p1out,pH,mq,retAns); apm1=cdot(retAns,cur2bplus); apm2=cdot(retAns,cur2bminus); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); apm3=cdot(retAns,cur2bplusFlip); apm4=cdot(retAns,cur2bminusFlip); } else { g_gH_HC(p1in,p1out,pH,mq,retAns); g_gH_HC(p1in,p1out,pH,mq2,retAnsb); app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus); app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb); app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip); app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip); // And non-conserving bits g_gH_HNC(p1in,p1out,pH,mq,retAns); g_gH_HNC(p1in,p1out,pH,mq2,retAnsb); apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus); apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb); apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip); apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip); } - return (abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1) - + abs2(apm2) + abs2(apm3) + abs2(apm4))*12./4./(3.*8.); - // factor = 12 (colour sum) /4 (hel avg) /(3*8) (col avg) + return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1) + + abs2(apm2) + abs2(apm3) + abs2(apm4); } #endif // RHEJ_BUILD_WITH_QCDLOOP double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta double s12,p1p,p2p; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p1p=p1.plus(); p2p=p2.plus(); } else { // opposite case p1p=p1.minus(); p2p=p2.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp)); temp=temp*conj(temp); return temp.real(); } double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.23) in hep-ph/0301013 double s12,php,p1p,phm; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 php=pH.plus(); phm=pH.minus(); p1p=p1.plus(); } else { // opposite case php=pH.minus(); phm=pH.plus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=-COM(0,1)*A/(2.*s12)*(conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)-pow(conj(p3perp)+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH))))); temp=temp*conj(temp); return temp.real(); } double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.22) in hep-ph/0301013 double s12,p2p,p1p; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p2p=p2.plus(); p1p=p1.plus(); } else { // opposite case p2p=p2.minus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=A/(2.*s12)*(sqrt(p2p/p1p)*p3perp*conj(p1perp)+sqrt(p1p/p2p)*p1perp*conj(p3perp)); temp=temp*conj(temp); return temp.real(); }