diff --git a/include/HEJ/Config.hh b/include/HEJ/Config.hh index d5aa15d..a2e1151 100644 --- a/include/HEJ/Config.hh +++ b/include/HEJ/Config.hh @@ -1,254 +1,264 @@ /** \file * \brief HEJ 2 configuration parameters * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include #include "fastjet/JetDefinition.hh" #include "yaml-cpp/yaml.h" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Fraction.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/optional.hh" #include "HEJ/output_formats.hh" namespace HEJ { //! Jet parameters struct JetParameters{ fastjet::JetDefinition def; /**< Jet Definition */ double min_pt{}; /**< Minimum Jet Transverse Momentum */ }; //! Settings for scale variation struct ScaleConfig{ //! Base scale choices std::vector base; //! Factors for multiplicative scale variation std::vector factors; //! Maximum ratio between renormalisation and factorisation scale double max_ratio{}; }; //! Settings for random number generator struct RNGConfig { //! Random number generator name std::string name; //! Optional initial seed optional seed; }; //! Settings for partial unweighting struct PartialUnweightConfig { //! Number of trials for training size_t trials; //! Maximum distance in standard deviations from mean logarithmic weight double max_dev; }; /**! Possible treatments for fixed-order input events. * * The program will decide on how to treat an event based on * the value of this enumeration. */ enum class EventTreatment{ reweight, /**< Perform resummation */ keep, /**< Keep the event */ discard, /**< Discard the event */ }; //! Container to store the treatments for various event types using EventTreatMap = std::map; //! Possible setting for the event weight enum class WeightType{ weighted, //!< weighted events unweighted_resum, //!< unweighted only resummation part partially_unweighted //!< mixed weighted and unweighted }; /**! Input parameters. * * This struct handles stores all configuration parameters * needed in a HEJ 2 run. * * \internal To add a new option: * 1. Add a member to the Config struct. * 2. Inside "src/YAMLreader.cc": * - Add the option name to the "supported" Node in * get_supported_options. * - Initialise the new Config member in to_Config. * The functions set_from_yaml (for mandatory options) and * set_from_yaml_if_defined (non-mandatory) may be helpful. * 3. Add a new entry (with short description) to config.yaml * 4. Update the user documentation in "doc/Sphinx/" */ struct Config { //! %Parameters for scale variation ScaleConfig scales; //! Resummation jet properties JetParameters resummation_jets; //! Fixed-order jet properties JetParameters fixed_order_jets; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. //! Use \ref soft_pt_regulator instead. double min_extparton_pt = 0.; //! \deprecated This is equivalent to\ref soft_pt_regulator //! and will be removed in future versions. optional> max_ext_soft_pt_fraction{}; //! @brief Maximum transverse momentum fraction from soft radiation in any //! tagging jet (e.g. extremal or qqbar jet) Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; //! Number of resummation configurations to generate per fixed-order event size_t trials{}; //! Maximal number of events optional max_events; //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Event output files names and formats std::vector output; //! Parameters for random number generation RNGConfig rng; //! Map to decide what to do for different event types EventTreatMap treat; //! %Parameters for custom analysis //! @deprecated use analyses_parameters instead YAML::Node analysis_parameters; //! %Parameters for custom analyses std::vector analyses_parameters; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; //! Type of event weight e.g. (un)weighted WeightType weight_type; //! Settings for partial unweighting optional unweight_config; }; //! Configuration options for the PhaseSpacePoint class struct PhaseSpacePointConfig { PhaseSpacePointConfig() = default; PhaseSpacePointConfig( JetParameters jet_param, double min_extparton_pt = 0., Fraction soft_pt_regulator = Fraction{DEFAULT_SOFT_PT_REGULATOR} ): jet_param{std::move(jet_param)}, min_extparton_pt{min_extparton_pt}, soft_pt_regulator{std::move(soft_pt_regulator)} {} //! Properties of resummation jets JetParameters jet_param; //! Minimum transverse momentum for extremal partons //! \deprecated This will be removed in future versions. //! Use \ref soft_pt_regulator instead. double min_extparton_pt = 0.; //! \deprecated This is equivalent to\ref soft_pt_regulator //! and will be removed in future versions. optional> max_ext_soft_pt_fraction{}; //! @brief Maximum transverse momentum fraction from soft radiation in any //! tagging jet (e.g. extremal or qqbar jet) Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; }; //! Configuration options for the MatrixElement class struct MatrixElementConfig { MatrixElementConfig() = default; MatrixElementConfig( bool log_correction, HiggsCouplingSettings Higgs_coupling, EWConstants ew_parameters, + Fraction soft_pt_regulator = Fraction{DEFAULT_SOFT_PT_REGULATOR}, double regulator_lambda = CLAMBDA ): log_correction{log_correction}, Higgs_coupling{std::move(Higgs_coupling)}, ew_parameters{std::move(ew_parameters)}, + soft_pt_regulator{soft_pt_regulator}, regulator_lambda{regulator_lambda} {} //! Whether to include the logarithmic correction from \f$\alpha_s\f$ running bool log_correction{}; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; //! elector weak parameters EWConstants ew_parameters; + //! @brief Maximum transverse momentum fraction from soft radiation in any + //! tagging jet (e.g. extremal or qqbar jet) + Fraction soft_pt_regulator{ DEFAULT_SOFT_PT_REGULATOR }; //! The regulator lambda for the subtraction terms double regulator_lambda = CLAMBDA; }; //! Configuration options for the EventReweighter class struct EventReweighterConfig { //! Settings for phase space point generation PhaseSpacePointConfig psp_config; //! Settings for matrix element calculation MatrixElementConfig ME_config; //! Access properties of resummation jets JetParameters & jet_param() { return psp_config.jet_param;} //! Access properties of resummation jets (const version) JetParameters const & jet_param() const { return psp_config.jet_param;} //! Treatment of the various event types EventTreatMap treat; }; /**! Extract PhaseSpacePointConfig from Config * * \internal We do not provide a PhaseSpacePointConfig constructor from Config * so that PhaseSpacePointConfig remains an aggregate. * This faciliates writing client code (e.g. the HEJ fixed-order generator) * that creates a PhaseSpacePointConfig *without* a Config object. * * @see to_MatrixElementConfig, to_EventReweighterConfig */ inline PhaseSpacePointConfig to_PhaseSpacePointConfig(Config const & conf) { return { conf.resummation_jets, conf.min_extparton_pt, conf.max_ext_soft_pt_fraction?*conf.max_ext_soft_pt_fraction :conf.soft_pt_regulator }; } /**! Extract MatrixElementConfig from Config * * @see to_PhaseSpacePointConfig, to_EventReweighterConfig */ inline MatrixElementConfig to_MatrixElementConfig(Config const & conf) { - return {conf.log_correction, conf.Higgs_coupling, - conf.ew_parameters, conf.regulator_lambda}; + return { + conf.log_correction, + conf.Higgs_coupling, + conf.ew_parameters, + conf.soft_pt_regulator, + conf.regulator_lambda + }; } /**! Extract EventReweighterConfig from Config * * @see to_PhaseSpacePointConfig, to_MatrixElementConfig */ inline EventReweighterConfig to_EventReweighterConfig(Config const & conf) { return { to_PhaseSpacePointConfig(conf), to_MatrixElementConfig(conf), conf.treat }; } } // namespace HEJ diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index 8e78b38..3a740a3 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,2135 +1,2135 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/MatrixElement.hh" #include #include #include #include #include #include #include #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/EWConstants.hh" #include "HEJ/Event.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/Hjets.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Particle.hh" #include "HEJ/Wjets.hh" #include "HEJ/Zjets.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/jets.hh" #include "HEJ/utility.hh" namespace HEJ { double MatrixElement::omega0( double alpha_s, double mur, fastjet::PseudoJet const & q_j ) const { const double lambda = param_.regulator_lambda; const double result = - alpha_s*N_C/M_PI*std::log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; return ( 1. + alpha_s/(4.*M_PI)*BETA0*std::log(mur*mur/(q_j.perp()*lambda)) )*result; } Weights MatrixElement::operator()(Event const & event) const { std::vector tree_kin_part=tree_kin(event); std::vector virtual_part=virtual_corrections(event); if(tree_kin_part.size() != virtual_part.size()) { throw std::logic_error("tree and virtuals have different sizes"); } Weights sum = Weights{0., std::vector(event.variations().size(), 0.)}; for(size_t i=0; i tree_kin_part=tree_kin(event); double sum = 0.; for(double i : tree_kin_part) { sum += i; } return tree_param(event)*sum; } Weights MatrixElement::tree_param(Event const & event) const { if(! is_resummable(event.type())) { return Weights{0., std::vector(event.variations().size(), 0.)}; } Weights result; // only compute once for each renormalisation scale std::unordered_map known; result.central = tree_param(event, event.central().mur); known.emplace(event.central().mur, result.central); for(auto const & var: event.variations()) { const auto ME_it = known.find(var.mur); if(ME_it == end(known)) { const double wt = tree_param(event, var.mur); result.variations.emplace_back(wt); known.emplace(var.mur, wt); } else { result.variations.emplace_back(ME_it->second); } } return result; } std::vector MatrixElement::virtual_corrections(Event const & event) const { if(! is_resummable(event.type())) { return {Weights{0., std::vector(event.variations().size(), 0.)}}; } // only compute once for each renormalisation scale std::unordered_map > known_vec; std::vector central_vec=virtual_corrections(event, event.central().mur); known_vec.emplace(event.central().mur, central_vec); for(auto const & var: event.variations()) { const auto ME_it = known_vec.find(var.mur); if(ME_it == end(known_vec)) { known_vec.emplace(var.mur, virtual_corrections(event, var.mur)); } } // At this stage known_vec contains one vector of virtual corrections for each mur value // Now put this into a vector of Weights std::vector result_vec; for(size_t i=0; isecond.at(i)); } result_vec.emplace_back(result); } return result_vec; } double MatrixElement::virtual_corrections_W( Event const & event, const double mur, Particle const & WBoson ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - partons[0].p; std::size_t first_idx = 0; std::size_t last_idx = partons.size() - 1; #ifndef NDEBUG bool wc = true; #endif bool wqq = false; // With extremal qqbar or unordered gluon outside the extremal // partons then it is not part of the FKL ladder and does not // contribute to the virtual corrections. W emitted from the // most backward leg must be taken into account in t-channel if (event.type() == event_type::unob) { q -= partons[1].p; ++first_idx; if (in[0].type != partons[1].type ){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } else if (event.type() == event_type::qqbar_exb) { q -= partons[1].p; ++first_idx; if (std::abs(partons[0].type) != std::abs(partons[1].type)){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } else { if(event.type() == event_type::unof || event.type() == event_type::qqbar_exf){ --last_idx; } if (in[0].type != partons[0].type ){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } std::size_t first_idx_qqbar = last_idx; std::size_t last_idx_qqbar = last_idx; //if qqbarMid event, virtual correction do not occur between qqbar pair. if(event.type() == event_type::qqbar_mid){ const auto backquark = std::find_if( begin(partons) + 1, end(partons) - 1 , [](Particle const & s){ return (s.type != pid::gluon); } ); if(backquark == end(partons) || (backquark+1)->type==pid::gluon) return 0; if(std::abs(backquark->type) != std::abs((backquark+1)->type)) { wqq=true; #ifndef NDEBUG wc=false; #endif } last_idx = std::distance(begin(partons), backquark); first_idx_qqbar = last_idx+1; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(std::size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q)*( partons[j+1].rapidity() - partons[j].rapidity() ); q -=partons[j+1].p; } // End Loop one if (last_idx != first_idx_qqbar) q -= partons[last_idx+1].p; if (wqq) q -= WBoson.p; for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){ exponent += omega0(alpha_s, mur, q)*( partons[j+1].rapidity() - partons[j].rapidity() ); q -= partons[j+1].p; } #ifndef NDEBUG if (wc) q -= WBoson.p; assert( nearby(q, -1*pb, norm) || is_AWZH_boson(partons.back().type) || event.type() == event_type::unof || event.type() == event_type::qqbar_exf ); #endif return std::exp(exponent); } std::vector MatrixElement::virtual_corrections_Z_qq( Event const & event, const double mur, Particle const & ZBoson ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q_t = pa - partons[0].p - ZBoson.p; fastjet::PseudoJet q_b = pa - partons[0].p; size_t first_idx = 0; size_t last_idx = partons.size() - 1; // Unordered gluon does not contribute to the virtual corrections if (event.type() == event_type::unob) { // Gluon is partons[0] and is already subtracted // partons[1] is the backward quark q_t -= partons[1].p; q_b -= partons[1].p; ++first_idx; } else if (event.type() == event_type::unof) { // End sum at forward quark --last_idx; } double sum_top=0.; double sum_bot=0.; double sum_mix=0.; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ const double dy = partons[j+1].rapidity() - partons[j].rapidity(); const double tmp_top = omega0(alpha_s, mur, q_t)*dy; const double tmp_bot = omega0(alpha_s, mur, q_b)*dy; sum_top += tmp_top; sum_bot += tmp_bot; sum_mix += (tmp_top + tmp_bot) / 2.; q_t -= partons[j+1].p; q_b -= partons[j+1].p; } return {exp(sum_top), exp(sum_bot), exp(sum_mix)}; } double MatrixElement::virtual_corrections_Z_qg( Event const & event, const double mur, Particle const & ZBoson, const bool is_gq_event ) const{ auto const & in = event.incoming(); const auto partons = filter_partons(event.outgoing()); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; #endif assert(std::is_sorted(partons.begin(), partons.end(), rapidity_less{})); assert(partons.size() >= 2); assert(pa.pz() < pb.pz()); // If this is a gq event, don't subtract the Z momentum from first q fastjet::PseudoJet q = (is_gq_event ? pa - partons[0].p : pa - partons[0].p - ZBoson.p); size_t first_idx = 0; size_t last_idx = partons.size() - 1; // Unordered gluon does not contribute to the virtual corrections if (event.type() == event_type::unob) { // Gluon is partons[0] and is already subtracted // partons[1] is the backward quark q -= partons[1].p; ++first_idx; } else if (event.type() == event_type::unof) { // End sum at forward quark --last_idx; } double sum=0.; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ sum += omega0(alpha_s, mur, q)*(partons[j+1].rapidity() - partons[j].rapidity()); q -= partons[j+1].p; } return exp(sum); } std::vector MatrixElement::virtual_corrections( Event const & event, const double mur ) const{ auto const & in = event.incoming(); auto const & out = event.outgoing(); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif const auto AWZH_boson = std::find_if( begin(out), end(out), [](Particle const & p){ return is_AWZH_boson(p); } ); if(AWZH_boson != end(out) && std::abs(AWZH_boson->type) == pid::Wp){ return {virtual_corrections_W(event, mur, *AWZH_boson)}; } if(AWZH_boson != end(out) && AWZH_boson->type == pid::Z_photon_mix){ if(is_gluon(in.back().type)){ // This is a qg event return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, false)}; } if(is_gluon(in.front().type)){ // This is a gq event return {virtual_corrections_Z_qg(event, mur, *AWZH_boson, true)}; } // This is a qq event return virtual_corrections_Z_qq(event, mur, *AWZH_boson); } assert(std::is_sorted(out.begin(), out.end(), rapidity_less{})); assert(out.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - out[0].p; std::size_t first_idx = 0; std::size_t last_idx = out.size() - 1; // if there is a Higgs boson, extremal qqbar or unordered gluon // outside the extremal partons then it is not part of the FKL // ladder and does not contribute to the virtual corrections if((out.front().type == pid::Higgs) || event.type() == event_type::unob || event.type() == event_type::qqbar_exb){ q -= out[1].p; ++first_idx; } if((out.back().type == pid::Higgs) || event.type() == event_type::unof || event.type() == event_type::qqbar_exf){ --last_idx; } std::size_t first_idx_qqbar = last_idx; std::size_t last_idx_qqbar = last_idx; //if central qqbar event, virtual correction do not occur between q-qbar. if(event.type() == event_type::qqbar_mid){ const auto backquark = std::find_if( begin(out) + 1, end(out) - 1 , [](Particle const & s){ return (s.type != pid::gluon && is_parton(s.type)); } ); if(backquark == end(out) || (backquark+1)->type==pid::gluon) return {0.}; last_idx = std::distance(begin(out), backquark); first_idx_qqbar = last_idx+1; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(std::size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } if (last_idx != first_idx_qqbar) q -= out[last_idx+1].p; for(std::size_t j = first_idx_qqbar; j < last_idx_qqbar; ++j){ exponent += omega0(alpha_s, mur, q)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } assert( nearby(q, -1*pb, norm) || out.back().type == pid::Higgs || event.type() == event_type::unof || event.type() == event_type::qqbar_exf ); return {std::exp(exponent)}; } namespace { //! Lipatov vertex for partons emitted into extremal jets CLHEP::HepLorentzVector CLipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ) { const CLHEP::HepLorentzVector p5 = qav-qbv; const CLHEP::HepLorentzVector CL = -(qav+qbv) + p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2)) - p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2)); return CL; } double C2Lipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ){ const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, p1, p2); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets double C2Lipatovots( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2())); const double kperp=(qav-qbv).perp(); if (kperp>lambda) return Cls; return Cls-4./(kperp*kperp); } double C2Lipatov_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2 ) { const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, p1, p2); const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, p1, p2); return -CL_t.dot(CL_b); } double C2Lipatovots_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, p1, p2) / sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2()); const double kperp = (qav_t - qbv_t).perp(); if (kperp > lambda){ return Cls; } return Cls - 4.0 / (kperp * kperp); } CLHEP::HepLorentzVector CLipatov( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ){ const CLHEP::HepLorentzVector p5 = qav-qbv; const CLHEP::HepLorentzVector CL = -(qav+qbv) + qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2. - qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2. + ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom)) + pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom)) - pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.; return CL; } //! Lipatov vertex double C2Lipatov( // B CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ){ const CLHEP::HepLorentzVector CL = CLipatov(qav, qbv, pim, pip, pom, pop); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction double C2Lipatovots( CLHEP::HepLorentzVector const & qav, CLHEP::HepLorentzVector const & qbv, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2())); const double kperp=(qav-qbv).perp(); if (kperp>lambda) return Cls; return Cls-4./(kperp*kperp); } double C2Lipatov_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & pim, CLHEP::HepLorentzVector const & pip, CLHEP::HepLorentzVector const & pom, CLHEP::HepLorentzVector const & pop ) { const CLHEP::HepLorentzVector CL_t = CLipatov(qav_t, qbv_t, pim, pip, pom, pop); const CLHEP::HepLorentzVector CL_b = CLipatov(qav_b, qbv_b, pim, pip, pom, pop); return -CL_t.dot(CL_b); } double C2Lipatovots_Mix( CLHEP::HepLorentzVector const & qav_t, CLHEP::HepLorentzVector const & qbv_t, CLHEP::HepLorentzVector const & qav_b, CLHEP::HepLorentzVector const & qbv_b, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & p2, const double lambda ) { const double Cls = C2Lipatov_Mix(qav_t, qbv_t, qav_b, qbv_b, pa, pb, p1, p2) / sqrt(qav_t.m2() * qbv_t.m2() * qav_b.m2() * qbv_b.m2()); const double kperp = (qav_t - qbv_t).perp(); if (kperp > lambda) { return Cls; } return Cls - 4.0 / (kperp * kperp); } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pg Unordered gluon momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering * * @note The unof contribution can be calculated by reversing the argument ordering. */ double ME_uno_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ using namespace currents; assert(aptype!=pid::gluon); // aptype cannot be gluon if (bptype==pid::gluon) { if (is_quark(aptype)) return ME_unob_qg(pg,p1,pa,pn,pb); return ME_unob_qbarg(pg,p1,pa,pn,pb); } if (is_antiquark(bptype)) { if (is_quark(aptype)) return ME_unob_qQbar(pg,p1,pa,pn,pb); return ME_unob_qbarQbar(pg,p1,pa,pn,pb); } //bptype == quark if (is_quark(aptype)) return ME_unob_qQ(pg,p1,pa,pn,pb); return ME_unob_qbarQ(pg,p1,pa,pn,pb); } /** Matrix element squared for tree-level current-current scattering * @param bptype Particle b PDG ID * @param pgin Incoming gluon momentum * @param pq Quark from splitting Momentum * @param pqbar Anti-quark from splitting Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal. * @returns ME Squared for Tree-Level Current-Current Scattering * * @note The forward qqbar contribution can be calculated by reversing the * argument ordering. */ double ME_qqbar_current( ParticleID bptype, CLHEP::HepLorentzVector const & pgin, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, bool const swap_qqbar ){ using namespace currents; if (bptype==pid::gluon) { if (swap_qqbar) // pq extremal return ME_Exqqbar_qqbarg(pgin,pq,pqbar,pn,pb); // pqbar extremal return ME_Exqqbar_qbarqg(pgin,pq,pqbar,pn,pb); } // b leg quark line if (swap_qqbar) //extremal pq return ME_Exqqbar_qqbarQ(pgin,pq,pqbar,pn,pb); return ME_Exqqbar_qbarqQ(pgin,pq,pqbar,pn,pb); } /* \brief Matrix element squared for central qqbar tree-level current-current * scattering * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param nabove Number of gluons emitted before central qqbarpair * @param nbelow Number of gluons emitted after central qqbarpair * @param pa Initial state a Momentum * @param pb Initial state b Momentum * @param pq Final state qbar Momentum * @param pqbar Final state q Momentum * @param partons Vector of all outgoing partons * @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering */ double ME_qqbar_mid_current( ParticleID aptype, ParticleID bptype, int nabove, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, std::vector const & partons ){ using namespace currents; // CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards) const bool swap_qqbar=pqbar.rapidity() < pq.rapidity(); double wt=1.; if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F; if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F; return wt*ME_Cenqqbar_qq(pa, pb, partons, is_antiquark(bptype), is_antiquark(aptype), swap_qqbar, nabove); } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ using namespace currents; if (aptype==pid::gluon && bptype==pid::gluon) { return ME_gg(pn,pb,p1,pa); } if (aptype==pid::gluon && bptype!=pid::gluon) { if (is_quark(bptype)) return ME_qg(pn,pb,p1,pa); return ME_qbarg(pn,pb,p1,pa); } if (bptype==pid::gluon && aptype!=pid::gluon) { if (is_quark(aptype)) return ME_qg(p1,pa,pn,pb); return ME_qbarg(p1,pa,pn,pb); } // they are both quark if (is_quark(bptype)) { if (is_quark(aptype)) return ME_qQ(pn,pb,p1,pa); return ME_qQbar(pn,pb,p1,pa); } if (is_quark(aptype)) return ME_qQbar(p1,pa,pn,pb); return ME_qbarQbar(pn,pb,p1,pa); } /** Matrix element squared for tree-level current-current scattering With W+Jets * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_W_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; // We know it cannot be gg incoming. assert(!(aptype==pid::gluon && bptype==pid::gluon)); if (aptype==pid::gluon && bptype!=pid::gluon) { if (is_quark(bptype)) return ME_W_qg(pn,plbar,pl,pb,p1,pa,Wprop); return ME_W_qbarg(pn,plbar,pl,pb,p1,pa,Wprop); } if (bptype==pid::gluon && aptype!=pid::gluon) { if (is_quark(aptype)) return ME_W_qg(p1,plbar,pl,pa,pn,pb,Wprop); return ME_W_qbarg(p1,plbar,pl,pa,pn,pb,Wprop); } // they are both quark if (wc){ // emission off b, (first argument pbout) if (is_quark(bptype)) { if (is_quark(aptype)) return ME_W_qQ(pn,plbar,pl,pb,p1,pa,Wprop); return ME_W_qQbar(pn,plbar,pl,pb,p1,pa,Wprop); } if (is_quark(aptype)) return ME_W_qbarQ(pn,plbar,pl,pb,p1,pa,Wprop); return ME_W_qbarQbar(pn,plbar,pl,pb,p1,pa,Wprop); } // emission off a, (first argument paout) if (is_quark(aptype)) { if (is_quark(bptype)) return ME_W_qQ(p1,plbar,pl,pa,pn,pb,Wprop); return ME_W_qQbar(p1,plbar,pl,pa,pn,pb,Wprop); } // a is anti-quark if (is_quark(bptype)) return ME_W_qbarQ(p1,plbar,pl,pa,pn,pb,Wprop); return ME_W_qbarQbar(p1,plbar,pl,pa,pn,pb,Wprop); } /** Matrix element squared for backwards uno tree-level current-current * scattering With W+Jets * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param pg Unordered gluon momentum * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for unob Tree-Level Current-Current Scattering * * @note The unof contribution can be calculated by reversing the argument ordering. */ double ME_W_uno_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; // we know they are not both gluons assert(bptype != pid::gluon || aptype != pid::gluon); if (bptype == pid::gluon && aptype != pid::gluon) { // b gluon => W emission off a if (is_quark(aptype)) return ME_Wuno_qg(p1,pa,pn,pb,pg,plbar,pl,Wprop); return ME_Wuno_qbarg(p1,pa,pn,pb,pg,plbar,pl,Wprop); } // they are both quark if (wc) {// emission off b, i.e. b is first current if (is_quark(bptype)){ if (is_quark(aptype)) return ME_W_unob_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop); return ME_W_unob_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop); } if (is_quark(aptype)) return ME_W_unob_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop); return ME_W_unob_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop); } // wc == false, emission off a, i.e. a is first current if (is_quark(aptype)) { if (is_quark(bptype)) //qq return ME_Wuno_qQ(p1,pa,pn,pb,pg,plbar,pl,Wprop); //qqbar return ME_Wuno_qQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop); } // a is anti-quark if (is_quark(bptype)) //qbarq return ME_Wuno_qbarQ(p1,pa,pn,pb,pg,plbar,pl,Wprop); //qbarqbar return ME_Wuno_qbarQbar(p1,pa,pn,pb,pg,plbar,pl,Wprop); } /** \brief Matrix element squared for backward qqbar tree-level current-current * scattering With W+Jets * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pa Initial state a Momentum * @param pb Initial state b Momentum * @param pq Final state q Momentum * @param pqbar Final state qbar Momentum * @param pn Final state n Momentum * @param plbar Final state anti-lepton momentum * @param pl Final state lepton momentum * @param swap_qqbar Ordering of qqbar pair. False: pqbar extremal. * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for qqbarb Tree-Level Current-Current Scattering * * @note calculate forwards qqbar contribution by reversing argument ordering. */ double ME_W_qqbar_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const swap_qqbar, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; // With qqbar we could have 2 incoming gluons and W Emission if (aptype==pid::gluon && bptype==pid::gluon) { //a gluon, b gluon gg->qqbarWg // This will be a wqqbar emission as there is no other possible W Emission // Site. if (swap_qqbar) return ME_WExqqbar_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop); return ME_WExqqbar_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop); } assert(aptype==pid::gluon && bptype!=pid::gluon ); //a gluon => W emission off b leg or qqbar if (!wc){ // W Emitted from backwards qqbar if (swap_qqbar) return ME_WExqqbar_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop); return ME_WExqqbar_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop); } // W Must be emitted from forwards leg. return ME_W_Exqqbar_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop); } /* \brief Matrix element squared for central qqbar tree-level current-current * scattering With W+Jets * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param nabove Number of gluons emitted before central qqbarpair * @param nbelow Number of gluons emitted after central qqbarpair * @param pa Initial state a Momentum * @param pb Initial state b Momentum\ * @param pq Final state qbar Momentum * @param pqbar Final state q Momentum * @param partons Vector of all outgoing partons * @param plbar Final state anti-lepton momentum * @param pl Final state lepton momentum * @param wqq Boolean. True siginfies W boson is emitted from Central qqbar * @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false. * @returns ME Squared for qqbar_mid Tree-Level Current-Current Scattering */ double ME_W_qqbar_mid_current( ParticleID aptype, ParticleID bptype, int nabove, int nbelow, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wqq, bool const wc, ParticleProperties const & Wprop ){ using namespace currents; // CAM factors for the qqbar amps, and qqbar ordering (default, pq backwards) const bool swap_qqbar=pqbar.rapidity() < pq.rapidity(); double wt=1.; if (aptype==pid::gluon) wt*=K_g(partons.front(),pa)/C_F; if (bptype==pid::gluon) wt*=K_g(partons.back(),pb)/C_F; if(wqq) return wt*ME_WCenqqbar_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype),is_antiquark(bptype), swap_qqbar, nabove, Wprop); return wt*ME_W_Cenqqbar_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype), is_antiquark(bptype), swap_qqbar, nabove, nbelow, wc, Wprop); } /** Matrix element squared for tree-level current-current scattering With Z+Jets * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param plbar Final state positron momentum * @param pl Final state electron momentum * @param Zprop Z properties * @param stw2 Value of sin(theta_w)^2 * @param ctw Value of cos(theta_w) * @returns ME Squared for Tree-Level Current-Current Scattering */ std::vector ME_Z_current( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, ParticleProperties const & Zprop, const double stw2, const double ctw ){ using namespace currents; // we know they are not both gluons assert(!is_gluon(aptype) || !is_gluon(bptype)); if(is_anyquark(aptype) && is_gluon(bptype)){ // This is a qg event return { ME_Z_qg(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) }; } if(is_gluon(aptype) && is_anyquark(bptype)){ // This is a gq event return { ME_Z_qg(pb,pa,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) }; } assert(is_anyquark(aptype) && is_anyquark(bptype)); // This is a qq event return ME_Z_qQ(pa,pb,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw); } /** Matrix element squared for backwards uno tree-level current-current * scattering With Z+Jets * * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param pg Unordered gluon momentum * @param plbar Final state positron momentum * @param pl Final state electron momentum * @param Zprop Z properties * @param stw2 Value of sin(theta_w)^2 * @param ctw Value of cos(theta_w) * @returns ME Squared for unob Tree-Level Current-Current Scattering * * @note The unof contribution can be calculated by reversing the argument ordering. */ std::vector ME_Z_uno_current( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, ParticleProperties const & Zprop, const double stw2, const double ctw ){ using namespace currents; // we know they are not both gluons assert(!is_gluon(aptype) || !is_gluon(bptype)); if (is_anyquark(aptype) && is_gluon(bptype)) { // This is a qg event return { ME_Zuno_qg(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw) }; } if (is_gluon(aptype) && is_anyquark(bptype)) { // This is a gq event return { ME_Zuno_qg(pb,pa,pg,pn,p1,plbar,pl,bptype,aptype,Zprop,stw2,ctw) }; } assert(is_anyquark(aptype) && is_anyquark(bptype)); // This is a qq event return ME_Zuno_qQ(pa,pb,pg,p1,pn,plbar,pl,aptype,bptype,Zprop,stw2,ctw); } /** \brief Matrix element squared for tree-level current-current scattering with Higgs * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared for Tree-Level Current-Current Scattering with Higgs */ double ME_Higgs_current( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb, double vev ){ using namespace currents; if (aptype==pid::gluon && bptype==pid::gluon) // gg initial state return ME_H_gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev); if (aptype==pid::gluon&&bptype!=pid::gluon) { if (is_quark(bptype)) return ME_H_qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.; return ME_H_qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4./9.; } if (bptype==pid::gluon && aptype!=pid::gluon) { if (is_quark(aptype)) return ME_H_qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.; return ME_H_qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev)*4./9.; } // they are both quark if (is_quark(bptype)) { if (is_quark(aptype)) return ME_H_qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.); return ME_H_qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.); } if (is_quark(aptype)) return ME_H_qbarQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.); return ME_H_qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb,vev)*4.*4./(9.*9.); } /** \brief Current matrix element squared with Higgs and unordered backward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param pg Unordered back Particle Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered backward emission * * @note This function assumes unordered gluon backwards from pa-p1 current. * For unof, reverse call order */ double ME_Higgs_current_uno( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pg, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb, double vev ){ using namespace currents; if (bptype==pid::gluon && aptype!=pid::gluon) { if (is_quark(aptype)) return ME_H_unob_gQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); return ME_H_unob_gQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); } // they are both quark if (is_quark(aptype)) { if (is_quark(bptype)) return ME_H_unob_qQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); return ME_H_unob_qbarQ(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); } if (is_quark(bptype)) return ME_H_unob_qQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); return ME_H_unob_qbarQbar(pg,p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb,vev); } void validate(MatrixElementConfig const & config) { #ifndef HEJ_BUILD_WITH_QCDLOOP if(!config.Higgs_coupling.use_impact_factors) { throw std::invalid_argument{ "Invalid Higgs coupling settings.\n" "HEJ without QCDloop support can only use impact factors.\n" "Set use_impact_factors to true or recompile HEJ.\n" }; } #endif if(config.Higgs_coupling.use_impact_factors && config.Higgs_coupling.mt != std::numeric_limits::infinity()) { throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } } // namespace MatrixElement::MatrixElement( std::function alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } std::vector MatrixElement::tree_kin( Event const & ev ) const { - if(! is_resummable(ev.type())) return {0.}; + if(!ev.valid_hej_state(param_.soft_pt_regulator)) return {0.}; auto AWZH_boson = std::find_if( begin(ev.outgoing()), end(ev.outgoing()), [](Particle const & p){return is_AWZH_boson(p);} ); if(AWZH_boson == end(ev.outgoing())) return {tree_kin_jets(ev)}; switch(AWZH_boson->type){ case pid::Higgs: return {tree_kin_Higgs(ev)}; case pid::Wp: case pid::Wm: return {tree_kin_W(ev)}; case pid::Z_photon_mix: return tree_kin_Z(ev); // TODO case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } namespace { constexpr int EXTREMAL_JET_IDX = 1; constexpr int NO_EXTREMAL_JET_IDX = 0; bool treat_as_extremal(Particle const & parton){ return parton.p.user_index() == EXTREMAL_JET_IDX; } template double FKL_ladder_weight( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn, double lambda ){ double wt = 1; auto qi = q0; for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){ assert(gluon_it->type == pid::gluon); const auto g = to_HepLorentzVector(*gluon_it); const auto qip1 = qi - g; if(treat_as_extremal(*gluon_it)){ wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A; } else{ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A; } qi = qip1; } return wt; } template std::vector FKL_ladder_weight_mix( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0_t, CLHEP::HepLorentzVector const & q0_b, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn, const double lambda ){ double wt_top = 1; double wt_bot = 1; double wt_mix = 1; auto qi_t = q0_t; auto qi_b = q0_b; for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){ assert(gluon_it->type == pid::gluon); const auto g = to_HepLorentzVector(*gluon_it); const auto qip1_t = qi_t - g; const auto qip1_b = qi_b - g; if(treat_as_extremal(*gluon_it)){ wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, lambda)*C_A; wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, lambda)*C_A; wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, lambda)*C_A; } else{ wt_top *= C2Lipatovots(qip1_t, qi_t, pa, pb, p1, pn, lambda)*C_A; wt_bot *= C2Lipatovots(qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A; wt_mix *= C2Lipatovots_Mix(qip1_t, qi_t, qip1_b, qi_b, pa, pb, p1, pn, lambda)*C_A; } qi_t = qip1_t; qi_b = qip1_b; } return {wt_top, wt_bot, wt_mix}; } std::vector tag_extremal_jet_partons( Event const & ev ){ auto out_partons = filter_partons(ev.outgoing()); if(out_partons.size() == ev.jets().size()){ // no additional emissions in extremal jets, don't need to tag anything for(auto & parton: out_partons){ parton.p.set_user_index(NO_EXTREMAL_JET_IDX); } return out_partons; } auto const & jets = ev.jets(); assert(jets.size() >= 2); auto most_backward = begin(jets); auto most_forward = end(jets) - 1; // skip jets caused by unordered emission or qqbar if(ev.type() == event_type::unob || ev.type() == event_type::qqbar_exb){ assert(jets.size() >= 3); ++most_backward; } else if(ev.type() == event_type::unof || ev.type() == event_type::qqbar_exf){ assert(jets.size() >= 3); --most_forward; } const auto extremal_jet_indices = ev.particle_jet_indices( {*most_backward, *most_forward} ); assert(extremal_jet_indices.size() == out_partons.size()); for(std::size_t i = 0; i < out_partons.size(); ++i){ assert(is_parton(out_partons[i])); const int idx = (extremal_jet_indices[i]>=0)? EXTREMAL_JET_IDX: NO_EXTREMAL_JET_IDX; out_partons[i].p.set_user_index(idx); } return out_partons; } double tree_kin_jets_qqbar_mid( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, double lambda ){ CLHEP::HepLorentzVector pq; CLHEP::HepLorentzVector pqbar; const auto backmidquark = std::find_if( begin(partons)+1, end(partons)-1, [](Particle const & s){ return s.type != pid::gluon; } ); assert(backmidquark!=end(partons)-1); if (is_quark(backmidquark->type)){ pq = to_HepLorentzVector(*backmidquark); pqbar = to_HepLorentzVector(*(backmidquark+1)); } else { pqbar = to_HepLorentzVector(*backmidquark); pq = to_HepLorentzVector(*(backmidquark+1)); } auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto q0 = pa - p1; // t-channel momentum after qqbar auto qqbart = q0; const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder_1 = (backmidquark); const auto begin_ladder_2 = (backmidquark+2); const auto end_ladder = cend(partons) - 1; for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){ qqbart -= to_HepLorentzVector(*parton_it); } const int nabove = std::distance(begin_ladder, backmidquark); std::vector partonsHLV; partonsHLV.reserve(partons.size()); for (std::size_t i = 0; i != partons.size(); ++i) { partonsHLV.push_back(to_HepLorentzVector(partons[i])); } const double current_factor = ME_qqbar_mid_current( aptype, bptype, nabove, pa, pb, pq, pqbar, partonsHLV ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder_1, q0, pa, pb, p1, pn, lambda )*FKL_ladder_weight( begin_ladder_2, end_ladder, qqbart, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_jets_qqbar(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, double lambda){ const bool swap_qqbar = is_quark(*BeginPart); const auto pgin = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1))); const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0))); const auto p1 = to_HepLorentzVector(*(BeginPart)); const auto pn = to_HepLorentzVector(*(EndPart-1)); assert((BeginIn)->type==pid::gluon); // Incoming a must be gluon. const double current_factor = ME_qqbar_current( (EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_qqbar )/(4.*(N_C*N_C - 1.)); const double ladder_factor = FKL_ladder_weight( (BeginPart+2), (EndPart-1), pgin-pq-pqbar, pgin, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_jets_uno(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, double lambda ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const double current_factor = ME_uno_current( (BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa ); const double ladder_factor = FKL_ladder_weight( (BeginPart+2), (EndPart-1), pa-p1-pg, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } } // namespace double MatrixElement::tree_kin_jets(Event const & ev) const { auto const & incoming = ev.incoming(); const auto partons = tag_extremal_jet_partons(ev); if (ev.type()==event_type::FKL){ const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); return ME_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa )/(4.*(N_C*N_C - 1.))*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, pa - p1, pa, pb, p1, pn, param_.regulator_lambda ); } if (ev.type()==event_type::unordered_backward){ return tree_kin_jets_uno(incoming.begin(), incoming.end(), partons.begin(), partons.end(), param_.regulator_lambda); } if (ev.type()==event_type::unordered_forward){ return tree_kin_jets_uno(incoming.rbegin(), incoming.rend(), partons.rbegin(), partons.rend(), param_.regulator_lambda); } if (ev.type()==event_type::extremal_qqbar_backward){ return tree_kin_jets_qqbar(incoming.begin(), incoming.end(), partons.begin(), partons.end(), param_.regulator_lambda); } if (ev.type()==event_type::extremal_qqbar_forward){ return tree_kin_jets_qqbar(incoming.rbegin(), incoming.rend(), partons.rbegin(), partons.rend(), param_.regulator_lambda); } if (ev.type()==event_type::central_qqbar){ return tree_kin_jets_qqbar_mid(incoming[0].type, incoming[1].type, to_HepLorentzVector(incoming[0]), to_HepLorentzVector(incoming[1]), partons, param_.regulator_lambda); } throw std::logic_error("Cannot reweight non-resummable processes in Pure Jets"); } namespace { double tree_kin_W_FKL( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, double lambda, ParticleProperties const & Wprop ){ auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder = cend(partons) - 1; bool wc = aptype==partons[0].type; //leg b emits w auto q0 = pa - p1; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_current( aptype, bptype, pn, pb, p1, pa, plbar, pl, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_W_uno(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, double lambda, ParticleProperties const & Wprop ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); bool wc = (BeginIn)->type==(BeginPart+1)->type; //leg b emits w auto q0 = pa - p1 - pg; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_uno_current( (BeginIn)->type, (BeginIn+1)->type, pn, pb, p1, pa, pg, plbar, pl, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_W_qqbar(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, double lambda, ParticleProperties const & Wprop ){ const bool swap_qqbar=is_quark(*BeginPart); const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pq = to_HepLorentzVector(*(BeginPart+(swap_qqbar?0:1))); const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_qqbar?1:0))); const auto p1 = to_HepLorentzVector(*(BeginPart)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const bool wc = (BeginIn+1)->type!=(EndPart-1)->type; //leg b emits w auto q0 = pa - pq - pqbar; if(!wc) q0 -= pl + plbar; const double current_factor = ME_W_qqbar_current( (BeginIn)->type, (BeginIn+1)->type, pa, pb, pq, pqbar, pn, plbar, pl, swap_qqbar, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } double tree_kin_W_qqbar_mid( ParticleID aptype, ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, double lambda, ParticleProperties const & Wprop ){ CLHEP::HepLorentzVector pq; CLHEP::HepLorentzVector pqbar; const auto backmidquark = std::find_if( begin(partons)+1, end(partons)-1, [](Particle const & s){ return s.type != pid::gluon; } ); assert(backmidquark!=end(partons)-1); if (is_quark(backmidquark->type)){ pq = to_HepLorentzVector(*backmidquark); pqbar = to_HepLorentzVector(*(backmidquark+1)); } else { pqbar = to_HepLorentzVector(*backmidquark); pq = to_HepLorentzVector(*(backmidquark+1)); } auto p1 = to_HepLorentzVector(partons.front()); auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; // t-channel momentum after qqbar auto qqbart = q0; bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqbar emit W bool wc = !wqq && (aptype==partons.front().type); //leg b emits w assert(!wqq || !wc); if(wqq){ // emission from qqbar qqbart -= pl + plbar; } else if(!wc) { // emission from leg a q0 -= pl + plbar; qqbart -= pl + plbar; } const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder_1 = (backmidquark); const auto begin_ladder_2 = (backmidquark+2); const auto end_ladder = cend(partons) - 1; for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){ qqbart -= to_HepLorentzVector(*parton_it); } const int nabove = std::distance(begin_ladder, backmidquark); const int nbelow = std::distance(begin_ladder_2, end_ladder); std::vector partonsHLV; partonsHLV.reserve(partons.size()); for (std::size_t i = 0; i != partons.size(); ++i) { partonsHLV.push_back(to_HepLorentzVector(partons[i])); } const double current_factor = ME_W_qqbar_mid_current( aptype, bptype, nabove, nbelow, pa, pb, pq, pqbar, partonsHLV, plbar, pl, wqq, wc, Wprop ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder_1, q0, pa, pb, p1, pn, lambda )*FKL_ladder_weight( begin_ladder_2, end_ladder, qqbart, pa, pb, p1, pn, lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } } // namespace double MatrixElement::tree_kin_W(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); #ifndef NDEBUG // assert that there is exactly one decay corresponding to the W assert(ev.decays().size() == 1); auto const & w_boson{ std::find_if(ev.outgoing().cbegin(), ev.outgoing().cend(), [] (Particle const & p) -> bool { return std::abs(p.type) == ParticleID::Wp; }) }; assert(w_boson != ev.outgoing().cend()); assert( static_cast(ev.decays().cbegin()->first) == std::distance(ev.outgoing().cbegin(), w_boson) ); #endif // find decay products of W auto const & decay{ ev.decays().cbegin()->second }; assert(decay.size() == 2); assert( ( is_anylepton(decay.at(0)) && is_anyneutrino(decay.at(1)) ) || ( is_anylepton(decay.at(1)) && is_anyneutrino(decay.at(0)) ) ); // get lepton & neutrino CLHEP::HepLorentzVector plbar; CLHEP::HepLorentzVector pl; if (decay.at(0).type < 0){ plbar = to_HepLorentzVector(decay.at(0)); pl = to_HepLorentzVector(decay.at(1)); } else{ pl = to_HepLorentzVector(decay.at(0)); plbar = to_HepLorentzVector(decay.at(1)); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto partons = tag_extremal_jet_partons(ev); if(ev.type() == FKL){ return tree_kin_W_FKL(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == unordered_backward){ return tree_kin_W_uno(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == unordered_forward){ return tree_kin_W_uno(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == extremal_qqbar_backward){ return tree_kin_W_qqbar(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == extremal_qqbar_forward){ return tree_kin_W_qqbar(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } assert(ev.type() == central_qqbar); return tree_kin_W_qqbar_mid(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } namespace { std::vector tree_kin_Z_FKL( const ParticleID aptype, const ParticleID bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, std::vector const & partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, const double lambda, ParticleProperties const & Zprop, const double stw2, const double ctw ){ const auto p1 = to_HepLorentzVector(partons[0]); const auto pn = to_HepLorentzVector(partons[partons.size() - 1]); const auto begin_ladder = cbegin(partons) + 1; const auto end_ladder = cend(partons) - 1; const std::vector current_factor = ME_Z_current( aptype, bptype, pn, pb, p1, pa, plbar, pl, Zprop, stw2, ctw ); std::vector ladder_factor; if(is_gluon(bptype)){ // This is a qg event const auto q0 = pa-p1-plbar-pl; ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda)); } else if(is_gluon(aptype)){ // This is a gq event const auto q0 = pa-p1; ladder_factor.push_back(FKL_ladder_weight(begin_ladder, end_ladder, q0, pa, pb, p1, pn, lambda)); } else { // This is a qq event const auto q0 = pa-p1-plbar-pl; const auto q1 = pa-p1; ladder_factor=FKL_ladder_weight_mix(begin_ladder, end_ladder, q0, q1, pa, pb, p1, pn, lambda); } std::vector result; for(size_t i=0; i std::vector tree_kin_Z_uno(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, const double lambda, ParticleProperties const & Zprop, const double stw2, const double ctw){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); const ParticleID aptype = (BeginIn)->type; const ParticleID bptype = (BeginIn+1)->type; const std::vector current_factor = ME_Z_uno_current( aptype, bptype, pn, pb, p1, pa, pg, plbar, pl, Zprop, stw2, ctw ); std::vector ladder_factor; if(is_gluon(bptype)){ // This is a qg event const auto q0 = pa-pg-p1-plbar-pl; ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda)); }else if(is_gluon(aptype)){ // This is a gq event const auto q0 = pa-pg-p1; ladder_factor.push_back(FKL_ladder_weight(BeginPart+2, EndPart-1, q0, pa, pb, p1, pn, lambda)); }else{ // This is a qq event const auto q0 = pa-pg-p1-plbar-pl; const auto q1 = pa-pg-p1; ladder_factor=FKL_ladder_weight_mix(BeginPart+2, EndPart-1, q0, q1, pa, pb, p1, pn, lambda); } std::vector result; for(size_t i=0; i MatrixElement::tree_kin_Z(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); // find decay products of Z auto const & decay{ ev.decays().cbegin()->second }; assert(decay.size() == 2); assert(is_anylepton(decay.at(0)) && !is_anyneutrino(decay.at(0)) && decay.at(0).type==-decay.at(1).type); // get leptons CLHEP::HepLorentzVector plbar; CLHEP::HepLorentzVector pl; if (decay.at(0).type < 0){ plbar = to_HepLorentzVector(decay.at(0)); pl = to_HepLorentzVector(decay.at(1)); } else{ pl = to_HepLorentzVector(decay.at(0)); plbar = to_HepLorentzVector(decay.at(1)); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto partons = tag_extremal_jet_partons(ev); const double stw2 = param_.ew_parameters.sin2_tw(); const double ctw = param_.ew_parameters.cos_tw(); if(ev.type() == FKL){ return tree_kin_Z_FKL(incoming[0].type, incoming[1].type, pa, pb, partons, plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } if(ev.type() == unordered_backward){ return tree_kin_Z_uno(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } if(ev.type() == unordered_forward){ return tree_kin_Z_uno(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Zprop(), stw2, ctw); } throw std::logic_error("Can only reweight FKL or uno processes in Z+Jets"); } double MatrixElement::tree_kin_Higgs(Event const & ev) const { if(is_uno(ev.type())){ return tree_kin_Higgs_between(ev); } if(ev.outgoing().front().type == pid::Higgs){ return tree_kin_Higgs_first(ev); } if(ev.outgoing().back().type == pid::Higgs){ return tree_kin_Higgs_last(ev); } return tree_kin_Higgs_between(ev); } namespace { // Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113 #ifdef HEJ_BUILD_WITH_QCDLOOP double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ){ if(type == pid::gluon) return currents::K_g(pout, pin); return C_F; } #endif // Colour factor in strict MRK limit double K_MRK(ParticleID type) { return (type == pid::gluon)?C_A:C_F; } } // namespace double MatrixElement::MH2_forwardH( CLHEP::HepLorentzVector const & p1out, CLHEP::HepLorentzVector const & p1in, ParticleID type2, CLHEP::HepLorentzVector const & p2out, CLHEP::HepLorentzVector const & p2in, CLHEP::HepLorentzVector const & pH, double t1, double t2 ) const{ using namespace currents; ignore(p2out, p2in); const double shat = p1in.invariantMass2(p2in); const double vev = param_.ew_parameters.vev(); // gluon case #ifdef HEJ_BUILD_WITH_QCDLOOP if(!param_.Higgs_coupling.use_impact_factors){ return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*ME_Houtside_gq( p1out, p1in, p2out, p2in, pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, vev )/(4*(N_C*N_C - 1)); } #endif return K_MRK(type2)/C_A*9./2.*shat*shat*( C2gHgp(p1in,p1out,pH,vev) + C2gHgm(p1in,p1out,pH,vev) )/(t1*t2); } double MatrixElement::tree_kin_Higgs_first(Event const & ev) const { auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); assert(outgoing.front().type == pid::Higgs); if(outgoing[1].type != pid::gluon) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(ev); } const auto pH = to_HepLorentzVector(outgoing.front()); const auto partons = tag_extremal_jet_partons( ev ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); const auto q0 = pa - p1 - pH; const double t1 = q0.m2(); const double t2 = (pn - pb).m2(); return MH2_forwardH( p1, pa, incoming[1].type, pn, pb, pH, t1, t2 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn, param_.regulator_lambda ); } double MatrixElement::tree_kin_Higgs_last(Event const & ev) const { auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); assert(outgoing.back().type == pid::Higgs); if(outgoing[outgoing.size()-2].type != pid::gluon) { assert(incoming.back().type == outgoing[outgoing.size()-2].type); return tree_kin_Higgs_between(ev); } const auto pH = to_HepLorentzVector(outgoing.back()); const auto partons = tag_extremal_jet_partons( ev ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; const double t1 = q0.m2(); const double t2 = (pn + pH - pb).m2(); return MH2_forwardH( pn, pb, incoming[0].type, p1, pa, pH, t2, t1 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn, param_.regulator_lambda ); } namespace { template double tree_kin_Higgs_uno(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, CLHEP::HepLorentzVector const & qH, CLHEP::HepLorentzVector const & qHp1, double mt, bool inc_bot, double mb, double vev ){ const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pg = to_HepLorentzVector(*BeginPart); const auto p1 = to_HepLorentzVector(*(BeginPart+1)); const auto pn = to_HepLorentzVector(*(EndPart-1)); return ME_Higgs_current_uno( (BeginIn)->type, (EndIn-1)->type, pg, pn, pb, p1, pa, qH, qHp1, mt, inc_bot, mb, vev ); } } // namespace double MatrixElement::tree_kin_Higgs_between(Event const & ev) const { using namespace event_type; auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); const auto the_Higgs = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return s.type == pid::Higgs; } ); assert(the_Higgs != end(outgoing)); const auto pH = to_HepLorentzVector(*the_Higgs); const auto partons = tag_extremal_jet_partons(ev); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector( partons[(ev.type() == unob)?1:0] ); auto pn = to_HepLorentzVector( partons[partons.size() - ((ev.type() == unof)?2:1)] ); auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing)); assert( (first_after_Higgs == end(partons) && ( (ev.type() == unob) || partons.back().type != pid::gluon )) || first_after_Higgs->rapidity() >= the_Higgs->rapidity() ); assert( (first_after_Higgs == begin(partons) && ( (ev.type() == unof) || partons.front().type != pid::gluon )) || (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity() ); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_Higgs == begin(partons)) ++first_after_Higgs; else if(first_after_Higgs == end(partons)) --first_after_Higgs; // t-channel momentum before Higgs auto qH = pa; for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){ qH -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor = NAN; if(ev.type() == FKL){ current_factor = ME_Higgs_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, param_.ew_parameters.vev() ); } else if(ev.type() == unob){ current_factor = C_A*C_A/2*tree_kin_Higgs_uno( begin(incoming), end(incoming), begin(partons), end(partons), qH, qH-pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, param_.ew_parameters.vev() ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(ev.type() == unof){ current_factor = C_A*C_A/2*tree_kin_Higgs_uno( rbegin(incoming), rend(incoming), rbegin(partons), rend(partons), qH-pH, qH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb, param_.ew_parameters.vev() ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ throw std::logic_error("Can only reweight FKL or uno processes in H+Jets"); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, q0, pa, pb, p1, pn, param_.regulator_lambda )*FKL_ladder_weight( first_after_Higgs, end_ladder, qH - pH, pa, pb, p1, pn, param_.regulator_lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } namespace { double get_AWZH_coupling(Event const & ev, double alpha_s, double alpha_w) { const auto AWZH_boson = std::find_if( begin(ev.outgoing()), end(ev.outgoing()), [](auto const & p){return is_AWZH_boson(p);} ); if(AWZH_boson == end(ev.outgoing())) return 1.; switch(AWZH_boson->type){ case pid::Higgs: return alpha_s*alpha_s; case pid::Wp: case pid::Wm: case pid::Z_photon_mix: return alpha_w*alpha_w; // TODO case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } } // namespace double MatrixElement::tree_param(Event const & ev, double mur) const { assert(is_resummable(ev.type())); const auto begin_partons = ev.begin_partons(); const auto end_partons = ev.end_partons(); const auto num_partons = std::distance(begin_partons, end_partons); const double alpha_s = alpha_s_(mur); const double gs2 = 4.*M_PI*alpha_s; double res = std::pow(gs2, num_partons); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(num_partons >= 2); const auto first_emission = std::next(begin_partons); const auto last_emission = std::prev(end_partons); for(auto parton = first_emission; parton != last_emission; ++parton){ res *= 1. + alpha_s/(2.*M_PI)*BETA0*std::log(mur/parton->perp()); } } return get_AWZH_coupling(ev, alpha_s, param_.ew_parameters.alpha_w())*res; } } // namespace HEJ diff --git a/t/ME_data/config_Z_ME.yml b/t/ME_data/config_Z_ME.yml index a219afc..76067c7 100644 --- a/t/ME_data/config_Z_ME.yml +++ b/t/ME_data/config_Z_ME.yml @@ -1,37 +1,39 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: reweight central qqbar: reweight non-resummable: discard scales: H_T/2 +soft pt regulator: 1 + log correction: false vev: 246.219650794137 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.1876 width: 2.495 random generator: name: mixmax seed: 1 diff --git a/t/ME_data/config_mt.yml b/t/ME_data/config_mt.yml index ae9abd1..d1a0f39 100644 --- a/t/ME_data/config_mt.yml +++ b/t/ME_data/config_mt.yml @@ -1,42 +1,44 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: discard central qqbar: discard non-resummable: discard scales: 125 +soft pt regulator: 1 + log correction: false random generator: name: mixmax seed: 1 vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 Higgs coupling: use impact factors: false mt: 174 include bottom: false diff --git a/t/ME_data/config_mtinf.yml b/t/ME_data/config_mtinf.yml index 1fb2269..d7a67b9 100644 --- a/t/ME_data/config_mtinf.yml +++ b/t/ME_data/config_mtinf.yml @@ -1,37 +1,39 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: discard central qqbar: discard non-resummable: discard scales: 125 +soft pt regulator: 1 + log correction: false vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 random generator: name: mixmax seed: 1 diff --git a/t/ME_data/config_mtmb.yml b/t/ME_data/config_mtmb.yml index 27a003c..6aa5b1d 100644 --- a/t/ME_data/config_mtmb.yml +++ b/t/ME_data/config_mtmb.yml @@ -1,43 +1,45 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: discard central qqbar: discard non-resummable: discard scales: 125 +soft pt regulator: 1 + log correction: false random generator: name: mixmax seed: 1 vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.385 width: 2.085 Z: mass: 91.187 width: 2.495 Higgs coupling: use impact factors: false mt: 174 include bottom: true mb: 4.7 diff --git a/t/ME_data/config_pure.yml b/t/ME_data/config_pure.yml index 327f679..1bdb9b2 100644 --- a/t/ME_data/config_pure.yml +++ b/t/ME_data/config_pure.yml @@ -1,37 +1,39 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: reweight central qqbar: discard non-resummable: discard scales: 125 +soft pt regulator: 1 + log correction: false vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.419 width: 2.0476 Z: mass: 91.187 width: 2.495 random generator: name: mixmax seed: 1 diff --git a/t/ME_data/config_w_ME.yml b/t/ME_data/config_w_ME.yml index fcf064d..7a24eb4 100644 --- a/t/ME_data/config_w_ME.yml +++ b/t/ME_data/config_w_ME.yml @@ -1,37 +1,39 @@ trials: 1 resummation jets: min pt: 30 algorithm: antikt R: 0.4 fixed order jets: min pt: 30 event treatment: FKL: reweight unordered: reweight extremal qqbar: reweight central qqbar: reweight non-resummable: discard scales: 125 +soft pt regulator: 1 + log correction: false vev: 246.2196508 particle properties: Higgs: mass: 125 width: 0.004165 W: mass: 80.419 width: 2.0476 Z: mass: 91.187 width: 2.495 random generator: name: mixmax seed: 1