diff --git a/include/RHEJ/Tensor.hh b/include/RHEJ/Tensor.hh index 4cba58b..e3867d2 100644 --- a/include/RHEJ/Tensor.hh +++ b/include/RHEJ/Tensor.hh @@ -1,199 +1,201 @@ /** \file * \brief Tensor Template Class declaration. * * This file contains the declaration of the Tensor Template class. This * is used to calculate some of the more complex currents within the * W+Jets implementation particularly. */ +#pragma once #include +///TODO remove function implementation from header template class Tensor { public: //Constructor Tensor(void); Tensor(COM x); //Destructor virtual ~Tensor(void); int rank(void){ return N; } int dim(void){ return D; } int len(){ return size; } COM at(int i){ return components[i]; } COM at(int i, int j) { return components[D*i +j]; } COM at(int i, int j, int k) { return components[D*(D*i + j)+ k]; } COM at(int i,int j, int k,int l) { return components[D*(D*(D*i +j) + k) + l]; } COM at(int i,int j, int k,int l, int m){ return components[D*(D*(D*(D*i + j) + k) + l) + m]; } bool isSet(void){ if(components.size()==0) return false; else return true; } void Fill(COM x){ components=x; } //Set component indexed by i,j,k,l,m void Set(int i,COM x){ components[i] = x; } void Set(int i, int j, COM x) { components[D*i +j] = x; } void Set(int i, int j, int k, COM x) { components[D*(D*i + j)+ k] = x; } void Set(int i,int j, int k,int l,COM x) { components[D*(D*(D*i +j) + k) + l] = x; } void Set(int i,int j, int k,int l, int m, COM x){ components[D*(D*(D*(D*i + j) + k) + l) + m] = x; } Tensor operator*(const double x){ Tensor newT; newT.components=components*COM(x,0); return newT; } Tensor operator*(const COM x){ Tensor newT; newT.components=components*x; return newT; } Tensor operator/(const double x){ Tensor newT; newT.components=components/COM(x,0); return newT; } Tensor operator/(const COM x){ Tensor newT; newT.components=components/x; return newT; } Tensor operator+(const Tensor T2){ Tensor newT; newT.components=components+T2.components; return newT; } Tensor operator-(const Tensor T2){ Tensor newT; newT.components=components-T2.components; return newT; } void operator+=(const Tensor T2){ components = components+T2.components; } void operator-=(const Tensor T2){ components=components-T2.components; } Tensor rightprod(const Tensor<1,D> T2){ Tensor newT; for(int i=0; i leftprod(const Tensor<1,D> T2){ Tensor newT; for(unsigned int j=0;j contract(const Tensor<1,D> T2, int k){ Tensor newT; for(int j=0; j components; private: int size; COM sign(unsigned int i){ if(i==0) return 1.; else return -1.; } }; template Tensor::Tensor(void) { size = pow(D,N); components.resize(size); } template Tensor::Tensor(COM x) { size = pow(D,N); components.resize(size, x); } template Tensor::~Tensor(void) {} // Tensor Functions: // Tensor<1,4> Sigma(int i, int j, bool hel); // Tensor<2,4> Metric(void); // int tensor2listindex(std::array indexlist); // int tensor2listindex(std::array indexlist); // void perms41(int same4, int diff, std::vector> * perms); // void perms32(int same3, int diff, std::vector> * perms); // void perms311(int same3, int diff1, int diff2, std::vector> * perms); // void perms221(int same2a, int same2b, int diff, std::vector> * perms); // void perms2111(int same2, int diff1,int diff2,int diff3, std::vector> * perms); // void perms21(int same, int diff, std::vector> * perms); // void perms111(int diff1, int diff2, int diff3, std::vector> * perms); Tensor<2,4> Metric(void); Tensor<1,4> TCurrent(CLHEP::HepLorentzVector p1, bool h1,CLHEP::HepLorentzVector p2, bool h2); Tensor<3,4> T3Current(CLHEP::HepLorentzVector p1, bool h1,CLHEP::HepLorentzVector p2, bool h2); Tensor<5,4> T5Current(CLHEP::HepLorentzVector p1, bool h1,CLHEP::HepLorentzVector p2, bool h2); Tensor<1,4> Construct1Tensor(CCurrent j); Tensor<1,4> Construct1Tensor(CLHEP::HepLorentzVector p); Tensor<1,4> eps(CLHEP::HepLorentzVector k, CLHEP::HepLorentzVector ref, bool pol); bool init_sigma_index(void); diff --git a/include/RHEJ/currents.hh b/include/RHEJ/currents.hh index ca8af5e..435650e 100644 --- a/include/RHEJ/currents.hh +++ b/include/RHEJ/currents.hh @@ -1,1322 +1,1324 @@ ////////////////////////////////////////////////// ////////////////////////////////////////////////// // 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 #include #include #include #include +///TODO add a namespace + typedef std::complex 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::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 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 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 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. */ CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin); //! 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 jio (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin); //! 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 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 jMZqQ (HLV, HLV, HLV, HLV, HLV, HLV, std::vector , std::vector < std::vector >, int, int, bool, bool); */ /* std::vector jMZqg (HLV, HLV, HLV, HLV, HLV, HLV, std::vector , std::vector < std::vector >, 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 b1b3e50..b136ebf 100644 --- a/src/Wjets.cc +++ b/src/Wjets.cc @@ -1,1968 +1,1973 @@ #include "RHEJ/currents.hh" #include "RHEJ/utility.hh" #include "RHEJ/Tensor.hh" +#include "RHEJ/Constants.hh" #include #include 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) + 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]; + COM cur[4]; - cur[0]=0.; - cur[1]=0.; - cur[2]=0.; - cur[3]=0.; - CCurrent sum(0.,0.,0.,0.); + 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; + is_sigma_index_set = true; else - return 0.; + 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; - double taW=(pa-pW).m2(); - double taW1=(pa-pW-p1).m2(); - double tb2=(pb-p2).m2(); - double tb2g=(pb-p2-pg).m2(); - double s1W=(p1+pW).m2(); - double s1gW=(p1+pW+pg).m2(); - double s1g=(p1+pg).m2(); - double tag=(pa-pg).m2(); - double taWg=(pa-pW-pg).m2(); - double ca =3.; - double cf =4./3.; + 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; /// 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> 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)); - } + 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. Note factor of 2 to account for sqrt(2) on each W vertex Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false)/2; //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 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 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. Note factor of 2 to account for sqrt(2) on each W vertex Tensor<1,4> ABCur = TCurrent(pl, false, plbar,false)/2; //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))); + 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 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 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. Note factor of 2 to account for sqrt(2) on each W vertex Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false)/2; //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))); + 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 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 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)) ); + 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 partons, bool hq, int nabove){ + Tensor <2,4> MCross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector partons, bool hq, int nabove){ - CLHEP::HepLorentzVector q1; - q1=pa; - for(int i=0;i Tq1 = Construct1Tensor(q1-pqbar); + Tensor<1,4> Tq1 = Construct1Tensor(q1-pqbar); - //Blank 3 gamma Current - Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq); + //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); + // 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.); + //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))); + for(int mu=0; mu<4;mu++){ + for(int nu=0;nu<4;nu++){ + Xcro.Set(mu,nu, (XCroCont.at(nu,mu))); + } } - } - return Xcro; -} + return Xcro; + } // Helper Functions Calculate the Uncrossed Contribution Tensor <2,4> MUncross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector partons, bool hq, int nabove){ CLHEP::HepLorentzVector q1; q1=pa; for(int i=0;i 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))); + 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 partons, bool hq, int nabove){ CLHEP::HepLorentzVector q1, q3; q1=pa; for(int i=0;i 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)) )); + 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 <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<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> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB; + Tensor<2,4> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB; - Tensor<2,4> T4bmMom(0.); + 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))); + 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))); + 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); + Tensor<1,4> T4bm = T4bmMom.contract(ABCurr,1); - return T4bm; -} + return T4bm; + } } // Anonymous Namespace helper functions //Functions which can be called elsewhere (declarations in currents.hh). // FKL W+Jets double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later - // Multiply by Cf^2 - return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); + return RHEJ::C_F*RHEJ::C_F*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later - // Multiply by Cf^2 - return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); + return RHEJ::C_F*RHEJ::C_F*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later - // Multiply by Cf^2 - return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); + return RHEJ::C_F*RHEJ::C_F*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later - // Multiply by Cf^2 - return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); + return RHEJ::C_F*RHEJ::C_F*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg->qenug scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj2p,mj2m; jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); 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 nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.; // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double sst = nonflipcolourmult*(a2Mmp+a2Mmm); // double sstsave=sst; // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg->qenug scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj2p,mj2m; jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); double ratio; // p2-/pb- in the notes // if (p2in.plus()>0) // if the gluon is the positive 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 nonflipcolourmult=(1.-1./9.)/2.*(ratio+1./ratio)+1./9.; // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double sst = nonflipcolourmult*(a2Mmp+a2Mmm); // double sstsave=sst; // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } // 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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(); - double cf=4./3.; + 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; } -//Naming scheme jM2-Wuno-g-({q/qbar}{Q/Qbar/g}) +///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; - double ca = 3.; - double cf = 4./3.; + const double ca = RHEJ::C_A; + const double cf = RHEJ::C_F; ///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; - double ca = 3.; - double cf = 4./3.; + const double ca = RHEJ::C_A; + const double cf = RHEJ::C_F; ///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; - double ca = 3.; - double cf = 4./3.; + const double ca = RHEJ::C_A; + const double cf = RHEJ::C_F; ///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; - double ca = 3.; - double cf = 4./3.; + const double ca = RHEJ::C_A; + const double cf = RHEJ::C_F; ///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; } // W+Jets qqxCentral double jM2WqqtoqQQq(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, std::vector 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))); + +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))); + +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))); + +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))); + +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 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))); + +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))); + +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))); + +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))); + +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