Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8724280
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
71 KB
Subscribers
None
View Options
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 3e92a2e..cf5ae44 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,2067 +1,2067 @@
#include "HEJ/currents.hh"
#include "HEJ/utility.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/Constants.hh"
#include <array>
#include <iostream>
namespace { // Helper Functions
// FKL W Helper Functions
void jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur)
{
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
current t65,vout,vin,temp2,temp3,temp5;
joo(pnu,helnu,pe,hele,t65);
vout[0]=pout.e();
vout[1]=pout.x();
vout[2]=pout.y();
vout[3]=pout.z();
vin[0]=pin.e();
vin[1]=pin.x();
vin[2]=pin.y();
vin[3]=pin.z();
COM brac615=cdot(t65,vout);
COM brac645=cdot(t65,vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
// noalias(temp)=prod(trans(CurrentOutOut(pout,helout,pnu,helout)),metric);
joo(pout,helout,pnu,helout,temp2);
// noalias(temp2)=prod(temp,ctemp);
COM prod1665=cdot(temp2,t65);
// noalias(temp)=prod(trans(Current(pe,helin,pin,helin)),metric);
// noalias(temp2)=prod(temp,ctemp);
joi(pe,helin,pin,helin,temp3);
COM prod5465=cdot(temp3,t65);
// noalias(temp)=prod(trans(Current(pnu,helin,pin,helin)),metric);
// noalias(temp2)=prod(temp,ctemp);
joo(pout,helout,pe,helout,temp2);
joi(pnu,helnu,pin,helin,temp3);
joi(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);
// cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb);
// cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb);
cadd(term1,term2,term3,sum);
// std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
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);
// cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb);
// cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb);
cadd(term1,term2,term3,sum);
// std::cout<<"term1: ("<<temp5[0]<<" "<<temp5[1]<<" "<<temp5[2]<<" "<<temp5[3]<<")"<<std::endl;
// std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n";
cur[0]=sum[0];
cur[1]=sum[1];
cur[2]=sum[2];
cur[3]=sum[3];
}
}
double WProp (const CLHEP::HepLorentzVector & plbar, const CLHEP::HepLorentzVector & pl){
COM propW = COM(0.,-1.)/((pl+plbar).m2() -HEJ::MW*HEJ::MW + COM(0.,1.)*HEJ::MW*HEJ::GammaW);
double PropFactor=(propW*conj(propW)).real();
return PropFactor;
}
CCurrent jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin)
{
COM cur[4];
cur[0]=0.;
cur[1]=0.;
cur[2]=0.;
cur[3]=0.;
CCurrent sum(0.,0.,0.,0.);
// NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5)
// Need to swap e and nu for events with W- --> e- nubar!
if (helin==helout && hele==helnu) {
CLHEP::HepLorentzVector qa=pout+pe+pnu;
CLHEP::HepLorentzVector qb=pin-pe-pnu;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pnu,helnu,pe,hele);
CCurrent vout(pout.e(),pout.x(),pout.y(),pout.z());
CCurrent vin(pin.e(),pin.x(),pin.y(),pin.z());
COM brac615=t65.dot(vout);
COM brac645=t65.dot(vin);
// prod1565 and prod6465 are zero for Ws (not Zs)!!
temp2 = joo(pout,helout,pnu,helout);
COM prod1665=temp2.dot(t65);
- temp3 = j(pe,helin,pin,helin);
+ temp3 = joi(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);
+ temp3=joi(pnu,helnu,pin,helin);
+ temp5=joi(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;
}
// Extremal quark current with W emission. Using Tensor class rather than CCurrent
Tensor <1,4> jW(HLV pin, HLV pout, HLV plbar, HLV pl, bool aqline){
// Build the external quark line W Emmision
Tensor<1,4> ABCurr = TCurrent(pl, false, plbar, false);
Tensor<1,4> Tp4W = Construct1Tensor((pout+pl+plbar));//p4+pw
Tensor<1,4> TpbW = Construct1Tensor((pin-pl-plbar));//pb-pw
Tensor<3,4> J4bBlank;
if (aqline){
J4bBlank = T3Current(pin,false,pout,false);
}
else{
J4bBlank = T3Current(pout,false,pin,false);
}
double t4AB = (pout+pl+plbar).m2();
double tbAB = (pin-pl-plbar).m2();
Tensor<2,4> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB;
Tensor<2,4> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB;
Tensor<2,4> T4bmMom(0.);
if (aqline){
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(0,-1)));
}
}
}
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))*(COM(0,1)));
}
}
}
Tensor<1,4> T4bm = T4bmMom.contract(ABCurr,1);
return T4bm;
}
// Relevant W+Jets Unordered Contribution Helper Functions
// W+Jets Uno
double jM2Wuno(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector plbar, CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector pa, bool h1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector pb, bool h2, bool pol)
{
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
//std::cout<<"Setting sigma_index...." << std::endl;
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
CLHEP::HepLorentzVector pW = pl+plbar;
CLHEP::HepLorentzVector q1g=pa-pW-p1-pg;
CLHEP::HepLorentzVector q1 = pa-p1-pW;
CLHEP::HepLorentzVector q2 = p2-pb;
const double taW = (pa-pW).m2();
const double taW1 = (pa-pW-p1).m2();
const double tb2 = (pb-p2).m2();
const double tb2g = (pb-p2-pg).m2();
const double s1W = (p1+pW).m2();
const double s1gW = (p1+pW+pg).m2();
const double s1g = (p1+pg).m2();
const double tag = (pa-pg).m2();
const double taWg = (pa-pW-pg).m2();
//use p1 as ref vec in pol tensor
Tensor<1,4> epsg = eps(pg,p2,pol);
Tensor<1,4> epsW = TCurrent(pl,false,plbar,false);
Tensor<1,4> j2b = TCurrent(p2,h2,pb,h2);
Tensor<1,4> Tq1q2 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+ p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<1,4> Tq1g = Construct1Tensor((-pg-q1))/taW1;
Tensor<1,4> Tq2g = Construct1Tensor((pg-q2));
Tensor<1,4> TqaW = Construct1Tensor((pa-pW));//pa-pw
Tensor<1,4> Tqag = Construct1Tensor((pa-pg));
Tensor<1,4> TqaWg = Construct1Tensor((pa-pg-pW));
Tensor<1,4> Tp1g = Construct1Tensor((p1+pg));
Tensor<1,4> Tp1W = Construct1Tensor((p1+pW));//p1+pw
Tensor<1,4> Tp1gW = Construct1Tensor((p1+pg+pW));//p1+pw+pg
Tensor<2,4> g=Metric();
Tensor<3,4> J31a = T3Current(p1, h1, pa, h1);
Tensor<2,4> J2_qaW =J31a.contract(TqaW/taW, 2);
Tensor<2,4> J2_p1W =J31a.contract(Tp1W/s1W, 2);
Tensor<3,4> L1a =J2_qaW.leftprod(Tq1q2);
Tensor<3,4> L1b =J2_p1W.leftprod(Tq1q2);
Tensor<3,4> L2a = J2_qaW.leftprod(Tq1g);
Tensor<3,4> L2b = J2_p1W.leftprod(Tq1g);
Tensor<3,4> L3 = (g.rightprod(J2_qaW.contract(Tq2g,1)+J2_p1W.contract(Tq2g,2)))/taW1;
Tensor<3,4> L(0.);
Tensor<5,4> J51a = T5Current(p1, h1, pa, h1);
Tensor<4,4> J_qaW = J51a.contract(TqaW,4);
Tensor<4,4> J_qag = J51a.contract(Tqag,4);
Tensor<4,4> J_p1gW = J51a.contract(Tp1gW,4);
Tensor<3,4> U1a = J_qaW.contract(Tp1g,2);
Tensor<3,4> U1b = J_p1gW.contract(Tp1g,2);
Tensor<3,4> U1c = J_p1gW.contract(Tp1W,2);
Tensor<3,4> U1(0.);
Tensor<3,4> U2a = J_qaW.contract(TqaWg,2);
Tensor<3,4> U2b = J_qag.contract(TqaWg,2);
Tensor<3,4> U2c = J_qag.contract(Tp1W,2);
Tensor<3,4> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L.Set(nu, mu, rho, L1a.at(nu,mu,rho) + L1b.at(nu,rho,mu)
+ L2a.at(mu,nu,rho) + L2b.at(mu,rho,nu) + L3.at(mu,nu,rho));
U1.Set(nu, mu, rho, U1a.at(nu, mu, rho) / (s1g*taW)
+ U1b.at(nu,rho,mu) / (s1g*s1gW) + U1c.at(rho,nu,mu) / (s1W*s1gW));
U2.Set(nu,mu,rho,U2a.at(mu,nu,rho) / (taWg*taW)
+ U2b.at(mu,rho,nu) / (taWg*tag) + U2c.at(rho,mu,nu) / (s1W*tag));
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1)).at(0);
double amp = HEJ::C_A*HEJ::C_F*HEJ::C_F/2.*(norm(X)+norm(Y)) - HEJ::C_F/2.*(X*conj(Y)).real();
double t1 = q1g.m2();
double t2 = q2.m2();
double WPropfact = WProp(plbar, pl);
//Divide by WProp
amp*=WPropfact;
//Divide by t-channels
amp/=(t1*t2);
//Average over initial states
amp/=(4.*HEJ::C_A*HEJ::C_A);
return amp;
}
// Relevant Wqqx Helper Functions.
//g->qxqlxl (Calculates gluon to qqx Current. See JV_\mu in WSubleading Notes)
Tensor <1,4> gtqqxW(CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar){
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
//blank 3 Gamma Current
Tensor<3,4> JV23 = T3Current(pq,false,pqbar,false);
// Components of g->qqW before W Contraction
Tensor<2,4> JV1 = JV23.contract((Tpq + TAB),2)/(s2AB);
Tensor<2,4> JV2 = JV23.contract((Tpqbar + TAB),2)/(s3AB);
// g->qqW Current. Note Minus between terms due to momentum flow.
// Also note: (-I)^2 from W vert. (I) from Quark prop.
Tensor<1,4> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2,4> MCrossW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double tcro1=(q3+pq).m2();
double tcro2=(q1-pqbar).m2();
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar,false);
//Blank 5 gamma Current
Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4,4> J_q3q = J523.contract((Tq3+Tpq),2);
Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
// Components of Crossed Vertex Contribution
Tensor<3,4> Xcro1 = J_q3q.contract((Tpqbar + TAB),3);
Tensor<3,4> Xcro2 = J_q3q.contract((Tq1-Tpqbar),3);
Tensor<3,4> Xcro3 = J_2AB.contract((Tq1-Tpqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2,4> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2,4> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2,4> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2,4> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro.Set(mu,nu, -(-Xcro1Cont.at(nu,mu)+Xcro2Cont.at(nu,mu)+Xcro3Cont.at(nu,mu)));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2,4> MUncrossW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0; i<nabove+1;i++){
q1=q1-partons.at(i);
}
q3 = q1 - pl - plbar - pq - pqbar;
double tunc1 = (q1-pq).m2();
double tunc2 = (q3+pqbar).m2();
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// Define llx current.
Tensor<1,4> ABCur = TCurrent(pl, false, plbar, false);
//Blank 5 gamma Current
Tensor<5,4> J523 = T5Current(pq,false,pqbar,false);
// 4 gamma currents (with 1 contraction already).
Tensor<4,4> J_2AB = J523.contract((Tpq+TAB),2);
Tensor<4,4> J_q1q = J523.contract((Tq1-Tpq),2);
// 2 Contractions taken care of.
Tensor<3,4> Xunc1 = J_2AB.contract((Tq3+Tpqbar),3);
Tensor<3,4> Xunc2 = J_q1q.contract((Tq3+Tpqbar),3);
Tensor<3,4> Xunc3 = J_q1q.contract((Tpqbar+TAB),3);
// Term Denominators Taken Care of at this stage
Tensor<2,4> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2,4> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2,4> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2,4> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc.Set(mu,nu,-(- Xunc1Cont.at(mu,nu)+Xunc2Cont.at(mu,nu) +Xunc3Cont.at(mu,nu)));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2,4> MSymW(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, int nabove){
double sa2=(pa+pq).m2();
double s12=(p1+pq).m2();
double sa3=(pa+pqbar).m2();
double s13=(p1+pqbar).m2();
double saA=(pa+pl).m2();
double s1A=(p1+pl).m2();
double saB=(pa+plbar).m2();
double s1B=(p1+plbar).m2();
double sb2=(pb+pq).m2();
double s42=(p4+pq).m2();
double sb3=(pb+pqbar).m2();
double s43=(p4+pqbar).m2();
double sbA=(pb+pl).m2();
double s4A=(p4+pl).m2();
double sbB=(pb+plbar).m2();
double s4B=(p4+plbar).m2();
double s23AB=(pl+plbar+pq+pqbar).m2();
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3=q1-pq-pqbar-plbar-pl;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Define Tensors to be used
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> TAB = Construct1Tensor(pl+plbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1,4> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13+s1A+s1B)) + Tpa*(t1/(sa2+sa3+saA+saB)));
Tensor<2,4> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43+s4A+s4B)) + Tpb*(t3/(sb2+sb3+sbA+sbB)));
Tensor<2,4> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar+TAB);
Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar-TAB);
Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
// Note the contraction of indices changes term by term
Tensor<2,4> X3g1Cont = X3g1.contract(JV,3);
Tensor<2,4> X3g2Cont = X3g2.contract(JV,2);
Tensor<2,4> X3g3Cont = X3g3.contract(JV,1);
// XSym is an amalgamation of x1a, X4b and X3g. Makes sense from a colour factor point of view.
Tensor<2,4>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym.Set(mu,nu, (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu) - X3g3Cont.at(nu,mu))
+ (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) );
}
}
return Xsym/s23AB;
}
Tensor <2,4> MCross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
Tensor<1,4> Tq1 = Construct1Tensor(q1-pqbar);
//Blank 3 gamma Current
Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2,4> XCroCont = J323.contract((Tq1),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2,4> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro.Set(mu,nu, (XCroCont.at(nu,mu)));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2,4> MUncross(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
Tensor<1,4> Tq1 = Construct1Tensor(q1-pq);
//Blank 3 gamma Current
Tensor<3,4> J323 = T3Current(pq,hq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2,4> XUncCont = J323.contract((Tq1),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2,4> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc.Set(mu,nu,-(XUncCont.at(mu,nu)));
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2,4> MSym(CLHEP::HepLorentzVector pa,CLHEP::HepLorentzVector p1,CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector p4, CLHEP::HepLorentzVector pq,CLHEP::HepLorentzVector pqbar, std::vector<HLV> partons, bool hq, int nabove){
CLHEP::HepLorentzVector q1, q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1-pq-pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
double s23 = (pq+pqbar).m2();
double sa2 = (pa+pq).m2();
double sa3 = (pa+pqbar).m2();
double s12 = (p1+pq).m2();
double s13 = (p1+pqbar).m2();
double sb2 = (pb+pq).m2();
double sb3 = (pb+pqbar).m2();
double s42 = (p4+pq).m2();
double s43 = (p4+pqbar).m2();
//Define Tensors to be used
Tensor<1,4> Tp1 = Construct1Tensor(p1);
Tensor<1,4> Tp4 = Construct1Tensor(p4);
Tensor<1,4> Tpa = Construct1Tensor(pa);
Tensor<1,4> Tpb = Construct1Tensor(pb);
Tensor<1,4> Tpq = Construct1Tensor(pq);
Tensor<1,4> Tpqbar = Construct1Tensor(pqbar);
Tensor<1,4> Tq1 = Construct1Tensor(q1);
Tensor<1,4> Tq3 = Construct1Tensor(q3);
Tensor<2,4> g=Metric();
Tensor<1,4> qqxCur = TCurrent(pq, hq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3,4> X1a = g.rightprod(Tp1*(t1/(s12+s13))+Tpa*(t1/(sa2+sa3)));
Tensor<2,4> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3,4> X4b = g.rightprod(Tp4*(t3/(s42+s43)) + Tpb*(t3/(sb2+sb3)));
Tensor<2,4> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3,4> X3g1 = g.leftprod(Tq1+Tpq+Tpqbar);
Tensor<3,4> X3g2 = g.leftprod(Tq3-Tpq-Tpqbar);
Tensor<3,4> X3g3 = g.leftprod((Tq1+Tq3));
// Note the contraction of indices changes term by term
Tensor<2,4> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2,4> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2,4> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2,4>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym.Set(mu, nu, COM(0,1) * ( (X3g1Cont.at(nu,mu) + X3g2Cont.at(mu,nu)
- X3g3Cont.at(nu,mu)) + (X1aCont.at(mu,nu) - X4bCont.at(mu,nu)) ) );
}
}
return Xsym/s23;
}
} // Anonymous Namespace helper functions
// W+Jets FKL Contributions
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);
joi(p2out,true,p2in,true,mj2p);
joi(p2out,false,p2in,false,mj2m);
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);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
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);
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);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
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);
joi(p2out,true,p2in,true,mj2p);
joi(p2out,false,p2in,false,mj2m);
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);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
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);
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);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf^2
return HEJ::C_F*HEJ::C_F*WPropfact*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
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);
joi(p2out,true,p2in,true,mj2p);
joi(p2out,false,p2in,false,mj2m);
// mj1m.mj2p
COM Mmp=cdot(mj1m,mj2p);
// mj1m.mj2m
COM Mmm=cdot(mj1m,mj2m);
const double K = K_g(p2out, p2in);
// sum of spinor strings ||^2
double a2Mmp=abs2(Mmp);
double a2Mmm=abs2(Mmm);
double sst = K/HEJ::C_A*(a2Mmp+a2Mmm);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf*Ca=4
return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
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);
joi(p2out,true,p2in,true,mj2p);
joi(p2out,false,p2in,false,mj2m);
// mj1m.mj2p
COM Mmp=cdot(mj1m,mj2p);
// mj1m.mj2m
COM Mmm=cdot(mj1m,mj2m);
const double K = K_g(p2out, p2in);
// sum of spinor strings ||^2
double a2Mmp=abs2(Mmp);
double a2Mmm=abs2(Mmm);
double sst = K/HEJ::C_A*(a2Mmp+a2Mmm);
double WPropfact = WProp(pe, pnu);
// Division by colour and Helicity average (Nc2-1)(4)
// Multiply by Cf*Ca=4
return HEJ::C_F*HEJ::C_A*WPropfact*sst/(q1.m2()*q2.m2()*(HEJ::N_C*HEJ::N_C - 1)*4);
}
// 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);
+ mj2p=joi(p2out,true,p2in,true);
+ mj2m=joi(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);
+ jgbp=joi(pg,true,p2in,true);
+ jgbm=joi(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 amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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 amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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);
+ mj2p=joi(p2out,true,p2in,true);
+ mj2m=joi(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);
+ jgbp=joi(pg,true,p2in,true);
+ jgbm=joi(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 amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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 amm,amp;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
amp=HEJ::C_F*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mp+U2mp);
double ampsq=-(amm+amp);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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);
+ mj1p=joi(p1out,true,p1in,true);
+ mj1m=joi(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);
+ jgap=joi(pg,true,p1in,true);
+ jgam=joi(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 amm,apm;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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);
+ mj1p=joi(p1out,true,p1in,true);
+ mj1m=joi(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);
+ jgap=joi(pg,true,p1in,true);
+ jgam=joi(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 amm,apm;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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 amm,apm;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// 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 amm,apm;
amm=HEJ::C_F*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1mm+U2mm);
apm=HEJ::C_F*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*HEJ::C_F*HEJ::C_F/3.*vabs2(U1pm+U2pm);
double ampsq=-(apm+amm);
//Divide by WProp
double WPropfact = WProp(pe, pnu);
ampsq*=WPropfact;
// Now add the t-channels
double th=q2.m2()*qg.m2();
ampsq/=th;
ampsq/=16.;
return ampsq;
}
///TODO make this comment more visible
/// Naming scheme jM2-Wuno-g-({q/qbar}{Q/Qbar/g})
///TODO Spit naming for more complicated functions?
/// e.g. jM2WqqtoqQQq -> jM2_Wqq_to_qQQq
double jM2WunogqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
//same as function above but actually obtaining the antiquark line by crossing symmetry, where p1out and p1in are expected to be negative.
//should give same result as jM2WunogqbarQ below (verified)
double jM2WunogqQ_crossqQ(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2WunogqQbar(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
return ME2;
}
double jM2Wunogqg(CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in)
{
//COM temp;
double ME2mpp=0.;
double ME2mpm=0.;
double ME2mmp=0.;
double ME2mmm=0.;
double ME2;
ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,true);
ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,true,false);
ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,true);
ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,false,p2out,p2in,false,false);
//Helicity sum
ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
double ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
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 ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
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 ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
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 ratio; // p2-/pb- in the notes
if (p2in.pz()>0.) // if the gluon is the positive
ratio=p2out.plus()/p2in.plus();
else // the gluon is the negative
ratio=p2out.minus()/p2in.minus();
double cam = ( (HEJ::C_A - 1/HEJ::C_A)*(ratio + 1./ratio)/2. + 1/HEJ::C_A)/HEJ::C_F;
ME2*=cam;
//Correct colour averaging
ME2*=(3.0/8.0);
return ME2;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1,4> qggm1(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
double t1 = (p3-pb)*(p3-pb);
Tensor<1,4> Tp3 = Construct1Tensor((p3));//p3
Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1,4> epsg = eps(pb, refmom, helg);
Tensor<3,4> qqCurBlank = T3Current(p2,hel2,p3,hel2);
Tensor<2,4> qqCur = qqCurBlank.contract(Tp3-Tpb,2);
Tensor<1,4> gqqCur = qqCur.contract(epsg,2)/t1;
return gqqCur*(-1);
}
//Function to calculate Term 2 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1,4> qggm2(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
double t1 = (p2-pb)*(p2-pb);
Tensor<1,4> Tp2 = Construct1Tensor((p2));//p2
Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1,4> epsg = eps(pb,refmom, helg);
Tensor<3,4> qqCurBlank = T3Current(p2,hel2,p3,hel2);
Tensor<2,4> qqCur = qqCurBlank.contract(Tp2-Tpb,2);
Tensor<1,4> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1,4> qggm3(CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool hel2, bool helg, CLHEP::HepLorentzVector refmom){
double s23 = (p2+p3)*(p2+p3);
Tensor<1,4> Tp2 = Construct1Tensor((p2));//p2
Tensor<1,4> Tp3 = Construct1Tensor((p3));//p3
Tensor<1,4> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1,4> epsg = eps(pb, refmom, helg);
Tensor<2,4> g=Metric();
Tensor<3,4> qqCurBlank1 = g.leftprod(Tp2+Tp3)/s23;
Tensor<3,4> qqCurBlank2 = g.leftprod(Tpb)/s23;
Tensor<1,4> Cur23 = TCurrent(p2,hel2, p3,hel2);
Tensor<2,4> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2,4> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2,4> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1,4> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double jM2WgqtoQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3,CLHEP::HepLorentzVector plbar,CLHEP::HepLorentzVector pl, bool aqlinepa){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;}
// 2 independent helicity choices (complex conjugation related).
Tensor<1,4> TMmmm1 = qggm1(pb,p2,p3,false,false, pa);
Tensor<1,4> TMmmm2 = qggm2(pb,p2,p3,false,false, pa);
Tensor<1,4> TMmmm3 = qggm3(pb,p2,p3,false,false, pa);
Tensor<1,4> TMpmm1 = qggm1(pb,p2,p3,false,true, pa);
Tensor<1,4> TMpmm2 = qggm2(pb,p2,p3,false,true, pa);
Tensor<1,4> TMpmm3 = qggm3(pb,p2,p3,false,true, pa);
// Build the external quark line W Emmision
Tensor<1,4> cur1a = jW(pa,p1,plbar,pl, aqlinepa);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1).at(0);
COM Mmmm2 = TMmmm2.contract(cur1a,1).at(0);
COM Mmmm3 = TMmmm3.contract(cur1a,1).at(0);
COM Mpmm1 = TMpmm1.contract(cur1a,1).at(0);
COM Mpmm2 = TMpmm2.contract(cur1a,1).at(0);
COM Mpmm3 = TMpmm3.contract(cur1a,1).at(0);
//Colour factors:
COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
cm1m1=8./3.;
cm2m2=8./3.;
cm3m3=6.;
cm1m2 =-1./3.;
cm1m3 = -3.*COM(0.,1.);
cm2m3 = 3.*COM(0.,1.);
//Sqaure and sum for each helicity config:
double Mmmm = real(cm1m1*pow(abs(Mmmm1),2)+cm2m2*pow(abs(Mmmm2),2)+cm3m3*pow(abs(Mmmm3),2)+2.*real(cm1m2*Mmmm1*conj(Mmmm2))+2.*real(cm1m3*Mmmm1*conj(Mmmm3))+2.*real(cm2m3*Mmmm2*conj(Mmmm3)));
double Mpmm = real(cm1m1*pow(abs(Mpmm1),2)+cm2m2*pow(abs(Mpmm2),2)+cm3m3*pow(abs(Mpmm3),2)+2.*real(cm1m2*Mpmm1*conj(Mpmm2))+2.*real(cm1m3*Mpmm1*conj(Mpmm3))+2.*real(cm2m3*Mpmm2*conj(Mpmm3)));
// Divide by WProp
double WPropfact = WProp(plbar, pl);
return (2*WPropfact*(Mmmm+Mpmm)/24./4.)/(pa-p1-pl-plbar).m2()/(p2+p3-pb).m2();
}
// W+Jets qqxCentral
double jM2WqqtoqQQq(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl, CLHEP::HepLorentzVector plbar, std::vector<HLV> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove)
{
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;}
HLV pq, pqbar, p1, p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1,4> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
if(!(aqlinepb)){
T4bp = TCurrent(p4, true, pb, true);
T4bm = TCurrent(p4, false, pb, false);}
else if(aqlinepb){
T4bp = TCurrent(pb, true, p4, true);
T4bm = TCurrent(pb, false, p4, false);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2,4> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2,4> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2,4> Xsym = MSymW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
// 4 Different Helicity Choices (Differs from Pure Jet Case, where there is also the choice in qqbar helicity.
// (- - hel choice)
COM M_mmUnc = (((Xunc).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1)).at(0);
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1)).at(0);
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1)).at(0);
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1)).at(0);
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1)).at(0);
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1)).at(0);
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1)).at(0);
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1)).at(0);
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar - pl - plbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
return amp;
}
// no wqq emission
double jM2WqqtoqQQqW(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb,CLHEP::HepLorentzVector pl,CLHEP::HepLorentzVector plbar, std::vector<CLHEP::HepLorentzVector> partons, bool aqlinepa, bool aqlinepb, bool qqxmarker, int nabove, int nbelow, bool forwards){
static bool is_sigma_index_set(false);
if(!is_sigma_index_set){
if(init_sigma_index())
is_sigma_index_set = true;
else
return 0.;
}
if (!forwards){ //If Emission from Leg a instead, flip process.
HLV dummymom = pa;
bool dummybool= aqlinepa;
int dummyint = nabove;
pa = pb;
pb = dummymom;
std::reverse(partons.begin(),partons.end());
qqxmarker = !(qqxmarker);
aqlinepa = aqlinepb;
aqlinepb = dummybool;
nabove = nbelow;
nbelow = dummyint;
}
HLV pq, pqbar, p1,p4;
if (qqxmarker){
pqbar = partons[nabove+1];
pq = partons[nabove+2];}
else{
pq = partons[nabove+1];
pqbar = partons[nabove+2];}
p1 = partons.front();
p4 = partons.back();
Tensor<1,4> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
T1ap = TCurrent(p1, true, pa, true);
T1am = TCurrent(p1, false, pa, false);}
else if(aqlinepa){
T1ap = TCurrent(pa, true, p1, true);
T1am = TCurrent(pa, false, p1, false);}
Tensor <1,4> T4bm = jW(pb, p4, plbar, pl, aqlinepb);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2,4> Xunc_m = MUncross(pa, pq, pqbar,partons, false, nabove);
Tensor<2,4> Xcro_m = MCross( pa, pq, pqbar,partons, false, nabove);
Tensor<2,4> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, false, nabove);
Tensor<2,4> Xunc_p = MUncross(pa, pq, pqbar,partons, true, nabove);
Tensor<2,4> Xcro_p = MCross( pa, pq, pqbar,partons, true, nabove);
Tensor<2,4> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, true, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1)).at(0);
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1)).at(0);
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1)).at(0);
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1)).at(0);
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1)).at(0);
//Colour factors:
COM cmsms,cmumu,cmcmc,cmsmu,cmsmc,cmumc;
cmsms=3.;
cmumu=4./3.;
cmcmc=4./3.;
cmsmu =3./2.*COM(0.,1.);
cmsmc = -3./2.*COM(0.,1.);
cmumc = -1./6.;
// Work Out Interference in each case of helicity:
double amp_mm = real(cmsms*pow(abs(M_mmSym),2)
+cmumu*pow(abs(M_mmUnc),2)
+cmcmc*pow(abs(M_mmCro),2)
+2.*real(cmsmu*M_mmSym*conj(M_mmUnc))
+2.*real(cmsmc*M_mmSym*conj(M_mmCro))
+2.*real(cmumc*M_mmUnc*conj(M_mmCro)));
double amp_mp = real(cmsms*pow(abs(M_mpSym),2)
+cmumu*pow(abs(M_mpUnc),2)
+cmcmc*pow(abs(M_mpCro),2)
+2.*real(cmsmu*M_mpSym*conj(M_mpUnc))
+2.*real(cmsmc*M_mpSym*conj(M_mpCro))
+2.*real(cmumc*M_mpUnc*conj(M_mpCro)));
double amp_pm = real(cmsms*pow(abs(M_pmSym),2)
+cmumu*pow(abs(M_pmUnc),2)
+cmcmc*pow(abs(M_pmCro),2)
+2.*real(cmsmu*M_pmSym*conj(M_pmUnc))
+2.*real(cmsmc*M_pmSym*conj(M_pmCro))
+2.*real(cmumc*M_pmUnc*conj(M_pmCro)));
double amp_pp = real(cmsms*pow(abs(M_ppSym),2)
+cmumu*pow(abs(M_ppUnc),2)
+cmcmc*pow(abs(M_ppCro),2)
+2.*real(cmsmu*M_ppSym*conj(M_ppUnc))
+2.*real(cmsmc*M_ppSym*conj(M_ppCro))
+2.*real(cmumc*M_ppUnc*conj(M_ppCro)));
double amp=((amp_mm+amp_mp+amp_pm+amp_pp)/(9.*4.));
CLHEP::HepLorentzVector q1,q3;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
q3 = q1 - pq - pqbar;
double t1 = (q1).m2();
double t3 = (q3).m2();
//Divide by t-channels
amp/=(t1*t1*t3*t3);
//Divide by WProp
double WPropfact = WProp(plbar, pl);
amp*=WPropfact;
return amp;
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 10:59 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4201333
Default Alt Text
(71 KB)
Attached To
rHEJ HEJ
Event Timeline
Log In to Comment