Page MenuHomeHEPForge
Files F10881898

No OneTemporary
Actions

View Options

diff --git a/include/RHEJ/currents.hh b/include/RHEJ/currents.hh
index 435650e..f426d73 100644
--- a/include/RHEJ/currents.hh
+++ b/include/RHEJ/currents.hh
@@ -1,1324 +1,1349 @@
//////////////////////////////////////////////////
//////////////////////////////////////////////////
// This source code is Copyright (2012) of //
// Jeppe R. Andersen and Jennifer M. Smillie //
// and is distributed under the //
// Gnu Public License version 2 //
// http://www.gnu.org/licenses/gpl-2.0.html //
// You are allowed to distribute and alter the //
// source under the conditions of the GPLv2 //
// as long as this copyright notice //
// is unaltered and distributed with the source //
// Any use should comply with the //
// MCNET GUIDELINES //
// for Event Generator Authors and Users //
// as distributed with this source code //
//////////////////////////////////////////////////
//////////////////////////////////////////////////
/** \file
* \brief Functions computing the square of current contractions.
*
* This file contains all the necessary functions to compute the current
* contractions for all valid HEJ processes. PJETS, H+JETS and W+JETS along with
* some unordered counterparts.
*/
#pragma once
#include <CLHEP/Vector/LorentzVector.h>
#include <complex>
#include <vector>
#include <valarray>
#include <limits>
///TODO add a namespace
typedef std::complex<double> COM;
typedef COM current[4];
typedef CLHEP::HepLorentzVector HLV;
//! The Higgs field vacuum expectation value in GeV
static constexpr double v = 246.;
constexpr double infinity = std::numeric_limits<double>::infinity();
constexpr double mb_default = 4.7;
void Setup_Currents(void);
//! Square of qQ->qenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qQ->qenuQ Scattering
*
* This returns the square of the current contractions in qQ->qenuQ scattering
* with an emission of a W Boson.
*/
double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qbarQ->qbarenuQ W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state quark
* @param p2in Momentum of intial state quark
* @returns Square of the current contractions for qbarQ->qbarenuQ Scattering
*
* This returns the square of the current contractions in qbarQ->qbarenuQ scattering
* with an emission of a W Boson.
*/
double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qQbar->qenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qQbar->qenuQbar Scattering
*
* This returns the square of the current contractions in qQbar->qenuQbar scattering
* with an emission of a W Boson.
*/
double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qbarQbar->qbarenuQbar W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state anti-quark
* @param p2in Momentum of intial state anti-quark
* @returns Square of the current contractions for qbarQbar->qbarenuQbar Scattering
*
* This returns the square of the current contractions in qbarQbar->qbarenuQbar scattering
* with an emission of a W Boson.
*/
double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qg->qenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qg->qenug Scattering
*
* This returns the square of the current contractions in qg->qenug scattering
* with an emission of a W Boson.
*/
double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qbarg->qbarenug W+Jets Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark
* @param p2out Momentum of final state gluon
* @param p2in Momentum of intial state gluon
* @returns Square of the current contractions for qbarg->qbarenug Scattering
*
* This returns the square of the current contractions in qbarg->qbarenug scattering
* with an emission of a W Boson.
*/
double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe,
CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
// W+Jets Unordered Functions
//! qQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in qQg->qQg scattering
* with an emission of a W Boson.
*/
double junobMWqQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg);
//! qbarQg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @param pg Momentum of final state unordered gluon
* @returns Square of the current contractions for qbarQ->qbarQg Scattering
*
* This returns the square of the current contractions in qbarQg->qbarQg scattering
* with an emission of a W Boson.
*/
double junobMWqbarQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg);
//! qQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @returns Square of the current contractions for qQbar->qQbarg Scattering
*
* This returns the square of the current contractions in qQbarg->qQbarg scattering
* with an emission of a W Boson.
*/
double junobMWqQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg);
//! qbarQbarg Wjets Unordered backwards opposite leg to W
/**
* @param p1out Momentum of final state anti-quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @param pg Momentum of final state unordered gluon
* @returns Square of the current contractions for qbarQbar->qbarQbarg Scattering
*
* This returns the square of the current contractions in qbarQbarg->qbarQbarg scattering
* with an emission of a W Boson.
*/
double junobMWqbarQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg);
//!Wjets Unordered forwards opposite leg to W
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @returns Square of the current contractions for qQ->gqQ Scattering
*
* This returns the square of the current contractions in qQg->gqQ scattering
* with an emission of a W Boson.
*/
double junofMWgqQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in);
//!Wjets Unordered forwards opposite leg to W
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state anti-quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This returns the square of the current contractions in qbarQg->gqbarQ scattering
* with an emission of a W Boson.
*/
double junofMWgqbarQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in);
//!Wjets Unordered forwards opposite leg to W
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This returns the square of the current contractions in qQbarg->gqQbar scattering
* with an emission of a W Boson.
*/
double junofMWgqQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in);
//!Wjets Unordered forwards opposite leg to W
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state anti-quark a
* @param pe Momentum of final state electron
* @param pnu Momentum of final state Neutrino
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This returns the square of the current contractions in qbarQbarg->gqbarQbar scattering
* with an emission of a W Boson.
*/
double junofMWgqbarQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in);
//!W+uno same leg
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in gqQ->gqQ scattering
* with an emission of a W Boson.
*/
double jM2WunogqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! TODO: What does this function do? Crossed contribution is Exqqx..?
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @returns Square of the current contractions for qQ->gqQ Scattering
*
* This returns the square of the current contractions in gqQ->gqQ scattering
* with an emission of a W Boson.
*/
double jM2WunogqQ_crossqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! W+uno same leg. quark anti-quark
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This returns the square of the current contractions in gqQbar->gqQbar scattering
* with an emission of a W Boson. (Unordered Same Leg)
*/
double jM2WunogqQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! W+uno same leg. quark gluon
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @returns Square of the current contractions for qg->gqg Scattering
*
* This returns the square of the current contractions in qg->gqg scattering
* with an emission of a W Boson.
*/
double jM2Wunogqg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! W+uno same leg. anti-quark quark
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state anti-quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state quark b
* @param p2in Momentum of intial state quark b
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This returns the square of the current contractions in qbarQ->gqbarQ scattering
* with an emission of a W Boson.
*/
double jM2WunogqbarQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! W+uno same leg. anti-quark anti-quark
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state anti-quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state anti-quark b
* @param p2in Momentum of intial state anti-quark b
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This returns the square of the current contractions in gqbarQbar->qbarQbar scattering
* with an emission of a W Boson.
*/
double jM2WunogqbarQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! W+uno same leg. anti-quark gluon
/**
* @param pg Momentum of final state unordered gluon
* @param p1out Momentum of final state anti-quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param p1in Momentum of initial state anti-quark a
* @param p2out Momentum of final state gluon b
* @param p2in Momentum of intial state gluon b
* @returns Square of the current contractions for ->gqbarg Scattering
*
* This returns the square of the current contractions in qbarg->gqbarg scattering
* with an emission of a W Boson.
*/
double jM2Wunogqbarg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//W+Jets qqxExtremal
//! W+Extremal qqx. qxqQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for ->qbarqQ Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2WgQtoqbarqQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//W+Jets qqxExtremal
//! W+Extremal qqx. qqxQ
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state anti-quark b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for ->qqbarQ Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2WgQtoqqbarQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//W+Jets qqxExtremal
//! W+Extremal qqx. gg->qxqg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon b
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for gg->qbarqg Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2Wggtoqbarqg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//W+Jets qqxExtremal
//! W+Extremal qqx. gg->qqxg
/**
* @param pgin Momentum of initial state gluon
* @param pqout Momentum of final state quark a
* @param plbar Momentum of final state anti-lepton
* @param pl Momentum of final state lepton
* @param pqbarout Momentum of final state anti-quark a
* @param p2out Momentum of initial state gluon a
* @param p2in Momentum of final state gluon b
* @returns Square of the current contractions for gg->qqbarg Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2Wggtoqqbarg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqbarout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
+//W+Jets qqxExtremal, W emission from opposite leg
+//! W+Extremal qqx. gg->qqxg. qqx on forwards leg, W emission backwards leg.
+/**
+ * @param pa Momentum of initial state (anti-)quark
+ * @param pb Momentum of initial state gluon
+ * @param p1 Momentum of final state (anti-)quark (after W emission)
+ * @param p2 Momentum of final state anti-quark
+ * @param p3 Momentum of final state quark
+ * @param plbar Momentum of final state anti-lepton
+ * @param pl Momentum of final state lepton
+ * @returns Square of the current contractions for gq->qqbarqW Scattering
+ *
+ * Calculates the square of the current contractions with extremal qqbar pair
+ * production. This is calculated via current contraction of existing currents.
+ * Assumes qqx split from forwards leg, W emission from backwards leg.
+ * Switch input (pa<->pb, p1<->p3) if calculating forwards qqx.
+ */
+double jM2WgqtoQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl);
+
//! W+Jets qqxCentral. qqx W emission.
/**
* @param pa Momentum of initial state particle a
* @param pb Momentum of initial state particle b
* @param pl Momentum of final state lepton
* @param plbar Momentum of final state anti-lepton
* @param partons Vector of outgoing parton momenta
* @param aqlinepa Bool: True= pa is anti-quark
* @param aqlinepb Bool: True= pb is anti-quark
* @param qqxmarker Bool: Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Number of lipatov vertices "above" qqbar pair
* @param nbelow Number of lipatov vertices "below" qqbar pair
* @returns Square of the current contractions for qq>qQQbarWq Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2WqqtoqQQq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove, int nbelow); //Doing
//emission from backwards leg
//! W+Jets qqxCentral. W emission from backwards leg.
/**
* @param ka HLV: Momentum of initial state particle a
* @param kb HLV: Momentum of initial state particle b
* @param pl HLV: Momentum of final state lepton
* @param plbar HLV: Momentum of final state anti-lepton
* @param partons Vector(HLV): outgoing parton momenta
* @param aqlinepa Bool: True= pa is anti-quark
* @param aqlinepb Bool: True= pb is anti-quark
* @param qqxmarker Bool: Ordering of the qqbar pair produced (qqx vs qxq)
* @param nabove Int: Number of lipatov vertices "above" qqbar pair
* @param nbelow Int: Number of lipatov vertices "below" qqbar pair
* @param forwards Bool: Swap to emission off front leg TODO:remove so args can be const
* @returns Square of the current contractions for qq>qQQbarWq Scattering
*
* Calculates the square of the current contractions with extremal qqbar pair
* production. This is calculated through the use of crossing symmetry.
*/
double jM2WqqtoqQQqW(HLV ka, HLV kb,HLV pl,HLV plbar, std::vector<HLV> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove, int nbelow, bool forwards); //Doing
//! Square of qQ->qQ Pure 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
* @returns Square of the current contractions for qQ->qQ Scattering
*
* This returns the square of the current contractions in qQ->qQ Pure Jet Scattering.
*/
double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qQbar->qQbar Pure 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
* @returns Square of the current contractions for qQbar->qQbar Scattering
*
* This returns the square of the current contractions in qQbar->qQbar Pure Jet Scattering.
* Note this can be used for qbarQ->qbarQ Scattering by inputting arguments appropriately.
*/
double jM2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qbarQbar->qbarQbar Pure 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
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*
* This returns the square of the current contractions in qbarQbar->qbarQbar Pure Jet Scattering.
*/
double jM2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qg->qg Pure 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
* @returns Square of the current contractions for qg->qg Scattering
*
* This returns the square of the current contractions in qg->qg Pure Jet Scattering.
* Note this can be used for gq->gq Scattering by inputting arguments appropriately.
*/
double jM2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of qbarg->qbarg Pure 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
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* This returns the square of the current contractions in qbarg->qbarg Pure Jet Scattering.
* Note this can be used for gqbar->gqbar Scattering by inputting arguments appropriately.
*/
double jM2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! Square of gg->gg Pure 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
* @returns Square of the current contractions for gg->gg Scattering
*
* This returns the square of the current contractions in gg->gg Pure Jet Scattering.
*/
double jM2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in);
//! 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 q1 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
* @returns Square of the current contractions for gg->gg Scattering
*
* This returns the square of the current contractions in gg->gg Higgs+Jet Scattering.
*
* g~p1 g~p2
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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
* @returns Square of the current contraction
*
*/
double MH2gq_outsideH(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector pH,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qg->qg Scattering
*
* This returns the square of the current contractions in qg->qg Higgs+Jet Scattering.
*
* q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qbarg->qbarg Scattering
*
* This returns the square of the current contractions in qbarg->qbarg Higgs+Jet Scattering.
*
* qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if g is backward, q1 is forward)
*/
double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qQ->qQ Scattering
*
* This returns the square of the current contractions in qQ->qQ Higgs+Jet Scattering.
*
* q~p1 Q~p2 (i.e. ALWAYS p1 for quark, p2 for quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Q is backward, q1 is forward)
*/
double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qQ->qQ Scattering
*
* This returns the square of the current contractions in qQbar->qQbar Higgs+Jet Scattering.
*
* q~p1 Qbar~p2 (i.e. ALWAYS p1 for quark, p2 for anti-quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Qbar is backward, q1 is forward)
*/
double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qbarQ->qbarQ Scattering
*
* This returns the square of the current contractions in qbarQ->qbarQ Higgs+Jet Scattering.
*
* qbar~p1 Q~p2 (i.e. ALWAYS p1 for anti-quark, p2 for quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Q is backward, q1 is forward)
*/
double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! 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 q1 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
* @returns Square of the current contractions for qbarQbar->qbarQbar Scattering
*
* This returns the square of the current contractions in qbarQbar->qbarQbar Higgs+Jet Scattering.
*
* qbar~p1 Qbar~p2 (i.e. ALWAYS p1 for anti-quark, p2 for anti-quark)
* should be called with q1 meant to be contracted with p2 in first part of vertex
* (i.e. if Qbar is backward, q1 is forward)
*/
double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in,
CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
// Unordered f
//! Square of qQ->gqQ Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered 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
* @returns Square of the current contractions for qQ->gqQ Scattering
*
* This returns the square of the current contractions in qQ->gqQ Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qQbar->gqQbar Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered 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
* @returns Square of the current contractions for qQbar->gqQbar Scattering
*
* This returns the square of the current contractions in qQbar->gqQbar Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qbarQ->gqbarQ Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered 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
* @returns Square of the current contractions for qbarQ->gqbarQ Scattering
*
* This returns the square of the current contractions in qbarQ->gqbarQ Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qbarQbar->gqbarQbar Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered 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
* @returns Square of the current contractions for qbarQbar->gqbarQbar Scattering
*
* This returns the square of the current contractions in qbarQbar->gqbarQbar Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qg->gqg Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @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
* @returns Square of the current contractions for qg->gqg Scattering
*
* This returns the square of the current contractions in qg->gqg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qbarg->gqbarg Higgs+Jets Unordered f Scattering Current
/**
* @param pg Momentum of unordered gluon
* @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
* @returns Square of the current contractions for qbarg->gbarg Scattering
*
* This returns the square of the current contractions in qbarg->gqbarg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: pg > p1out >> p2out
*/
double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,
CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//Unordered b
//! Square of qbarQ->qbarQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for qbarQ->qbarQg Scattering
*
* This returns the square of the current contractions in qbarQ->qbarQg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qQ->qQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for qQ->qQg Scattering
*
* This returns the square of the current contractions in qQ->qQg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qQbar->qQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state quark
* @param p1in Momentum of initial state quark
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for qQbar->qQbarg Scattering
*
* This returns the square of the current contractions in qQbar->qQbarg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of qbarQbar->qbarQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state anti-quark
* @param p1in Momentum of initial state anti-quark
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for qbarQbar->qbarQbarg Scattering
*
* This returns the square of the current contractions in qbarQbar->qbarQbarg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of gQbar->gQbarg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for gQbar->gQbarg Scattering
*
* This returns the square of the current contractions in gQbar->gQbarg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
//! Square of gQ->gQg Higgs+Jets Unordered b Scattering Current
/**
* @param p1out Momentum of final state gluon
* @param p1in Momentum of initial state gluon
* @param pg Momentum of unordered b 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
* @returns Square of the current contractions for gQ->gQg Scattering
*
* This returns the square of the current contractions in gQ->gQg Higgs+Jet Scattering.
*
* This construction is taking rapidity order: p1out >> p2out > pg
*/
double jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in,
CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out,
CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1,
CLHEP::HepLorentzVector qH2,
double mt = infinity,
bool include_bottom = false, double mb = mb_default);
// impact factors for Higgs + jet
//! Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.22) in Hep-ph/0301013 with modifications
*
* This gives the impact factor. First it determines first whether this is the case
* p1p\sim php>>p3p or the opposite
*/
double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector pH);
//! Implements Eq. (4.23) in hep-ph/0301013 with modifications to incoming plus momenta
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.23) in Hep-ph/0301013
*
* This gives the impact factor. First it determines first whether this is the case
* p1p\sim php>>p3p or the opposite
*/
double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector pH);
//! Implements Eq. (4.22) in hep-ph/0301013
/**
* @param p2 Momentum of Particle 2
* @param p1 Momentum of Particle 1
* @param pH Momentum of Higgs
* @returns Value of Eq. (4.22) in Hep-ph/0301013
*
* This gives the impact factor. First it determines first whether this is the case
* p1p\sim php>>p3p or the opposite
*/
double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1,
CLHEP::HepLorentzVector pH);
/** \class CCurrent currents.hh "include/RHEJ/currents.hh"
* \brief This is the a new class structure for currents.
*/
class CCurrent
{
public:
CCurrent(COM sc0, COM sc1, COM sc2, COM sc3)
:c0(sc0),c1(sc1),c2(sc2),c3(sc3)
{};
CCurrent(const CLHEP::HepLorentzVector p)
{
c0=p.e();
c1=p.px();
c2=p.py();
c3=p.pz();
};
CCurrent()
{};
CCurrent operator+(const CCurrent& other);
CCurrent operator-(const CCurrent& other);
CCurrent operator*(const double x);
CCurrent operator*(const COM x);
CCurrent operator/(const double x);
CCurrent operator/(const COM x);
friend std::ostream& operator<<(std::ostream& os, const CCurrent& cur);
COM dot(CLHEP::HepLorentzVector p1);
COM dot(CCurrent p1);
COM c0,c1,c2,c3;
private:
};
/* std::ostream& operator <<(std::ostream& os, const CCurrent& cur); */
CCurrent operator * ( double x, CCurrent& m);
CCurrent operator * ( COM x, CCurrent& m);
CCurrent operator / ( double x, CCurrent& m);
CCurrent operator / ( COM x, CCurrent& m);
//! Current ???
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
void j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin,current &cur);
//! Current ???
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur);
//! Current ???
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur);
+//! Current ???
+/**
+ * These functions are a mess. There are many more defined in the source file than declared in the
+ * header - and the arguments are mislabelled in some cases. Need to investigate.
+ */
+void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur);
//! Current ???
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);
//! Current <incoming state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
CCurrent jio (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);
//! Current <outgoing state | mu | outgoing state>
/**
* These functions are a mess. There are many more defined in the source file than declared in the
* header - and the arguments are mislabelled in some cases. Need to investigate.
*/
CCurrent joo (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin);
/* // Coupling values */
/* const double stw2 = 0.2222; */
/* const double ctw = sqrt(1.0 - stw2); */
/* const double gs = 1.217716; */
/* const double gw = 0.653232911; */
/* const double Zem = (-1.0 / 2.0 + stw2) / ctw; */
/* const double Zep = stw2 / ctw; */
/* const double Zum = ( 1.0 / 2.0 - 2.0 * stw2 / 3.0) / ctw; */
/* const double Zup = - 2.0 * stw2 / 3.0 / ctw; */
/* const double Zdm = (-1.0 / 2.0 + 1.0 / 3.0 * stw2) / ctw; */
/* const double Zdp = stw2 / 3.0 / ctw; */
/* const double RWeak = -pow(gw, 2.0); */
/* const double Strong = pow(gs, 4.0); */
/* const double ee = pow(gw, 2.0) * stw2; */
/* std::vector <double> jMZqQ (HLV, HLV, HLV, HLV, HLV, HLV, std::vector <double>, std::vector < std::vector <double> >, int, int, bool, bool); */
/* std::vector <double> jMZqg (HLV, HLV, HLV, HLV, HLV, HLV, std::vector <double>, std::vector < std::vector <double> >, int, int, bool, bool); */
/* void jZ (HLV, HLV, HLV, HLV, bool, bool, current); */
/* void jZbar (HLV, HLV, HLV, HLV, bool, bool, current); */
/* COM PZ(double); */
/* double Zq (int, bool); */
/* double Gq (int); */
inline COM cdot(const current & j1, const current & j2)
{
return j1[0]*j2[0]-j1[1]*j2[1]-j1[2]*j2[2]-j1[3]*j2[3];
}
inline COM cdot(const HLV & p, const current & j1) {
return j1[0]*p.e()-j1[1]*p.x()-j1[2]*p.y()-j1[3]*p.z();
}
inline void cmult(const COM & factor, const current & j1, current &cur)
{
cur[0]=factor*j1[0];
cur[1]=factor*j1[1];
cur[2]=factor*j1[2];
cur[3]=factor*j1[3];
}
// WHY!?!
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, const current & j5, current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0]+j4[0]+j5[0];
sum[1]=j1[1]+j2[1]+j3[1]+j4[1]+j5[1];
sum[2]=j1[2]+j2[2]+j3[2]+j4[2]+j5[2];
sum[3]=j1[3]+j2[3]+j3[3]+j4[3]+j5[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
const current & j4, current &sum) {
sum[0] = j1[0] + j2[0] + j3[0] + j4[0];
sum[1] = j1[1] + j2[1] + j3[1] + j4[1];
sum[2] = j1[2] + j2[2] + j3[2] + j4[2];
sum[3] = j1[3] + j2[3] + j3[3] + j4[3];
}
inline void cadd(const current & j1, const current & j2, const current & j3,
current &sum)
{
sum[0]=j1[0]+j2[0]+j3[0];
sum[1]=j1[1]+j2[1]+j3[1];
sum[2]=j1[2]+j2[2]+j3[2];
sum[3]=j1[3]+j2[3]+j3[3];
}
inline void cadd(const current & j1, const current & j2, current &sum)
{
sum[0]=j1[0]+j2[0];
sum[1]=j1[1]+j2[1];
sum[2]=j1[2]+j2[2];
sum[3]=j1[3]+j2[3];
}
inline double abs2(const COM & a)
{
return (a*conj(a)).real();
}
inline double vabs2(const CCurrent & cur)
{
return abs2(cur.c0)-abs2(cur.c1)-abs2(cur.c2)-abs2(cur.c3);
}
inline double vre(const CCurrent & a, const CCurrent & b)
{
return real(a.c0*conj(b.c0)-a.c1*conj(b.c1)-a.c2*conj(b.c2)-a.c3*conj(b.c3));
}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 7e02575..8d75651 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1769 +1,2019 @@
#include "RHEJ/currents.hh"
#include "RHEJ/utility.hh"
#include "RHEJ/Tensor.hh"
#include "RHEJ/Constants.hh"
#include <array>
#include <iostream>
namespace { // Helper Functions
// FKL W Helper Functions
void jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe,
bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin,
bool helin, current cur)
{
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
current t65,vout,vin,temp2,temp3,temp5;
joo(pnu,helnu,pe,hele,t65);
vout[0]=pout.e();
vout[1]=pout.x();
vout[2]=pout.y();
vout[3]=pout.z();
vin[0]=pin.e();
vin[1]=pin.x();
vin[2]=pin.y();
vin[3]=pin.z();
COM brac615=cdot(t65,vout);
COM brac645=cdot(t65,vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
joo(pout,helout,pnu,helout,temp2);
COM prod1665=cdot(temp2,t65);
j(pe,helin,pin,helin,temp3);
COM prod5465=cdot(temp3,t65);
joo(pout,helout,pe,helout,temp2);
j(pnu,helnu,pin,helin,temp3);
j(pout,helout,pin,helin,temp5);
current term1,term2,term3,sum;
cmult(2.*brac615/ta+2.*brac645/tb,temp5,term1);
cmult(prod1665/ta,temp3,term2);
cmult(-prod5465/tb,temp2,term3);
cadd(term1,term2,term3,sum);
cur[0]=sum[0];
cur[1]=sum[1];
cur[2]=sum[2];
cur[3]=sum[3];
}
}
void jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
{
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
current t65,vout,vin,temp2,temp3,temp5;
joo(pnu,helnu,pe,hele,t65);
vout[0]=pout.e();
vout[1]=pout.x();
vout[2]=pout.y();
vout[3]=pout.z();
vin[0]=pin.e();
vin[1]=pin.x();
vin[2]=pin.y();
vin[3]=pin.z();
COM brac615=cdot(t65,vout);
COM brac645=cdot(t65,vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
joo(pe,helout,pout,helout,temp2); // temp2 is <5|alpha|1>
COM prod5165=cdot(temp2,t65);
jio(pin,helin,pnu,helin,temp3); // temp3 is <4|alpha|6>
COM prod4665=cdot(temp3,t65);
joo(pnu,helout,pout,helout,temp2); // temp2 is now <6|mu|1>
jio(pin,helin,pe,helin,temp3); // temp3 is now <4|mu|5>
jio(pin,helin,pout,helout,temp5); // temp5 is <4|mu|1>
current term1,term2,term3,sum;
cmult(-2.*brac615/ta-2.*brac645/tb,temp5,term1);
cmult(-prod5165/ta,temp3,term2);
cmult(prod4665/tb,temp2,term3);
cadd(term1,term2,term3,sum);
cur[0]=sum[0];
cur[1]=sum[1];
cur[2]=sum[2];
cur[3]=sum[3];
}
}
CCurrent jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin)
{
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pnu,helnu,pe,hele);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(pout,helout,pnu,helout);
COM prod1665=temp2.dot(t65);
temp3 = j(pe,helin,pin,helin);
COM prod5465=temp3.dot(t65);
temp2=joo(pout,helout,pe,helout);
temp3=j(pnu,helnu,pin,helin);
temp5=j(pout,helout,pin,helin);
CCurrent term1,term2,term3;
term1=(2.*brac615/ta+2.*brac645/tb)*temp5;
term2=(prod1665/ta)*temp3;
term3=(-prod5465/tb)*temp2;
sum=term1+term2+term3;
}
return sum;
}
CCurrent jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin)
{
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pnu,helnu,pe,hele);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(pe,helout,pout,helout); // temp2 is <5|alpha|1>
COM prod5165=temp2.dot(t65);
temp3 = jio(pin,helin,pnu,helin); // temp3 is <4|alpha|6>
COM prod4665=temp3.dot(t65);
temp2=joo(pnu,helout,pout,helout); // temp2 is now <6|mu|1>
temp3=jio(pin,helin,pe,helin); // temp3 is now <4|mu|5>
temp5=jio(pin,helin,pout,helout); // temp5 is <4|mu|1>
CCurrent term1,term2,term3;
term1 =(-2.*brac615/ta-2.*brac645/tb)*temp5;
term2 =(-prod5165/ta)*temp3;
term3 =(prod4665/tb)*temp2;
sum = term1 + term2 + term3;
}
return sum;
}
// Relevant W+Jets Unordered Contribution Helper Functions
// W+Jets Uno
double jM2Wuno(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pa, bool h1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector pb, bool h2, bool pol)
{
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
//std::cout<<"Setting sigma_index...." << std::endl;
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
CLHEP::HepLorentzVector pW = pl+plbar;
CLHEP::HepLorentzVector q1g=pa-pW-p1-pg;
CLHEP::HepLorentzVector q1 = pa-p1-pW;
CLHEP::HepLorentzVector q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double tb2 = (pb-p2).m2();
const double tb2g = (pb-p2-pg).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
const double ca = RHEJ::C_A;
const double cf = RHEJ::C_F; ///<TODO directly use RHEJ constants
//use p1 as ref vec in pol tensor
Tensor<1,4> epsg = eps(pg,p2,pol);
Tensor<1,4> epsW = TCurrent(pl,false,plbar,false);
Tensor<1,4> j2b = TCurrent(p2,h2,pb,h2);
Tensor<1,4> Tq1q2 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+ p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<1,4> Tq1g = Construct1Tensor((-pg-q1))/taW1;
Tensor<1,4> Tq2g = Construct1Tensor((pg-q2));
Tensor<1,4> TqaW = Construct1Tensor((pa-pW));//pa-pw
Tensor<1,4> Tqag = Construct1Tensor((pa-pg));
Tensor<1,4> TqaWg = Construct1Tensor((pa-pg-pW));
Tensor<1,4> Tp1g = Construct1Tensor((p1+pg));
Tensor<1,4> Tp1W = Construct1Tensor((p1+pW));//p1+pw
Tensor<1,4> Tp1gW = Construct1Tensor((p1+pg+pW));//p1+pw+pg
Tensor<2,4> g=Metric();
Tensor<3,4> J31a = T3Current(p1, h1, pa, h1);
Tensor<2,4> J2_qaW =J31a.contract(TqaW/taW, 2);
Tensor<2,4> J2_p1W =J31a.contract(Tp1W/s1W, 2);
Tensor<3,4> L1a =J2_qaW.leftprod(Tq1q2);
Tensor<3,4> L1b =J2_p1W.leftprod(Tq1q2);
Tensor<3,4> L2a = J2_qaW.leftprod(Tq1g);
Tensor<3,4> L2b = J2_p1W.leftprod(Tq1g);
Tensor<3,4> L3 = (g.rightprod(J2_qaW.contract(Tq2g,1)+J2_p1W.contract(Tq2g,2)))/taW1;
Tensor<3,4> L(0.);
Tensor<5,4> J51a = T5Current(p1, h1, pa, h1);
Tensor<4,4> J_qaW = J51a.contract(TqaW,4);
Tensor<4,4> J_qag = J51a.contract(Tqag,4);
Tensor<4,4> J_p1gW = J51a.contract(Tp1gW,4);
Tensor<3,4> U1a = J_qaW.contract(Tp1g,2);
Tensor<3,4> U1b = J_p1gW.contract(Tp1g,2);
Tensor<3,4> U1c = J_p1gW.contract(Tp1W,2);
Tensor<3,4> U1(0.);
Tensor<3,4> U2a = J_qaW.contract(TqaWg,2);
Tensor<3,4> U2b = J_qag.contract(TqaWg,2);
Tensor<3,4> U2c = J_qag.contract(Tp1W,2);
Tensor<3,4> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L.Set(nu, mu, rho, L1a.at(nu,mu,rho) + L1b.at(nu,rho,mu)
+ L2a.at(mu,nu,rho) + L2b.at(mu,rho,nu) + L3.at(mu,nu,rho));
U1.Set(nu, mu, rho, U1a.at(nu, mu, rho) / (s1g*taW)
+ U1b.at(nu,rho,mu) / (s1g*s1gW) + U1c.at(rho,nu,mu) / (s1W*s1gW));
U2.Set(nu,mu,rho,U2a.at(mu,nu,rho) / (taWg*taW)
+ U2b.at(mu,rho,nu) / (taWg*tag) + U2c.at(rho,mu,nu) / (s1W*tag));
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
double amp = ca*cf*cf/2.*(norm(X)+norm(Y)) - cf/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
//Divide by t-channels
amp/=(t1*t2);
//Average over initial states
amp/=(4.*ca*ca);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1,4> gtqqxW(CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar){
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
//blank 3 Gamma Current
Tensor<3,4> JV23 = T3Current(pq,false,pqbar,false);
// Components of g->qqW before W Contraction
Tensor<2,4> JV1 = JV23.contract((Tpq + TAB),2)/(s2AB);
Tensor<2,4> JV2 = JV23.contract((Tpqbar + TAB),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1,4> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2,4> MCrossW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar,false);
//Blank 5 gamma Current
Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4,4> J_q3q = J523.contract((Tq3+Tpq),2);
Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
// Components of Crossed Vertex Contribution
Tensor<3,4> Xcro1 = J_q3q.contract((Tpqbar + TAB),3);
Tensor<3,4> Xcro2 = J_q3q.contract((Tq1-Tpqbar),3);
Tensor<3,4> Xcro3 = J_2AB.contract((Tq1-Tpqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2,4> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2,4> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2,4> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2,4> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro.Set(mu,nu, -(-Xcro1Cont.at(nu,mu)+Xcro2Cont.at(nu,mu)+Xcro3Cont.at(nu,mu)));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2,4> MUncrossW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
//Blank 5 gamma Current
Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
Tensor<4,4> J_q1q = J523.contract((Tq1-Tpq),2);
// 2 Contractions taken care of.
Tensor<3,4> Xunc1 = J_2AB.contract((Tq3+Tpqbar),3);
Tensor<3,4> Xunc2 = J_q1q.contract((Tq3+Tpqbar),3);
Tensor<3,4> Xunc3 = J_q1q.contract((Tpqbar+TAB),3);
// Term Denominators Taken Care of at this stage
Tensor<2,4> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2,4> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2,4> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2,4> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc.Set(mu,nu,-(- Xunc1Cont.at(mu,nu)+Xunc2Cont.at(mu,nu) +Xunc3Cont.at(mu,nu)));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2,4> MSymW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Define Tensors to be used
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1,4> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13+s1A+s1B)) + Tpa*(t1/(sa2+sa3+saA+saB)));
Tensor<2,4> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43+s4A+s4B)) + Tpb*(t3/(sb2+sb3+sbA+sbB)));
Tensor<2,4> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar+TAB);
Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar-TAB);
Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
// Note the contraction of indices changes term by term
Tensor<2,4> X3g1Cont = X3g1.contract(JV,3);
Tensor<2,4> X3g2Cont = X3g2.contract(JV,2);
Tensor<2,4> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g. Makes sense from a colour factor point of view.
Tensor<2,4>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym.Set(mu,nu, (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu) - X3g3Cont.at(nu,mu))
+ (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) );
}
}
return Xsym/s23AB;
}
Tensor <2,4> MCross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
Tensor<1,4> Tq1 = Construct1Tensor(q1-pqbar);
//Blank 3 gamma Current
Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2,4> XCroCont = J323.contract((Tq1),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2,4> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro.Set(mu,nu, (XCroCont.at(nu,mu)));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2,4> MUncross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
Tensor<1,4> Tq1 = Construct1Tensor(q1-pq);
//Blank 3 gamma Current
Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2,4> XUncCont = J323.contract((Tq1),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2,4> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc.Set(mu,nu,-(XUncCont.at(mu,nu)));
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2,4> MSym(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
//Define Tensors to be used
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
Tensor<1,4> qqxCur = TCurrent(pq, hq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13))+Tpa*(t1/(sa2+sa3)));
Tensor<2,4> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43)) + Tpb*(t3/(sb2+sb3)));
Tensor<2,4> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar);
Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar);
Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
// Note the contraction of indices changes term by term
Tensor<2,4> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2,4> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2,4> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2,4>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym.Set(mu, nu, COM(0,1) * ( (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu)
- X3g3Cont.at(nu,mu)) + (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) ) );
}
}
return Xsym/s23;
}
Tensor <1,4> jW4bEmit(HLV pb, HLV p4, HLV pl, HLV plbar, bool aqlinepb){
// Build the external quark line W Emmision
Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false)/2;
Tensor<1,4> Tp4W = Construct1Tensor((p4+pl+plbar));//p4+pw
Tensor<1,4> TpbW = Construct1Tensor((pb-pl-plbar));//pb-pw
Tensor<3,4> J4bBlank;
if (aqlinepb){
J4bBlank = T3Current(pb,false,p4,false);
}
else{
J4bBlank = T3Current(p4,false,pb,false);
}
double t4AB = (p4+pl+plbar).m2();
double tbAB = (pb-pl-plbar).m2();
Tensor<2,4> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB;
Tensor<2,4> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB;
Tensor<2,4> T4bmMom(0.);
if (aqlinepb){
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom.Set(mu,nu, (J4b1.at(nu,mu) + J4b2.at(mu,nu))*(COM(-1,0)));
}
}
}
else{
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom.Set(nu,mu, (J4b1.at(nu,mu) + J4b2.at(mu,nu)));
}
}
}
Tensor<1,4> T4bm = T4bmMom.contract(ABCurr,1);
return T4bm;
}
} // Anonymous Namespace helper functions
//Functions which can be called elsewhere (declarations in currents.hh).
// W+Jets Unordered Contributions
//qQ->qQWg_unob
double junobMWqQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj1m,mj2p,mj2m;
CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
CLHEP::HepLorentzVector q3=-(p2in-p2out);
mj1m=jW(p1out,false,pe,false,pnu,false,p1in,false);
mj2p=j(p2out,true,p2in,true);
mj2m=j(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgbm,jgbp,j2gm,j2gp;
j2gp=joo(p2out,true,pg,true);
j2gm=joo(p2out,false,pg,false);
jgbp=j(pg,true,p2in,true);
jgbm=j(pg,false,p2in,false);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
const double cf=RHEJ::C_F;
double amm,amp;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
// Now add the t-channels
double th=q2.m2()*q1.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qQbar->qQbarWg_unob
double junobMWqQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj1m,mj2p,mj2m;
CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
CLHEP::HepLorentzVector q3=-(p2in-p2out);
mj1m=jW(p1out,false,pe,false,pnu,false,p1in,false);
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgbm,jgbp,j2gm,j2gp;
j2gp=joo(pg,true,p2out,true);
j2gm=joo(pg,false,p2out,false);
jgbp=jio(p2in,true,pg,true);
jgbm=jio(p2in,false,pg,false);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
const double cf=RHEJ::C_F;
double amm,amp;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
// Now add the t-channels
double th=q2.m2()*q1.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qbarQ->qbarQWg_unob
double junobMWqbarQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj1m,mj2p,mj2m;
CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
CLHEP::HepLorentzVector q3=-(p2in-p2out);
mj1m=jWbar(p1out,false,pe,false,pnu,false,p1in,false);
mj2p=j(p2out,true,p2in,true);
mj2m=j(p2out,false,p2in,false);
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgbm,jgbp,j2gm,j2gp;
j2gp=joo(p2out,true,pg,true);
j2gm=joo(p2out,false,pg,false);
jgbp=j(pg,true,p2in,true);
jgbm=j(pg,false,p2in,false);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
const double cf=RHEJ::C_F;
double amm,amp;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
// Now add the t-channels
double th=q2.m2()*q1.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qbarQbar->qbarQbarWg_unob
double junobMWqbarQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pg)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj1m,mj2p,mj2m;
CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pg);
CLHEP::HepLorentzVector q3=-(p2in-p2out);
mj1m=jWbar(p1out,false,pe,false,pnu,false,p1in,false);
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
// Dot products of these which occur again and again
COM MWmp=mj1m.dot(mj2p); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgbm,jgbp,j2gm,j2gp;
j2gp=joo(pg,true,p2out,true);
j2gm=joo(pg,false,p2out,false);
jgbp=jio(p2in,true,pg,true);
jgbm=jio(p2in,false,pg,false);
CCurrent qsum(q2+q3);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in);
CCurrent p2o(p2out);
CCurrent p2i(p2in);
Lmm=((-1.)*qsum*(MWmm) + (-2.*mj1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmm/2.))/q3.m2();
Lmp=((-1.)*qsum*(MWmp) + (-2.*mj1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mj1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MWmp/2.))/q3.m2();
U1mm=(jgbm.dot(mj1m)*j2gm+2.*p2o*MWmm)/(p2out+pg).m2();
U1mp=(jgbp.dot(mj1m)*j2gp+2.*p2o*MWmp)/(p2out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj1m)*jgbm+2.*p2i*MWmm)/(p2in-pg).m2();
U2mp=((-1.)*j2gp.dot(mj1m)*jgbp+2.*p2i*MWmp)/(p2in-pg).m2();
const double cf=RHEJ::C_F;
double amm,amp;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
// Now add the t-channels
double th=q2.m2()*q1.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
////////////////////////////////////////////////////////////////////
//qQ->qQWg_unof
double junofMWgqQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj2m,mj1p,mj1m;
CLHEP::HepLorentzVector q1=p1in-p1out;
CLHEP::HepLorentzVector qg=p1in-p1out-pg;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
mj2m=jW(p2out,false,pe,false,pnu,false,p2in,false);
mj1p=j(p1out,true,p1in,true);
mj1m=j(p1out,false,p1in,false);
// Dot products of these which occur again and again
COM MWpm=mj1p.dot(mj2m); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgam,jgap,j2gm,j2gp;
j2gp=joo(p1out,true,pg,true);
j2gm=joo(p1out,false,pg,false);
jgap=j(pg,true,p1in,true);
jgam=j(pg,false,p1in,false);
CCurrent qsum(q1+qg);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
CCurrent p1o(p1out);
CCurrent p1i(p1in);
Lmm=(qsum*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
const double cf=RHEJ::C_F;
double amm,apm;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
// Now add the t-channels
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qQbar->qQbarWg_unof
double junofMWgqQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj2m,mj1p,mj1m;
CLHEP::HepLorentzVector q1=p1in-p1out;
CLHEP::HepLorentzVector qg=p1in-p1out-pg;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
mj2m=jWbar(p2out,false,pe,false,pnu,false,p2in,false);
mj1p=j(p1out,true,p1in,true);
mj1m=j(p1out,false,p1in,false);
// Dot products of these which occur again and again
COM MWpm=mj1p.dot(mj2m); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgam,jgap,j2gm,j2gp;
j2gp=joo(p1out,true,pg,true);
j2gm=joo(p1out,false,pg,false);
jgap=j(pg,true,p1in,true);
jgam=j(pg,false,p1in,false);
CCurrent qsum(q1+qg);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
CCurrent p1o(p1out);
CCurrent p1i(p1in);
Lmm=(qsum*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
const double cf=RHEJ::C_F;
double amm,apm;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
// Now add the t-channels
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qbarQ->qbarQWg_unof
double junofMWgqbarQ (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj2m,mj1p,mj1m;
CLHEP::HepLorentzVector q1=p1in-p1out;
CLHEP::HepLorentzVector qg=p1in-p1out-pg;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
mj2m=jW(p2out,false,pe,false,pnu,false,p2in,false);
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
// Dot products of these which occur again and again
COM MWpm=mj1p.dot(mj2m); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgam,jgap,j2gm,j2gp;
j2gp=joo(pg,true,p1out,true);
j2gm=joo(pg,false,p1out,false);
jgap=jio(p1in,true,pg,true);
jgam=jio(p1in,false,pg,false);
CCurrent qsum(q1+qg);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
CCurrent p1o(p1out);
CCurrent p1i(p1in);
Lmm=(qsum*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
const double cf=RHEJ::C_F;
double amm,apm;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
// Now add the t-channels
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
//qbarQbar->qbarQbarWg_unof
double junofMWgqbarQbar (CLHEP::HepLorentzVector pg,CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu, CLHEP::HepLorentzVector p2in)
// Calculates the square of the current contractions for qQ->qenuQ scattering
// p1: quark (with W emittance)
// p2: Quark
{
CCurrent mj2m,mj1p,mj1m;
CLHEP::HepLorentzVector q1=p1in-p1out;
CLHEP::HepLorentzVector qg=p1in-p1out-pg;
CLHEP::HepLorentzVector q2=-(p2in-p2out-pe-pnu);
mj2m=jWbar(p2out,false,pe,false,pnu,false,p2in,false);
mj1p=jio(p1in,true,p1out,true);
mj1m=jio(p1in,false,p1out,false);
// Dot products of these which occur again and again
COM MWpm=mj1p.dot(mj2m); // And now for the Higgs ones
COM MWmm=mj1m.dot(mj2m);
CCurrent jgam,jgap,j2gm,j2gp;
j2gp=joo(pg,true,p1out,true);
j2gm=joo(pg,false,p1out,false);
jgap=jio(p1in,true,pg,true);
jgam=jio(p1in,false,pg,false);
CCurrent qsum(q1+qg);
CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in);
CCurrent p1o(p1out);
CCurrent p1i(p1in);
Lmm=(qsum*(MWmm) + (-2.*mj2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWmm/2.))/q1.m2();
Lpm=(qsum*(MWpm) + (-2.*mj2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mj2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MWpm/2.))/q1.m2();
U1mm=(jgam.dot(mj2m)*j2gm+2.*p1o*MWmm)/(p1out+pg).m2();
U1pm=(jgap.dot(mj2m)*j2gp+2.*p1o*MWpm)/(p1out+pg).m2();
U2mm=((-1.)*j2gm.dot(mj2m)*jgam+2.*p1i*MWmm)/(p1in-pg).m2();
U2pm=((-1.)*j2gp.dot(mj2m)*jgap+2.*p1i*MWpm)/(p1in-pg).m2();
const double cf=RHEJ::C_F;
double amm,apm;
amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm);
apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
// Now add the t-channels
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
///TODO make this comment more visible
/// Naming scheme jM2-Wuno-g-({q/qbar}{Q/Qbar/g})
///TODO Spit naming for more complicated functions?
/// e.g. jM2WqqtoqQQq -> jM2_Wqq_to_qQQq
double jM2WunogqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
//same as function above but actually obtaining the antiquark line by crossing symmetry, where p1out and p1in are expected to be negative.
//should give same result as jM2WunogqbarQ below (verified)
double jM2WunogqQ_crossqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2WunogqQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2Wunogqg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
const double ca = RHEJ::C_A;
const double cf = RHEJ::C_F; ///<TODO directly use RHEJ constants
double ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (ca - 1/ca)*(ratio + 1./ratio)/2. + 1/ca)/cf;
ME2*=cam;
return ME2;
}
double jM2WunogqbarQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2WunogqbarQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2Wunogqbarg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,true,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
const double ca = RHEJ::C_A;
const double cf = RHEJ::C_F; ///<TODO directly use RHEJ constants
double ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (ca - 1/ca)*(ratio + 1./ratio)/2. + 1/ca)/cf;
ME2*=cam;
return ME2;
}
// W+Jets qqxExtremal
// W+Jets qqxExtremal Currents - wqq emission
double jM2WgQtoqbarqQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
//Correct colour averaging
ME2*=(3.0/8.0);
return ME2;
}
double jM2WgQtoqqbarQ(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqbarout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in){
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,true);
ME2mpm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,false);
ME2mmp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,true);
ME2mmm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
//Correct colour averaging
ME2*=(3.0/8.0);
return ME2;
}
double jM2Wggtoqbarqg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqbarout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
const double ca = RHEJ::C_A;
const double cf = RHEJ::C_F; ///<TODO directly use RHEJ constants
double ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (ca - 1/ca)*(ratio + 1./ratio)/2. + 1/ca)/cf;
ME2*=cam;
//Correct colour averaging
ME2*=(3.0/8.0);
return ME2;
}
double jM2Wggtoqqbarg(CLHEP::HepLorentzVector pgin, CLHEP::HepLorentzVector pqbarout,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pqout, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in){
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,true);
ME2mpm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,true,false);
ME2mmp = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,true);
ME2mmm = jM2Wuno(-pgin, pqbarout,plbar,pl,-pqout,true,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
const double ca = RHEJ::C_A;
const double cf = RHEJ::C_F; ///<TODO directly use RHEJ constants
double ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (ca - 1/ca)*(ratio + 1./ratio)/2. + 1/ca)/cf;
ME2*=cam;
//Correct colour averaging
ME2*=(3.0/8.0);
return ME2;
}
+namespace {
+//First, a function for generating polarisation tensors. Output as 'current'.
+void eps(CLHEP::HepLorentzVector refmom, CLHEP::HepLorentzVector kb, bool hel, current &ep){
+ current curm,curp;
+ //Recall - positive helicity eps has negative helicity choices for spinors and vice versa
+ j(refmom,true,kb,true,curm);
+ j(refmom,false,kb,false,curp);
+ double norm=1.;
+ if(kb.z()<0.)
+ norm *= sqrt(2.*refmom.plus()*kb.minus());
+ if(kb.z()>0.)
+ norm = sqrt(2.*refmom.minus()*kb.plus());
+ if(hel==false){
+ ep[0] = curm[0]/norm;
+ ep[1] = curm[1]/norm;
+ ep[2] = curm[2]/norm;
+ ep[3] = curm[3]/norm;
+ }
+ if(hel==true){
+ ep[0] = curp[0]/norm;
+ ep[1] = curp[1]/norm;
+ ep[2] = curp[2]/norm;
+ ep[3] = curp[3]/norm;
+ }
+}
+
+
+//Now build up each part of the squared amplitude
+
+COM qWggm1(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom){
+ current cur33, cur23, curb3, cur2b, ep;
+ joo(p3, helchain, p3, helchain,cur33);
+ joo(p2,helchain,p3,helchain,cur23);
+ jio(pb,helchain,p3,helchain,curb3);
+ joi(p2,helchain,pb,helchain,cur2b);
+
+ // Build the external quark line W Emmision
+ Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false);
+ Tensor<1,4> Tp1W = Construct1Tensor((p1+pl+plbar));//p1+pw
+ Tensor<1,4> TqaW = Construct1Tensor((pa-pl-plbar));//pa-pw
+ Tensor<3,4> J1aBlank = T3Current(p1,false,pa,false);
+ double t1AB = (p1+pl+plbar).m2();
+ double taAB = (pa-pl-plbar).m2();
+ Tensor<2,4> J1a1 = (J1aBlank.contract(Tp1W,2))/t1AB;
+ Tensor<2,4> J1a2 = (J1aBlank.contract(TqaW,2))/taAB;
+ Tensor<1,4> cur1a = J1a1.contract(ABCurr,1) + J1a2.contract(ABCurr,2);
+
+ double t2 = (p3-pb)*(p3-pb);
+
+ //Create vertex
+ COM v1[4][4];
+ for(int u=0; u<4;u++)
+ {
+ for(int v=0; v<4; v++)
+ {
+ v1[u][v]=(cur23[u]*cur33[v]-cur2b[u]*curb3[v])/t2*(-1.);
+ }
+ }
+ //Dot in current and eps
+ //Metric tensor
+ double eta[4][4]={};
+ eta[0][0]=1.;
+ eta[1][1]=-1.;
+ eta[2][2]=-1.;
+ eta[3][3]=-1.;
+ //eps
+ eps(refmom,pb,helb, ep);
+ COM M1=0.;
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ for(int k=0; k<4; k++){
+ for(int l=0; l<4;l++){
+ M1+= eta[i][k]*cur1a.at(k)*(v1[i][j])*ep[l]*eta[l][j];
+ }
+ }
+ }
+ }
+ return M1;
+}
+
+COM qWggm2(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom){
+ current cur22, cur23, curb3, cur2b, ep;
+ joo(p2, helchain, p2, helchain,cur22);
+ joo(p2,helchain,p3,helchain,cur23);
+ jio(pb,helchain,p3,helchain,curb3);
+ joi(p2,helchain,pb,helchain,cur2b);
+
+ // Build the external quark line W Emmision
+ Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false);
+ Tensor<1,4> Tp1W = Construct1Tensor((p1+pl+plbar));//p1+pw
+ Tensor<1,4> TqaW = Construct1Tensor((pa-pl-plbar));//pa-pw
+
+ Tensor<3,4> J1aBlank = T3Current(p1,false,pa,false);
+
+ double t1AB = (p1+pl+plbar).m2();
+ double taAB = (pa-pl-plbar).m2();
+
+ Tensor<2,4> J1a1 = (J1aBlank.contract(Tp1W,2))/t1AB;
+ Tensor<2,4> J1a2 = (J1aBlank.contract(TqaW,2))/taAB;
+
+ Tensor<1,4> cur1a = J1a1.contract(ABCurr,1) + J1a2.contract(ABCurr,2);
+
+ double t2t = (p2-pb)*(p2-pb);
+ //Create vertex
+ COM v2[4][4]={};
+ for(int u=0; u<4;u++)
+ {
+ for(int v=0; v<4; v++)
+ {
+ v2[u][v]=(cur22[v]*cur23[u]-cur2b[v]*curb3[u])/t2t;
+ }
+ }
+ //Dot in current and eps
+ //Metric tensor
+ double eta[4][4]={};
+ eta[0][0]=1.;
+ eta[1][1]=-1.;
+ eta[2][2]=-1.;
+ eta[3][3]=-1.;
+ //eps
+ eps(refmom,pb,helb, ep);
+ COM M2=0.;
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ for(int k=0; k<4; k++){
+ for(int l=0; l<4;l++){
+ M2+= eta[i][k]*cur1a.at(k)*(v2[i][j])*ep[l]*eta[l][j];
+ }
+ }
+ }
+ }
+ return M2;
+}
+
+COM qWggm3(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom){
+ //3 gluon vertex bit
+ double eta[4][4]={};
+ eta[0][0]=1.;
+ eta[1][1]=-1.;
+ eta[2][2]=-1.;
+ eta[3][3]=-1.;
+ current spincur,ep;
+ double s23 = (p2+p3)*(p2+p3);
+ joo(p2,helchain,p3,helchain,spincur);
+
+ // Build the external quark line W Emmision
+ Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false);
+ Tensor<1,4> Tp1W = Construct1Tensor((p1+pl+plbar));//p1+pw
+ Tensor<1,4> TqaW = Construct1Tensor((pa-pl-plbar));//pa-pw
+
+ Tensor<3,4> J1aBlank = T3Current(p1,false,pa,false);
+
+ double t1AB = (p1+pl+plbar).m2();
+ double taAB = (pa-pl-plbar).m2();
+
+ Tensor<2,4> J1a1 = (J1aBlank.contract(Tp1W,2))/t1AB;
+ Tensor<2,4> J1a2 = (J1aBlank.contract(TqaW,2))/taAB;
+
+ Tensor<1,4> cur1a = J1a1.contract(ABCurr,1) + J1a2.contract(ABCurr,2);
+
+ //Redefine relevant momenta as currents - for ease of calling correct part of vector
+ current ka,k2,k3,kb;
+ kb[0]=pb.e();
+ kb[1]=pb.x();
+ kb[2]=pb.y();
+ kb[3]=pb.z();
+ k2[0]=p2.e();
+ k2[1]=p2.x();
+ k2[2]=p2.y();
+ k2[3]=p2.z();
+ k3[0]=p3.e();
+ k3[1]=p3.x();
+ k3[2]=p3.y();
+ k3[3]=p3.z();
+ ka[0]=pa.e();
+ ka[1]=pa.x();
+ ka[2]=pa.y();
+ ka[3]=pa.z();
+ COM V3g[4][4]={};
+ for(int u=0;u<4;u++){
+ for(int v=0;v<4;v++){
+ for(int p=0;p<4;p++){
+ for(int r=0; r<4;r++){
+ V3g[u][v] += COM(0.,1.)*(((2.*k2[v]+2.*k3[v])*eta[u][p] - (2.*kb[u])*eta[p][v]+2.*kb[p]*eta[u][v])*spincur[r]*eta[r][p])/s23;
+ }
+ }
+ }
+ }
+ COM diffextrabit[4][4]={};
+ for(int u=0;u<4;u++) {
+ for(int v=0;v<4;v++) {
+ diffextrabit[u][v] = 0.;
+ }
+ }
+
+ //Dot in current and eps
+ //eps
+ eps(refmom,pb,helb, ep);
+ COM M3=0.;
+ for(int i=0;i<4;i++){
+ for(int j=0;j<4;j++){
+ for(int k=0; k<4; k++){
+ for(int l=0; l<4;l++){
+ M3+= eta[i][k]*cur1a.at(k)*(V3g[i][j]+diffextrabit[i][j])*ep[l]*eta[l][j];
+ }
+ }
+ }
+ }
+ return M3;
+}
+}
+
+// no wqq emission
+double jM2WgqtoQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl){
+
+ // 4 indepedent helicity choices (complex conjugation related).
+
+ //Need to evalute each independent hel configuration and store that result somewhere
+ COM Mmmm1 = qWggm1(pa,pb,p1,p2,p3,pl,plbar,false,false,false, pa);
+ COM Mmmm2 = qWggm2(pa,pb,p1,p2,p3,pl,plbar,false,false,false, pa);
+ COM Mmmm3 = qWggm3(pa,pb,p1,p2,p3,pl,plbar,false,false,false, pa);
+ COM Mmmp1 = qWggm1(pa,pb,p1,p2,p3,pl,plbar,false,true,false, pa);
+ COM Mmmp2 = qWggm2(pa,pb,p1,p2,p3,pl,plbar,false,true,false, pa);
+ COM Mmmp3 = qWggm3(pa,pb,p1,p2,p3,pl,plbar,false,true,false, pa);
+ COM Mpmm1 = qWggm1(pa,pb,p1,p2,p3,pl,plbar,false,false,true, pa);
+ COM Mpmm2 = qWggm2(pa,pb,p1,p2,p3,pl,plbar,false,false,true, pa);
+ COM Mpmm3 = qWggm3(pa,pb,p1,p2,p3,pl,plbar,false,false,true, pa);
+ COM Mpmp1 = qWggm1(pa,pb,p1,p2,p3,pl,plbar,false,true,true, pa);
+ COM Mpmp2 = qWggm2(pa,pb,p1,p2,p3,pl,plbar,false,true,true, pa);
+ COM Mpmp3 = qWggm3(pa,pb,p1,p2,p3,pl,plbar,false,true,true, pa);
+
+ //Colour factors:
+ COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
+ cm1m1=8./3.;
+ cm2m2=8./3.;
+ cm3m3=6.;
+ cm1m2 =-1./3.;
+ cm1m3 = -3.*COM(0.,1.);
+ cm2m3 = 3.*COM(0.,1.);
+
+ //Sqaure and sum for each helicity config:
+ double Mmmm,Mmmp,Mpmm,Mpmp;
+ Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
+ Mmmp = real(cm1m1*pow(abs(Mmmp1),2)+cm2m2*pow(abs(Mmmp2),2)+cm3m3*pow(abs(Mmmp3),2)+2.*real(cm1m2*Mmmp1*conj(Mmmp2))+2.*real(cm1m3*Mmmp1*conj(Mmmp3))+2.*real(cm2m3*Mmmp2*conj(Mmmp3)));
+ Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
+ Mpmp = real(cm1m1*pow(abs(Mpmp1),2)+cm2m2*pow(abs(Mpmp2),2)+cm3m3*pow(abs(Mpmp3),2)+2.*real(cm1m2*Mpmp1*conj(Mpmp2))+2.*real(cm1m3*Mpmp1*conj(Mpmp3))+2.*real(cm2m3*Mpmp2*conj(Mpmp3)));
+
+ return ((Mmmm+Mmmp+Mpmm+Mpmp)/24./4.)/(pa-p1).m2()/(p2+p3-pb).m2();
+}
+
// W+Jets qqxCentral
double jM2WqqtoqQQq(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove)
{
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;}
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1,4> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
if(!(aqlinepb)){
T4bp = TCurrent(p4, true, pb, true);
T4bm = TCurrent(p4, false, pb, false);}
else if(aqlinepb){
T4bp = TCurrent(pb, true, p4, true);
T4bm = TCurrent(pb, false, p4, false);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2,4> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2,4> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2,4> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1)).at(0);
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1)).at(0);
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1)).at(0);
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1)).at(0);
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1)).at(0);
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1)).at(0);
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1)).at(0);
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1)).at(0);
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
return amp;
}
// no wqq emission
double jM2WqqtoqQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<CLHEP::HepLorentzVector> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove, int nbelow, bool forwards){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
if (!forwards){ //If Emission from Leg a instead, flip process.
HLV dummymom = pa;
bool dummybool= aqlinepa;
int dummyint = nabove;
pa = pb;
pb = dummymom;
std::reverse(partons.begin(),partons.end());
qqxmarker = !(qqxmarker);
aqlinepa = aqlinepb;
aqlinepb = dummybool;
nabove = nbelow;
nbelow = dummyint;
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1,4> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
Tensor <1,4> T4bm = jW4bEmit(pb, p4, pl, plbar, aqlinepb);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2,4> Xunc_m = MUncross(pa, pq, pqbar,partons, false, nabove);
Tensor<2,4> Xcro_m = MCross( pa, pq, pqbar,partons, false, nabove);
Tensor<2,4> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, false, nabove);
Tensor<2,4> Xunc_p = MUncross(pa, pq, pqbar,partons, true, nabove);
Tensor<2,4> Xcro_p = MCross( pa, pq, pqbar,partons, true, nabove);
Tensor<2,4> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, true, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1)).at(0);
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1)).at(0);
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
return amp;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:55 AM (6 h, 56 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983197
Default Alt Text
(135 KB)

Event Timeline