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 <H1DD>,...,<H12DD>;
+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<double> 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<double> 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 063f06f..544d1b7 100644
--- a/src/Hjets.cc
+++ b/src/Hjets.cc
@@ -1,1006 +1,478 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/jets.hh"
#include "HEJ/Hjets.hh"
#include <array>
#include <cassert>
#include <cmath>
#include <complex>
#include <limits>
#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<double>::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<std::complex<double>> result(3);
static auto ql_B0 = [](){
ql::Bubble<std::complex<double>,double,double> ql_B0;
ql_B0.setCacheSize(100);
return ql_B0;
}();
static std::vector<double> masses(2);
static std::vector<double> 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<std::complex<double>> result(3);
static auto ql_C0 = [](){
ql::Triangle<std::complex<double>,double,double> ql_C0;
ql_C0.setCacheSize(100);
return ql_C0;
}();
static std::vector<double> masses(3);
static std::vector<double> 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 const & q3, double mq)
- {
- static std::vector<std::complex<double>> result(3);
- static auto ql_D0 = [](){
- ql::Box<std::complex<double>,double,double> ql_D0;
- ql_D0.setCacheSize(100);
- return ql_D0;
- }();
- static std::vector<double> masses(4);
- static std::vector<double> 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<COM, 2> 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<COM, 2> res = {qH1.dot(qH2), 1.};
for(auto & tt: res) tt /= (3.*M_PI*vev);
return res;
}
std::array<COM, 2> 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<minus, minus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_mp = HEJ::j_h_j<minus, plus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_pm = HEJ::j_h_j<plus, minus>(
p1out, p1in, p2out, p2in, qH11, qH12, qH21, qH22, T_pref[0], T_pref[1]
);
const COM amp_pp = HEJ::j_h_j<plus, plus>(
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<minus, minus>(
pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
);
const COM amp_mp = HEJ::jgh_j<minus, plus>(
pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
);
const COM amp_pm = HEJ::jgh_j<plus, minus>(
pg, pb, pn, ph1, ph2, T_pref[0], T_pref[1]
);
const COM amp_pp = HEJ::jgh_j<plus, plus>(
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<Helicity h1, Helicity h2, Helicity hg>
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<COM, 2> const & T_pref
) {
// TODO: code duplication with Wjets and pure jets
const COM u1 = U1_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
const COM u2 = U2_h_j<h1, h2, hg>(pa,p1,pb,p2,pg,qH11,qH12,qH21,qH22,T_pref[0],T_pref[1]);
const COM l = L_h_j<h1, h2, hg>(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<double, 2, 2, 2> amp{};
// NOLINTNEXTLINE
#define ASSIGN_HEL(RES, J, H1, H2, HG) \
RES[H1][H2][HG] = J<H1, H2, HG>( \
p1in, p2in, pg, p1out, p2out, qH11, qH12, qH21, qH22, T_pref \
)
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 const q2=-(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const Delta = s12*s34 - S1*S2;
- double const 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 const q2 = -(k1+k2+kh);
- double const S1 = 2.*k1.dot(q2);
- double const S2 = 2.*k2.dot(q2);
- double const s12 = 2.*k1.dot(k2);
- double const s34 = q2.m2();
- double const 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<Helicity hout, Helicity h2>
- 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<hout, h2>(pa, p1, pb, p2, ph1, ph2, mq):
- jh_j_backward<hout, h2>(pa, p1, pb, p2, ph1, ph2, mq)
- ;
- }
-
- template<Helicity hout, Helicity h2>
- 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<hout, h2>(pa, p1, pb, p2, ph1, ph2, mq);
- if(include_bottom) {
- res += 4.*mq2*mq2/vev*jh_j<hout, h2>(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<plus, plus>(
- pa, p1, pb, p2, ph1, ph2,
- mq, include_bottom, mq2, vev
- );
- const COM appm = amp_jh_j<plus, minus>(
- pa, p1, pb, p2, ph1, ph2,
- mq, include_bottom, mq2, vev
- );
- const COM apmp = amp_jh_j<minus, plus>(
- pa, p1, pb, p2, ph1, ph2,
- mq, include_bottom, mq2, vev
- );
- const COM apmm = amp_jh_j<minus, minus>(
- 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 mt, const bool include_bottom, const double mb, 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,
- mt, include_bottom, mb, 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),
- mt, include_bottom, mb, 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