diff --git a/current_generator/jgh_j.frm b/current_generator/jgh_j.frm new file mode 100644 index 0000000..6a4dddb --- /dev/null +++ b/current_generator/jgh_j.frm @@ -0,0 +1,39 @@ +*/** +* \brief Contraction of g -> H current with FKL current +* +* \authors The HEJ collaboration (see AUTHORS for details) +* \date 2020 +* \copyright GPLv2 or later +* +* The Higgs boson momentum is decomposed as pH = pH1 + pH2 +* such that pH1, pH2 are both lightlike with positive energy +*/ +#include- include/helspin.frm +#include- include/write.frm + +v pg,pb,pn,pH,pH1,pH2,pr; +i mu,nu; + +#do h1={+,-} + #do h2={+,-} + l [jgh_j `h1'`h2'] = ( + Eps(`h1'1, mu, pg, pr) + *VH(mu, nu, pg, pg-pH) + *Current(`h2'1, pn, nu, pb) + ); + #enddo +#enddo + +multiply replace_( + pH, pH1 + pH2, + pr, pb +); + +#call ContractCurrents +.sort +format O4; +format c; +#call WriteHeader(`OUTPUT') +#call WriteOptimised(`OUTPUT',jgh_j,2,pg,pb,pn,pH1,pH2,COMPLEX:,T1,T2) +#call WriteFooter(`OUTPUT') +.end diff --git a/include/HEJ/Hjets.hh b/include/HEJ/Hjets.hh index bc87e04..9049807 100644 --- a/include/HEJ/Hjets.hh +++ b/include/HEJ/Hjets.hh @@ -1,412 +1,431 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ /** \file * \brief Functions computing the square of current contractions in H+Jets. * * This file contains all the H+Jet specific components to compute * the current contractions for valid HEJ processes, to form a full * H+Jets ME, currently one would have to use functions from the * jets.hh header also. We have FKL and also unordered components for * H+Jets. */ #pragma once #include "CLHEP/Vector/LorentzVector.h" namespace HEJ { namespace currents { using HLV = CLHEP::HepLorentzVector; //! Square of gg->gg Higgs+Jets Scattering Current /** * @param p1out Momentum of final state gluon * @param p1in Momentum of initial state gluon * @param p2out Momentum of final state gluon * @param p2in Momentum of intial state gluon * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for gg->gg * * g~p1 g~p2 * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if g is backward, qH1 is forward) */ double ME_H_gg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); + /** + * @brief Square of gq->Hq current contraction + * @param ph Outgoing Higgs boson momentum + * @param pg Incoming gluon momentum + * @param pn Momentum of outgoing particle in FKL current + * @param pb Momentum of incoming particle in FKL current + * @param mt top mass (inf or value) + * @param inc_bottom whether to include bottom mass effects (true) or not (false) + * @param mb bottom mass (value) + * @param vev Higgs vacuum expectation value + * + * Calculates || j_{gH}^\mu j_\mu ||^2 + */ + double ME_jgH_j( + HLV const & ph, HLV const & pg, + HLV const & pn, HLV const & pb, + const double mt, const bool inc_bottom, const double mb, const double vev + ); + //! Square of gq->gq Higgs+Jets Scattering Current with Higgs before Gluon /** * @param p1out Momentum of final state gluon * @param p1in Momentum of initial state gluon * @param p2out Momentum of final state gluon * @param p2in Momentum of intial state gluon * @param pH Momentum of Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contraction */ double ME_Houtside_gq(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & pH, double mt, bool include_bottom, double mb, double vev); //! Square of qg->qg Higgs+Jets Scattering Current /** * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state gluon * @param p2in Momentum of intial state gluon * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qg->qg * * q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if g is backward, qH1 is forward) */ double ME_H_qg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qbarg->qbarg Higgs+Jets Scattering Current /** * @param p1out Momentum of final state anti-quark * @param p1in Momentum of initial state anti-quark * @param p2out Momentum of final state gluon * @param p2in Momentum of intial state gluon * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qbarg->qbarg * * qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if g is backward, qH1 is forward) */ double ME_H_qbarg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qQ->qQ Higgs+Jets Scattering Current /** * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qQ->qQ * * q~p1 Q~p2 (i.e. ALWAYS p1 for quark, p2 for quark) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if Q is backward, qH1 is forward) */ double ME_H_qQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qQbar->qQbar Higgs+Jets Scattering Current /** * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state anti-quark * @param p2in Momentum of intial state anti-quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qQ->qQ * * q~p1 Qbar~p2 (i.e. ALWAYS p1 for quark, p2 for anti-quark) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if Qbar is backward, qH1 is forward) */ double ME_H_qQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qbarQ->qbarQ Higgs+Jets Scattering Current /** * @param p1out Momentum of final state anti-quark * @param p1in Momentum of initial state anti-quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qbarQ->qbarQ * * qbar~p1 Q~p2 (i.e. ALWAYS p1 for anti-quark, p2 for quark) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if Q is backward, qH1 is forward) */ double ME_H_qbarQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qbarQbar->qbarQbar Higgs+Jets Scattering Current /** * @param p1out Momentum of final state anti-quark * @param p1in Momentum of initial state anti-quark * @param p2out Momentum of final state anti-quark * @param p2in Momentum of intial state anti-quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for * qbarQbar->qbarQbar * * qbar~p1 Qbar~p2 (i.e. ALWAYS p1 for anti-quark, p2 for anti-quark) * should be called with qH1 meant to be contracted with p2 in first part of * vertex (i.e. if Qbar is backward, qH1 is forward) */ double ME_H_qbarQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! @name Unordered backwards //! @{ //! Square of qbarQ->qbarQg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state anti-quark * @param p1in Momentum of initial state anti-quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for * qbarQ->qbarQg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qQ->qQg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for qQ->qQg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qQbar->qQbarg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state quark * @param p1in Momentum of initial state quark * @param p2out Momentum of final state anti-quark * @param p2in Momentum of intial state anti-quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for * qQbar->qQbarg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of qbarQbar->qbarQbarg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state anti-quark * @param p1in Momentum of initial state anti-quark * @param p2out Momentum of final state anti-quark * @param p2in Momentum of intial state anti-quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for * qbarQbar->qbarQbarg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of gQbar->gQbarg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state gluon * @param p1in Momentum of initial state gluon * @param p2out Momentum of final state anti-quark * @param p2in Momentum of intial state anti-quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for * gQbar->gQbarg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_gQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! Square of gQ->gQg Higgs+Jets Unordered b Scattering Current /** * @param pg Momentum of unordered b gluon * @param p1out Momentum of final state gluon * @param p1in Momentum of initial state gluon * @param p2out Momentum of final state quark * @param p2in Momentum of intial state quark * @param qH1 Momentum of t-channel propagator before Higgs * @param qH2 Momentum of t-channel propagator after Higgs * @param mt Top quark mass * @param include_bottom Specifies whether bottom corrections are included * @param mb Bottom quark mass * @param vev Vacuum expectation value * @returns Square of the current contractions for gQ->gQg * * This construction is taking rapidity order: p1out >> p2out > pg */ double ME_H_unob_gQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev); //! @} //! @name impact factors for Higgs + jet //! @{ //! @brief Implements Eq. (4.22) in \cite DelDuca:2003ba with modifications to //! incoming plus momenta /** * @param p2 Momentum of Particle 2 * @param p1 Momentum of Particle 1 * @param pH Momentum of Higgs * @param vev Vacuum expectation value * @returns impact factor * * This gives the impact factor. First it determines whether this is the * case \f$p1p\sim php\gg p3p\f$ or the opposite */ double C2gHgm(HLV const & p2, HLV const & p1, HLV const & pH, double vev); //! @brief Implements Eq. (4.23) in \cite DelDuca:2003ba with modifications to //! incoming plus momenta /** * @param p2 Momentum of Particle 2 * @param p1 Momentum of Particle 1 * @param pH Momentum of Higgs * @param vev Vacuum expectation value * @returns impact factor * * This gives the impact factor. First it determines whether this is the * case \f$p1p\sim php\gg p3p\f$ or the opposite */ double C2gHgp(HLV const & p2, HLV const & p1, HLV const & pH, double vev); //! Implements Eq. (4.21) in \cite DelDuca:2003ba /** * @param p2 Momentum of Particle 2 * @param p1 Momentum of Particle 1 * @param pH Momentum of Higgs * @param vev Vacuum expectation value * @returns impact factor * * This gives the impact factor. First it determines whether this is the * case \f$p1p\sim php\gg p3p\f$ or the opposite * * @TODO remove this function is not used */ double C2qHqm(HLV const & p2, HLV const & p1, HLV const & pH, double vev); //! @} } // namespace currents } // namespace HEJ diff --git a/include/HEJ/MatrixElement.hh b/include/HEJ/MatrixElement.hh index ddc25d6..d99e637 100644 --- a/include/HEJ/MatrixElement.hh +++ b/include/HEJ/MatrixElement.hh @@ -1,204 +1,202 @@ /** \file * \brief Contains the MatrixElement Class * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include "fastjet/PseudoJet.hh" #include "HEJ/Config.hh" #include "HEJ/PDG_codes.hh" #include "HEJ/Parameters.hh" namespace CLHEP { class HepLorentzVector; } namespace HEJ { class Event; struct 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. * * If there is more than one entry in the returned vector, each entry * corresponds to the contribution from the interference of two * channels. The order of these entries matches the one returned by * the tree_kin member function, but is otherwise unspecified. * * 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 */ std::vector 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. * * If there is more than one entry in the returned vector, each entry * corresponds to the contribution from the interference of two * channels. The order of these entries matches the one returned by * the virtual_corrections member function, but is otherwise unspecified. * * @see tree, tree_param */ std::vector 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; std::vector virtual_corrections_Z_qq( Event const & event, double mur, Particle const & ZBoson ) const; double virtual_corrections_Z_qg( Event const & event, double mur, Particle const & ZBoson, bool is_gq_event ) const; std::vector 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 ) const; double tree_kin_jets( Event const & ev ) const; double tree_kin_W( Event const & ev ) const; std::vector tree_kin_Z( 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; double MH2_forwardH( - CLHEP::HepLorentzVector const & p1out, - CLHEP::HepLorentzVector const & p1in, - pid::ParticleID type2, - CLHEP::HepLorentzVector const & p2out, - CLHEP::HepLorentzVector const & p2in, - CLHEP::HepLorentzVector const & pH, - double t1, double t2 + ParticleID type2, + CLHEP::HepLorentzVector const & pg, + CLHEP::HepLorentzVector const & pb, + CLHEP::HepLorentzVector const & pH, + CLHEP::HepLorentzVector const & pn ) const; std::function alpha_s_; MatrixElementConfig param_; }; } // namespace HEJ diff --git a/src/Hjets.cc b/src/Hjets.cc index babae9f..1036644 100644 --- a/src/Hjets.cc +++ b/src/Hjets.cc @@ -1,972 +1,1013 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "HEJ/jets.hh" #include "HEJ/Hjets.hh" #include #include #include #include #include #include "HEJ/ConfigFlags.hh" #include "HEJ/Constants.hh" #include "HEJ/LorentzVector.hh" #include "HEJ/utility.hh" // generated headers #include "HEJ/currents/j_h_j.hh" #include "HEJ/currents/juno_h_j.hh" +#include "HEJ/currents/jgh_j.hh" #ifdef HEJ_BUILD_WITH_QCDLOOP #include "qcdloop/qcdloop.h" #include "HEJ/currents/jh_j.hh" #else #include "HEJ/exceptions.hh" #endif namespace HEJ { namespace currents { namespace { // short hand for math functions using std::norm; using std::abs; using std::conj; using std::pow; using std::sqrt; constexpr double infinity = std::numeric_limits::infinity(); // NOLINT // Loop integrals #ifdef HEJ_BUILD_WITH_QCDLOOP const COM LOOPRWFACTOR = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4); COM B0DD(HLV const & q, double mq) { static std::vector> result(3); static auto ql_B0 = [](){ ql::Bubble,double,double> ql_B0; ql_B0.setCacheSize(100); return ql_B0; }(); static std::vector masses(2); static std::vector momenta(1); for(auto & m: masses) m = mq*mq; momenta.front() = q.m2(); ql_B0.integral(result, 1, masses, momenta); return result[0]; } COM C0DD(HLV const & q1, HLV const & q2, double mq) { static std::vector> result(3); static auto ql_C0 = [](){ ql::Triangle,double,double> ql_C0; ql_C0.setCacheSize(100); return ql_C0; }(); static std::vector masses(3); static std::vector momenta(3); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = (q1+q2).m2(); ql_C0.integral(result, 1, masses, momenta); return result[0]; } COM D0DD(HLV const & q1, HLV const & q2, HLV q3, double mq) { static std::vector> result(3); static auto ql_D0 = [](){ ql::Box,double,double> ql_D0; ql_D0.setCacheSize(100); return ql_D0; }(); static std::vector masses(4); static std::vector momenta(6); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = q3.m2(); momenta[3] = (q1+q2+q3).m2(); momenta[4] = (q1+q2).m2(); momenta[5] = (q2+q3).m2(); ql_D0.integral(result, 1, masses, momenta); return result[0]; } // Kallen lambda functions, see eq:lambda in developer manual double lambda(const double s1, const double s2, const double s3) { return s1*s1 + s2*s2 + s3*s3 - 2*s1*s2 - 2*s1*s3 - 2*s2*s3; } // eq: T_1 in developer manual COM T1(HLV const & q1, HLV const & q2, const double m) { const double q12 = q1.m2(); const double q22 = q2.m2(); const HLV ph = q1 - q2; const double ph2 = ph.m2(); const double lam = lambda(q12, q22, ph2); assert(m > 0.); const double m2 = m*m; return - C0DD(q1, -q2, m)*(2.*m2 + 1./2.*(q12 + q22 - ph2) + 2.*q12*q22*ph2/lam) - (B0DD(q2, m) - B0DD(ph, m))*(q22 - q12 - ph2)*q22/lam - (B0DD(q1, m) - B0DD(ph, m))*(q12 - q22 - ph2)*q12/lam - 1.; } // eq: T_2 in developer manual COM T2(HLV const & q1, HLV const & q2, const double m) { const double q12 = q1.m2(); const double q22 = q2.m2(); const HLV ph = q1 - q2; const double ph2 = ph.m2(); const double lam = lambda(q12, q22, ph2); assert(m > 0.); const double m2 = m*m; return C0DD(q1, -q2, m)*( 4.*m2/lam*(ph2 - q12 - q22) - 1. - 4.*q12*q22/lam*( 1 + 3.*ph2*(q12 + q22 - ph2)/lam ) ) - (B0DD(q2, m) - B0DD(ph, m))*(1. + 6.*q12/lam*(q22 - q12 + ph2))*2.*q22/lam - (B0DD(q1, m) - B0DD(ph, m))*(1. + 6.*q22/lam*(q12 - q22 + ph2))*2.*q12/lam - 2.*(q12 + q22 - ph2)/lam; } #else // no QCDloop COM T1(HLV const & /*q1*/, HLV const & /*q2*/, double /*mt*/){ throw std::logic_error{"T1 called without QCDloop support"}; } COM T2(HLV const & /*q1*/, HLV const & /*q2*/, double /*mt*/){ throw std::logic_error{"T2 called without QCDloop support"}; } #endif // prefactors of g^{\mu \nu} and q_2^\mu q_1^\nu in Higgs boson emission vertex // see eq:VH in developer manual, but *without* global factor \alpha_s std::array TT( HLV const & qH1, HLV const & qH2, const double mt, const bool inc_bottom, const double mb, const double vev ) { if(mt == infinity) { std::array res = {qH1.dot(qH2), 1.}; for(auto & tt: res) tt /= (3.*M_PI*vev); return res; } std::array res = {T1(qH1, qH2, mt), T2(qH1, qH2, mt)}; for(auto & tt: res) tt *= mt*mt; if(inc_bottom) { res[0] += mb*mb*T1(qH1, qH2, mb); res[1] += mb*mb*T2(qH1, qH2, mb); } for(auto & tt: res) tt /= M_PI*vev; return res; } /** * @brief Higgs+Jets FKL Contributions, function to handle all incoming types. * @param p1out Outgoing Particle 1 * @param p1in Incoming Particle 1 * @param p2out Outgoing Particle 2 * @param p2in Incoming Particle 2 * @param qH1 t-channel momenta into higgs vertex * @param qH2 t-channel momenta out of higgs vertex * @param mt top mass (inf or value) * @param inc_bottom whether to include bottom mass effects (true) or not (false) * @param mb bottom mass (value) * @param vev Higgs vacuum expectation value * * Calculates j^\mu H j_\mu. FKL with higgs vertex somewhere in the FKL chain. * Handles all possible incoming states. */ double j_h_j( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, const double mt, const bool inc_bottom, const double mb, const double vev ){ using helicity::plus; using helicity::minus; const auto qqH1 = split_into_lightlike(qH1); const HLV qH11 = qqH1.first; const HLV qH12 = -qqH1.second; const auto qqH2 = split_into_lightlike(qH2); const HLV qH21 = qqH2.first; const HLV qH22 = -qqH2.second; // since qH1.m2(), qH2.m2() < 0 the following assertions are always true assert(qH11.e() > 0); assert(qH12.e() > 0); assert(qH21.e() > 0); assert(qH22.e() > 0); const auto T_pref = TT(qH1, qH2, mt, inc_bottom, mb, vev); const COM amp_mm = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_mp = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_pm = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); const COM amp_pp = HEJ::j_h_j( p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1] ); static constexpr double num_hel = 4.; // square of amplitudes, averaged over helicities const double amp2 = (norm(amp_mm) + norm(amp_mp) + norm(amp_pm) + norm(amp_pp))/num_hel; return amp2/((p1in-p1out).m2()*(p2in-p2out).m2()*qH1.m2()*qH2.m2()); } // } } // namespace double ME_H_qQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_qQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_qbarQ(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_qbarQbar(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_qg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev) * K_g(p2out,p2in)/C_A; } double ME_H_qbarg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev) * K_g(p2out,p2in)/C_A; } double ME_H_gg(HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return j_h_j(p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev) * K_g(p2out,p2in)/C_A * K_g(p1out,p1in)/C_A; } //@} + + double ME_jgH_j( + HLV const & ph, HLV const & pg, + HLV const & pn, HLV const & pb, + const double mt, const bool inc_bottom, const double mb, const double vev + ){ + using helicity::plus; + using helicity::minus; + + const auto pH = split_into_lightlike(ph); + const HLV ph1 = pH.first; + const HLV ph2 = pH.second; + // since pH.m2() > 0 the following assertions are always true + assert(ph1.e() > 0); + assert(ph2.e() > 0); + + const auto T_pref = TT(pg, pg-ph, mt, inc_bottom, mb, vev); + + const COM amp_mm = HEJ::jgh_j( + pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1] + ); + const COM amp_mp = HEJ::jgh_j( + pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1] + ); + const COM amp_pm = HEJ::jgh_j( + pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1] + ); + const COM amp_pp = HEJ::jgh_j( + pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1] + ); + + static constexpr double num_hel = 4.; + + // square of amplitudes, averaged over helicities + const double amp2 = (norm(amp_mm) + norm(amp_mp) + norm(amp_pm) + norm(amp_pp))/num_hel; + + return amp2/((pg-ph).m2()*(pb-pn).m2()); + } + + namespace { template double amp_juno_h_j( HLV const & pa, HLV const & pb, HLV const & pg, HLV const & p1, HLV const & p2, HLV const & qH11, HLV const & qH12, HLV const & qH21, HLV const & qH22, std::array const & T_pref ) { // TODO: code duplication with Wjets and pure jets const COM u1 = U1_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); const COM u2 = U2_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); const COM l = L_h_j(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]); return 2.*C_F*std::real((l-u1)*std::conj(l+u2)) + 2.*C_F*C_F/3.*std::norm(u1+u2) ; } //@{ /** * @brief Higgs+Jets Unordered Contributions, function to handle all incoming types. * @param pg Unordered Gluon momenta * @param p1out Outgoing Particle 1 * @param p1in Incoming Particle 1 * @param p2out Outgoing Particle 2 * @param p2in Incoming Particle 2 * @param qH1 t-channel momenta into higgs vertex * @param qH2 t-channel momenta out of higgs vertex * @param mt top mass (inf or value) * @param inc_bottom whether to include bottom mass effects (true) or not (false) * @param mb bottom mass (value) * * Calculates j_{uno}^\mu H j_\mu. Unordered with higgs vertex * somewhere in the FKL chain. Handles all possible incoming states. */ double juno_h_j( HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, const double mt, const bool incBot, const double mb, const double vev ){ using helicity::plus; using helicity::minus; const auto qqH1 = split_into_lightlike(qH1); const HLV qH11 = qqH1.first; const HLV qH12 = -qqH1.second; const auto qqH2 = split_into_lightlike(qH2); const HLV qH21 = qqH2.first; const HLV qH22 = -qqH2.second; // since qH1.m2(), qH2.m2() < 0 the following assertions are always true assert(qH11.e() > 0); assert(qH12.e() > 0); assert(qH21.e() > 0); assert(qH22.e() > 0); const auto T_pref = TT(qH1, qH2, mt, incBot, mb, vev); // only 4 out of the 8 helicity amplitudes are independent // we still compute all of them for better numerical stability (mirror check) MultiArray amp; #define ASSIGN_HEL(RES, J, H1, H2, HG) \ RES[H1][H2][HG] = J( \ p1in, p2in, pg, p1out, p2out, qH11, qH12, qH21, qH22, T_pref \ ) // NOLINT ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, minus); ASSIGN_HEL(amp, amp_juno_h_j, minus, minus, plus); ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, minus); ASSIGN_HEL(amp, amp_juno_h_j, minus, plus, plus); ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, minus); ASSIGN_HEL(amp, amp_juno_h_j, plus, minus, plus); ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, minus); ASSIGN_HEL(amp, amp_juno_h_j, plus, plus, plus); #undef ASSIGN_HEL const HLV q1 = p1in-p1out; // Top End const HLV q2 = p2out-p2in; // Bottom End const HLV qg = p1in-p1out-pg; // Extra bit post-gluon double ampsq = 0.; for(Helicity h1: {minus, plus}) { for(Helicity h2: {minus, plus}) { for(Helicity hg: {minus, plus}) { ampsq += amp[h1][h2][hg]; } } } ampsq /= 16.*qg.m2()*qH1.m2()*qH2.m2()*q2.m2(); // Factor of (Cf/Ca) for each quark to match ME_H_qQ. ampsq*=C_F*C_F/C_A/C_A; return ampsq; } } // namespace double ME_H_unob_qQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_unob_qbarQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_unob_qQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_unob_qbarQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev); } double ME_H_unob_gQ(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev) *K_g(p2out,p2in)/C_F; } double ME_H_unob_gQbar(HLV const & pg, HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & qH1, HLV const & qH2, double mt, bool include_bottom, double mb, double vev ){ return juno_h_j(pg, p1out, p1in, p2out, p2in, qH1, qH2, mt, include_bottom, mb, vev) *K_g(p2out,p2in)/C_F; } //@} // Begin finite mass stuff #ifdef HEJ_BUILD_WITH_QCDLOOP namespace { // All the stuff needed for the box functions in qg->qgH now... COM E1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2=-(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) + S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) + 2.*(s34 + S1)*(s34 + S1)/Delta + S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + 2.*(s34 + S1)/Delta + S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) - C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1 + k2, q2, mq)*2.*s34/ S1 - (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 + S1)/(S1*Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(2.*s34*(s34 + S1)*(S1 - S2)/(Delta*Sigma) + 2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S1)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } COM F1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR*(-S2*D0DD(k1, k2, q2, mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) - s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 + S2)*(S2 - S1)/(Delta*Sigma) + (B0DD( k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S2)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } COM G1(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return LOOPRWFACTOR*(S2*D0DD(k1, q2, k2, mq)*(Delta/s12/s12 - 4.*mq*mq/s12) - S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) - S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) - C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) - C0DD(k1, q2, mq))*4.*mq*mq/ s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./ s12 + (B0DD(k1 + q2, mq) - B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1))); } COM E4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR* (-s12*D0DD(k2, k1, q2, mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) - s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S1),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta + 2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) + s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq)) *((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma))); } COM F4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR* (-s12*D0DD(k1, k2, q2, mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) + s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S2),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/Delta + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq)) *(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma))); } COM G4(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return LOOPRWFACTOR* (-D0DD(k1, q2, k2, mq)*(Delta/s12 + (s12 + S1)/2. - 4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD( k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD( k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*2./(s12 + S2)); } COM E10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta + 12.*mq*mq*S1*(s34 + S1)/Delta/Delta - 4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) + C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) - 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) + 8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s34*(4.*s12 + 3.*S1 + S2)/(Delta*Sigma) + 8.*s12*s34*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*S1*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } COM F10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return LOOPRWFACTOR* (s12*D0DD(k1, k2, q2, mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta + 12.*mq*mq*s34*(s12 + S1)/Delta/Delta - 4.*s12*pow((s34 + S2),2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*(s34 + S2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) - C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) + 4.*s12*s34*(S2 - S1)/(Delta*Sigma) + 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) - (B0DD( k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) + 8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(-12*s34*(2*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s12*s34*s34/(S2*Delta*Delta) + 4.*s34*S1/(Delta*Sigma) - 4.*s34*(s12*s34*(2.*s12 + S2) - S1*S1*(2.*s12 + S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*(2.*s34 + S1)/(Delta*Sigma) - 8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } COM G10(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ HLV q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return LOOPRWFACTOR* (-D0DD(k1, q2, k2, mq)*(1. + 4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*4.*(s34 + S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) - B0DD(k2 + q2, mq))*4.*s34/(Delta*S2)); } COM H1DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq); } COM H4DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq); } COM H10DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq); } COM H2DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return -1.*H1DD(k2,k1,kh,mq); } COM H5DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return -1.*H4DD(k2,k1,kh,mq); } COM H12DD(HLV const & k1, HLV const & k2, HLV const & kh, double mq){ return -1.*H10DD(k2,k1,kh,mq); } HLV parity_flip(HLV const & p){ HLV flippedVector; flippedVector.setE(p.e()); flippedVector.setX(-p.x()); flippedVector.setY(-p.y()); flippedVector.setZ(-p.z()); return flippedVector; } template COM jh_j( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p2, HLV const & ph1, HLV const & ph2, double const mq ) { return (pa.z() > 0)? jh_j_forward(pa, p1, pb, p2, ph1, ph2, mq): jh_j_backward(pa, p1, pb, p2, ph1, ph2, mq) ; } template COM amp_jh_j( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p2, HLV const & ph1, HLV const & ph2, double const mq, bool const include_bottom, double const mq2, double const vev ) { COM res = 4.*mq*mq/vev*jh_j(pa, p1, pb, p2, ph1, ph2, mq); if(include_bottom) { res += 4.*mq2*mq2/vev*jh_j(pa, p1, pb, p2, ph1, ph2, mq2); } return res; } // sum over jh_j helicity amplitudes squared with + incoming gluon double ampsq_sum_jh_j( HLV const & pa, HLV const & p1, HLV const & pb, HLV const & p2, HLV const & ph1, HLV const & ph2, double const mq, bool const include_bottom, double const mq2, double const vev ) { using helicity::plus; using helicity::minus; using std::norm; const COM appp = amp_jh_j( pa, p1, pb, p2, ph1, ph2, mq, include_bottom, mq2, vev ); const COM appm = amp_jh_j( pa, p1, pb, p2, ph1, ph2, mq, include_bottom, mq2, vev ); const COM apmp = amp_jh_j( pa, p1, pb, p2, ph1, ph2, mq, include_bottom, mq2, vev ); const COM apmm = amp_jh_j( pa, p1, pb, p2, ph1, ph2, mq, include_bottom, mq2, vev ); return norm(appp) + norm(appm) + norm(apmp) + norm(apmm); } } // namespace // Higgs emitted close to gluon with full mt effects. double ME_Houtside_gq( HLV const & p1out, HLV const & p1in, HLV const & p2out, HLV const & p2in, HLV const & pH, const double mq, const bool include_bottom, const double mq2, const double vev ){ using helicity::plus; using helicity::minus; const auto ph = split_into_lightlike(pH); assert(ph.first.e() > 0); assert(ph.second.e() > 0); // incoming gluon with + helicity const double ME_plus = ampsq_sum_jh_j( p1in, p1out, p2in, p2out, ph.first, ph.second, mq, include_bottom, mq2, vev ); // incoming gluon with - helicity const double ME_minus = ampsq_sum_jh_j( parity_flip(p1in), parity_flip(p1out), parity_flip(p2in), parity_flip(p2out), parity_flip(ph.first), parity_flip(ph.second), mq, include_bottom, mq2, vev ); const double prop = m2(p1in - p1out - pH); return (ME_plus + ME_minus)/(prop*prop); } #endif // HEJ_BUILD_WITH_QCDLOOP double C2gHgm(HLV const & p2, HLV const & p1, HLV const & pH, double vev){ const double A=1./(3.*M_PI*vev); // Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta double s12 = NAN; double p1p = NAN; double p2p = NAN; COM p1perp; COM p3perp; COM phperp; // Determine first whether this is the case p1p\sim php>>p3p or the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p1p=p1.plus(); p2p=p2.plus(); } else { // opposite case p1p=p1.minus(); p2p=p2.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp +p1p/p2p*p1perp*conj(p3perp)); temp=temp*conj(temp); return temp.real(); } double C2gHgp(HLV const & p2, HLV const & p1, HLV const & pH, double vev){ const double A=1./(3.*M_PI*vev); // Implements Eq. (4.23) in hep-ph/0301013 double s12 = NAN; double php = NAN; double p1p = NAN; double phm = NAN; COM p1perp; COM p3perp; COM phperp; // Determine first whether this is the case p1p\sim php>>p3p or the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 php=pH.plus(); phm=pH.minus(); p1p=p1.plus(); } else { // opposite case php=pH.minus(); phm=pH.plus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=-COM(0,1)*A/(2.*s12)*( conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p) +s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm) -pow(conj(p3perp)+(1.+php/p1p)*conj(p1perp),2) /((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH)))) ); temp=temp*conj(temp); return temp.real(); } double C2qHqm(HLV const & p2, HLV const & p1, HLV const & pH, double vev){ const double A=1./(3.*M_PI*vev); // Implements Eq. (4.21) in hep-ph/0301013 double s12 = NAN; double p2p = NAN; double p1p = NAN; COM p1perp; COM p3perp; COM phperp; // Determine first whether this is the case p1p\sim php>>p3p or the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p2p=p2.plus(); p1p=p1.plus(); } else { // opposite case p2p=p2.minus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=A/(2.*s12)*( sqrt(p2p/p1p)*p3perp*conj(p1perp) +sqrt(p1p/p2p)*p1perp*conj(p3perp) ); temp=temp*conj(temp); return temp.real(); } } // namespace currents } // namespace HEJ diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index 594bf6a..8a969a3 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,2141 +1,2107 @@ /** * \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 "CLHEP/Vector/LorentzVector.h" #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/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 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::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::qqxexb) { 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::qqxexf){ --last_idx; } if (in[0].type != partons[0].type ){ q -= WBoson.p; #ifndef NDEBUG wc=false; #endif } } std::size_t first_idx_qqx = last_idx; std::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(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_qqx = 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_qqx) q -= partons[last_idx+1].p; if (wqq) q -= WBoson.p; for(std::size_t j = first_idx_qqx; j < last_idx_qqx; ++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::qqxexf ); #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 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; } std::size_t first_idx_qqx = last_idx; std::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(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_qqx) q -= out[last_idx+1].p; for(std::size_t j = first_idx_qqx; j < last_idx_qqx; ++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::qqxexf ); 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_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal. * @returns ME Squared for Tree-Level Current-Current Scattering * * @note The qqxf contribution can be calculated by reversing the argument ordering. */ double ME_qqx_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_q_qx ){ using namespace currents; if (bptype==pid::gluon) { if (swap_q_qx) // pq extremal return ME_Exqqx_qqbarg(pgin,pq,pqbar,pn,pb); // pqbar extremal return ME_Exqqx_qbarqg(pgin,pq,pqbar,pn,pb); } // b leg quark line if (swap_q_qx) //extremal pq return ME_Exqqx_qqbarQ(pgin,pq,pqbar,pn,pb); return ME_Exqqx_qbarqQ(pgin,pq,pqbar,pn,pb); } /* \brief Matrix element squared for central qqx 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 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 * @returns ME Squared for qqxmid Tree-Level Current-Current Scattering */ double ME_qqxmid_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 qqx amps, and qqbar ordering (default, pq backwards) const bool swap_q_qx=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_Cenqqx_qq(pa, pb, partons, is_antiquark(bptype), is_antiquark(aptype), swap_q_qx, 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 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 swap_q_qx Boolean. Ordering of qqbar pair. False: pqbar extremal. * @param wc Boolean. True->W Emitted from b. Else; emitted from leg a * @returns ME Squared for qqxb Tree-Level Current-Current Scattering * * @note calculate forwards qqx contribution by reversing argument ordering. */ double ME_W_qqx_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_q_qx, 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 wqqx emission as there is no other possible W Emission // Site. if (swap_q_qx) return ME_WExqqx_qqbarg(pa, pqbar, plbar, pl, pq, pn, pb, Wprop); return ME_WExqqx_qbarqg(pa, pq, plbar, pl, pqbar, pn, pb, Wprop); } assert(aptype==pid::gluon && bptype!=pid::gluon ); //a gluon => W emission off b leg or qqx if (!wc){ // W Emitted from backwards qqx if (swap_q_qx) return ME_WExqqx_qqbarQ(pa, pqbar, plbar, pl, pq, pn, pb, Wprop); return ME_WExqqx_qbarqQ(pa, pq, plbar, pl, pqbar, pn, pb, Wprop); } // W Must be emitted from forwards leg. return ME_W_Exqqx_QQq(pb, pa, pn, pq, pqbar, plbar, pl, is_antiquark(bptype), Wprop); } /* \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( 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 qqx amps, and qqbar ordering (default, pq backwards) const bool swap_q_qx=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_WCenqqx_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype),is_antiquark(bptype), swap_q_qx, nabove, Wprop); return wt*ME_W_Cenqqx_qq(pa, pb, pl, plbar, partons, is_antiquark(aptype), is_antiquark(bptype), swap_q_qx, 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); } CLHEP::HepLorentzVector to_HepLorentzVector(Particle const & particle){ return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()}; } 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.}; 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 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 = 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_qqxmid( 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 qqx auto qqxt = 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){ qqxt -= 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_qqxmid_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, qqxt, pa, pb, p1, pn, lambda ); return current_factor*ladder_factor; } template double tree_kin_jets_qqx(InIter BeginIn, InIter EndIn, partIter BeginPart, partIter EndPart, double lambda){ const bool swap_q_qx = is_quark(*BeginPart); const auto pgin = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(EndIn-1)); const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1))); const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?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_qqx_current( (EndIn-1)->type, pgin, pq, pqbar, pn, pb, swap_q_qx )/(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_qqxb){ return tree_kin_jets_qqx(incoming.begin(), incoming.end(), partons.begin(), partons.end(), param_.regulator_lambda); } if (ev.type()==event_type::extremal_qqxf){ return tree_kin_jets_qqx(incoming.rbegin(), incoming.rend(), partons.rbegin(), partons.rend(), param_.regulator_lambda); } if (ev.type()==event_type::central_qqx){ return tree_kin_jets_qqxmid(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_qqx(InIter BeginIn, partIter BeginPart, partIter EndPart, const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl, double lambda, ParticleProperties const & Wprop ){ const bool swap_q_qx=is_quark(*BeginPart); const auto pa = to_HepLorentzVector(*BeginIn); const auto pb = to_HepLorentzVector(*(BeginIn+1)); const auto pq = to_HepLorentzVector(*(BeginPart+(swap_q_qx?0:1))); const auto pqbar = to_HepLorentzVector(*(BeginPart+(swap_q_qx?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_qqx_current( (BeginIn)->type, (BeginIn+1)->type, pa, pb, pq, pqbar, pn, plbar, pl, swap_q_qx, 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_qqxmid( 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 qqx auto qqxt = q0; bool wqq = backmidquark->type != -(backmidquark+1)->type; // qqx emit W bool wc = !wqq && (aptype==partons.front().type); //leg b emits w assert(!wqq || !wc); if(wqq){ // emission from qqx qqxt -= pl + plbar; } else if(!wc) { // emission from leg a q0 -= pl + plbar; qqxt -= 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){ qqxt -= 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_qqxmid_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, qqxt, 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_qqxb){ return tree_kin_W_qqx(cbegin(incoming), cbegin(partons), cend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } if(ev.type() == extremal_qqxf){ return tree_kin_W_qqx(crbegin(incoming), crbegin(partons), crend(partons), plbar, pl, param_.regulator_lambda, param_.ew_parameters.Wprop()); } assert(ev.type() == central_qqx); return tree_kin_W_qqxmid(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{ + ParticleID type2, + CLHEP::HepLorentzVector const & pg, + CLHEP::HepLorentzVector const & pb, + CLHEP::HepLorentzVector const & pH, + CLHEP::HepLorentzVector const & pn + ) 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); + return K(type2, pn, pb)*ME_jgH_j( + pH, pg, pn, pb, + param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, + param_.Higgs_coupling.mb, vev + )/(N_C*N_C - 1); } 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) { + const auto partons = tag_extremal_jet_partons(ev); + + if(is_anyquark(incoming.front())) { 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 auto pg = to_HepLorentzVector(incoming.front()); + const auto pb = to_HepLorentzVector(incoming.back()); - const double t1 = q0.m2(); - const double t2 = (pn - pb).m2(); + const auto pH = to_HepLorentzVector(outgoing.front()); + const auto pn = to_HepLorentzVector(outgoing.back()); return MH2_forwardH( - p1, pa, incoming[1].type, pn, pb, pH, - t1, t2 + incoming.back().type, + pg, pb, pH, pn )*FKL_ladder_weight( - begin(partons) + 1, end(partons) - 1, - q0, pa, pb, p1, pn, + begin(partons), end(partons) - 1, + pg - pH, pg, pb, pH, 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) { + const auto partons = tag_extremal_jet_partons(ev); + + if(is_anyquark(incoming.back())) { 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.front()); + const auto pg = to_HepLorentzVector(incoming.back()); - 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(); + const auto p1 = to_HepLorentzVector(partons.front()); + const auto pH = to_HepLorentzVector(outgoing.back()); return MH2_forwardH( - pn, pb, incoming[0].type, p1, pa, pH, - t2, t1 + incoming.front().type, + pg, pa, pH, p1 )*FKL_ladder_weight( - begin(partons) + 1, end(partons) - 1, - q0, pa, pb, p1, pn, + begin(partons) + 1, end(partons), + pa - p1, pa, pg, p1, pH, 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/ME_h_mt_tree.dat b/t/ME_data/ME_h_mt_tree.dat index 2df83f3..3c409d0 100644 --- a/t/ME_data/ME_h_mt_tree.dat +++ b/t/ME_data/ME_h_mt_tree.dat @@ -1,807 +1,807 @@ 3.992094874e-08 -4.185821773e-05 +8.196772518e-06 3.886272139e-05 1.677832055e-02 4.050725953e-06 -4.055999567e-07 +4.056594514e-08 6.400462689e-07 -7.808091068e-07 +8.225909805e-08 6.053699559e-07 -2.319585446e-05 +3.114948650e-06 4.729758692e-05 1.621148815e-06 -5.231361193e-04 +1.928174823e-04 1.533703249e-05 9.573732998e-06 3.385087104e-02 8.137786341e-07 6.075300603e+02 5.951445813e-06 -1.348604339e-02 +1.923897843e-03 1.928186046e+00 5.710179361e-04 7.986990576e-04 -3.294994143e-05 -7.238377879e-04 +4.277256679e-05 +9.043190439e-04 7.024748657e-04 3.785614363e-02 -1.030019421e-06 +5.468590073e-07 2.008211091e+00 8.826462516e-08 4.627520004e-02 1.062169630e-01 9.340895791e-07 -7.480214167e-03 +1.531699393e-03 2.268339315e-02 3.600279122e-05 6.713118017e-07 -5.887169405e-03 +1.328329939e-03 7.567735364e-05 -5.345132291e-03 +5.786088949e-04 6.078088479e-04 8.320774902e-03 -3.346589526e-03 -3.866215769e-05 +4.545712675e-03 +2.958686194e-05 1.512660797e-04 5.522082262e-03 -3.877123396e-02 +4.737988069e-03 2.308611056e-07 5.469487827e-05 8.465334170e+00 -1.859344209e-04 -7.070711820e-06 +1.665315530e-04 +4.568237154e-06 5.082189599e-06 -1.433966763e-04 -2.108127322e-01 +8.661377654e-04 +1.619932721e-01 2.385636219e-02 4.438425107e+00 -2.701708040e-04 +4.661442265e-05 1.264773773e+00 5.041102815e-03 -2.763710407e-06 -6.964145214e-02 +3.716226630e-07 +1.239416175e-02 2.802069699e-03 -7.001209819e+00 +3.343368015e-01 7.203729133e-08 -4.770672076e-06 -1.891164022e-05 +1.189573607e-06 +4.521427902e-06 1.519941271e-05 1.114780681e-06 9.522309887e-04 9.263417548e-07 2.281568524e+00 1.328903015e-03 1.582906949e-03 2.692415322e+04 -1.205472913e-06 +7.173623704e-07 6.882327922e-04 1.063345487e-07 1.965024170e-02 -1.948146595e-06 +1.971623819e-07 4.393282876e-05 2.557397536e-02 4.961734615e+04 5.373795656e-06 1.069130607e-06 -2.803644104e-06 +1.895419629e-06 1.064996211e-04 -2.235242597e-04 +2.286364134e-05 7.837271762e-08 -2.699241349e-03 -5.812480053e-04 +3.565256274e-04 +3.238039599e-04 7.537519714e-02 -1.327685553e-03 -6.402315259e-05 +1.448194489e-03 +1.777422392e-05 1.147351214e-01 1.830765877e-01 1.136708763e-04 -1.828136154e-05 -2.655898063e-02 +1.953958263e-05 +3.144277863e-03 1.074103594e-03 5.057074674e-06 -1.092011261e-01 -9.236810553e-05 -9.185737640e+01 +1.058011396e-03 +1.462240721e-04 +5.698656425e+00 1.941656711e-02 8.942285136e-05 -3.601711816e-06 -3.259191724e-04 +2.063754311e-06 +5.671138604e-04 1.297220016e-04 1.168241824e-06 1.325832643e-04 3.029510688e+03 1.376379474e+00 -8.648106025e-02 -3.613453340e+01 +1.039863249e-01 +1.871097159e+01 1.782427563e-06 -1.016584963e-05 +8.185184700e-06 1.612906427e-03 -1.229884842e+00 +3.154486811e-02 3.636108827e-04 -2.703392541e-06 +1.486885165e-06 4.711472902e-06 1.707180334e-03 3.177104373e-01 -9.400835278e-04 -9.106525885e-02 +1.047618895e-03 +9.578091675e-02 2.624330979e-01 2.743869739e-02 -3.024143211e-06 -2.792768205e-03 +3.647400439e-07 +9.042081225e-04 1.890057728e-06 -1.770224622e-06 +8.498183932e-07 2.287175741e-05 -1.313886754e-05 +1.259372143e-05 2.859899271e-02 -1.518969937e-04 +5.352785148e-04 5.981959910e+00 4.271356353e+08 2.903377608e-04 1.798370637e+02 1.147588563e+00 -4.865983903e-06 +8.158198162e-07 2.932938348e-02 -1.238851181e+01 +2.579055974e+00 5.857377557e-04 1.233660800e-06 -5.854387384e-05 +4.663209870e-05 7.215231888e-01 -3.632504266e-04 +1.935911035e-05 1.098584101e-07 7.556062318e-04 3.450815469e-06 1.899482017e-03 -1.763070185e+00 +3.249063119e-01 6.624702819e-04 -4.349779721e+02 +1.127573782e+03 1.465647177e-05 -3.885460673e-03 +2.172213593e-04 5.365838158e-03 -9.878213852e-06 -2.629911519e-03 +5.541309490e-06 +1.323333114e-04 4.238909442e-02 1.808140555e-05 -8.517809026e-04 -1.965513473e+00 -2.184291473e-06 +2.065719509e-04 +1.740842154e+01 +2.317453897e-07 1.937148606e-04 6.421521178e-05 4.354545939e+02 8.570152334e-05 3.166251700e-03 -3.139728923e-07 +1.074557759e-07 1.068967309e-01 1.132391851e-07 5.424886276e-06 6.669088071e-06 -2.897766828e-03 -4.934832650e-03 +1.915154764e-02 +7.082251904e-04 2.417562879e-05 3.776703888e-08 5.135619353e-06 -1.042254650e-04 -1.363662187e+04 +5.966374165e-06 +2.336514675e+03 1.169042834e+02 -3.869061515e+02 -2.268279738e-01 +2.464038527e+02 +2.629013486e-01 2.331313778e-06 -1.020465182e-04 +7.553585245e-05 4.678573801e-03 4.271243283e-01 -1.084748834e-04 +7.683190450e-05 1.982331822e-06 2.248892692e-07 2.940369951e-05 -3.575535659e-01 -8.032884075e-06 +1.624502189e-01 +1.991493821e-06 4.885213404e-03 5.315825514e-04 -1.058261576e-07 +1.740468263e-08 1.116237563e-05 2.199703410e-06 -4.835456896e-06 -7.387759362e-03 +2.090107678e-06 +1.649862554e-03 5.758980869e-05 1.895971551e-06 -2.287620980e-03 +1.561753804e-03 7.491129176e-03 -9.445787374e-02 -7.659990515e-02 -1.380427275e-03 -1.315056721e-04 +1.518788183e-02 +4.161585757e-02 +7.658687799e-04 +8.666042251e-05 3.120838079e-07 -6.665033847e-07 +2.009776042e-07 2.243972471e-07 1.392182883e-07 3.111525022e-04 3.119219908e-03 3.190913642e-07 5.136773321e-04 -1.193962463e+00 -1.712479583e-05 +2.367557977e-01 +8.327397584e-06 2.523382412e-05 -4.021133825e-06 +2.113099081e-06 4.475792015e-07 3.590719160e-05 6.502851811e-07 2.670456447e-03 1.697306266e-01 -1.508902869e-02 -2.470827689e-05 -1.031014784e-03 -2.053297095e-04 +1.126961557e-02 +5.886943530e-05 +1.129649217e-03 +5.624466197e-05 3.063400862e-06 1.174532311e+03 -9.361653347e-06 +4.890087057e-06 7.091644248e-07 4.533040396e-05 3.947065691e-06 3.265370023e-03 1.286330059e-02 1.755757132e-02 -1.038603024e-04 -8.572677272e-02 +5.764809635e-06 +4.980128265e-02 7.479933532e-07 3.028129925e-02 2.521276297e-07 1.073442658e-06 5.899456289e-04 -3.973301390e-06 -6.759815353e+01 -8.904777259e-06 -3.989890077e-05 -3.390793260e-06 -1.175480487e-01 +2.589323889e-06 +1.682309143e+01 +4.083970000e-06 +1.037182523e-05 +2.599598125e-06 +6.417212980e-03 1.360433647e-07 7.177289496e-03 3.373827012e-06 1.333420649e-08 8.081973985e-04 8.003311035e-08 5.465880989e-07 9.880715609e-07 1.095536228e-04 4.997618552e-03 1.326013609e-02 3.701464255e-10 1.770693543e-09 6.184290168e-10 1.372839961e-08 -1.149320956e-05 +4.370563494e-07 1.085732929e-10 7.599253295e-06 1.628454581e-06 7.182423584e-08 8.292324554e-09 2.257110951e-08 -3.069783839e-07 +2.231223523e-07 7.666821934e-08 1.173780382e-07 6.007280547e-04 1.464744990e-05 1.348042340e-07 -3.843469421e-07 -7.753772161e-07 +2.463453284e-07 +3.749235270e-06 1.703489267e-02 1.789180400e-04 3.066687326e-06 2.222684120e+04 6.717770841e-02 -1.144741442e+00 +5.341756033e-01 7.090400404e-07 8.686008641e-04 3.450945082e-01 1.636463113e-03 6.817527667e-07 1.028383115e+02 1.721842018e-06 -6.431190489e-03 +1.467250438e-02 2.918858849e-04 6.129798687e-03 -3.008832924e-07 -1.065045579e-05 +2.486666846e-07 +5.748231156e-05 5.809373613e-07 2.927954692e-03 1.388370713e-08 1.356586223e-06 2.813984721e-02 1.297746410e-09 7.786568980e+01 3.264184500e-06 2.836681306e-07 -4.398094311e-07 +3.852958379e-08 3.598802506e-04 1.868220397e-06 2.849410612e-07 2.596087476e-02 4.592205330e-04 -1.459850163e-03 +8.197478301e-05 2.484254753e-08 2.391190808e-11 4.240951797e-07 1.387291457e-05 8.034780869e-07 2.124451168e-07 8.390942376e-07 2.361998231e-06 6.247996695e-02 1.828058984e-07 2.162487170e-08 6.312203852e-10 1.337381780e-07 -4.874341409e-06 +6.062997927e-07 1.702327829e-05 1.011609745e-04 1.533468085e-02 1.667951887e+00 5.030267420e-10 1.139451323e-05 3.183032762e-07 3.235639818e-10 2.430181002e-04 1.381401338e-06 4.979010339e-01 9.020006681e-05 2.759486784e-02 1.339412310e-02 7.515614514e-07 2.253136074e-04 3.711042475e-04 1.854708579e-08 5.705200359e-07 2.473009659e-09 1.514970229e-04 2.322274317e-11 1.726476060e-05 7.341296496e-09 -9.197227698e-06 +5.411694965e-06 5.279672173e-04 2.475544869e-02 6.201410725e-06 1.711624401e-08 7.902721792e-09 9.324230363e-03 4.075955349e-07 4.688015189e-03 6.749004939e-10 2.999224438e-03 6.019268002e-05 8.787361497e-08 6.945252864e-08 9.595782277e-08 1.625437720e-05 5.133881400e-07 8.769140427e-08 1.022666334e-08 3.132242163e-07 4.193068850e+01 4.631926277e-03 -1.144469109e-05 -5.664388534e-08 +1.218191114e-04 +1.189850982e-07 3.375398402e-06 2.598800453e-07 -2.337574645e-03 +9.990353196e-04 1.274531455e-05 6.215715396e-02 1.775364321e-08 2.374516724e-07 9.753985675e-05 3.869868336e-01 4.390578321e-07 3.000736624e-08 3.867882672e-05 -1.050104466e-05 +6.739472512e-06 1.399824935e-06 4.687501738e-06 2.378767375e-07 2.412995792e-06 2.704096959e-04 1.381808608e-02 1.707843925e-07 1.214986888e+00 1.767481603e+02 3.854751270e-02 7.420331844e-05 2.259667315e-04 2.942925975e-03 8.757517651e-03 5.344369331e-08 3.593653536e+02 6.467206229e-07 2.745454045e-07 8.574175060e-08 4.255704047e-07 2.002284348e-05 5.493924956e-08 1.802539443e-03 1.601027440e+01 7.907947992e+00 5.027886303e-12 1.274815065e-06 -3.367942241e-03 +5.691083917e-04 1.740566241e-05 -1.774837005e-06 +1.341515010e-06 2.208876407e-04 9.243270257e-05 8.524429842e-05 2.419623547e-06 1.230884542e-04 -8.907205979e-09 +3.222989050e-08 2.566123942e-07 8.539518312e-04 -3.043047889e-01 +3.329624812e-02 6.429865611e-04 1.201072347e-07 4.750256525e-07 2.576400162e+02 2.763586288e-09 1.267891193e-11 1.458368254e-09 3.481231126e+01 8.491610679e-05 -4.663481890e-06 +3.411851905e-07 1.904436262e-05 3.709893359e-07 1.885579069e-08 -2.676142195e-02 -1.523362214e+01 +1.490033152e-02 +2.306733771e+01 2.868038539e-10 1.357926076e-06 9.510763538e-07 3.078854471e-04 9.821307420e-04 1.300034923e-05 4.615435789e-03 1.691264974e-02 9.625087635e-07 8.484139033e-06 8.754745805e-06 1.274545153e-03 6.656659217e-07 5.641878480e-11 2.383536690e-01 1.621172643e-07 2.434094753e-04 2.309462437e-03 7.523563623e-09 7.235599966e-12 2.156490604e-02 6.454117346e-09 7.417210735e-08 6.106701982e-04 2.258166597e-07 8.795188382e-07 1.685348984e-07 -3.386721515e-06 +7.652551375e-07 5.471138589e-08 4.593203864e-04 3.358979136e-10 8.951701347e-09 1.343141523e-01 2.488084203e-06 9.204318421e-06 7.372957808e-06 3.592101975e+00 8.031445428e-06 4.722279636e-07 2.188554117e-07 2.191013639e-06 3.125129233e-07 1.431464700e-05 2.604187492e-04 2.318813104e+01 -1.775352611e-08 +1.374687235e-09 4.153948605e-06 1.759323469e-03 4.247828764e-09 3.314663642e+06 9.546118067e-04 3.628713279e-02 5.151960839e-01 -4.787966943e-03 +5.340843199e-03 2.306063053e-06 5.992993419e-02 4.335225445e-04 2.006054641e-08 1.805085571e-04 3.245279608e-04 -6.198119690e-07 +5.490962073e-06 2.877595699e-08 -4.213932767e-02 +7.930134255e-03 8.509180682e-07 1.898404662e+01 3.391894018e-04 1.572894562e-07 4.829147344e-04 4.124974312e-07 5.920628259e-09 1.133085150e-02 5.866745018e-05 1.477890940e-08 3.373201885e-01 1.783285709e-06 1.019277350e-07 1.077236434e-08 -1.229835899e-08 +6.876864321e-09 2.972773815e-03 4.579972943e-06 6.660223604e-11 3.255274106e-06 6.049147519e-02 4.838831696e-09 1.103007279e-06 -3.877771858e-04 -6.589496231e-04 +1.587132068e-04 +3.189738930e-04 2.252402531e-04 4.286469542e-06 1.183819823e-07 1.041053006e-03 8.023777540e-06 1.752202725e-10 6.760656342e-08 5.978357995e+02 1.570112402e-07 4.259084825e-02 2.734813579e-08 5.126265230e-04 1.262174839e-10 2.604124552e-07 1.853779761e-06 1.073157373e-10 1.186846920e-07 3.373467406e-02 6.462294348e-04 3.617927133e-09 1.235675676e-05 1.014307383e-05 9.670573618e-06 2.734079755e-04 9.022790348e-09 -7.773060762e-07 +1.753850348e-07 3.192976198e-05 1.605410097e-11 3.347506568e-03 6.096357632e-08 3.761983860e-04 2.002368022e-04 3.663856060e-11 9.014051583e-09 2.039080244e-08 2.382698390e-08 4.826461958e-09 3.679235457e-11 1.064541901e-07 2.647936408e-07 -3.642977680e-11 +3.225572113e-11 2.472948608e-02 3.511106044e-05 3.074534950e-09 1.016903387e-09 2.865192597e-10 8.729591282e-06 5.263810095e-09 3.162734995e-11 1.553836336e-03 -6.724165748e-07 -2.487746874e-09 +2.358540404e-08 +8.419002281e-10 7.601811027e-10 1.002426515e-02 1.099404522e-06 3.911025646e-06 3.707210974e-09 6.845537328e-09 1.288663850e-09 6.213720922e-12 3.370298227e-06 -4.073205483e-09 +9.673533191e-10 7.619743854e-05 6.153142911e-05 1.627750340e-06 7.071597014e-09 4.073669344e-10 2.323555275e+05 8.054141816e-03 8.651093356e-10 2.242451061e-10 3.216554842e-08 1.393140866e-06 2.192548352e-08 1.589958782e-08 4.319282434e-05 5.033873119e-04 1.969565661e-03 3.026962238e-05 5.847143807e-04 7.159197474e-06 6.852087190e-06 3.993958708e-06 5.152611742e-03 2.371431793e-09 9.761043438e-06 1.179390494e-08 1.500987015e-08 3.198124316e-11 3.669563179e-01 1.279595204e-09 4.407960907e-04 -1.201825090e-07 +1.541770012e-08 8.758634293e-07 7.197319124e-11 1.456265421e-06 2.118209999e-02 1.166460184e-02 3.975708970e-09 8.585133916e-09 1.431700406e-08 8.083628402e-10 7.710293690e-04 2.486389630e-12 2.209639132e-06 -3.015485921e-01 -1.011048497e-05 +1.219143263e-01 +7.821533720e-06 2.717537412e-09 3.117310948e-06 1.196435320e-06 -4.481379516e-09 +1.490502179e-09 2.669049021e-08 3.465183725e-09 5.746377675e-06 6.827640940e-05 2.205905568e-07 8.277136503e-08 5.233327912e-07 6.831993026e-11 1.722970086e-06 9.209494681e-07 1.853861729e-04 2.597674405e-03 6.127253802e-11 1.815659011e-04 7.719959906e-14 1.355959480e-08 1.068097702e-06 4.571150240e-07 7.835485051e-07 2.443967251e+00 7.521608515e-11 3.772716632e-06 1.492486651e-07 3.947354347e-09 8.108470196e-09 1.226200398e-04 2.265314247e-11 1.841207960e-03 1.113595985e-02 -2.535644631e-06 +3.333353103e-07 4.652540757e-11 4.191956716e-10 4.620562396e-05 4.971924276e-04 7.584566548e-09 4.136067428e-13 2.639942371e-03 2.560094895e-08 6.494205457e-08 3.813228849e-11 2.394276536e-01 9.383928172e-08 2.354446386e-03 1.169788005e-11 1.635056993e-09 8.927013480e-08 3.401987129e-10 -6.359659489e-07 +9.632532344e-07 8.846280299e-07 8.063601622e-13 -1.272867559e-07 +8.804357796e-08 5.056101098e-08 3.717053057e-05 1.584125694e-06 4.622564202e-08 3.862253555e-05 2.369045122e-04 -8.076951359e-08 +2.214677222e-07 2.028201866e-10 1.490632954e-05 2.809039801e-06 4.286284920e-10 1.225531203e-09 9.651458117e-10 5.073720873e-06 1.168279305e-03 -3.125393506e-07 +2.920249868e-07 2.924547371e-06 3.901702794e-06 -1.812373580e-08 +7.800021966e-09 9.208991156e-06 4.592993319e-03 1.243474714e-11 1.620460995e-06 4.067811815e-04 6.212555249e-02 2.225663359e-12 1.546065372e-10 2.095252325e-01 2.734979770e-02 1.919282224e-12 1.583212739e-03 -1.842092784e+00 +2.828972283e-01 8.594625100e-01 9.586561470e-02 5.893601451e-11 9.940877292e-12 5.380308774e-02 1.865793667e-04 6.607021968e-08 3.838609063e-07 9.049824655e-07 3.888723805e-04 2.777900492e-11 2.684360703e-04 5.913695711e-06 5.737197844e-04 4.110113744e-09 5.107891323e-08 3.390290129e-09 8.095451133e-12 4.321661893e-08 4.873283944e-07 1.274003933e-11 9.104409500e-12 7.179607836e-07 3.906896479e-11 1.948313091e-08 1.164976881e-01 1.112886046e-06 -1.022484093e-06 +4.361734352e-07 2.591738179e-07 1.275159404e-04 4.499451474e-04 5.923928986e-10 9.029183686e-06 1.165979189e-08 4.397445726e-07 1.101331326e-15 7.092692390e-08 6.766025857e-05 2.582987643e-06 2.626015584e-07 2.252897510e-10 1.582238936e-05 -1.110207514e-04 +2.124085547e-05 1.437364367e-09 3.309200992e-03 4.409409406e-05 3.202533864e-09 7.209532734e-10 3.001628290e-04 4.616089749e-05 8.022866888e-10 6.706777161e-12 3.828630640e-05 7.926932401e-05 4.488498744e-05 1.348692501e-06 9.161443367e-10 8.678114985e-07 1.926070493e-06 4.765469595e-09 1.893168874e-08 1.549578157e-08 9.229013716e-06 1.035062224e-04 6.357487377e-07 4.034696570e-06 1.504164503e-08 9.700520388e-12 2.293483524e-05 2.338097661e-02 3.059032675e-03 5.427962171e-08 1.401989839e-03 1.579634046e-12 1.072621482e+01 6.864694660e-04 3.075076448e-07 4.970636771e-09 2.187827570e-08 diff --git a/t/ME_data/ME_h_mtinf_tree.dat b/t/ME_data/ME_h_mtinf_tree.dat index 0a1434f..7b8a4d3 100644 --- a/t/ME_data/ME_h_mtinf_tree.dat +++ b/t/ME_data/ME_h_mtinf_tree.dat @@ -1,807 +1,807 @@ 5.572804568e-08 -2.820771959e-05 +7.733052565e-06 3.888389298e-05 2.042806864e-02 4.680070730e-06 -1.526564425e-07 +3.986132965e-08 6.409292912e-07 -5.895267757e-07 +7.965976214e-08 6.467450370e-07 -1.483353367e-05 +2.966246407e-06 4.871451037e-05 1.941740277e-06 -4.430985498e-04 +1.888801668e-04 1.597702620e-05 1.229649514e-05 3.422004512e-02 8.538221861e-07 5.785658479e+02 6.315861962e-06 -1.308773607e-02 +1.993364479e-03 2.249050403e+00 1.159555900e-03 7.916035504e-04 -3.014257182e-05 -1.601499960e-04 +5.358231482e-05 +1.132256948e-03 9.764963622e-04 3.871461664e-02 -3.207564898e-07 +6.311772268e-07 1.938663548e+00 8.986523335e-08 4.606848096e-02 1.664459283e-01 1.059082305e-06 -4.397799541e-03 +1.509034797e-03 2.194313999e-02 3.746404383e-05 6.881839613e-07 -6.326635550e-03 +1.249421212e-03 2.323756323e-04 -4.620157100e-03 +5.675699760e-04 5.863695633e-04 1.141179377e-02 -2.782529791e-03 -2.283025018e-05 +4.580570661e-03 +2.837300381e-05 1.815046198e-04 9.223059152e-03 -3.244681476e-02 +4.861319433e-03 3.089668561e-07 5.373616917e-05 1.201443762e+01 -1.866292793e-04 -4.281174307e-06 +1.931944033e-04 +8.350132849e-06 9.885414773e-05 -1.131417601e-04 -1.411312700e-01 +8.818949784e-04 +1.555676074e-01 2.370571321e-02 5.430257614e+00 -2.196112759e-04 +4.522187081e-05 1.255259399e+00 5.534475149e-03 -1.490132281e-06 -4.094139359e-02 +3.591615295e-07 +1.390652027e-02 3.105378110e-03 -4.886950565e+00 +3.359717020e-01 9.549405497e-08 -4.289881031e-06 -1.261969556e-05 +1.139757378e-06 +4.451372181e-06 1.504176295e-05 4.024167287e-06 9.287030781e-04 3.575383717e-06 2.205209754e+00 1.345140269e-03 1.944571866e-03 2.590234738e+04 -4.789731617e-07 +7.661844320e-07 1.265615683e-03 2.194569070e-07 2.564179368e-02 -1.735139373e-06 +1.917476016e-07 5.374038239e-05 3.608593751e-02 4.828245194e+04 6.210373583e-06 8.257528074e-06 -2.432389755e-06 +1.940573073e-06 8.192357630e-04 -1.953389911e-04 +2.388705403e-05 8.449199830e-08 -2.494709041e-03 -5.259481531e-04 +3.470747463e-04 +3.277775251e-04 8.386707002e-02 -1.107702869e-03 -1.673540132e-05 +1.481848847e-03 +1.949945883e-05 1.208500717e-01 2.018148288e-01 1.466103051e-04 -1.495761867e-05 -1.793952356e-02 +2.191712369e-05 +3.065352390e-03 1.468013763e-03 6.097640805e-06 -1.145453769e-01 -6.033545136e-05 -8.309124414e+01 +1.015360716e-03 +1.455096615e-04 +5.889686948e+00 1.949243022e-02 1.087303162e-04 -5.946680972e-07 -2.078428330e-04 +2.049386967e-06 +5.509716144e-04 1.537887580e-04 1.356242586e-06 1.974324073e-04 2.963883364e+03 1.431916564e+00 -7.883498059e-02 -3.372884128e+01 +1.078616845e-01 +1.796690978e+01 1.940705323e-06 -8.496241084e-06 +1.253258061e-05 1.668054846e-03 -1.104399324e+00 +3.522247815e-02 4.749137893e-04 -2.250590696e-06 +1.585422170e-06 5.016010297e-06 2.680527032e-03 3.199000005e-01 -6.817976032e-04 -6.150184679e-02 +1.016549548e-03 +9.014438386e-02 3.785330696e-01 3.056569994e-02 -1.281744929e-06 -2.332077392e-03 +3.677057648e-07 +9.914719085e-04 1.183459824e-05 -1.706974631e-06 +8.537388750e-07 2.249345228e-05 -1.297721679e-05 +1.569046912e-05 2.934136148e-02 -1.570407685e-04 +8.851338772e-04 6.389068161e+00 4.837678972e+08 2.910973086e-04 1.716201719e+02 1.631035646e+00 -3.743085237e-06 +8.829866494e-07 3.066208474e-02 -1.304137762e+01 +2.688394911e+00 6.836375721e-04 2.114055578e-06 -2.780369101e-05 +4.485509365e-05 1.195553498e+00 -2.234658305e-04 +1.848548301e-05 1.485581537e-07 7.396640117e-04 6.598939949e-06 2.119548369e-03 -1.561680586e+00 +3.116415021e-01 6.658022304e-04 -4.724282754e+02 +1.458066815e+03 1.430389956e-05 -3.772331655e-03 +2.625824174e-04 5.369080172e-03 -4.438842246e-06 -2.336918233e-03 +2.370008214e-05 +1.411742942e-04 1.690138110e-01 1.829564964e-05 -6.883322589e-04 -1.804051793e+00 -1.736607582e-06 +2.051207741e-04 +1.872050738e+01 +2.550452687e-07 2.255678771e-04 6.524599620e-05 4.372467535e+02 8.387224326e-05 3.343506359e-03 -1.890263288e-07 +2.854101208e-07 1.110841737e-01 2.100158475e-07 6.050673169e-06 6.731118778e-06 -2.104439733e-03 -3.011438923e-03 +1.802124455e-02 +7.036560727e-04 2.674458314e-05 4.144657087e-08 6.897914197e-06 -4.406776189e-05 -1.315696681e+04 +5.771319258e-06 +2.279022969e+03 1.121949844e+02 -4.072498680e+02 -1.944309706e-01 +2.850696815e+02 +2.590926771e-01 3.627260888e-06 -8.933360379e-05 +7.634301399e-05 4.530032022e-03 4.366576508e-01 -9.763726361e-05 +7.607425884e-05 2.773608322e-06 1.698189395e-06 3.062031893e-05 -3.566098292e-01 -4.856826323e-06 +1.698921313e-01 +1.898966670e-06 5.045283380e-03 5.299933444e-04 -2.781944269e-08 +3.017540386e-08 1.275786883e-05 9.664916790e-06 -2.533942575e-06 -6.437171912e-03 +2.077047482e-06 +1.822359375e-03 6.595509942e-05 1.892603441e-06 -1.839316783e-03 +1.511182400e-03 7.097187691e-02 -9.275503249e-02 -8.054671830e-02 -7.827390196e-04 -1.299825904e-04 +1.630340232e-02 +4.616139930e-02 +7.280060014e-04 +9.002810899e-05 6.948078299e-07 -4.624337301e-07 +1.981505848e-07 2.351049919e-07 1.218637919e-06 3.374685455e-04 3.125835543e-03 6.317642047e-07 5.304174029e-04 -1.037834687e+00 -1.574520165e-05 +2.251947692e-01 +8.085900174e-06 4.025622878e-05 -2.138630832e-06 +3.579791076e-06 4.624321565e-07 3.557996624e-05 7.632352346e-07 4.087274884e-03 1.707433480e-01 -9.388651019e-03 -1.084988719e-05 -9.918063878e-04 -1.280830068e-04 +1.094536785e-02 +5.812909858e-05 +1.130435999e-03 +5.935609372e-05 3.013380422e-06 1.140558348e+03 -6.138041700e-06 +4.746179921e-06 5.936317197e-06 4.476488634e-05 4.866820693e-06 3.364268242e-03 1.345732182e-02 1.691848296e-02 -8.648323648e-05 -7.719194047e-02 +6.469953760e-06 +5.004395458e-02 2.261298664e-06 3.032245343e-02 6.391553476e-07 2.699173814e-06 6.756168369e-04 -2.416507126e-06 -6.467006973e+01 -8.329641063e-06 -3.183693467e-05 -1.709549257e-06 -9.650189020e-02 +2.755355122e-06 +1.718242742e+01 +4.657914496e-06 +9.996247392e-06 +2.629072903e-06 +7.102602193e-03 1.391403765e-07 7.263556652e-03 4.999370887e-06 1.625061232e-08 8.509561784e-04 1.202205227e-07 7.144871653e-07 9.471851526e-07 1.191493923e-04 5.840966853e-03 1.370040570e-02 5.713297919e-10 3.607906530e-08 7.307412331e-10 1.451344169e-08 -1.724590008e-05 +4.677166985e-07 1.522273099e-10 7.303492170e-06 2.394106739e-06 9.439029291e-08 1.197178286e-08 2.852370502e-08 -2.918666188e-07 +2.402523519e-07 1.633081697e-07 4.999725784e-07 4.364098403e-03 1.530966988e-05 5.067352329e-07 -3.099247619e-07 -2.158143194e-06 +2.520162972e-07 +3.743151285e-06 2.052969833e-02 1.958876337e-04 9.931888956e-06 2.504301598e+04 7.117111486e-02 -1.129363603e+00 +5.435560334e-01 8.519028290e-07 8.611561913e-04 3.313925797e-01 1.685927042e-03 7.737216703e-07 1.356040671e+02 1.715670931e-06 -1.765374179e-02 +1.394834207e-02 3.045426205e-04 9.039006896e-03 -4.862966762e-07 -9.910844952e-07 +2.394536597e-07 +5.543261305e-05 6.749068743e-07 2.759944309e-03 4.200025166e-07 1.836184773e-06 2.906214889e-02 1.426547480e-09 1.236256079e+02 3.871989502e-06 2.818399887e-07 -2.178107596e-07 +3.742818303e-08 3.622199797e-04 1.953122542e-06 2.920307346e-07 2.735652104e-02 4.826860535e-04 -1.292659386e-03 +8.429485887e-05 4.952619212e-08 9.362955197e-11 4.786061319e-07 1.721659726e-05 1.051237231e-06 1.003149697e-06 1.453591818e-06 6.002229724e-06 7.422256382e-02 4.977610976e-07 2.615599662e-08 6.718184084e-10 1.459950517e-07 -7.829096020e-06 +7.273135157e-07 2.230570626e-05 1.132328395e-04 1.544043912e-02 1.732457391e+00 1.154545164e-08 1.291848361e-05 3.342745364e-07 3.485361981e-09 2.706474718e-04 5.574804238e-06 8.990146655e-01 1.001326662e-04 2.233523312e-02 1.335995585e-02 7.467887166e-07 2.231001789e-04 3.699152099e-04 1.894965441e-08 1.654294075e-06 4.086577889e-09 1.755105427e-04 8.416393596e-11 5.823380307e-04 7.687882493e-09 -9.965822218e-06 +5.276775857e-06 7.023227535e-04 2.470065095e-02 6.259343015e-06 3.462534698e-08 1.721326261e-08 9.100990839e-03 6.068046546e-07 5.073613628e-03 8.104795790e-10 4.489318495e-03 7.059335059e-05 1.116936224e-07 9.447042548e-08 1.142507724e-07 2.495452814e-05 1.102725863e-06 9.170351616e-08 1.116173034e-08 3.160808071e-07 4.719294278e+01 5.609765930e-03 -2.640897278e-05 -3.211672292e-08 +1.271062114e-04 +1.193166957e-07 3.933913946e-06 1.705704126e-06 -2.687621190e-03 +1.029819588e-03 1.353055173e-05 6.037711935e-02 4.178013966e-08 2.754380072e-07 1.284315046e-04 4.613545669e-01 8.243835252e-07 1.031488143e-06 4.518197192e-05 -2.299612194e-05 +6.428165315e-06 1.469929984e-06 4.854750586e-06 2.231002492e-07 5.180713437e-06 2.793839809e-04 1.366295653e-02 1.616729204e-06 1.076928081e+01 1.777645145e+02 4.165016415e-02 7.327803251e-05 2.312992445e-04 3.278010006e-03 1.857660707e-02 4.395091584e-07 8.377603390e+02 7.221567849e-07 7.394039845e-07 9.034821635e-08 6.938495910e-07 1.962974288e-05 7.831554318e-08 2.417915157e-03 1.681184118e+01 8.358803602e+00 2.506211004e-11 1.351591326e-06 -3.909223946e-03 +5.653784700e-04 1.874994001e-05 -1.059999529e-06 +1.273948784e-06 2.518063912e-04 2.797398009e-04 8.269820486e-05 2.741694330e-06 1.362989830e-04 -4.102693857e-09 +3.236163537e-08 3.180061341e-07 1.076406164e-03 -2.895143691e-01 +3.310137746e-02 7.799632930e-04 1.245152763e-07 6.645743195e-07 2.622509110e+02 2.169016271e-07 2.388605910e-10 1.819082046e-09 3.475786517e+01 1.057002104e-04 -3.534833550e-06 +3.418764175e-07 1.874181702e-05 4.841351554e-07 2.302607917e-07 -2.925781695e-02 -6.513680610e+01 +1.696858652e-02 +2.620683902e+01 4.055168507e-10 1.438034814e-06 1.064547563e-06 3.082750857e-04 1.073504174e-03 1.300681473e-05 6.459601030e-03 1.895987606e-02 1.060829403e-06 8.238797134e-06 9.317793462e-06 1.278356229e-03 2.279610719e-06 6.248891603e-10 2.530038522e-01 1.627768885e-07 3.293036130e-04 2.529327920e-03 1.911162552e-08 8.591395213e-12 3.257413794e-02 7.190258956e-09 8.073987581e-08 6.258596186e-04 2.468811944e-07 9.951052105e-07 1.743583334e-07 -4.744668344e-06 +7.892382885e-07 1.326313019e-07 5.072538045e-04 3.985459568e-10 1.104303535e-08 1.692316356e-01 3.763376717e-06 1.038624337e-05 8.048797417e-06 3.622786876e+00 2.101538644e-05 5.040009061e-07 2.335240322e-07 2.498424766e-06 4.202295601e-07 4.830236598e-05 2.774322272e-04 2.330464677e+01 -2.127545991e-08 +2.256028846e-09 8.724210005e-06 1.781787412e-03 4.502368871e-09 3.186046302e+06 9.687848173e-04 4.668030680e-02 5.164858709e-01 -5.049608334e-03 +5.778084949e-03 2.351688446e-06 7.460050065e-02 6.143272317e-04 3.188169157e-08 2.788262516e-04 3.793487471e-04 -8.181505685e-07 +5.337517599e-06 2.933159250e-08 -1.641859923e-02 +7.844566698e-03 8.517988386e-07 1.843228168e+01 4.363104434e-04 1.649030026e-07 5.022484285e-04 5.146033092e-07 6.127315124e-09 1.194048916e-02 7.296523578e-05 2.261779127e-08 3.550279519e-01 4.665801541e-06 1.133805742e-07 1.077545331e-08 -4.317767277e-09 +7.131867193e-09 3.024651909e-03 4.900211468e-06 8.425164568e-10 3.260771789e-06 7.229895794e-02 9.030803015e-09 2.800132225e-06 -5.852057905e-04 -3.863934208e-04 +1.589382128e-04 +3.832403987e-04 2.302342216e-04 7.639352281e-06 1.200219161e-07 1.046101258e-03 8.161065404e-06 4.074095042e-10 7.474610092e-08 8.139530652e+02 2.860319094e-07 5.402717193e-02 1.094737399e-07 7.218336860e-04 1.363337213e-10 2.897166555e-07 5.466416822e-06 1.046528666e-09 1.214075429e-07 3.517308930e-02 6.876570099e-04 4.381704653e-09 1.425064448e-05 1.092939333e-05 1.013613284e-05 4.161445140e-04 3.929296641e-08 -1.861345650e-06 +2.029432413e-07 3.350149307e-05 1.942422023e-11 3.921387549e-03 6.092617772e-08 4.006639955e-04 2.193198014e-04 4.588823322e-11 2.247656490e-08 4.235884774e-07 2.699312268e-08 7.118759240e-09 4.661647359e-11 1.283410386e-07 3.251622127e-07 -2.281843258e-10 +8.943926432e-11 2.504071095e-02 5.445082047e-05 3.290489223e-09 1.115318901e-09 3.484841326e-10 1.191475982e-05 9.490624444e-09 3.852807809e-11 1.668368723e-03 -1.354227113e-06 -3.547185971e-09 +3.149765928e-08 +8.480556366e-10 1.128183821e-09 1.205968563e-02 2.973950143e-06 7.161620283e-06 9.571931819e-09 1.477120229e-08 1.360026180e-09 6.536532052e-12 3.756762190e-06 -3.519031871e-09 +1.003492002e-09 7.705635137e-05 6.527963639e-05 1.900806678e-06 9.604137497e-09 4.832543450e-10 2.461535131e+05 8.317793495e-03 1.288549452e-09 2.984905834e-10 3.442753033e-08 1.446019240e-06 2.554651308e-08 1.716746196e-08 4.218315871e-05 6.263553513e-04 2.045309149e-03 4.000660563e-05 5.947272843e-04 1.599749747e-05 1.072982129e-05 4.537764019e-06 5.320579544e-03 3.835647646e-09 1.154139540e-05 1.347260164e-08 1.754674352e-08 3.517542867e-09 3.604393244e-01 1.896710657e-09 4.711292756e-04 -1.566983764e-07 +2.519062555e-08 1.025894195e-06 1.607487456e-10 2.074443931e-06 2.716520917e-02 1.306997660e-02 7.745194101e-09 1.008059614e-08 5.325487074e-08 8.890288216e-10 8.313145689e-04 4.728469346e-12 2.657420217e-06 -3.166893117e-01 -1.834172030e-05 +1.203719769e-01 +8.070335574e-06 3.333556889e-09 4.293987770e-06 1.491532386e-06 -1.410273888e-08 +3.346466355e-09 2.890119270e-08 4.758693888e-09 6.648097567e-06 7.550647813e-05 7.063391619e-07 2.802247921e-07 5.253163265e-07 7.788619114e-11 1.927421367e-06 9.100970575e-07 6.011167546e-04 3.591861909e-03 7.601695154e-11 3.437030489e-04 2.831401705e-13 1.700113295e-08 1.120067516e-06 5.630141622e-07 8.019872998e-07 2.370711819e+00 9.945530280e-10 4.398816915e-05 2.658536988e-07 4.070765645e-09 1.348964882e-08 1.346100783e-04 2.397829940e-11 2.216977235e-03 1.189249860e-02 -2.450245944e-06 +3.514778160e-07 8.735058436e-11 5.128984744e-10 5.483813146e-05 5.291911467e-04 1.678041194e-08 5.087404421e-13 3.470363247e-03 2.618358550e-08 6.675635132e-08 3.972641717e-11 2.395412045e-01 1.018953051e-07 2.618325163e-03 1.491391997e-11 2.263240793e-09 1.092156400e-07 4.240317799e-10 -8.177875085e-07 +9.763229490e-07 9.662018142e-07 1.204063649e-12 -1.737324984e-07 +9.132831222e-08 9.262311981e-08 4.417568826e-05 2.309013645e-06 1.033175136e-07 5.304100308e-05 5.585876907e-04 -1.268736665e-07 +2.702525918e-07 4.226412120e-10 1.496595206e-05 4.128704267e-06 1.069473089e-08 1.288077396e-09 1.976544269e-09 7.208953037e-06 1.198806306e-03 -5.010388968e-07 +4.205581376e-07 2.893094105e-06 3.846974156e-06 -1.917703115e-08 +9.196702346e-09 9.122224995e-06 5.017367317e-03 1.467921358e-11 1.712696979e-06 1.576706134e-03 6.218478551e-02 2.210917575e-11 1.833388455e-10 2.076192027e-01 2.666210155e-02 3.173840670e-12 1.736922913e-03 -1.270319507e+00 +2.846260134e-01 1.248995403e+00 1.027514498e-01 1.303631291e-10 9.717345904e-12 5.470246326e-02 2.344214153e-04 8.645766604e-08 4.546502119e-07 9.821436471e-07 4.543371326e-04 2.734115483e-11 2.831850423e-04 1.396680491e-05 6.540453499e-04 6.962926737e-09 5.502582362e-08 3.370076910e-09 1.449987618e-11 4.711965187e-08 4.729467396e-07 1.929793020e-11 8.523475458e-11 7.387954245e-07 1.080040278e-10 2.281110242e-08 1.166677915e-01 1.409613611e-06 -1.147187181e-06 +4.514266010e-07 2.672310396e-07 1.255071010e-04 6.229838431e-04 6.419598350e-10 1.064850629e-05 1.156380518e-08 4.852815617e-07 6.863856159e-14 7.264596037e-08 8.284444693e-05 3.088627657e-05 6.018435532e-07 2.429391340e-10 1.743400033e-05 -7.000642750e-05 +2.123828380e-05 1.477860697e-09 3.426271735e-03 4.850294661e-05 3.693472519e-09 8.188840272e-10 3.134182902e-04 4.963596239e-05 1.459290217e-09 2.466456916e-11 5.230009936e-05 1.211786947e-04 5.437715515e-05 1.641328630e-06 1.042738834e-09 9.306578287e-07 2.233023059e-06 6.994767103e-09 5.103840066e-08 2.180291056e-08 1.054043488e-05 1.337487784e-04 7.502939118e-07 5.169726504e-06 1.781680794e-08 1.074536464e-10 2.730551306e-05 2.550106620e-02 3.354202915e-03 9.945180494e-08 3.884971018e-03 1.801945180e-12 1.090804491e+01 6.883835961e-04 5.290759535e-07 7.522574713e-09 3.280439980e-08 diff --git a/t/ME_data/ME_h_mtmb_tree.dat b/t/ME_data/ME_h_mtmb_tree.dat index 0ba3603..ad6bd71 100644 --- a/t/ME_data/ME_h_mtmb_tree.dat +++ b/t/ME_data/ME_h_mtmb_tree.dat @@ -1,807 +1,807 @@ 4.036423148e-08 -3.940695462e-05 +7.736848255e-06 3.971238751e-05 1.704916390e-02 4.126163153e-06 -4.120221319e-07 +4.133766165e-08 6.515738458e-07 -7.794000157e-07 +8.291597961e-08 6.202166544e-07 -2.273189686e-05 +3.055343331e-06 4.920813253e-05 1.645597849e-06 -5.308682891e-04 +1.961234164e-04 1.548455272e-05 9.714753846e-06 3.502266406e-02 8.370742698e-07 6.069969581e+02 6.064159817e-06 -1.379402693e-02 +1.985298245e-03 1.960086246e+00 5.767474795e-04 8.255987630e-04 -3.390595917e-05 -7.438763677e-04 +4.405908649e-05 +9.315333388e-04 7.109077582e-04 3.938946000e-02 -1.058566813e-06 +5.645949929e-07 1.999418757e+00 9.091732441e-08 4.795920875e-02 1.076748080e-01 9.486837249e-07 -7.638363821e-03 +1.562860509e-03 2.337001022e-02 3.677679730e-05 6.849081872e-07 -5.390885072e-03 +1.218640337e-03 7.633470907e-05 -5.345400607e-03 +5.890114065e-04 6.189558681e-04 8.432994191e-03 -3.442253987e-03 -3.843951823e-05 +4.672976133e-03 +2.944749147e-05 1.543375743e-04 5.581577133e-03 -3.950854564e-02 +4.884634521e-03 2.346338154e-07 5.658368514e-05 8.586756438e+00 -1.916769456e-04 -7.218548577e-06 +1.719104149e-04 +4.667893099e-06 5.110672307e-06 -1.477663167e-04 -2.103313165e-01 +8.920555618e-04 +1.615939616e-01 2.423820202e-02 4.531340391e+00 -2.719470313e-04 +4.706751300e-05 1.312005381e+00 5.118439829e-03 -2.732686015e-06 -7.187306762e-02 +3.738025342e-07 +1.280402982e-02 2.872395003e-03 -7.193746078e+00 +3.434829822e-01 7.307477806e-08 -4.677914766e-06 -1.913225677e-05 +1.182185461e-06 +4.611796306e-06 1.558401995e-05 1.125457825e-06 9.938467298e-04 9.305719108e-07 2.351577836e+00 1.351848038e-03 1.604570791e-03 2.725148428e+04 -1.234759714e-06 +7.412335698e-07 6.959827491e-04 1.069546678e-07 2.002109665e-02 -1.958000049e-06 +1.995123777e-07 4.452818467e-05 2.597028907e-02 5.063411137e+04 5.471108411e-06 1.074665005e-06 -2.864245369e-06 +1.953754359e-06 1.070605830e-04 -2.289704455e-04 +2.360615790e-05 8.007807626e-08 -2.693200621e-03 -5.984302131e-04 +3.610629401e-04 +3.331715455e-04 7.711293504e-02 -1.369089840e-03 -6.589295759e-05 +1.492536307e-03 +1.836757682e-05 1.177180255e-01 1.871326492e-01 1.153295778e-04 -1.879598427e-05 -2.690894990e-02 +2.018590409e-05 +3.187839556e-03 1.085933997e-03 5.144887011e-06 -1.054704620e-01 -9.467940750e-05 -9.451233173e+01 +1.054304440e-03 +1.498081308e-04 +5.879291316e+00 2.029765471e-02 9.076488466e-05 -3.660807690e-06 -3.293606172e-04 +2.112793665e-06 +5.733726407e-04 1.320311226e-04 1.191386384e-06 1.341747592e-04 3.090047756e+03 1.403839816e+00 -8.913504370e-02 -3.586958366e+01 +1.073140504e-01 +1.866226673e+01 1.822509069e-06 -1.040930061e-05 +8.389617599e-06 1.663315754e-03 -1.266898109e+00 +3.259078056e-02 3.675662279e-04 -2.785066676e-06 +1.536320239e-06 4.832687364e-06 1.726784601e-03 3.288072256e-01 -9.498002614e-04 -8.414617436e-02 +1.058019447e-03 +8.851033210e-02 2.659712131e-01 2.787024900e-02 -3.083583490e-06 -2.873679676e-03 +3.749902643e-07 +9.343818265e-04 1.905511571e-06 -1.790088684e-06 +8.730061394e-07 2.337835933e-05 -1.351489964e-05 +1.297453696e-05 2.913100451e-02 -1.547540316e-04 +5.478435495e-04 6.089047514e+00 4.379840484e+08 3.001194264e-04 1.817496329e+02 1.164476362e+00 -5.011329237e-06 +8.430765400e-07 2.985758357e-02 -1.276801201e+01 +2.662480089e+00 5.979372364e-04 1.243000641e-06 -5.844423344e-05 +4.662882332e-05 7.282185572e-01 -3.584110205e-04 +1.911236170e-05 1.111648768e-07 7.823789154e-04 3.486592726e-06 1.934888323e-03 -1.738847775e+00 +3.234911367e-01 6.783546728e-04 -4.470619232e+02 +1.160451138e+03 1.500005936e-05 -3.990822914e-03 +2.239842727e-04 5.562357380e-03 -1.002281688e-05 -2.692339517e-03 +5.621101723e-06 +1.367336953e-04 4.261028205e-02 1.835360553e-05 -8.668223263e-04 -2.014791060e+00 -2.248997227e-06 +2.114757082e-04 +1.798949264e+01 +2.394741875e-07 1.986324999e-04 6.578850492e-05 4.502179572e+02 8.855059478e-05 3.237727416e-03 -3.182235485e-07 +1.093420259e-07 1.106011705e-01 1.143715836e-07 5.522698510e-06 6.890923645e-06 -2.672771798e-03 -5.045078327e-03 +1.765912134e-02 +7.251891194e-04 2.462746837e-05 3.856346551e-08 5.216536125e-06 -1.050167951e-04 -1.367775165e+04 +6.007013755e-06 +2.369775354e+03 1.186670853e+02 -3.979706551e+02 -2.281946477e-01 +2.543801550e+02 +2.682901038e-01 2.351156237e-06 -1.048852232e-04 +7.769821656e-05 4.845925968e-03 4.408732612e-01 -1.112548897e-04 +7.856962492e-05 2.011543913e-06 2.266763426e-07 3.059571023e-05 -3.689019494e-01 -7.844255620e-06 +1.677396858e-01 +1.959921970e-06 5.025542164e-03 5.486520955e-04 -1.075092826e-07 +1.779895054e-08 1.132860013e-05 2.213263103e-06 -4.938232858e-06 -7.616499860e-03 +2.140323057e-06 +1.704817163e-03 5.859877168e-05 1.934315184e-06 -2.304443733e-03 +1.572922025e-03 7.553653114e-03 -9.696050391e-02 -7.884528318e-02 -1.348357805e-03 -1.357842010e-04 +1.569442565e-02 +4.299987069e-02 +7.474807797e-04 +8.944381958e-05 3.141317968e-07 -6.709252930e-07 +2.051394799e-07 2.303965173e-07 1.399944656e-07 3.172292798e-04 3.244924593e-03 3.215239406e-07 5.351823628e-04 -1.165409009e+00 -1.698720114e-05 +2.315013359e-01 +8.415215133e-06 2.551101888e-05 -4.114120001e-06 +2.161795380e-06 4.567066262e-07 3.716863062e-05 6.613977572e-07 2.708492813e-03 1.754571252e-01 -1.525562358e-02 -2.519367322e-05 -1.040323131e-03 -2.111103983e-04 +1.139082627e-02 +6.013044846e-05 +1.159312026e-03 +5.809826011e-05 3.163729413e-06 1.195048391e+03 -9.442817960e-06 +4.939698613e-06 7.139106015e-07 4.701054611e-05 4.016940958e-06 3.392971071e-03 1.324572094e-02 1.809116765e-02 -1.068491608e-04 -8.801093474e-02 +5.955431988e-06 +5.116326723e-02 7.543629353e-07 3.144888884e-02 2.541618224e-07 1.081287052e-06 6.019085856e-04 -4.102207206e-06 -6.950122208e+01 -9.158826406e-06 -3.969242776e-05 -3.461103196e-06 -1.215029690e-01 +2.675294230e-06 +1.733420387e+01 +4.217619832e-06 +1.039839924e-05 +2.674339279e-06 +6.630844226e-03 1.388249177e-07 7.461768963e-03 3.402909074e-06 1.353552111e-08 8.196497731e-04 8.070767847e-08 5.530827384e-07 9.776453894e-07 1.114655452e-04 5.077992189e-03 1.369512011e-02 3.740658941e-10 1.781418456e-09 6.274192504e-10 1.401209043e-08 -1.183617756e-05 +4.516147958e-07 1.097744907e-10 7.681004360e-06 1.642521973e-06 7.268695991e-08 8.422665758e-09 2.289208384e-08 -3.164564910e-07 +2.305737007e-07 7.734269324e-08 1.182053682e-07 6.038027033e-04 1.493672352e-05 1.355625923e-07 -3.951040786e-07 -7.931076516e-07 +2.538809292e-07 +3.845119784e-06 1.730384601e-02 1.844042121e-04 3.085693244e-06 2.260606209e+04 6.813153964e-02 -1.179740446e+00 +5.501077576e-01 7.181881004e-07 8.815710180e-04 3.504735286e-01 1.700680109e-03 6.910992987e-07 1.043157794e+02 1.814172816e-06 -6.263269091e-03 +1.432388062e-02 2.976015172e-04 6.200943181e-03 -2.965807565e-07 -1.073076862e-05 +2.490212037e-07 +5.767096764e-05 5.990564946e-07 2.921459438e-03 1.391998767e-08 1.376391218e-06 2.868305907e-02 1.321869635e-09 7.854832125e+01 3.317601049e-06 2.896302554e-07 -4.433606537e-07 +3.895052741e-08 3.755558047e-04 1.885141400e-06 2.979456575e-07 2.644454811e-02 4.727532404e-04 -1.505876650e-03 +8.453223851e-05 2.507321344e-08 2.401908115e-11 4.336051545e-07 1.411161765e-05 8.172635648e-07 2.133943837e-07 8.472465938e-07 2.378645368e-06 6.373870068e-02 1.839733051e-07 2.196405434e-08 6.488703749e-10 1.367098502e-07 -5.025676102e-06 +6.253120021e-07 1.726155379e-05 1.028043800e-04 1.581189528e-02 1.714873385e+00 5.042793140e-10 1.167632099e-05 3.269051929e-07 3.251831260e-10 2.489763266e-04 1.392964403e-06 5.023658342e-01 9.167250740e-05 2.647450934e-02 1.348673997e-02 7.683645278e-07 2.341732570e-04 3.779103872e-04 1.907922842e-08 5.731777278e-07 2.491851220e-09 1.535566443e-04 2.338525694e-11 1.730607356e-05 7.427438437e-09 -9.321434871e-06 +5.487384543e-06 5.356705978e-04 2.577249280e-02 6.339805093e-06 1.722225805e-08 7.949601147e-09 9.467943594e-03 4.124552284e-07 4.752378199e-03 6.855892451e-10 3.030134010e-03 6.111501060e-05 8.896968839e-08 7.036626305e-08 9.767086779e-08 1.641504920e-05 5.167665362e-07 8.929765297e-08 1.041787025e-08 3.214441781e-07 4.278749951e+01 4.713220366e-03 -1.171095095e-05 -5.812379355e-08 +1.257667096e-04 +1.221758281e-07 3.435883898e-06 2.608975632e-07 -2.404695727e-03 +1.030456592e-03 1.308690410e-05 6.442728868e-02 1.785133042e-08 2.426793708e-07 9.869946867e-05 3.931244362e-01 4.427447162e-07 3.012213417e-08 3.937011867e-05 -1.030781975e-05 +6.637002401e-06 1.477893935e-06 4.793521076e-06 2.274482489e-07 2.426949341e-06 2.790263661e-04 1.411695733e-02 1.718655950e-07 1.222952974e+00 1.841651984e+02 3.962193099e-02 7.575329958e-05 2.302985029e-04 2.986165595e-03 8.808106987e-03 5.367867416e-08 3.629488229e+02 6.565633407e-07 2.763279312e-07 8.792666186e-08 4.312444225e-07 2.084090097e-05 5.553313339e-08 1.824845466e-03 1.628190320e+01 8.032792019e+00 5.045965101e-12 1.309316910e-06 -3.447366362e-03 +5.828461023e-04 1.771992742e-05 -1.732654970e-06 +1.305278849e-06 2.253917357e-04 9.285755927e-05 8.707117109e-05 2.458948201e-06 1.264146894e-04 -9.259005878e-09 +3.310491833e-08 2.612269543e-07 8.656036107e-04 -3.092116235e-01 +3.410094197e-02 6.541188422e-04 1.248762684e-07 4.828081916e-07 2.666096671e+02 2.779624128e-09 1.275654743e-11 1.622124438e-09 3.535882708e+01 8.637645168e-05 -4.779015648e-06 +3.502689223e-07 1.985078302e-05 3.769679004e-07 1.893519630e-08 -2.764248045e-02 -1.569774109e+01 +1.538847359e-02 +2.382434533e+01 2.898043591e-10 1.381188536e-06 9.778652394e-07 3.129079728e-04 1.006792969e-03 1.324187996e-05 4.670316575e-03 1.714481366e-02 9.746261540e-07 8.754150483e-06 8.871860946e-06 1.302177141e-03 6.687989253e-07 5.658176860e-11 2.428070629e-01 1.652839606e-07 2.464020929e-04 2.376283608e-03 7.577069770e-09 7.369154255e-12 2.187410343e-02 6.583332200e-09 7.586000565e-08 6.230704627e-04 2.301894399e-07 8.933274543e-07 1.742752087e-07 -3.464376873e-06 +7.893638481e-07 5.505250354e-08 4.664050600e-04 3.406810468e-10 9.055849938e-09 1.358011703e-01 2.511809081e-06 9.330572939e-06 7.547511867e-06 3.643813774e+00 8.070569932e-06 4.812451470e-07 2.258393748e-07 2.233315236e-06 3.159599249e-07 1.438419527e-05 2.656523503e-04 2.401748292e+01 -1.813967533e-08 +1.407146054e-09 4.182507375e-06 1.817768494e-03 4.358306112e-09 3.366166577e+06 9.748286443e-04 3.665246802e-02 5.257456411e-01 -4.936459997e-03 +5.519266935e-03 2.345660154e-06 6.085901852e-02 4.397967866e-04 2.027126088e-08 1.824801845e-04 3.284634641e-04 -6.255318133e-07 +5.554104325e-06 2.940093142e-08 -4.307129122e-02 +8.106348513e-03 8.679944530e-07 1.929082941e+01 3.427794650e-04 1.620709084e-07 4.976512597e-04 4.199111596e-07 6.038191476e-09 1.165738541e-02 5.973500915e-05 1.493338629e-08 3.492347166e-01 1.795227293e-06 1.037966149e-07 1.116776921e-08 -1.263951412e-08 +7.096793504e-09 3.063922251e-03 4.659880760e-06 6.685144842e-11 3.344874212e-06 6.174930345e-02 4.874739894e-09 1.108852165e-06 -3.950363917e-04 -6.789868863e-04 +1.629092705e-04 +3.289616837e-04 2.342719974e-04 4.323244263e-06 1.217044592e-07 1.076613707e-03 8.266921906e-06 1.764292648e-10 6.887713583e-08 6.045824798e+02 1.581575863e-07 4.321419033e-02 2.747231122e-08 5.176464067e-04 1.283947898e-10 2.635146077e-07 1.866099754e-06 1.077625828e-10 1.227034128e-07 3.462488088e-02 6.572734893e-04 3.682533741e-09 1.258102652e-05 1.039023360e-05 9.860921117e-06 2.764566092e-04 9.083830083e-09 -8.019616684e-07 +1.810615225e-07 3.309101491e-05 1.630414134e-11 3.408138743e-03 6.234054635e-08 3.798460345e-04 2.056820068e-04 3.723185048e-11 9.116442122e-09 2.047558303e-08 2.432333844e-08 4.876844945e-09 3.729175624e-11 1.086734095e-07 2.678246345e-07 -3.699977967e-11 +3.281030349e-11 2.518544572e-02 3.543478653e-05 3.144209706e-09 1.038332373e-09 2.905763262e-10 8.845631144e-06 5.299197118e-09 3.212858955e-11 1.574326626e-03 -6.902833399e-07 -2.553781901e-09 +2.425194383e-08 +8.654050737e-10 7.677924584e-10 1.019360586e-02 1.103600437e-06 3.940498178e-06 3.727582277e-09 6.887332609e-09 1.302011294e-09 6.527101708e-12 3.443844857e-06 -4.187145645e-09 +9.983167336e-10 7.948118615e-05 6.284425219e-05 1.649790491e-06 7.173724647e-09 4.153972683e-10 2.405566428e+05 8.355925992e-03 8.745341043e-10 2.269492014e-10 3.309506110e-08 1.413945418e-06 2.232475177e-08 1.631049240e-08 4.510105717e-05 5.086183012e-04 1.988765291e-03 3.062340138e-05 5.993816281e-04 7.223618099e-06 6.919829161e-06 4.084309097e-06 5.317909021e-03 2.390820360e-09 9.926494617e-06 1.203124209e-08 1.587085982e-08 3.209965319e-11 3.730564190e-01 1.294323217e-09 4.525994100e-04 -1.229663713e-07 +1.578300749e-08 8.883022051e-07 7.242134641e-11 1.469712006e-06 2.146381110e-02 1.190600444e-02 4.002908719e-09 8.716239815e-09 1.439376814e-08 8.214348516e-10 7.841363574e-04 2.502302822e-12 2.237540399e-06 -3.060990927e-01 -1.042379315e-05 +1.245214956e-01 +8.068275432e-06 2.747503032e-09 3.152499662e-06 1.213038341e-06 -4.561045660e-09 +1.519089765e-09 2.785775036e-08 3.505026909e-09 5.852871676e-06 6.930894599e-05 2.218608906e-07 8.336025248e-08 5.321324864e-07 6.957006774e-11 1.756920847e-06 9.201475866e-07 1.862585228e-04 2.623905961e-03 6.200680271e-11 1.830818882e-04 7.754292544e-14 1.375087808e-08 1.089787443e-06 4.669668495e-07 7.929650811e-07 2.541277492e+00 7.548645527e-11 3.792553481e-06 1.507840761e-07 3.976718923e-09 8.170369827e-09 1.244907359e-04 2.305047052e-11 1.868748852e-03 1.151371320e-02 -2.616767864e-06 +3.443112060e-07 4.685416671e-11 4.257552002e-10 4.676237105e-05 5.067525573e-04 7.631196519e-09 4.186397884e-13 2.672149954e-03 2.649317492e-08 6.788068043e-08 3.881131771e-11 2.442520598e-01 9.611315772e-08 2.532528157e-03 1.187572603e-11 1.653186265e-09 9.024678356e-08 3.482565421e-10 -6.545561410e-07 +9.913457595e-07 9.067490575e-07 8.135923216e-13 -1.312975615e-07 +9.086122702e-08 5.089372500e-08 3.779620535e-05 1.597856376e-06 4.649833682e-08 3.899419379e-05 2.385486151e-04 -8.272888822e-08 +2.283018442e-07 2.048011836e-10 1.598484276e-05 2.842204256e-06 4.300810814e-10 1.236712825e-09 9.724480523e-10 5.137729447e-06 1.206863564e-03 -3.205098174e-07 +2.997179873e-07 2.968956835e-06 3.892655757e-06 -1.871135757e-08 +8.048486118e-09 9.315475683e-06 4.685745597e-03 1.264393372e-11 1.651517618e-06 4.093555503e-04 6.291217780e-02 2.232241817e-12 1.581588279e-10 2.176076635e-01 2.778841978e-02 1.933575772e-12 1.606840276e-03 -1.889515674e+00 +2.907179979e-01 8.696453002e-01 9.849266934e-02 5.923627073e-11 1.000520546e-11 5.484410290e-02 1.900079973e-04 6.710428134e-08 4.023506025e-07 9.262433958e-07 3.971033910e-04 2.798384077e-11 2.733321071e-04 5.949377984e-06 5.804131920e-04 4.149708488e-09 5.189626516e-08 3.459693749e-09 8.152801222e-12 4.442083122e-08 4.885241039e-07 1.286023795e-11 9.146948664e-12 7.438541958e-07 3.933356062e-11 1.992948079e-08 1.233914164e-01 1.179935767e-06 -1.052253323e-06 +4.500519258e-07 2.614263483e-07 1.305318834e-04 4.563775945e-04 6.061435715e-10 9.124754450e-06 1.181449830e-08 4.457412126e-07 1.105658535e-15 7.403561460e-08 6.857923072e-05 2.595980444e-06 2.642703583e-07 2.275429709e-10 1.608232344e-05 -1.139157860e-04 +2.179365199e-05 1.471558450e-09 3.413025936e-03 4.519315643e-05 3.265308485e-09 7.421214523e-10 3.097172719e-04 4.702159026e-05 8.124970186e-10 6.737184939e-12 3.875912321e-05 8.002679747e-05 4.570485914e-05 1.372447479e-06 9.295972807e-10 9.091336789e-07 2.002021098e-06 4.844911625e-09 1.902794787e-08 1.566263784e-08 9.357485889e-06 1.048780982e-04 6.452709776e-07 4.085398049e-06 1.535693061e-08 9.742658204e-12 2.622006774e-05 2.375471653e-02 3.141422539e-03 5.474762658e-08 1.414376429e-03 1.609237488e-12 1.107110125e+01 7.034088451e-04 3.105924621e-07 5.029108582e-09 2.207617143e-08