No OneTemporary
Actions

Size

23 KB

Subscribers

None

View Options

	diff --git a/src/jets.cc b/src/jets.cc
	index 48dcffb..05572a5 100644
	--- a/src/jets.cc
	+++ b/src/jets.cc
	@@ -1,388 +1,704 @@
	/**
	* \authors The HEJ collaboration (see AUTHORS for details)
	* \date 2019
	* \copyright GPLv2 or later
	*/
	#include "HEJ/jets.hh"

	#include "HEJ/Constants.hh"

	// Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113
	// @TODO: this is not a current and should be moved somewhere else
	double K_g(double p1minus, double paminus) {
	return 1./2.(p1minus/paminus + paminus/p1minus)(HEJ::C_A - 1./HEJ::C_A) + 1./HEJ::C_A;
	}
	double K_g(
	HLV const & pout,
	HLV const & pin
	) {
	if(pin.z() > 0) return K_g(pout.plus(), pin.plus());
	return K_g(pout.minus(), pin.minus());
	}

	CCurrent CCurrent::operator+(const CCurrent& other)
	{
	COM result_c0=c0 + other.c0;
	COM result_c1=c1 + other.c1;
	COM result_c2=c2 + other.c2;
	COM result_c3=c3 + other.c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator-(const CCurrent& other)
	{
	COM result_c0=c0 - other.c0;
	COM result_c1=c1 - other.c1;
	COM result_c2=c2 - other.c2;
	COM result_c3=c3 - other.c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator*(const double x)
	{
	COM result_c0=x*CCurrent::c0;
	COM result_c1=x*CCurrent::c1;
	COM result_c2=x*CCurrent::c2;
	COM result_c3=x*CCurrent::c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator/(const double x)
	{
	COM result_c0=CCurrent::c0/x;
	COM result_c1=CCurrent::c1/x;
	COM result_c2=CCurrent::c2/x;
	COM result_c3=CCurrent::c3/x;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator*(const COM x)
	{
	COM result_c0=x*CCurrent::c0;
	COM result_c1=x*CCurrent::c1;
	COM result_c2=x*CCurrent::c2;
	COM result_c3=x*CCurrent::c3;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	CCurrent CCurrent::operator/(const COM x)
	{
	COM result_c0=(CCurrent::c0)/x;
	COM result_c1=(CCurrent::c1)/x;
	COM result_c2=(CCurrent::c2)/x;
	COM result_c3=(CCurrent::c3)/x;

	return CCurrent(result_c0,result_c1,result_c2,result_c3);
	}

	std::ostream& operator <<(std::ostream& os, const CCurrent& cur)
	{
	os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")";
	return os;
	}

	CCurrent operator * ( double x, CCurrent& m)
	{
	return m*x;
	}

	CCurrent operator * ( COM x, CCurrent& m)
	{
	return m*x;
	}

	CCurrent operator / ( double x, CCurrent& m)
	{
	return m/x;
	}

	CCurrent operator / ( COM x, CCurrent& m)
	{
	return m/x;
	}

	COM CCurrent::dot(HLV p1)
	{
	// Current goes (E,px,py,pz)
	// Vector goes (px,py,pz,E)
	return p1[3]c0-p1[0]c1-p1[1]c2-p1[2]c3;
	}

	COM CCurrent::dot(CCurrent p1)
	{
	return p1.c0c0-p1.c1c1-p1.c2c2-p1.c3c3;
	}

	//Current Functions
	void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) {
	cur[0]=0.;
	cur[1]=0.;
	cur[2]=0.;
	cur[3]=0.;

	const double sqpop = sqrt(pout.plus());
	const double sqpom = sqrt(pout.minus());
	const COM poperp = pout.x() + COM(0, 1) * pout.y();

	if (helout != helin) {
	throw std::invalid_argument{"Non-matching helicities"};
	} else if (helout == false) { // negative helicity
	if (pin.plus() > pin.minus()) { // if forward
	const double sqpip = sqrt(pin.plus());
	cur[0] = sqpop * sqpip;
	cur[1] = sqpom * sqpip * poperp / abs(poperp);
	cur[2] = -COM(0,1) * cur[1];
	cur[3] = cur[0];
	} else { // if backward
	const double sqpim = sqrt(pin.minus());
	cur[0] = -sqpom * sqpim * poperp / abs(poperp);
	cur[1] = -sqpim * sqpop;
	cur[2] = COM(0,1) * cur[1];
	cur[3] = -cur[0];
	}
	} else { // positive helicity
	if (pin.plus() > pin.minus()) { // if forward
	const double sqpip = sqrt(pin.plus());
	cur[0] = sqpop * sqpip;
	cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp);
	cur[2] = COM(0,1) * cur[1];
	cur[3] = cur[0];
	} else { // if backward
	const double sqpim = sqrt(pin.minus());
	cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp);
	cur[1] = -sqpim * sqpop;
	cur[2] = -COM(0,1) * cur[1];
	cur[3] = -cur[0];
	}
	}
	}

	CCurrent joi (HLV pout, bool helout, HLV pin, bool helin)
	{
	current cur;
	joi(pout, helout, pin, helin, cur);
	return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	}

	void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) {
	joi(pout, !helout, pin, !helin, cur);
	}

	CCurrent jio (HLV pin, bool helin, HLV pout, bool helout)
	{
	current cur;
	jio(pin, helin, pout, helout, cur);
	return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	}

	void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) {

	// Zero our current
	cur[0] = 0.0;
	cur[1] = 0.0;
	cur[2] = 0.0;
	cur[3] = 0.0;
	if (heli!=helj) {
	throw std::invalid_argument{"Non-matching helicities"};
	} else if ( heli == true ) { // If positive helicity swap momenta
	std::swap(pi,pj);
	}

	const double sqpjp = sqrt(pj.plus());
	const double sqpjm = sqrt(pj.minus());
	const double sqpip = sqrt(pi.plus());
	const double sqpim = sqrt(pi.minus());

	const COM piperp = pi.x() + COM(0,1) * pi.y();
	const COM pjperp = pj.x() + COM(0,1) * pj.y();
	const COM phasei = piperp / abs(piperp);
	const COM phasej = pjperp / abs(pjperp);

	cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
	cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej);
	cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej));
	cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp;
	}

	CCurrent joo (HLV pi, bool heli, HLV pj, bool helj)
	{
	current cur;
	joo(pi, heli, pj, helj, cur);
	return CCurrent(cur[0],cur[1],cur[2],cur[3]);
	}
	namespace{
	//@{
	/**
	* @brief Pure Jet FKL Contributions, function to handle all incoming types.
	* @param p1out Outgoing Particle 1.
	* @param p1in Incoming Particle 1.
	* @param p2out Outgoing Particle 2
	* @param p2in Incoming Particle 2
	*
	* Calculates j_\mu j^\mu.
	* Handles all possible incoming states. Helicity doesn't matter since we sum
	* over all of them.
	*/
	double j_j(HLV const & p1out, HLV const & p1in,
	HLV const & p2out, HLV const & p2in
	){
	HLV const q1=p1in-p1out;
	HLV const q2=-(p2in-p2out);
	current mj1m,mj1p,mj2m,mj2p;

	// Note need to flip helicities in anti-quark case.
	joi(p1out, false, p1in, false, mj1p);
	joi(p1out, true, p1in, true, mj1m);
	joi(p2out, false, p2in, false, mj2p);
	joi(p2out, true, p2in, true, mj2m);

	COM const Mmp=cdot(mj1m,mj2p);
	COM const Mmm=cdot(mj1m,mj2m);
	COM const Mpp=cdot(mj1p,mj2p);
	COM const Mpm=cdot(mj1p,mj2m);

	double const sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm);

	// Multiply by Cf^2
	return HEJ::C_FHEJ::C_F(sst)/(q1.m2()*q2.m2());
	}
	} //anonymous namespace
	double ME_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in);
	}

	double ME_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in);
	}

	double ME_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in);
	}

	double ME_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/HEJ::C_F;
	}

	double ME_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in)*K_g(p2out, p2in)/(HEJ::C_F);
	}

	double ME_gg(HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return j_j(p1out, p1in, p2out, p2in)K_g(p1out, p1in)K_g(p2out, p2in)/(HEJ::C_F*HEJ::C_F);
	}
	//@}

	namespace{
	double juno_j(HLV const & pg, HLV const & p1out,
	HLV const & p1in, HLV const & p2out, HLV const & p2in
	){
	// This construction is taking rapidity order: pg > p1out >> p2out
	HLV q1=p1in-p1out; // Top End
	HLV q2=-(p2in-p2out); // Bottom End
	HLV qg=p1in-p1out-pg; // Extra bit post-gluon

	// Note <p1\|eps\|pa> current split into two by gauge choice.
	// See James C's Thesis (p72). <p1\|eps\|pa> -> <p1\|pg><pg\|pa>
	CCurrent mj1p=joi(p1out, false, p1in, false);
	CCurrent mj1m=joi(p1out, true, p1in, true);
	CCurrent jgap=joi(pg, false, p1in, false);
	CCurrent jgam=joi(pg, true, p1in, true);

	// Note for function joo(): <p1+\|pg+> = <pg-\|p1->.
	CCurrent j2gp=joo(p1out, false, pg, false);
	CCurrent j2gm=joo(p1out, true, pg, true);

	CCurrent mj2p=joi(p2out, false, p2in, false);
	CCurrent mj2m=joi(p2out, true, p2in, true);

	// Dot products of these which occur again and again
	COM Mmp=mj1m.dot(mj2p);
	COM Mmm=mj1m.dot(mj2m);
	COM Mpp=mj1p.dot(mj2p);
	COM Mpm=mj1p.dot(mj2m);

	CCurrent p1o(p1out),p2o(p2out),p2i(p2in),qsum(q1+qg),p1i(p1in);

	CCurrent Lmm=(qsum(Mmm)+(-2.mj2m.dot(pg))mj1m+2.mj1m.dot(pg)*mj2m
	+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()Mmm/2.))/q1.m2();
	CCurrent Lmp=(qsum(Mmp) + (-2.mj2p.dot(pg))mj1m+2.mj1m.dot(pg)*mj2p
	+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()Mmp/2.))/q1.m2();
	CCurrent Lpm=(qsum(Mpm) + (-2.mj2m.dot(pg))mj1p+2.mj1p.dot(pg)*mj2m
	+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()Mpm/2.))/q1.m2();
	CCurrent Lpp=(qsum(Mpp) + (-2.mj2p.dot(pg))mj1p+2.mj1p.dot(pg)*mj2p
	+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))(qg.m2()Mpp/2.))/q1.m2();

	CCurrent U1mm=(jgam.dot(mj2m)j2gm+2.p1o*Mmm)/(p1out+pg).m2();
	CCurrent U1mp=(jgam.dot(mj2p)j2gm+2.p1o*Mmp)/(p1out+pg).m2();
	CCurrent U1pm=(jgap.dot(mj2m)j2gp+2.p1o*Mpm)/(p1out+pg).m2();
	CCurrent U1pp=(jgap.dot(mj2p)j2gp+2.p1o*Mpp)/(p1out+pg).m2();
	CCurrent U2mm=((-1.)j2gm.dot(mj2m)jgam+2.p1iMmm)/(p1in-pg).m2();
	CCurrent U2mp=((-1.)j2gm.dot(mj2p)jgam+2.p1iMmp)/(p1in-pg).m2();
	CCurrent U2pm=((-1.)j2gp.dot(mj2m)jgap+2.p1iMpm)/(p1in-pg).m2();
	CCurrent U2pp=((-1.)j2gp.dot(mj2p)jgap+2.p1iMpp)/(p1in-pg).m2();

	constexpr double cf=HEJ::C_F;

	double amm=cf(2.vre(Lmm-U1mm,Lmm+U2mm))+2.cfcf/3.*vabs2(U1mm+U2mm);
	double amp=cf(2.vre(Lmp-U1mp,Lmp+U2mp))+2.cfcf/3.*vabs2(U1mp+U2mp);
	double apm=cf(2.vre(Lpm-U1pm,Lpm+U2pm))+2.cfcf/3.*vabs2(U1pm+U2pm);
	double app=cf(2.vre(Lpp-U1pp,Lpp+U2pp))+2.cfcf/3.*vabs2(U1pp+U2pp);
	double ampsq=-(amm+amp+apm+app);

	//Divide by t-channels
	ampsq/=q2.m2()*qg.m2();
	ampsq/=16.;

	// Factor of (Cf/Ca) for each quark to match j_j.
	ampsq=(HEJ::C_FHEJ::C_F)/(HEJ::C_A*HEJ::C_A);

	return ampsq;

	}
	}

	//Unordered bits for pure jet
	double ME_unob_qQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in);
	}

	double ME_unob_qbarQ (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in);
	}

	double ME_unob_qQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in);
	}

	double ME_unob_qbarQbar (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in);
	}

	double ME_unob_qg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
	}

	double ME_unob_qbarg (HLV pg, HLV p1out, HLV p1in, HLV p2out, HLV p2in){
	return juno_j(pg, p1out, p1in, p2out, p2in)*K_g(p2out,p2in)/HEJ::C_F;
	}

	+
	+// Extremal qqx
	double ME_Exqqx_qbarqQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
	return juno_j(-pgin, pqout, -pqbarout, p2out, p2in)*(3.0/8.0);
	}

	double ME_Exqqx_qqbarQ(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
	return juno_j(-pgin, pqbarout, -pqout, p2out, p2in)*(3.0/8.0);
	}

	double ME_Exqqx_qbarqg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
	return juno_j(-pgin, pqout, -pqbarout, p2out, p2in)(3.0/8.0)K_g(p2out,p2in)/HEJ::C_F;
	}

	double ME_Exqqx_qqbarg(HLV pgin, HLV pqout, HLV pqbarout, HLV p2out, HLV p2in){
	return juno_j(-pgin, pqbarout, -pqout, p2out, p2in)(3.0/8.0)K_g(p2out,p2in)/HEJ::C_F;
	}
	+
	+// qg -> qQQ~
	+
	+//First, a function for generating polarisation tensors. Output as 'current'. Should be general for any refmom now that I've added a bit to the j function.
	+
	+void eps(CLHEP::HepLorentzVector refmom, CLHEP::HepLorentzVector kb, bool hel, current &ep){
	+ current curm,curp;
	+ //Recall - positive helicity eps has negative helicity choices for spinors and vice versa
	+ joi(refmom,true,kb,true,curm);
	+ joi(refmom,false,kb,false,curp);
	+ double norm=1.;
	+ if(kb.z()<0.)
	+ norm = sqrt(2.refmom.plus()*kb.minus());
	+ if(kb.z()>0.)
	+ norm = sqrt(2.refmom.minus()kb.plus());
	+//cout << norm << endl;
	+ if(hel==false){
	+ ep[0] = curm[0]/norm;
	+ ep[1] = curm[1]/norm;
	+ ep[2] = curm[2]/norm;
	+ ep[3] = curm[3]/norm;
	+ }
	+ if(hel==true){
	+ ep[0] = curp[0]/norm;
	+ ep[1] = curp[1]/norm;
	+ ep[2] = curp[2]/norm;
	+ ep[3] = curp[3]/norm;
	+ }
	+ //cout << ep[0] << " " << ep[1] << " " << ep[2] <<" " << ep[3] << endl;
	+}
	+
	+//Now build up each part of the sqaured amplitude
	+
	+COM qggm1(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom, bool aqline){
	+ //Since everything is defined with currents, need to use compeleness relation to expand p slash. i.e. pslash = \|p><p\|. Only one helicity 'survives' as defined by the helicities of the spinors at the end of the chain.
	+ current cur33, cur23, curb3, cur2b, cur1a, ep;
	+ joo(p3, helchain, p3, helchain,cur33);
	+ joo(p2,helchain,p3,helchain,cur23);
	+ jio(pb,helchain,p3,helchain,curb3);
	+ joi(p2,helchain,pb,helchain,cur2b);
	+ if(aqline==true)
	+ jio(pa, heltop, p1, heltop,cur1a);
	+ if(aqline==false)
	+ joi(p1, heltop, pa, heltop,cur1a);
	+
	+ double t2 = (p3-pb)*(p3-pb);
	+ //Create vertex
	+ COM v1[4][4];
	+ for(int u=0; u<4;u++)
	+ {
	+ for(int v=0; v<4; v++)
	+ {
	+ v1[u][v]=(cur23[u]cur33[v]-cur2b[u]curb3[v])/t2*(-1.);
	+ }
	+ }
	+ //Dot in current and eps
	+ //Metric tensor
	+ double eta[4][4]={};
	+ eta[0][0]=1.;
	+ eta[1][1]=-1.;
	+ eta[2][2]=-1.;
	+ eta[3][3]=-1.;
	+ //eps
	+ eps(refmom,pb,helb, ep);
	+ COM M1=0.;
	+ for(int i=0;i<4;i++){
	+ for(int j=0;j<4;j++){
	+ for(int k=0; k<4; k++){
	+ for(int l=0; l<4;l++){
	+ M1+= eta[i][k]cur1a[k](v1[i][j])ep[l]eta[l][j];
	+ }
	+ }
	+ }
	+ }
	+ return M1;
	+ //return v1[0][0];
	+}
	+
	+COM qggm2(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom, bool aqline){
	+ //Since everything is defined with currents, need to use compeleness relation to expand p slash. i.e. pslash = \|p><p\|. Only one helicity 'survives' as defined by the helicities of the spinors at the end of the chain.
	+ current cur22, cur23, curb3, cur2b, cur1a, ep;
	+ joo(p2, helchain, p2, helchain,cur22);
	+ joo(p2,helchain,p3,helchain,cur23);
	+ jio(pb,helchain,p3,helchain,curb3);
	+ joi(p2,helchain,pb,helchain,cur2b);
	+
	+ if(aqline==true)
	+ jio(pa, heltop, p1, heltop,cur1a);
	+ if(aqline==false)
	+ joi(p1, heltop, pa, heltop,cur1a);
	+
	+ double t2t = (p2-pb)*(p2-pb);
	+ //Create vertex
	+ COM v2[4][4]={};
	+ for(int u=0; u<4;u++)
	+ {
	+ for(int v=0; v<4; v++)
	+ {
	+ v2[u][v]=(cur22[v]cur23[u]-cur2b[v]curb3[u])/t2t;
	+ }
	+ }
	+ //Dot in current and eps
	+ //Metric tensor
	+ double eta[4][4]={};
	+ eta[0][0]=1.;
	+ eta[1][1]=-1.;
	+ eta[2][2]=-1.;
	+ eta[3][3]=-1.;
	+ //eps
	+ eps(refmom,pb,helb, ep);
	+ COM M2=0.;
	+ for(int i=0;i<4;i++){
	+ for(int j=0;j<4;j++){
	+ for(int k=0; k<4; k++){
	+ for(int l=0; l<4;l++){
	+ M2+= eta[i][k]cur1a[k](v2[i][j])ep[l]eta[l][j];
	+ }
	+ }
	+ }
	+ }
	+ return M2;
	+ //return v2[0][0];
	+}
	+
	+COM qggm3(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool helchain, bool heltop, bool helb,CLHEP::HepLorentzVector refmom, bool aqline){
	+ //3 gluon vertex bit
	+ //Metric tensor
	+ double eta[4][4]={};
	+ eta[0][0]=1.;
	+ eta[1][1]=-1.;
	+ eta[2][2]=-1.;
	+ eta[3][3]=-1.;
	+ current spincur,ep,cur1a;
	+ double s23 = (p2+p3)*(p2+p3);
	+ //double t1 = (p2+p3-pb)*(p2+p3-pb);
	+ //double sab= (pa+pb)*(pa+pb);
	+ joo(p2,helchain,p3,helchain,spincur);
	+ if(aqline==true)
	+ jio(pa, heltop, p1, heltop,cur1a);
	+ if(aqline==false)
	+ joi(p1, heltop, pa, heltop,cur1a);
	+ //Redefine relevant momenta as currents - for ease of calling correct part of vector
	+ current ka,k2,k3,kb;
	+ kb[0]=pb.e();
	+ kb[1]=pb.x();
	+ kb[2]=pb.y();
	+ kb[3]=pb.z();
	+ k2[0]=p2.e();
	+ k2[1]=p2.x();
	+ k2[2]=p2.y();
	+ k2[3]=p2.z();
	+ k3[0]=p3.e();
	+ k3[1]=p3.x();
	+ k3[2]=p3.y();
	+ k3[3]=p3.z();
	+ ka[0]=pa.e();
	+ ka[1]=pa.x();
	+ ka[2]=pa.y();
	+ ka[3]=pa.z();
	+ COM V3g[4][4]={};
	+ for(int u=0;u<4;u++){
	+ for(int v=0;v<4;v++){
	+ for(int p=0;p<4;p++){
	+ for(int r=0; r<4;r++){
	+ V3g[u][v] += COM(0.,1.)(((2.k2[v]+2.k3[v])eta[u][p] - (2.kb[u])eta[p][v]+2.kb[p]eta[u][v])spincur[r]eta[r][p])/s23;
	+ }
	+ }
	+ }
	+ }
	+ //Needed to add an extra eta to get the correct index strucutre
	+ //Extra bit for s and u diagrams
	+ /*COM extrafact;
	+ extrafact = 2.((sqrt(p2.minus()p3.minus())(p2.perp()conj(p3.perp()))/abs(p2.perp())/abs(p3.perp()))/pb.minus())*t1/s23;
	+ if(helchain==false){
	+ extrafact = conj(extrafact);
	+ }
	+ COM extrabit[4][4]={};
	+ for(int u=0;u<4;u++) {
	+ for(int v=0;v<4;v++) {
	+ extrabit[u][v] = COM(0.,1.)extrafacteta[u][v];
	+ }
	+ }*/
	+ //Alternative extra bit - I will also choose the gauge such that this part vanishes
	+ COM diffextrabit[4][4]={};
	+ //If pa is used as a reference vector, this extra bit doesn't contriubtue. I will always use this choice.
	+ //if(refmom==pa)
	+ //{
	+ for(int u=0;u<4;u++) {
	+ for(int v=0;v<4;v++) {
	+ diffextrabit[u][v] = 0.;
	+ }
	+ }
	+ //}
	+ //Otherwise
	+ //else{
	+ //for(int u=0;u<4;u++) {
	+ // for(int v=0;v<4;v++) {
	+ // diffextrabit[u][v] = COM(0.,1.)(spincur[u]2.ka[v])t1/s23/sab;
	+ // }
	+ //}
	+ //}
	+
	+ //Dot in current and eps
	+ //eps
	+ eps(refmom,pb,helb, ep);
	+ COM M3=0.;
	+ for(int i=0;i<4;i++){
	+ for(int j=0;j<4;j++){
	+ for(int k=0; k<4; k++){
	+ for(int l=0; l<4;l++){
	+ M3+= eta[i][k]cur1a[k](V3g[i][j]+diffextrabit[i][j])ep[l]eta[l][j];
	+ }
	+ }
	+ }
	+ }
	+ return M3;
	+ //return V3g[0][0];
	+}
	+
	+//Now the function to give helicity/colour sum/average
	+double MqgtqQQ(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p3, bool aqline, bool qqxmarker)
	+{
	+
	+
	+ //If qqxmarker is true, switch the order of quark and anti-quark
	+ CLHEP::HepLorentzVector ka,kb,k1,k2,k3;
	+
	+ if(qqxmarker==true){
	+ ka=pa;
	+ kb=pb;
	+ k1=p1;
	+ k2=p3;
	+ k3=p2;
	+ } else {
	+ ka=pa;
	+ kb=pb;
	+ k1=p1;
	+ k2=p2;
	+ k3=p3;
	+ }
	+
	+ // 4 indepedent helicity choices (complex conjugation related).
	+
	+ //Need to evalute each independent hel configuration and store that result somewhere
	+
	+ COM Mmmm1 = qggm1(ka,kb,k1,k2,k3,false,false,false, ka, aqline);
	+ COM Mmmm2 = qggm2(ka,kb,k1,k2,k3,false,false,false, ka, aqline);
	+ COM Mmmm3 = qggm3(ka,kb,k1,k2,k3,false,false,false, ka, aqline);
	+ COM Mmmp1 = qggm1(ka,kb,k1,k2,k3,false,true,false, ka, aqline);
	+ COM Mmmp2 = qggm2(ka,kb,k1,k2,k3,false,true,false, ka, aqline);
	+ COM Mmmp3 = qggm3(ka,kb,k1,k2,k3,false,true,false, ka, aqline);
	+ COM Mpmm1 = qggm1(ka,kb,k1,k2,k3,false,false,true, ka, aqline);
	+ COM Mpmm2 = qggm2(ka,kb,k1,k2,k3,false,false,true, ka, aqline);
	+ COM Mpmm3 = qggm3(ka,kb,k1,k2,k3,false,false,true, ka, aqline);
	+ COM Mpmp1 = qggm1(ka,kb,k1,k2,k3,false,true,true, ka, aqline);
	+ COM Mpmp2 = qggm2(ka,kb,k1,k2,k3,false,true,true, ka, aqline);
	+ COM Mpmp3 = qggm3(ka,kb,k1,k2,k3,false,true,true, ka, aqline);
	+
	+ //Colour factors:
	+ COM cm1m1,cm2m2,cm3m3,cm1m2,cm1m3,cm2m3;
	+ cm1m1=8./3.;
	+ cm2m2=8./3.;
	+ cm3m3=6.;
	+ cm1m2 =-1./3.;
	+ cm1m3 = -3.*COM(0.,1.);
	+ cm2m3 = 3.*COM(0.,1.);
	+
	+ //Sqaure and sum for each helicity config:
	+ double Mmmm,Mmmp,Mpmm,Mpmp;
	+ Mmmm = real(cm1m1pow(abs(Mmmm1),2)+cm2m2pow(abs(Mmmm2),2)+cm3m3pow(abs(Mmmm3),2)+2.real(cm1m2Mmmm1conj(Mmmm2))+2.real(cm1m3Mmmm1conj(Mmmm3))+2.real(cm2m3Mmmm2conj(Mmmm3)));
	+ Mmmp = real(cm1m1pow(abs(Mmmp1),2)+cm2m2pow(abs(Mmmp2),2)+cm3m3pow(abs(Mmmp3),2)+2.real(cm1m2Mmmp1conj(Mmmp2))+2.real(cm1m3Mmmp1conj(Mmmp3))+2.real(cm2m3Mmmp2conj(Mmmp3)));
	+ Mpmm = real(cm1m1pow(abs(Mpmm1),2)+cm2m2pow(abs(Mpmm2),2)+cm3m3pow(abs(Mpmm3),2)+2.real(cm1m2Mpmm1conj(Mpmm2))+2.real(cm1m3Mpmm1conj(Mpmm3))+2.real(cm2m3Mpmm2conj(Mpmm3)));
	+ Mpmp = real(cm1m1pow(abs(Mpmp1),2)+cm2m2pow(abs(Mpmp2),2)+cm3m3pow(abs(Mpmp3),2)+2.real(cm1m2Mpmp1conj(Mpmp2))+2.real(cm1m3Mpmp1conj(Mpmp3))+2.real(cm2m3Mpmp2conj(Mpmp3)));
	+
	+ //cout << "Each part: " <<endl;
	+ //cout << cm1m1*pow(abs(Mmmm1),2) <<endl;
	+ //cout<< cm2m2*pow(abs(Mmmm2),2) <<endl;
	+ //cout << cm3m3*pow(abs(Mmmm3),2) <<endl;
	+ //cout << 2.real(cm1m2Mmmm1*conj(Mmmm2))<<endl;
	+ //cout << 2.real(cm1m3Mmmm1*conj(Mmmm3)) <<endl;
	+ //cout << 2.real(cm2m3Mmmm2*conj(Mmmm3)) <<endl;
	+
	+ //cout << cm1m1*pow(abs(Mmmp1),2) <<endl;
	+ //cout<< cm2m2*pow(abs(Mmmp2),2) <<endl;
	+ //cout << cm3m3*pow(abs(Mmmp3),2) <<endl;
	+ //cout << 2.real(cm1m2Mmmp1*conj(Mmmp2))<<endl;
	+ //cout << 2.real(cm1m3Mmmp1*conj(Mmmp3)) <<endl;
	+ //cout << 2.real(cm2m3Mmmp2*conj(Mmmp3)) <<endl;
	+
	+ //cout << cm1m1*pow(abs(Mpmm1),2) <<endl;
	+ //cout<< cm2m2*pow(abs(Mpmm2),2) <<endl;
	+ //cout << cm3m3*pow(abs(Mpmm3),2) <<endl;
	+ //cout << 2.real(cm1m2Mpmm1*conj(Mpmm2))<<endl;
	+ //cout << 2.real(cm1m3Mpmm1*conj(Mpmm3)) <<endl;
	+ //cout << 2.real(cm2m3Mpmm2*conj(Mpmm3)) <<endl;
	+
	+ //cout << cm1m1*pow(abs(Mpmp1),2) <<endl;
	+ //cout<< cm2m2*pow(abs(Mpmp2),2) <<endl;
	+ //cout << cm3m3*pow(abs(Mpmp3),2) <<endl;
	+ //cout << 2.real(cm1m2Mpmp1*conj(Mpmp2))<<endl;
	+ //cout << 2.real(cm1m3Mpmp1*conj(Mpmp3)) <<endl;
	+ //cout << 2.real(cm2m3Mpmp2*conj(Mpmp3)) <<endl;
	+
	+ //cout << "Each hel config seperately: " <<endl;
	+ //cout <<Mmmm << endl;
	+ //cout << Mmmp << endl;
	+ //cout << Mpmm << endl;
	+ //cout << Mpmp << endl;
	+
	+ //Result (averaged, without coupling). Factor of 2 for the helicity configurations we didn't need to evalute explicity
	+
	+ return (2.*(Mmmm+Mmmp+Mpmm+Mpmp)/24./4.)/(ka-k1).m2()/(k2+k3-kb).m2();
	+ //return 2.*(Mmmm)/24./4.;
	+}

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 5:39 AM (5 h, 36 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4982788
Default Alt Text: (23 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions