diff --git a/include/HEJ/Constants.hh b/include/HEJ/Constants.hh index 4e2c052..f73bb41 100644 --- a/include/HEJ/Constants.hh +++ b/include/HEJ/Constants.hh @@ -1,33 +1,33 @@ /** \file * \brief Header file defining all global constants used for HEJ * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once namespace HEJ{ /// @name QCD parameters //@{ constexpr double N_C = 3.; //!< number of Colours constexpr double C_A = N_C; //!< \f$C_A\f$ constexpr double C_F = (N_C*N_C - 1.)/(2.*N_C); //!< \f$C_F\f$ constexpr double t_f = 0.5; //!< \f$t_f\f$ constexpr double n_f = 5.; //!< number light flavours constexpr double beta0 = 11./3.*C_A - 4./3.*t_f*n_f; //!< \f$\beta_0\f$ //@} /// @name QFT parameters //@{ constexpr double vev = 246.2196508; //!< Higgs vacuum expectation value in GeV constexpr double gw = 0.653233; constexpr double MW = 80.419; // The W mass in GeV/c^2 constexpr double GammaW = 2.0476; // the W width in GeV/c^2 //@} /// @name Generation Parameters //@{ constexpr double CLAMBDA = 0.2; //!< Scale for virtual correction, \f$\lambda\f$ cf. eq. (20) in \cite Andersen:2011hs - constexpr double CMINPT = CLAMBDA; //!< minimal \f$p_t\f$ of all partons + constexpr double CMINPT = 0.2; //!< minimal \f$p_t\f$ of all partons //@} } diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh index ca5a1ee..8df28b9 100644 --- a/include/HEJ/MatrixElement.hh +++ b/include/HEJ/MatrixElement.hh @@ -1,191 +1,191 @@ /** \file * \brief Contains the MatrixElement Class * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Weights.hh" #include "HEJ/config.hh" namespace CLHEP { class HepLorentzVector; } namespace HEJ{ class Event; class Particle; //! Class to calculate the squares of matrix elements class MatrixElement{ public: /** \brief MatrixElement Constructor * @param alpha_s Function taking the renormalisation scale * and returning the strong coupling constant * @param conf General matrix element settings */ MatrixElement( std::function alpha_s, MatrixElementConfig conf ); /** * \brief squares of regulated HEJ matrix elements * @param event The event for which to calculate matrix elements * @returns The squares of HEJ matrix elements including virtual corrections * * This function returns one value for the central parameter choice * and one additional value for each entry in \ref Event.variations(). * See eq. (22) in \cite Andersen:2011hs for the definition of the squared * matrix element. * * \internal Relation to standard HEJ Met2: MatrixElement = Met2*shat^2/(pdfta*pdftb) */ Weights operator()(Event const & event) const; //! Squares of HEJ tree-level matrix elements /** * @param event The event for which to calculate matrix elements * @returns The squares of HEJ matrix elements without virtual corrections * * cf. eq. (22) in \cite Andersen:2011hs */ Weights tree(Event const & event) const; /** * \brief Virtual corrections to matrix element squares * @param event The event for which to calculate matrix elements * @returns The virtual corrections to the squares of the matrix elements * * The all order virtual corrections to LL in the MRK limit is * given by replacing 1/t in the scattering amplitude according to the * lipatov ansatz. * * cf. second-to-last line of eq. (22) in \cite Andersen:2011hs * note that indices are off by one, i.e. out[0].p corresponds to p_1 */ Weights virtual_corrections(Event const & event) const; /** * \brief Scale-dependent part of tree-level matrix element squares * @param event The event for which to calculate matrix elements * @returns The scale-dependent part of the squares of the * tree-level matrix elements * * The tree-level matrix elements factorises into a renormalisation-scale * dependent part, given by the strong coupling to some power, and a * scale-independent remainder. This function only returns the former parts * for the central scale choice and all \ref Event.variations(). * * @see tree, tree_kin */ Weights tree_param( Event const & event ) const; /** * \brief Kinematic part of tree-level matrix element squares * @param event The event for which to calculate matrix elements * @returns The kinematic part of the squares of the * tree-level matrix elements * * The tree-level matrix elements factorises into a renormalisation-scale * dependent part, given by the strong coupling to some power, and a * scale-independent remainder. This function only returns the latter part. * Since it does not depend on the parameter variations, only a single value * is returned. * * @see tree, tree_param */ double tree_kin(Event const & event) const; private: double tree_param( Event const & event, double mur ) const; double virtual_corrections_W( Event const & event, double mur, Particle const & WBoson ) const; double virtual_corrections( Event const & event, double mur ) const; //! \internal cf. last line of eq. (22) in \cite Andersen:2011hs double omega0( double alpha_s, double mur, - fastjet::PseudoJet const & q_j, double lambda + fastjet::PseudoJet const & q_j ) const; double tree_kin_jets( Event const & ev ) const; double tree_kin_W( Event const & ev ) const; double tree_kin_Higgs( Event const & ev ) const; double tree_kin_Higgs_first( Event const & ev ) const; double tree_kin_Higgs_last( Event const & ev ) const; /** * \internal * \brief Higgs inbetween extremal partons. * * Note that in the case of unordered emission, the Higgs is *always* * treated as if in between the extremal (FKL) partons, even if its * rapidity is outside the extremal parton rapidities */ double tree_kin_Higgs_between( Event const & ev ) const; double tree_param_partons( double alpha_s, double mur, std::vector const & partons ) const; std::vector in_extremal_jet_indices( std::vector const & partons ) const; std::vector tag_extremal_jet_partons( Event const & ev ) const; double MH2_forwardH( CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, pid::ParticleID type2, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double t1, double t2 ) const; std::function alpha_s_; MatrixElementConfig param_; }; } diff --git a/include/HEJ/config.hh b/include/HEJ/config.hh index 2b9d9f6..51f973c 100644 --- a/include/HEJ/config.hh +++ b/include/HEJ/config.hh @@ -1,175 +1,191 @@ /** \file * \brief HEJ 2 configuration parameters * * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #pragma once #include #include "fastjet/JetDefinition.hh" #include "yaml-cpp/yaml.h" +#include "HEJ/Constants.hh" #include "HEJ/event_types.hh" #include "HEJ/HiggsCouplingSettings.hh" #include "HEJ/optional.hh" #include "HEJ/output_formats.hh" #include "HEJ/ScaleFunction.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; }; /**! 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; /**! 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 double min_extparton_pt; //! Maximum transverse momentum fraction from soft radiation in extremal jets double max_ext_soft_pt_fraction; + //! The regulator lambda for the subtraction terms + double regulator_lambda = CLAMBDA; //! Number of resummation configurations to generate per fixed-order event int trials; //! 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 analyses YAML::Node analysis_parameters; //! Settings for effective Higgs-gluon coupling HiggsCouplingSettings Higgs_coupling; }; //! Configuration options for the PhaseSpacePoint class struct PhaseSpacePointConfig { //! Properties of resummation jets JetParameters jet_param; //! Minimum transverse momentum for extremal partons double min_extparton_pt; //! Maximum transverse momentum fraction from soft radiation in extremal jets double max_ext_soft_pt_fraction; }; //! Configuration options for the MatrixElement class struct MatrixElementConfig { + MatrixElementConfig() = default; + MatrixElementConfig( + bool log_correction, + HiggsCouplingSettings Higgs_coupling, + double regulator_lambda = CLAMBDA + ): + log_correction(log_correction), + Higgs_coupling(Higgs_coupling), + 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; + //! 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; //! Properties of resummation jets JetParameters 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 }; } /**! Extract MatrixElementConfig from Config * * @see to_PhaseSpacePointConfig, to_EventReweighterConfig */ inline MatrixElementConfig to_MatrixElementConfig(Config const & conf) { - return {conf.log_correction, conf.Higgs_coupling}; + return {conf.log_correction, conf.Higgs_coupling, 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.resummation_jets, conf.treat }; } } // namespace HEJ diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index e153997..17322b6 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,1722 +1,1757 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/MatrixElement.hh" #include #include #include #include #include #include #include #include "CLHEP/Vector/LorentzVector.h" #include "fastjet/ClusterSequence.hh" #include "HEJ/Constants.hh" #include "HEJ/currents.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/event_types.hh" #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" #include "HEJ/Particle.hh" #include "HEJ/utility.hh" namespace HEJ{ //cf. last line of eq. (22) in \ref Andersen:2011hs double MatrixElement::omega0( double alpha_s, double mur, - fastjet::PseudoJet const & q_j, double lambda + fastjet::PseudoJet const & q_j ) const { + const double lambda = param_.regulator_lambda; const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; // use alpha_s(sqrt(q_j*lambda)), evolved to mur return ( 1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda)) )*result; } Weights MatrixElement::operator()( Event const & event ) const { return tree(event)*virtual_corrections(event); } Weights MatrixElement::tree( Event const & event ) const { return tree_param(event)*tree_kin(event); } Weights MatrixElement::tree_param( Event const & event ) const { if(! is_HEJ(event.type())) { return Weights{0., std::vector(event.variations().size(), 0.)}; } Weights result; // only compute once for each renormalisation scale std::unordered_map known; result.central = tree_param(event, event.central().mur); known.emplace(event.central().mur, result.central); for(auto const & var: event.variations()) { const auto ME_it = known.find(var.mur); if(ME_it == end(known)) { const double wt = tree_param(event, var.mur); result.variations.emplace_back(wt); known.emplace(var.mur, wt); } else { result.variations.emplace_back(ME_it->second); } } return result; } Weights MatrixElement::virtual_corrections( Event const & event ) const { if(! is_HEJ(event.type())) { return Weights{0., std::vector(event.variations().size(), 0.)}; } Weights result; // only compute once for each renormalisation scale std::unordered_map known; result.central = virtual_corrections(event, event.central().mur); known.emplace(event.central().mur, result.central); for(auto const & var: event.variations()) { const auto ME_it = known.find(var.mur); if(ME_it == end(known)) { const double wt = virtual_corrections(event, var.mur); result.variations.emplace_back(wt); known.emplace(var.mur, wt); } else { result.variations.emplace_back(ME_it->second); } } return result; } double MatrixElement::virtual_corrections_W( Event const & event, 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; size_t first_idx = 0; size_t last_idx = partons.size() - 1; bool wc = true; bool wqq = false; // With extremal qqx 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::FKL) { if (in[0].type != partons[0].type ){ q -= WBoson.p; wc = false; } } else if (event.type() == event_type::unob) { q -= partons[1].p; ++first_idx; if (in[0].type != partons[1].type ){ q -= WBoson.p; wc = false; } } else if (event.type() == event_type::qqxexb) { q -= partons[1].p; ++first_idx; if (abs(partons[0].type) != abs(partons[1].type)){ q -= WBoson.p; wc = false; } } if(event.type() == event_type::unof || event.type() == event_type::qqxexf){ --last_idx; } size_t first_idx_qqx = last_idx; size_t last_idx_qqx = last_idx; //if qqxMid event, virtual correction do not occur between //qqx pair. if(event.type() == event_type::qqxmid){ 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(abs(backquark->type) != abs((backquark+1)->type)) { wqq=true; wc=false; } last_idx = std::distance(begin(partons), backquark); first_idx_qqx = last_idx+1; } 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)*( + 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_qqx) q -= partons[last_idx+1].p; if (wqq) q -= WBoson.p; for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){ - exponent += omega0(alpha_s, mur, q, CLAMBDA)*( + exponent += omega0(alpha_s, mur, q)*( partons[j+1].rapidity() - partons[j].rapidity() ); q -= partons[j+1].p; } 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::qqxexf ); return exp(exponent); } double MatrixElement::virtual_corrections( Event const & event, double mur ) const{ auto const & in = event.incoming(); auto const & out = event.outgoing(); fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif const auto AWZH_boson = std::find_if( begin(out), end(out), [](Particle const & p){ return is_AWZH_boson(p); } ); if(AWZH_boson != end(out) && abs(AWZH_boson->type) == pid::Wp){ return virtual_corrections_W(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; size_t first_idx = 0; size_t last_idx = out.size() - 1; // if there is a Higgs boson, extremal qqx 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::qqxexb){ q -= out[1].p; ++first_idx; } if((out.back().type == pid::Higgs) || event.type() == event_type::unof || event.type() == event_type::qqxexf){ --last_idx; } size_t first_idx_qqx = last_idx; size_t last_idx_qqx = last_idx; //if qqxMid event, virtual correction do not occur between //qqx pair. if(event.type() == event_type::qqxmid){ 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_qqx = last_idx+1; } 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)*( + exponent += omega0(alpha_s, mur, q)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } if (last_idx != first_idx_qqx) q -= out[last_idx+1].p; for(size_t j = first_idx_qqx; j < last_idx_qqx; ++j){ - exponent += omega0(alpha_s, mur, q, CLAMBDA)*( + 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::qqxexf ); return exp(exponent); } } // namespace HEJ namespace { //! Lipatov vertex for partons emitted into extremal jets double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2)) - p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2)); return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets - double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, - CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) - { + double C2Lipatovots( + CLHEP::HepLorentzVector qav, + CLHEP::HepLorentzVector qbv, + CLHEP::HepLorentzVector p1, + CLHEP::HepLorentzVector p2, + double lambda + ) { double kperp=(qav-qbv).perp(); - if (kperp>HEJ::CLAMBDA) + if (kperp>lambda) 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 C2Lipatovots( + CLHEP::HepLorentzVector qav, + CLHEP::HepLorentzVector qbv, + CLHEP::HepLorentzVector pa, + CLHEP::HepLorentzVector pb, + CLHEP::HepLorentzVector p1, + CLHEP::HepLorentzVector p2, + double lambda + ) { double kperp=(qav-qbv).perp(); - if (kperp>HEJ::CLAMBDA) + if (kperp>lambda) 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 * @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( 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, bool const wc ){ // We know it cannot be gg incoming. assert(!(aptype==21 && bptype==21)); if (aptype==21&&bptype!=21) { if (bptype > 0) return jMWqg(pn,plbar,pl,pb,p1,pa); else return jMWqbarg(pn,plbar,pl,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jMWqg(p1,plbar,pl,pa,pn,pb); else return jMWqbarg(p1,plbar,pl,pa,pn,pb); } else { // they are both quark if (wc==true){ // emission off b, (first argument pbout) if (bptype>0) { if (aptype>0) return jMWqQ(pn,plbar,pl,pb,p1,pa); else return jMWqQbar(pn,plbar,pl,pb,p1,pa); } else { if (aptype>0) return jMWqbarQ(pn,plbar,pl,pb,p1,pa); else return jMWqbarQbar(pn,plbar,pl,pb,p1,pa); } } else{ // emission off a, (first argument paout) if (aptype > 0) { if (bptype > 0) return jMWqQ(p1,plbar,pl,pa,pn,pb); else return jMWqQbar(p1,plbar,pl,pa,pn,pb); } else { // a is anti-quark if (bptype > 0) return jMWqbarQ(p1,plbar,pl,pa,pn,pb); else return jMWqbarQbar(p1,plbar,pl,pa,pn,pb); } } } throw std::logic_error("unknown particle types"); } /** 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 */ double ME_W_unob_current( int aptype, int 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 ){ // we know they are not both gluons if (bptype == 21 && aptype != 21) { // b gluon => W emission off a if (aptype > 0) return jM2Wunogqg(pg,p1,plbar,pl,pa,pn,pb); else return jM2Wunogqbarg(pg,p1,plbar,pl,pa,pn,pb); } else { // they are both quark if (wc==true) {// emission off b, i.e. b is first current if (bptype>0){ if (aptype>0) return junobMWqQg(pn,plbar,pl,pb,p1,pa,pg); else return junobMWqQbarg(pn,plbar,pl,pb,p1,pa,pg); } else{ if (aptype>0) return junobMWqbarQg(pn,plbar,pl,pb,p1,pa,pg); else return junobMWqbarQbarg(pn,plbar,pl,pb,p1,pa,pg); } } else {// wc == false, emission off a, i.e. a is first current if (aptype > 0) { if (bptype > 0) //qq return jM2WunogqQ(pg,p1,plbar,pl,pa,pn,pb); else //qqbar return jM2WunogqQbar(pg,p1,plbar,pl,pa,pn,pb); } else { // a is anti-quark if (bptype > 0) //qbarq return jM2WunogqbarQ(pg,p1,plbar,pl,pa,pn,pb); else //qbarqbar return jM2WunogqbarQbar(pg,p1,plbar,pl,pa,pn,pb); } } } } /** Matrix element squared for uno forward 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 unof Tree-Level Current-Current Scattering */ double ME_W_unof_current( int aptype, int 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 ){ // we know they are not both gluons if (aptype==21 && bptype!=21) {//a gluon => W emission off b if (bptype > 0) return jM2Wunogqg(pg, pn,plbar, pl, pb, p1, pa); else return jM2Wunogqbarg(pg, pn,plbar, pl, pb, p1, pa); } else { // they are both quark if (wc==true) {// emission off b, i.e. b is first current if (bptype>0){ if (aptype>0) return jM2WunogqQ(pg,pn,plbar,pl,pb,p1,pa); else return jM2WunogqQbar(pg,pn,plbar,pl,pb,p1,pa); } else{ if (aptype>0) return jM2WunogqbarQ(pg,pn,plbar,pl,pb,p1,pa); else return jM2WunogqbarQbar(pg,pn,plbar,pl,pb,p1,pa); } } else {// wc == false, emission off a, i.e. a is first current if (aptype > 0) { if (bptype > 0) //qq return junofMWgqQ(pg,pn,pb,p1,plbar,pl,pa); else //qqbar return junofMWgqQbar(pg,pn,pb,p1,plbar,pl,pa); } else { // a is anti-quark if (bptype > 0) //qbarq return junofMWgqbarQ(pg,pn,pb,p1,plbar,pl,pa); else //qbarqbar return junofMWgqbarQbar(pg,pn,pb,p1,plbar,pl,pa); } } } } /** \brief Matrix element squared for backward qqx 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 wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for qqxb Tree-Level Current-Current Scattering */ double ME_W_qqxb_current( int aptype, int 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 wc ){ // CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal) bool swapQuarkAntiquark=false; double CFbackward; if (pqbar.rapidity() > pq.rapidity()){ swapQuarkAntiquark=true; CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pq.minus())+(pq.minus())/pa.minus())+1./3.)*3./4.; } else{ CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(pqbar.minus())+(pqbar.minus())/pa.minus())+1./3.)*3./4.; } // With qqbar we could have 2 incoming gluons and W Emission if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg // This will be a wqqx emission as there is no other possible W Emission Site. if (swapQuarkAntiquark){ return jM2Wggtoqqbarg(pa, pqbar, plbar, pl, pq, pn,pb)*CFbackward;} else { return jM2Wggtoqbarqg(pa, pq, plbar, pl, pqbar, pn,pb)*CFbackward;} } else if (aptype==21&&bptype!=21 ) {//a gluon => W emission off b leg or qqx if (wc!=1){ // W Emitted from backwards qqx if (swapQuarkAntiquark){ return jM2WgQtoqqbarQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;} else{ return jM2WgQtoqbarqQ(pa, pq, plbar, pl, pqbar, pn, pb)*CFbackward;} } else { // W Must be emitted from forwards leg. if(bptype > 0){ if (swapQuarkAntiquark){ return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, false)*CFbackward;} else{ return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, false)*CFbackward;} } else { if (swapQuarkAntiquark){ return jM2WgqtoQQqW(pb, pa, pn, pqbar, pq, plbar, pl, true)*CFbackward;} else{ return jM2WgqtoQQqW(pb, pa, pn, pq, pqbar, plbar, pl, true)*CFbackward;} } } } else{ throw std::logic_error("Incompatible incoming particle types with qqxb"); } } /* \brief Matrix element squared for forward qqx 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 p1 Final state 1 Momentum * @param plbar Final state anti-lepton momentum * @param pl Final state lepton momentum * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for qqxf Tree-Level Current-Current Scattering */ double ME_W_qqxf_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wc ){ // CAM factors for the qqx amps, and qqbar ordering (default, qbar extremal) bool swapQuarkAntiquark=false; double CFforward; if (pqbar.rapidity() < pq.rapidity()){ swapQuarkAntiquark=true; CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pq.plus())+(pq.plus())/pb.plus())+1./3.)*3./4.; } else{ CFforward = (0.5*(3.-1./3.)*(pb.plus()/(pqbar.plus())+(pqbar.plus())/pb.plus())+1./3.)*3./4.; } // With qqbar we could have 2 incoming gluons and W Emission if (aptype==21&&bptype==21) {//a gluon, b gluon gg->qqbarWg // This will be a wqqx emission as there is no other possible W Emission Site. if (swapQuarkAntiquark){ return jM2Wggtoqqbarg(pb, pqbar, plbar, pl, pq, p1,pa)*CFforward;} else { return jM2Wggtoqbarqg(pb, pq, plbar, pl, pqbar, p1,pa)*CFforward;} } else if (bptype==21&&aptype!=21) {// b gluon => W emission off a or qqx if (wc==1){ // W Emitted from forwards qqx if (swapQuarkAntiquark){ return jM2WgQtoqbarqQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;} else { return jM2WgQtoqqbarQ(pb, pq, plbar,pl, pqbar, p1, pa)*CFforward;} } // W Must be emitted from backwards leg. if (aptype > 0){ if (swapQuarkAntiquark){ return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, false)*CFforward;} else{ return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, false)*CFforward;} } else { if (swapQuarkAntiquark){ return jM2WgqtoQQqW(pa,pb, p1, pqbar, pq, plbar, pl, true)*CFforward;} else{ return jM2WgqtoQQqW(pa,pb, p1, pq, pqbar, plbar, pl, true)*CFforward;} } } else{ throw std::logic_error("Incompatible incoming particle types with qqxf"); } } /* \brief Matrix element squared for central qqx 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 qqxpair * @param nbelow Number of gluons emitted after central qqxpair * @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 qqx * @param wc Boolean. wc=true signifies w boson emitted from leg b; if wqq=false. * @returns ME Squared for qqxmid Tree-Level Current-Current Scattering */ double ME_W_qqxmid_current( int aptype, int bptype, int nabove, int nbelow, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & pq, CLHEP::HepLorentzVector const & pqbar, std::vector partons, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl, bool const wqq, bool const wc ){ // CAM factors for the qqx amps, and qqbar ordering (default, pq backwards) bool swapQuarkAntiquark=false; if (pqbar.rapidity() < pq.rapidity()){ swapQuarkAntiquark=true; } double CFforward = (0.5*(3.-1./3.)*(pb.plus()/(partons[partons.size()-1].plus())+(partons[partons.size()-1].plus())/pb.plus())+1./3.)*3./4.; double CFbackward = (0.5*(3.-1./3.)*(pa.minus()/(partons[0].minus())+(partons[0].minus())/pa.minus())+1./3.)*3./4.; double wt=1.; if (aptype==21) wt*=CFbackward; if (bptype==21) wt*=CFforward; if (aptype <=0 && bptype <=0){ // Both External AntiQuark if (wqq==1){//emission from central qqbar return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons,true,true, swapQuarkAntiquark, nabove); } else if (wc==1){//emission from b leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, true); } else { // emission from a leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true,true, swapQuarkAntiquark, nabove, nbelow, false); } } // end both antiquark else if (aptype<=0){ // a is antiquark if (wqq==1){//emission from central qqbar return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove); } else if (wc==1){//emission from b leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons,false,true, swapQuarkAntiquark, nabove, nbelow, true); } else { // emission from a leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, true, swapQuarkAntiquark, nabove, nbelow, false); } } // end a is antiquark else if (bptype<=0){ // b is antiquark if (wqq==1){//emission from central qqbar return wt*jM2WqqtoqQQq(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove); } else if (wc==1){//emission from b leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, true); } else { // emission from a leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, true, false, swapQuarkAntiquark, nabove, nbelow, false); } } //end b is antiquark else{ //Both Quark or gluon if (wqq==1){//emission from central qqbar return wt*jM2WqqtoqQQq(pa, pb, pl, plbar, partons, false, false, swapQuarkAntiquark, nabove);} else if (wc==1){//emission from b leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, true); } else { // emission from a leg return wt*jM2WqqtoqQQqW(pa, pb, pl,plbar, partons, false, false, swapQuarkAntiquark, nabove, nbelow, false); } } } /** \brief Matrix element squared for tree-level current-current scattering with Higgs * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared for Tree-Level Current-Current Scattering with Higgs */ double ME_Higgs_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype==21) // gg initial state return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else if (aptype==21&&bptype!=21) { if (bptype > 0) return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; else return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; } else if (bptype==21&&aptype!=21) { if (aptype > 0) return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; else return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; } else { // they are both quark if (bptype>0) { if (aptype>0) return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } else { if (aptype>0) return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered forward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param punof Unordered Particle Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered forward emission */ double ME_Higgs_current_unof( int aptype, int bptype, CLHEP::HepLorentzVector const & punof, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (aptype>0) return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered backward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param punob Unordered back Particle Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered backward emission */ double ME_Higgs_current_unob( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & punob, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (bptype==21&&aptype!=21) { if (aptype > 0) return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (aptype>0) { if (bptype>0) return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (bptype>0) return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } CLHEP::HepLorentzVector to_HepLorentzVector(HEJ::Particle const & particle){ return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()}; } void validate(HEJ::MatrixElementConfig const & config) { #ifndef HEJ_BUILD_WITH_QCDLOOP if(!config.Higgs_coupling.use_impact_factors) { throw std::invalid_argument{ "Invalid Higgs coupling settings.\n" "HEJ without QCDloop support can only use impact factors.\n" "Set use_impact_factors to true or recompile HEJ.\n" }; } #endif if(config.Higgs_coupling.use_impact_factors && config.Higgs_coupling.mt != std::numeric_limits::infinity()) { throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } } // namespace anonymous namespace HEJ{ MatrixElement::MatrixElement( std::function alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } double MatrixElement::tree_kin( Event const & ev ) const { if(! is_HEJ(ev.type())) return 0.; auto AWZH_boson = std::find_if( begin(ev.outgoing()), end(ev.outgoing()), [](Particle const & p){return is_AWZH_boson(p);} ); if(AWZH_boson == end(ev.outgoing())){ return tree_kin_jets(ev); } switch(AWZH_boson->type){ case pid::Higgs: { return tree_kin_Higgs(ev); } case pid::Wp: case pid::Wm: return tree_kin_W(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 + 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)*C_A; + wt *= C2Lipatovots(qip1, qi, pa, pb, lambda)*C_A; } else{ - wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A; + wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn, lambda)*C_A; } qi = qip1; } return wt; } } // namespace anonymous std::vector MatrixElement::tag_extremal_jet_partons( Event const & ev ) const{ auto out_partons = filter_partons(ev.outgoing()); if(out_partons.size() == ev.jets().size()){ // no additional emissions in extremal jets, don't need to tag anything for(auto & parton: out_partons){ parton.p.set_user_index(no_extremal_jet_idx); } return out_partons; } // TODO: avoid reclustering fastjet::ClusterSequence cs(to_PseudoJet(out_partons), ev.jet_def()); const auto jets = sorted_by_rapidity(cs.inclusive_jets(ev.min_jet_pt())); assert(jets.size() >= 2); auto most_backward = begin(jets); auto most_forward = end(jets) - 1; // skip jets caused by unordered emission or qqx if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){ assert(jets.size() >= 3); ++most_backward; } else if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){ assert(jets.size() >= 3); --most_forward; } const auto extremal_jet_indices = cs.particle_jet_indices( {*most_backward, *most_forward} ); assert(extremal_jet_indices.size() == out_partons.size()); for(size_t i = 0; i < out_partons.size(); ++i){ assert(HEJ::is_parton(out_partons[i])); const int idx = (extremal_jet_indices[i]>=0)? extremal_jet_idx: no_extremal_jet_idx; out_partons[i].p.set_user_index(idx); } return out_partons; } double MatrixElement::tree_kin_jets( Event const & ev ) const { auto const & incoming = ev.incoming(); const auto partons = tag_extremal_jet_partons(ev); if(is_uno(ev.type())){ throw not_implemented("unordered emission not implemented for pure jets"); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); return ME_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa )/(4*(N_C*N_C - 1))*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, - pa - p1, pa, pb, p1, pn + pa - p1, pa, pb, p1, pn, + param_.regulator_lambda ); } namespace{ double tree_kin_W_FKL( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; bool wc = true; auto q0 = pa - p1; if (aptype!=partons[0].type) { //leg a emits w wc = false; q0 -=pl + plbar; } const double current_factor = ME_W_current( aptype, bptype, pn, pb, p1, pa, plbar, pl, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda ); return current_factor*ladder_factor; } double tree_kin_W_unob( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { auto pg = to_HepLorentzVector(partons[0]); auto p1 = to_HepLorentzVector(partons[1]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto begin_ladder = begin(partons) + 2; auto end_ladder = end(partons) - 1; bool wc = true; auto q0 = pa - p1 -pg; if (aptype!=partons[1].type) { //leg a emits w wc = false; q0 -=pl + plbar; } const double current_factor = ME_W_unob_current( aptype, bptype, pn, pb, p1, pa, pg, plbar, pl, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double tree_kin_W_unof( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 2]); auto pg = to_HepLorentzVector(partons[partons.size() - 1]); auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 2; bool wc = true; auto q0 = pa - p1; if (aptype!=partons[0].type) { //leg a emits w wc = false; q0 -=pl + plbar; } const double current_factor = ME_W_unof_current( aptype, bptype, pn, pb, p1, pa, pg, plbar, pl, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double tree_kin_W_qqxb( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { HLV pq,pqbar; if(is_quark(partons[0])){ pq = to_HepLorentzVector(partons[0]); pqbar = to_HepLorentzVector(partons[1]); } else{ pq = to_HepLorentzVector(partons[1]); pqbar = to_HepLorentzVector(partons[0]); } auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto begin_ladder = begin(partons) + 2; auto end_ladder = end(partons) - 1; bool wc = true; auto q0 = pa - pq - pqbar; if (partons[1].type!=partons[0].type) { //leg a emits w wc = false; q0 -=pl + plbar; } const double current_factor = ME_W_qqxb_current( aptype, bptype, pa, pb, pq, pqbar, pn, plbar, pl, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double tree_kin_W_qqxf( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { HLV pq,pqbar; if(is_quark(partons[partons.size() - 1])){ pq = to_HepLorentzVector(partons[partons.size() - 1]); pqbar = to_HepLorentzVector(partons[partons.size() - 2]); } else{ pq = to_HepLorentzVector(partons[partons.size() - 2]); pqbar = to_HepLorentzVector(partons[partons.size() - 1]); } auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 2; bool wc = true; auto q0 = pa - p1; if (aptype!=partons[0].type) { //leg a emits w wc = false; q0 -=pl + plbar; } const double current_factor = ME_W_qqxf_current( aptype, bptype, pa, pb, pq, pqbar, p1, plbar, pl, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double tree_kin_W_qqxmid( int aptype, int bptype, HLV pa, HLV pb, std::vector const & partons, - HLV plbar, HLV pl + HLV plbar, HLV pl, + double lambda ) { HLV pq,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 qqx auto qqxt = q0; bool wc, wqq; if (backmidquark->type == -(backmidquark+1)->type){ // Central qqx does not emit wqq=false; if (aptype==partons[0].type) { wc = true; } else{ wc = false; q0-=pl+plbar; } } else{ wqq = true; wc = false; qqxt-=pl+plbar; } auto begin_ladder = begin(partons) + 1; auto end_ladder_1 = (backmidquark); auto begin_ladder_2 = (backmidquark+2); auto end_ladder = end(partons) - 1; for(auto parton_it = begin_ladder; parton_it < begin_ladder_2; ++parton_it){ qqxt -= to_HepLorentzVector(*parton_it); } int nabove = std::distance(begin_ladder, backmidquark); int nbelow = std::distance(begin_ladder_2, end_ladder); std::vector partonsHLV; partonsHLV.reserve(partons.size()); for (size_t i = 0; i != partons.size(); ++i) { partonsHLV.push_back(to_HepLorentzVector(partons[i])); } const double current_factor = ME_W_qqxmid_current( aptype, bptype, nabove, nbelow, pa, pb, pq, pqbar, partonsHLV, plbar, pl, wqq, wc ); const double ladder_factor = FKL_ladder_weight( begin_ladder, end_ladder_1, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + lambda )*FKL_ladder_weight( begin_ladder_2, end_ladder, - qqxt, pa, pb, p1, pn + qqxt, pa, pb, p1, pn, + lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } } // namespace anonymous double MatrixElement::tree_kin_W(Event const & ev) const { using namespace event_type; auto const & incoming(ev.incoming()); auto const & decays(ev.decays()); 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 pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto partons = tag_extremal_jet_partons(ev); if(ev.type() == unordered_backward){ return tree_kin_W_unob(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } if(ev.type() == unordered_forward){ return tree_kin_W_unof(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } if(ev.type() == extremal_qqxb){ return tree_kin_W_qqxb(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } if(ev.type() == extremal_qqxf){ return tree_kin_W_qqxf(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } if(ev.type() == central_qqx){ return tree_kin_W_qqxmid(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } return tree_kin_W_FKL(incoming[0].type, incoming[1].type, - pa, pb, partons, plbar, pl); + pa, pb, partons, plbar, pl, + param_.regulator_lambda); } double MatrixElement::tree_kin_Higgs( Event const & ev ) const { if(is_uno(ev.type())){ return tree_kin_Higgs_between(ev); } if(ev.outgoing().front().type == pid::Higgs){ return tree_kin_Higgs_first(ev); } if(ev.outgoing().back().type == pid::Higgs){ return tree_kin_Higgs_last(ev); } return tree_kin_Higgs_between(ev); } namespace { // Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113 #ifdef HEJ_BUILD_WITH_QCDLOOP // TODO: code duplication with currents.cc double K_g(double p1minus, double paminus) { return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; } double K_g( CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(type == ParticleID::gluon) return K_g(pout, pin); return C_F; } #endif // Colour factor in strict MRK limit double K_MRK(ParticleID type) { return (type == ParticleID::gluon)?C_A:C_F; } } double MatrixElement::MH2_forwardH( CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, ParticleID type2, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double t1, double t2 ) const{ ignore(p2out, p2in); const double shat = p1in.invariantMass2(p2in); // gluon case #ifdef HEJ_BUILD_WITH_QCDLOOP if(!param_.Higgs_coupling.use_impact_factors){ return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH( p1out, p1in, p2out, p2in, pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb )/(4*(N_C*N_C - 1)); } #endif return K_MRK(type2)/C_A*9./2.*shat*shat*( C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH) )/(t1*t2); } double MatrixElement::tree_kin_Higgs_first( Event const & ev ) const { auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); assert(outgoing.front().type == pid::Higgs); if(outgoing[1].type != pid::gluon) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(ev); } const auto pH = to_HepLorentzVector(outgoing.front()); const auto partons = tag_extremal_jet_partons( ev ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); const auto q0 = pa - p1 - pH; const double t1 = q0.m2(); const double t2 = (pn - pb).m2(); return MH2_forwardH( p1, pa, incoming[1].type, pn, pb, pH, t1, t2 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, - q0, pa, pb, p1, pn + 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 + q0, pa, pb, p1, pn, + param_.regulator_lambda ); } double MatrixElement::tree_kin_Higgs_between( Event const & ev ) const { using namespace event_type; auto const & incoming = ev.incoming(); auto const & outgoing = ev.outgoing(); const auto the_Higgs = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return s.type == pid::Higgs; } ); assert(the_Higgs != end(outgoing)); const auto pH = to_HepLorentzVector(*the_Higgs); const auto partons = tag_extremal_jet_partons(ev); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector( partons[(ev.type() == unob)?1:0] ); auto pn = to_HepLorentzVector( partons[partons.size() - ((ev.type() == unof)?2:1)] ); auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing)); assert( (first_after_Higgs == end(partons) && ( (ev.type() == unob) || partons.back().type != pid::gluon )) || first_after_Higgs->rapidity() >= the_Higgs->rapidity() ); assert( (first_after_Higgs == begin(partons) && ( (ev.type() == unof) || partons.front().type != pid::gluon )) || (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity() ); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_Higgs == begin(partons)) ++first_after_Higgs; else if(first_after_Higgs == end(partons)) --first_after_Higgs; // t-channel momentum before Higgs auto qH = pa; for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){ qH -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if(ev.type() == unob){ current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(ev.type() == unof){ current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ current_factor = ME_Higgs_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, - q0, pa, pb, p1, pn + q0, pa, pb, p1, pn, + param_.regulator_lambda )*FKL_ladder_weight( first_after_Higgs, end_ladder, - qH - pH, pa, pb, p1, pn + qH - pH, pa, pb, p1, pn, + param_.regulator_lambda ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double MatrixElement::tree_param_partons( double alpha_s, double mur, std::vector const & partons ) const{ const double gs2 = 4.*M_PI*alpha_s; double wt = std::pow(gs2, partons.size()); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(partons.size() >= 2); for(size_t i = 1; i < partons.size()-1; ++i){ wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp()); } } return wt; } double MatrixElement::tree_param( Event const & ev, double mur ) const{ assert(is_HEJ(ev.type())); const auto & out = ev.outgoing(); const double alpha_s = alpha_s_(mur); auto AWZH_boson = std::find_if( begin(out), end(out), [](auto const & p){return is_AWZH_boson(p);} ); double AWZH_coupling = 1.; if(AWZH_boson != end(out)){ switch(AWZH_boson->type){ case pid::Higgs: { AWZH_coupling = alpha_s*alpha_s; break; } // TODO case pid::Wp:{ AWZH_coupling = gw*gw*gw*gw/4; break; } case pid::Wm:{ AWZH_coupling = gw*gw*gw*gw/4; break; } case pid::photon: case pid::Z: default: throw not_implemented("Emission of boson of unsupported type"); } } if(ev.type() == event_type::unob || ev.type() == event_type::qqxexb){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(out) + 1, end(out)}) ); } if(ev.type() == event_type::unof || ev.type() == event_type::qqxexf){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(out), end(out) - 1}) ); } return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(out)); } } // namespace HEJ diff --git a/src/YAMLreader.cc b/src/YAMLreader.cc index 72919e4..1d51bd4 100644 --- a/src/YAMLreader.cc +++ b/src/YAMLreader.cc @@ -1,471 +1,475 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/YAMLreader.hh" #include #include #include #include #include #include #include #include #include "HEJ/ScaleFunction.hh" #include "HEJ/event_types.hh" #include "HEJ/output_formats.hh" +#include "HEJ/Constants.hh" namespace HEJ{ class Event; namespace{ //! Get YAML tree of supported options /** * The configuration file is checked against this tree of options * in assert_all_options_known. */ YAML::Node const & get_supported_options(){ const static YAML::Node supported = [](){ YAML::Node supported; static const auto opts = { "trials", "min extparton pt", "max ext soft pt fraction", "FKL", "unordered", "qqx", "non-HEJ", "scales", "scale factors", "max scale ratio", "import scales", - "log correction", "event output", "analysis" + "log correction", "event output", "analysis", "regulator parameter" }; // add subnodes to "supported" - the assigned value is irrelevant for(auto && opt: opts) supported[opt] = ""; for(auto && jet_opt: {"min pt", "algorithm", "R"}){ supported["resummation jets"][jet_opt] = ""; supported["fixed order jets"][jet_opt] = ""; } for(auto && opt: {"mt", "use impact factors", "include bottom", "mb"}){ supported["Higgs coupling"][opt] = ""; } for(auto && opt: {"name", "seed"}){ supported["random generator"][opt] = ""; } return supported; }(); return supported; } fastjet::JetAlgorithm to_JetAlgorithm(std::string const & algo){ using namespace fastjet; static const std::map known = { {"kt", kt_algorithm}, {"cambridge", cambridge_algorithm}, {"antikt", antikt_algorithm}, {"cambridge for passive", cambridge_for_passive_algorithm}, {"plugin", plugin_algorithm} }; const auto res = known.find(algo); if(res == known.end()){ throw std::invalid_argument("Unknown jet algorithm " + algo); } return res->second; } EventTreatment to_EventTreatment(std::string const & name){ static const std::map known = { {"reweight", EventTreatment::reweight}, {"keep", EventTreatment::keep}, {"discard", EventTreatment::discard} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown event treatment " + name); } return res->second; } } // namespace anonymous namespace detail{ void set_from_yaml(fastjet::JetAlgorithm & setting, YAML::Node const & yaml){ setting = to_JetAlgorithm(yaml.as()); } void set_from_yaml(EventTreatment & setting, YAML::Node const & yaml){ setting = to_EventTreatment(yaml.as()); } void set_from_yaml(ParticleID & setting, YAML::Node const & yaml){ setting = to_ParticleID(yaml.as()); } } // namespace detail JetParameters get_jet_parameters( YAML::Node const & node, std::string const & entry ){ assert(node); JetParameters result; fastjet::JetAlgorithm jet_algo = fastjet::antikt_algorithm; double R; set_from_yaml_if_defined(jet_algo, node, entry, "algorithm"); set_from_yaml(R, node, entry, "R"); result.def = fastjet::JetDefinition{jet_algo, R}; set_from_yaml(result.min_pt, node, entry, "min pt"); return result; } RNGConfig to_RNGConfig( YAML::Node const & node, std::string const & entry ){ assert(node); RNGConfig result; set_from_yaml(result.name, node, entry, "name"); set_from_yaml_if_defined(result.seed, node, entry, "seed"); return result; } HiggsCouplingSettings get_Higgs_coupling( YAML::Node const & node, std::string const & entry ){ assert(node); static constexpr double mt_max = 2e4; #ifndef HEJ_BUILD_WITH_QCDLOOP if(node[entry]){ throw std::invalid_argument{ "Higgs coupling settings require building HEJ 2 " "with QCDloop support" }; } #endif HiggsCouplingSettings settings; set_from_yaml_if_defined(settings.mt, node, entry, "mt"); set_from_yaml_if_defined(settings.mb, node, entry, "mb"); set_from_yaml_if_defined(settings.include_bottom, node, entry, "include bottom"); set_from_yaml_if_defined(settings.use_impact_factors, node, entry, "use impact factors"); if(settings.use_impact_factors){ if(settings.mt != std::numeric_limits::infinity()){ throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } else{ // huge values of the top mass are numerically unstable settings.mt = std::min(settings.mt, mt_max); } return settings; } FileFormat to_FileFormat(std::string const & name){ static const std::map known = { {"Les Houches", FileFormat::Les_Houches}, {"HepMC", FileFormat::HepMC} }; const auto res = known.find(name); if(res == known.end()){ throw std::invalid_argument("Unknown file format " + name); } return res->second; } std::string extract_suffix(std::string const & filename){ size_t separator = filename.rfind('.'); if(separator == filename.npos) return {}; return filename.substr(separator + 1); } FileFormat format_from_suffix(std::string const & filename){ const std::string suffix = extract_suffix(filename); if(suffix == "lhe") return FileFormat::Les_Houches; if(suffix == "hepmc") return FileFormat::HepMC; throw std::invalid_argument{ "Can't determine format for output file " + filename }; } void assert_all_options_known( YAML::Node const & conf, YAML::Node const & supported ){ if(!conf.IsMap()) return; if(!supported.IsMap()) throw invalid_type{"must not have sub-entries"}; for(auto const & entry: conf){ const auto name = entry.first.as(); if(! supported[name]) throw unknown_option{name}; /* check sub-options, e.g. 'resummation jets: min pt' * we don't check analysis sub-options * those depend on the analysis being used and should be checked there * similar for "import scales" */ if(name != "analysis" && name != "import scales"){ try{ assert_all_options_known(conf[name], supported[name]); } catch(unknown_option const & ex){ throw unknown_option{name + ": " + ex.what()}; } catch(invalid_type const & ex){ throw invalid_type{name + ": " + ex.what()}; } } } } } // namespace HEJ namespace YAML { Node convert::encode(HEJ::OutputFile const & outfile) { Node node; node[to_string(outfile.format)] = outfile.name; return node; }; bool convert::decode(Node const & node, HEJ::OutputFile & out) { switch(node.Type()){ case NodeType::Map: { YAML::const_iterator it = node.begin(); out.format = HEJ::to_FileFormat(it->first.as()); out.name = it->second.as(); return true; } case NodeType::Scalar: out.name = node.as(); out.format = HEJ::format_from_suffix(out.name); return true; default: return false; } } } // namespace YAML namespace HEJ{ namespace detail{ void set_from_yaml(OutputFile & setting, YAML::Node const & yaml){ setting = yaml.as(); } } namespace{ void update_fixed_order_jet_parameters( JetParameters & fixed_order_jets, YAML::Node const & yaml ){ if(!yaml["fixed order jets"]) return; set_from_yaml_if_defined( fixed_order_jets.min_pt, yaml, "fixed order jets", "min pt" ); fastjet::JetAlgorithm algo = fixed_order_jets.def.jet_algorithm(); set_from_yaml_if_defined(algo, yaml, "fixed order jets", "algorithm"); double R = fixed_order_jets.def.R(); set_from_yaml_if_defined(R, yaml, "fixed order jets", "R"); fixed_order_jets.def = fastjet::JetDefinition{algo, R}; } // like std::stod, but throw if not the whole string can be converted double to_double(std::string const & str){ std::size_t pos; const double result = std::stod(str, &pos); if(pos < str.size()){ throw std::invalid_argument(str + " is not a valid double value"); } return result; } using EventScale = double (*)(Event const &); void import_scale_functions( std::string const & file, std::vector const & scale_names, std::unordered_map & known ) { auto handle = dlopen(file.c_str(), RTLD_NOW); char * error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; for(auto const & scale: scale_names) { void * sym = dlsym(handle, scale.c_str()); error = dlerror(); if(error != nullptr) throw std::runtime_error{error}; known.emplace(scale, reinterpret_cast(sym)); } } auto get_scale_map( YAML::Node const & yaml ) { std::unordered_map scale_map; scale_map.emplace("H_T", H_T); scale_map.emplace("max jet pperp", max_jet_pt); scale_map.emplace("jet invariant mass", jet_invariant_mass); scale_map.emplace("m_j1j2", m_j1j2); if(yaml["import scales"]) { if(! yaml["import scales"].IsMap()) { throw invalid_type{"Entry 'import scales' is not a map"}; } for(auto const & import: yaml["import scales"]) { const auto file = import.first.as(); const auto scale_names = import.second.IsSequence() ?import.second.as>() :std::vector{import.second.as()}; import_scale_functions(file, scale_names, scale_map); } } return scale_map; } // simple (as in non-composite) scale functions /** * An example for a simple scale function would be H_T, * H_T/2 is then composite (take H_T and then divide by 2) */ ScaleFunction parse_simple_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ) { assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); const auto it = known.find(scale_fun); if(it != end(known)) return {it->first, it->second}; try{ const double scale = to_double(scale_fun); return {scale_fun, FixedScale{scale}}; } catch(std::invalid_argument const &){} throw std::invalid_argument{"Unknown scale choice: " + scale_fun}; } std::string trim_front(std::string const & str){ const auto new_begin = std::find_if( begin(str), end(str), [](char c){ return ! std::isspace(c); } ); return std::string(new_begin, end(str)); } std::string trim_back(std::string str){ size_t pos = str.size() - 1; // use guaranteed wrap-around behaviour to check whether we have // traversed the whole string for(; pos < str.size() && std::isspace(str[pos]); --pos) {} str.resize(pos + 1); // note that pos + 1 can be 0 return str; } ScaleFunction parse_ScaleFunction( std::string const & scale_fun, std::unordered_map const & known ){ assert( scale_fun.empty() || (!std::isspace(scale_fun.front()) && !std::isspace(scale_fun.back())) ); // parse from right to left => a/b/c gives (a/b)/c const size_t delim = scale_fun.find_last_of("*/"); if(delim == scale_fun.npos){ return parse_simple_ScaleFunction(scale_fun, known); } const std::string first = trim_back(std::string{scale_fun, 0, delim}); const std::string second = trim_front(std::string{scale_fun, delim+1}); if(scale_fun[delim] == '/'){ return parse_ScaleFunction(first, known) / parse_ScaleFunction(second, known); } else{ assert(scale_fun[delim] == '*'); return parse_ScaleFunction(first, known) * parse_ScaleFunction(second, known); } } EventTreatMap get_event_treatment( YAML::Node const & yaml ){ using namespace event_type; EventTreatMap treat { {no_2_jets, EventTreatment::discard}, {bad_final_state, EventTreatment::discard}, {FKL, EventTreatment::reweight}, {unob, EventTreatment::keep}, {unof, EventTreatment::keep}, {qqxexb, EventTreatment::keep}, {qqxexf, EventTreatment::keep}, {qqxmid, EventTreatment::keep}, {nonHEJ, EventTreatment::keep} }; set_from_yaml(treat.at(FKL), yaml, "FKL"); set_from_yaml(treat.at(unob), yaml, "unordered"); treat.at(unof) = treat.at(unob); set_from_yaml(treat.at(qqxexb), yaml, "qqx"); set_from_yaml(treat.at(qqxexf), yaml, "qqx"); set_from_yaml(treat.at(qqxmid), yaml, "qqx"); set_from_yaml(treat.at(nonHEJ), yaml, "non-HEJ"); if(treat[nonHEJ] == EventTreatment::reweight){ throw std::invalid_argument{"Cannot reweight non-HEJ events"}; } return treat; } Config to_Config(YAML::Node const & yaml){ try{ assert_all_options_known(yaml, get_supported_options()); } catch(unknown_option const & ex){ throw unknown_option{std::string{"Unknown option '"} + ex.what() + "'"}; } Config config; config.resummation_jets = get_jet_parameters(yaml, "resummation jets"); config.fixed_order_jets = config.resummation_jets; update_fixed_order_jet_parameters(config.fixed_order_jets, yaml); set_from_yaml(config.min_extparton_pt, yaml, "min extparton pt"); + // Sets the standard value, then changes this if defined + config.regulator_lambda=CLAMBDA; + set_from_yaml_if_defined(config.regulator_lambda, yaml, "regulator parameter"); config.max_ext_soft_pt_fraction = std::numeric_limits::infinity(); set_from_yaml_if_defined( config.max_ext_soft_pt_fraction, yaml, "max ext soft pt fraction" ); set_from_yaml(config.trials, yaml, "trials"); set_from_yaml(config.log_correction, yaml, "log correction"); config.treat = get_event_treatment(yaml); set_from_yaml_if_defined(config.output, yaml, "event output"); config.rng = to_RNGConfig(yaml, "random generator"); set_from_yaml_if_defined(config.analysis_parameters, yaml, "analysis"); config.scales = to_ScaleConfig(yaml); config.Higgs_coupling = get_Higgs_coupling(yaml, "Higgs coupling"); return config; } } // namespace anonymous ScaleConfig to_ScaleConfig(YAML::Node const & yaml){ ScaleConfig config; auto scale_funs = get_scale_map(yaml); std::vector scales; set_from_yaml(scales, yaml, "scales"); config.base.reserve(scales.size()); std::transform( begin(scales), end(scales), std::back_inserter(config.base), [scale_funs](auto const & entry){ return parse_ScaleFunction(entry, scale_funs); } ); set_from_yaml_if_defined(config.factors, yaml, "scale factors"); config.max_ratio = std::numeric_limits::infinity(); set_from_yaml_if_defined(config.max_ratio, yaml, "max scale ratio"); return config; } Config load_config(std::string const & config_file){ try{ return to_Config(YAML::LoadFile(config_file)); } catch(...){ std::cerr << "Error reading " << config_file << ":\n "; throw; } } } // namespace HEJ