diff --git a/current_generator/include/currents.frm b/current_generator/include/currents.frm index 9f9a68d..c2cfd6c 100644 --- a/current_generator/include/currents.frm +++ b/current_generator/include/currents.frm @@ -1,134 +1,66 @@ */** * \brief Predefined current functions and polarisation vectors * * \authors The HEJ collaboration (see AUTHORS for details) * \date 2020 * \copyright GPLv2 or later */ #ifndef `$HEJCurrentsIncluded' #$HEJCurrentsIncluded = 1; -* Extremal Higgs boson emission vertex -cfunction JH; * FKL current for vector boson cfunction JV; * Polarisation vectors cfunction Eps; * Polarisation vectors contracted with Higgs emission vertex cfunction EPSH; * Higgs emission vertex cfunction VH; cfunction m2,plus,minus,perphat,sqrt,conj,dot; -cf T1,T2; -cf ,...,; +s T1,T2; * internal symbols are all uppercase symbols H, H1, HL, M; indices MU1,...,MU50; vectors PA,P1,PH,PLBAR,PL,PG,PR,Q1,Q2; #procedure InsSpecialCurrents #call InsJV - #call InsJH #call InsEps #call InsVH #endprocedure * Replace vector boson FKL current #procedure InsJV() * we use a loop here to introduce a new dummy index for each occurence repeat; once JV(H1?, HL?, MU1?, PA?, P1?, PLBAR?, PL?) = Current(HL, PL, MU2, PLBAR)*( + Current(H1, P1, MU2, P1+PLBAR+PL, MU1, PA)/m2(P1+PLBAR+PL) + Current(H1, P1, MU1, PA-PLBAR-PL, MU2, PA)/m2(PA-PLBAR-PL) ); sum MU2; endrepeat; #endprocedure * Replace polarisation vectors by currents #procedure InsEps() id Conjugate(Eps(H1?, MU1?, PG?, PR?)) = - Eps(-H1, MU1, PG, PR); id Eps(-1, MU1?, PG?, PR?) = sqrt(2)/2*SpinorChain(PR, MU1, PG)/AngleChain(PG,PR); id Eps(+1, MU1?, PG?, PR?) = sqrt(2)/2*SpinorChain(PG, MU1, PR)/SquareChain(PG,PR); #endprocedure * Higgs emission vertex * see eq:VH in developer manual, but without overall prefactor #procedure InsVH() id VH(MU1?, MU2?, Q1?, Q2?) = d_(MU1, MU2)*T1 - Q2(MU1)*Q1(MU2)*T2; #endprocedure -* Higgs currents -* see eq:jH_same_helicity and jH_helicity_flip in developer manual -* the overall prefactor is omitted -#procedure InsJH() - - id JH(+1, +1, MU1?, P1?, PA?, PH?, M?) = - conj(perphat(P1))*m2(PA - P1 - PH)*( - - sqrt(2*minus(P1)/minus(PA))/square(PA, P1)*Conjugate(EPSH(+1, MU1, P1, PH, PA, M)) - + sqrt(2*minus(PA)/minus(P1))/angle(P1, PA)*EPSH(+1, MU1, PA, PH, P1, M) - ) - + sqrt(2)*Current(+1, P1, PH, PA)*( - + EPSH(+1, MU1, PA, PH, P1, M)/angle(P1, PA) - + Conjugate(EPSH(+1, MU1, P1, PH, PA, M))/square(PA, P1) - - angle(P1, PA)/(2*m2(PA - PH))*T2(PA, PA-PH, M)*Conjugate(Eps(+1, MU1, P1, PA)) - - square(PA, P1)/(2*m2(P1 + PH))*T2(P1+PH, P1, M)*Eps(+1, MU1, PA, P1) - - 4*pi_^2*i_*H4DD(P1, -PA, PH, M)/square(PA, P1)*Conjugate(Eps(+1, MU1, P1, PA)) - + 4*pi_^2*i_*H5DD(P1, -PA, PH, M)/angle(P1, PA)*Eps(+1, MU1, PA, P1) - ) - - 4*pi_^2*i_*Current(+1, P1, PH, PA)^2/m2(PA-P1)*( - PA(MU1)*H10DD(P1, -PA, PH, M) - P1(MU1)*H12DD(P1, -PA, PH, M) - ) - ; - -* TODO: -* signs of negative-helicity polarisation vectors -* and of the terms proportional to H10, H12 -* differ from developer manual and arXiv:1812.08072 - id JH(+1, -1, MU1?, P1?, PA?, PH?, M?) = - m2(PA - P1 - PH)/square(PA, P1)*( - + sqrt(2*minus(P1)/minus(PA))*conj(perphat(P1))*Conjugate(EPSH(-1, MU1, P1, PH, PA, M)) - + sqrt(2*minus(PA)/minus(P1))*perphat(P1)*EPSH(+1, MU1, PA, PH, P1, M) - ) - + sqrt(2)*Current(+1, P1, PH, PA)*( - - Conjugate(EPSH(-1, MU1, P1, PH, PA, M))/square(PA, P1) - + angle(P1, PA)/(2*m2(PA - PH))*T2(PA, PA-PH, M)*Conjugate(Eps(-1, MU1, P1, PA)) - + 4*pi_^2*i_*H4DD(P1, -PA, PH, M)/square(PA, P1)*Conjugate(Eps(-1, MU1, P1, PA)) - ) - + sqrt(2)*Current(+1, PA, PH, P1)*( - + EPSH(+1, MU1, PA, PH, P1, M)/square(PA, P1) - - angle(P1, PA)/(2*m2(P1 + PH))*T2(P1+PH, P1, M)*Eps(+1, MU1, PA, P1) - + 4*pi_^2*i_*H5DD(P1, -PA, PH, M)/square(PA, P1)*Eps(+1, MU1, PA, P1) - ) - + 4*pi_^2*i_*Current(+1, P1, PH, PA)*Current(+1, PA, PH, P1)/square(PA, P1)^2*( - PA(MU1)*H10DD(P1, -PA, PH, M) - P1(MU1)*H12DD(P1, -PA, PH, M) - ) - + angle(P1, PA)/square(PA, P1)*( - + 8*pi_^2*i_*H1DD(P1, -PA, PH, M)*P1(MU1) - - 8*pi_^2*i_*H2DD(P1, -PA, PH, M)*PA(MU1) - + 2*dot(P1, PH)*T2(P1+PH, P1, M)/m2(P1+PH)*PA(MU1) - - 2*dot(PA, PH)*T2(PA, PA-PH, M)/m2(PA-PH)*P1(MU1) - + T1(PA-P1, PA-P1-PH, M)/m2(PA-P1)*(P1(MU1) + PA(MU1)) - - dot(P1+PA, PH)/m2(PA-P1)*T2(PA-P1, PA-P1-PH, M)*(P1(MU1) - PA(MU1)) - ); - -* eq:eps_H in developer manual - id EPSH(H1?, MU1?, PA?, PH?, P1?, M?) = T2(PA, PA-PH, M)/m2(PA - PH)*( - dot(PA,PH)*Eps(H1, MU1, PA, P1) - Eps(H1, PH, PA, P1)*PA(MU1) - ); - id Conjugate(EPSH(H1?, MU1?, P1?, PH?, PA?, M?)) = -T2(P1+PH, P1, M)/m2(P1 + PH)*( - dot(P1,PH)*Conjugate(Eps(H1, MU1, P1, PA)) - Conjugate(Eps(H1, PH, P1, PA))*P1(MU1) - ); - -#endprocedure - #endif diff --git a/current_generator/jh_j.frm b/current_generator/jh_j.frm deleted file mode 100644 index caf2252..0000000 --- a/current_generator/jh_j.frm +++ /dev/null @@ -1,77 +0,0 @@ -*/** -* \brief Contraction of extremal Higgs boson emission current with FKL current -* -* \authors The HEJ collaboration (see AUTHORS for details) -* \date 2020 -* \copyright GPLv2 or later -* -*/ -#include- include/helspin.frm -#include- include/write.frm - -v pb,p2,pa,p1,ph,ph1,ph2; -i mu; -s mq; - -* backward (p1.z < 0) -* the first helicity corresponds to the outgoing gluon, -* where - corresponds to a helicity flip -* h2 is the helicity of the current without Higgs emission -#do h2={+,-} - l [jh_j_backward +`h2'] = ( - JH(+1, +1, mu, p1, pa, ph, mq) - *Current(`h2'1, p2, mu, pb) - ); - l [jh_j_backward -`h2'] = ( - JH(+1, -1, mu, p1, pa, ph, mq) - *Current(`h2'1, p2, mu, pb) - ); -#enddo - -multiply replace_( - ph, ph1 + ph2 -); - -#call ContractCurrents -.sort -skip; - -* forward (p1.z > 0) -#do hout={+,-} - #do h2={+,-} - l [jh_j_forward `hout'`h2'] = [jh_j_backward `hout'`h2']; - #enddo -#enddo - -multiply replace_(minus,plus); -id perphat(p1) = -1; -id conj(perphat(p1)) = -1; -.sort - -format O4; -format c; -#call WriteHeader(`OUTPUT') -#write <`OUTPUT'> " namespace currents {" -#write <`OUTPUT'> " namespace {" -#do i=1,2 - #write <`OUTPUT'> " std::complex T`i'(" - #write <`OUTPUT'> " CLHEP::HepLorentzVector const & q1," - #write <`OUTPUT'> " CLHEP::HepLorentzVector const & q2," - #write <`OUTPUT'> " double m" - #write <`OUTPUT'> " );" -#enddo -#do i={1,2,4,5,10,12} - #write <`OUTPUT'> " std::complex H`i'DD(" - #write <`OUTPUT'> " CLHEP::HepLorentzVector const & k1," - #write <`OUTPUT'> " CLHEP::HepLorentzVector const & k2," - #write <`OUTPUT'> " CLHEP::HepLorentzVector const & kh," - #write <`OUTPUT'> " double mq" - #write <`OUTPUT'> " );" -#enddo -#write <`OUTPUT'> " }" -#write <`OUTPUT'> " }" -#do DIR={backward,forward} - #call WriteOptimised(`OUTPUT',jh_j_`DIR',2,pa,p1,pb,p2,ph1,ph2,DOUBLE:,mq) -#enddo -#call WriteFooter(`OUTPUT') -.end diff --git a/include/HEJ/Hjets.hh b/include/HEJ/Hjets.hh index 9049807..676d3ef 100644 --- a/include/HEJ/Hjets.hh +++ b/include/HEJ/Hjets.hh @@ -1,431 +1,360 @@ /** * \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/src/Hjets.cc b/src/Hjets.cc index 1036644..9513739 100644 --- a/src/Hjets.cc +++ b/src/Hjets.cc @@ -1,1013 +1,477 @@ /** * \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