Page MenuHomeHEPForge
Files F7878947

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
index c55f9d0..f1d7867 100644
--- a/include/HEJ/Tensor.hh
+++ b/include/HEJ/Tensor.hh
@@ -1,210 +1,210 @@
/** \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.
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include <array>
#include <complex>
#include <valarray>
namespace CLHEP {
class HepLorentzVector;
}
class CCurrent;
namespace HEJ {
typedef std::complex<double> COM;
///@TODO put in some namespace
namespace detail {
static constexpr std::size_t d = 4;
//! \internal Compute integer powers at compile time
inline
constexpr std::size_t power(unsigned base, unsigned exp) {
if(exp == 0) return 1;
// use exponentiation by squaring
// there are potentially faster implementations using bit shifts
// but we don't really care because everything's done at compile time
// and this is easier to understand
const std::size_t sqrt = power(base, exp/2);
return (exp % 2)?base*sqrt*sqrt:sqrt*sqrt;
}
}
template <unsigned int N>
class Tensor{
public:
static constexpr std::size_t d = detail::d;
//! Constructor
Tensor() = default;
explicit Tensor(COM x);
//! Rank of Tensor
constexpr unsigned rank() const{
return N;
};
//! total number of entries
std::size_t size() const {
return components.size();
};
//! Tensor element with given indices
template<typename... Indices>
COM const & operator()(Indices... i) const;
//! Tensor element with given indices
template<typename... Indices>
COM & operator()(Indices... rest);
//! implicit conversion to complex number for rank 0 tensors (scalars)
operator COM() const;
Tensor<N> & operator*=(COM const & x);
Tensor<N> & operator/=(COM const & x);
Tensor<N> & operator+=(Tensor<N> const & T2);
Tensor<N> & operator-=(Tensor<N> const & T2);
template<unsigned M>
friend bool operator==(Tensor<M> const & T1, Tensor<M> const & T2);
//! Outer product of two tensors
template<unsigned M>
Tensor<N+M> outer(Tensor<M> const & T2) const;
//! Outer product of two tensors
template<unsigned K, unsigned M>
friend Tensor<K+M> outer(Tensor<K> const & T1, Tensor<M> const & T2);
/**
* \brief T^(mu1...mk..mN)T2_(muk) contract kth index, where k member of [1,N]
* @param T2 Tensor of Rank 1 to contract with.
* @param k int to contract Tensor T2 with from original Tensor.
* @returns T1.contract(T2,1) -> T1^(mu,nu,rho...)T2_mu
*/
Tensor<N-1> contract(Tensor<1> const & T2, int k);
//! Complex conjugate tensor
Tensor<N> & complex_conj();
std::array<COM, detail::power(d, N)> components;
private:
};
template<unsigned N>
Tensor<N> operator*(Tensor<N> t, COM const & x);
template<unsigned N>
Tensor<N> operator/(Tensor<N> t, COM const & x);
template<unsigned N>
Tensor<N> operator+(Tensor<N> T1, Tensor<N> const & T2);
template<unsigned N>
Tensor<N> operator-(Tensor<N> T1, Tensor<N> const & T2);
template<unsigned N>
bool operator==(Tensor<N> const & T1, Tensor<N> const & T2);
template<unsigned N>
bool operator!=(Tensor<N> const & T1, Tensor<N> const & T2);
namespace helicity {
enum Helicity {
minus, plus
};
}
using helicity::Helicity;
inline
constexpr Helicity flip(Helicity h) {
return (h == helicity::plus)?helicity::minus:helicity::plus;
}
/**
* \brief Returns diag(+---) metric
* @returns metric {(1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,-1)}
*/
Tensor<2> metric();
/**
* \brief Calculates current <p1|mu|p2>
* @param p1 Momentum of Particle 1
* @param p2 Momentum of Particle 2
* @param h Helicity of the Spinors
* @returns Tensor T^mu = <p1|mu|p2>
*
* The vector current is helicity conserving, so we only need to state
* one helicity.
*
* @note in/out configuration considered in calculation
*/
- Tensor<1> TCurrent(
+ Tensor<1> current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
);
/**
* \brief Calculates current <p1|mu nu rho|p2>
* @param p1 Momentum of Particle 1
* @param h1 Helicity of Particle 1
* @param p2 Momentum of Particle 2
* @param h2 Helicity of Particle 2
* @returns Tensor T^mu^nu^rho = <p1|mu nu rho|p2>
*
* @note in/out configuration considered in calculation
*/
- Tensor<3> T3Current(CLHEP::HepLorentzVector p1, Helicity h1,
+ Tensor<3> rank3_current(CLHEP::HepLorentzVector p1, Helicity h1,
CLHEP::HepLorentzVector p2, Helicity h2);
/**
* \brief Calculates current <p1|mu nu rho tau sigma|p2>
* @param p1 Momentum of Particle 1
* @param h1 Helicity of Particle 1
* @param p2 Momentum of Particle 2
* @param h2 Helicity of Particle 2
* @returns Tensor T^mu^nu^rho = <p1|mu nu rho tau sigma|p2>
*
* @note in/out configuration considered in calculation
*/
- Tensor<5> T5Current(CLHEP::HepLorentzVector p1, Helicity h1,
+ Tensor<5> rank5_current(CLHEP::HepLorentzVector p1, Helicity h1,
CLHEP::HepLorentzVector p2, Helicity h2);
/**
* \brief Convert from CCurrent class
* @param j Current in CCurrent format
* @returns Current in Tensor Format
*/
Tensor<1> Construct1Tensor(CCurrent j);
/**
* \brief Convert from HLV class
* @param p Current in HLV format
* @returns Current in Tensor Format
*/
Tensor<1> Construct1Tensor(CLHEP::HepLorentzVector p);
/**
* \brief Construct Epsilon (Polarisation) Tensor
* @param k Momentum of incoming/outgoing boson
* @param ref Reference momentum for calculation
* @param pol Polarisation of boson
* @returns Polarisation Tensor E^mu
*/
Tensor<1> eps(CLHEP::HepLorentzVector k, CLHEP::HepLorentzVector ref, Helicity pol);
//! Initialises Tensor values by iterating over permutations of gamma matrices.
bool init_sigma_index();
}
// implementation of template functions
#include "HEJ/detail/Tensor_impl.hh"
diff --git a/src/Tensor.cc b/src/Tensor.cc
index 92251e1..a228668 100644
--- a/src/Tensor.cc
+++ b/src/Tensor.cc
@@ -1,827 +1,827 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/Tensor.hh"
#include <array>
#include <iostream>
#include <valarray>
#include <CLHEP/Vector/LorentzVector.h>
#include "HEJ/currents.hh" // only for CCurrent (?)
namespace HEJ {
namespace{
// Tensor Template definitions
short int sigma_index5[1024];
short int sigma_index3[64];
std::valarray<COM> permfactor5;
std::valarray<COM> permfactor3;
short int helfactor5[2][1024];
short int helfactor3[2][64];
// 2D sigma matrices
const COM sigma0[2][2] = { {1.,0.}, {0., 1.} };
const COM sigma1[2][2] = { {0.,1.}, {1., 0.} };
const COM sigma2[2][2] = { {0,-1.*COM(0,1)}, {1.*COM(0,1), 0.} };
const COM sigma3[2][2] = { {1.,0.}, {0., -1.} };
Tensor<1> Sigma(int i, int j, Helicity hel){
Tensor<1> newT;
if(hel==helicity::plus){
newT.components[0]=sigma0[i][j];
newT.components[1]=sigma1[i][j];
newT.components[2]=sigma2[i][j];
newT.components[3]=sigma3[i][j];
} else {
newT.components[0]= sigma0[i][j];
newT.components[1]=-sigma1[i][j];
newT.components[2]=-sigma2[i][j];
newT.components[3]=-sigma3[i][j];
}
return newT;
}
// map from a list of 5 tensor lorentz indices onto a single index 0<=i<1024
// in 4 dimensional spacetime
int tensor2listindex(std::array<int,5> indexlist){
int mu=indexlist[0];
int nu=indexlist[1];
int sigma=indexlist[2];
int tau=indexlist[3];
int rho=indexlist[4];
int myindex;
myindex = 256*mu+64*nu+16*sigma+4*tau+rho;
if(myindex<0||myindex>1023){
std::cerr<<"bad index in tensor2listindex "<<std::endl;
return 1024;
}
return myindex;
}
// map from a list of 3 tensor lorentz indices onto a single index 0<=i<64 in
// 4 dimensional spacetime
int tensor2listindex(std::array<int,3> indexlist){
int mu=indexlist[0];
int nu=indexlist[1];
int sigma=indexlist[2];
int myindex;
myindex = 16*mu+4*nu+sigma;
if(myindex<0||myindex>64){
std::cerr<<"bad index in tensor2listindex "<<std::endl;
return 64;
}
return myindex;
}
// generate all unique perms of vectors {a,a,a,a,b}, return in perms
// set_permfactor is a bool which encodes the anticommutation relations of the
// sigma matrices namely if we have sigma0, set_permfactor= false because it
// commutes with all others otherwise we need to assign a minus sign to odd
// permutations, set in permfactor
// note, inital perm is always even
void perms41(int same4, int diff, std::vector<std::array<int,5>> & perms){
bool set_permfactor(true);
if(same4==0||diff==0)
set_permfactor=false;
for(int diffpos=0;diffpos<5;diffpos++){
std::array<int,5> tempvec={same4,same4,same4,same4,same4};
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor){
if(diffpos%2==1)
permfactor5[tensor2listindex(tempvec)]=-1.;
}
}
}
// generate all unique perms of vectors {a,a,a,b,b}, return in perms
// note, inital perm is always even
void perms32(int same3, int diff, std::vector<std::array<int,5>> & perms){
bool set_permfactor(true);
if(same3==0||diff==0)
set_permfactor=false;
for(int diffpos1=0;diffpos1<5;diffpos1++){
for(int diffpos2=diffpos1+1;diffpos2<5;diffpos2++){
std::array<int,5> tempvec={same3,same3,same3,same3,same3};
tempvec[diffpos1]=diff;
tempvec[diffpos2]=diff;
perms.push_back(tempvec);
if(set_permfactor){
if((diffpos2-diffpos1)%2==0)
permfactor5[tensor2listindex(tempvec)]=-1.;
}
}
}
}
// generate all unique perms of vectors {a,a,a,b,c}, return in perms
// we have two bools since there are three distinct type of sigma matrices to
// permute, so need to test if the 3xrepeated = sigma0 or if one of the
// singles is sigma0
// if diff1/diff2 are distinct, flipping them results in a change of perm,
// otherwise it'll be symmetric under flip -> encode this in set_permfactor2
// as long as diff2!=0 can ensure inital perm is even
// if diff2=0 then it is not possible to ensure inital perm even -> encode in
// bool signflip
void perms311(int same3, int diff1, int diff2,
std::vector<std::array<int,5>> & perms
){
bool set_permfactor2(true);
bool same0(false);
bool diff0(false);
bool signflip(false); // if true, inital perm is odd
if(same3==0) // is the repeated matrix sigma0?
same0 = true;
else if(diff2==0){ // is one of the single matrices sigma0
diff0=true;
if((diff1%3-same3)!=-1)
signflip=true;
} else if(diff1==0){
std::cerr<<"Note, only first and last argument may be zero"<<std::endl;
return;
}
// possible outcomes: tt, ft, tf
for(int diffpos1=0;diffpos1<5;diffpos1++){
for(int diffpos2=diffpos1+1;diffpos2<5;diffpos2++){
std::array<int,5> tempvec={same3,same3,same3,same3,same3};
tempvec[diffpos1]=diff1;
tempvec[diffpos2]=diff2;
perms.push_back(tempvec);
if(!same0 && !diff0){
// full antisymmetric under exchange of same3,diff1,diff2
if((diffpos2-diffpos1)%2==0){
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd perm
// if this perm is odd, swapping diff1<->diff2 automatically even
set_permfactor2=false;
} else {
permfactor5[tensor2listindex(tempvec)]=COM(0,1); // even perm
// if this perm is even, swapping diff1<->diff2 automatically odd
set_permfactor2=true;
}
} else if(diff0){// one of the single matrices is sigma0
if(signflip){ // initial config is odd!
if(diffpos1%2==1){
permfactor5[tensor2listindex(tempvec)]=COM(0,1); // even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2=false;
} else {
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2=true;
}
} else {// diff0=true, initial config is even
if(diffpos1%2==1){
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd perm
// initally symmetric under diff1<->diff2 =>
// if this perm is odd, automatically odd for first diffpos2
set_permfactor2=true;
} else {
permfactor5[tensor2listindex(tempvec)]=COM(0,1); // even perm
// initally symmetric under diff1<->diff2 =>
// if this perm is even, automatically even for first diffpos2
set_permfactor2=false;
}
}
if((diffpos2-diffpos1-1)%2==1)
set_permfactor2=!set_permfactor2; // change to account for diffpos2
} else if(same0){
// the repeated matrix is sigma0
// => only relative positions of diff1, diff2 matter.
// always even before flip because diff1pos<diff2pos
permfactor5[tensor2listindex(tempvec)]=COM(0,1);
// if this perm is odd, swapping diff1<->diff2 automatically odd
set_permfactor2=true;
}
tempvec[diffpos1]=diff2;
tempvec[diffpos2]=diff1;
perms.push_back(tempvec);
if(set_permfactor2)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor2listindex(tempvec)]=COM(0,1);
}
}
} // end perms311
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
void perms221(int same2a, int same2b, int diff,
std::vector<std::array<int,5>> & perms
){
bool set_permfactor1(true);
bool set_permfactor2(true);
if(same2b==0){
std::cerr<<"second entry in perms221() shouldn't be zero" <<std::endl;
return;
} else if(same2a==0)
set_permfactor1=false;
else if(diff==0)
set_permfactor2=false;
for(int samepos=0;samepos<5;samepos++){
int permcount = 0;
for(int samepos2=samepos+1;samepos2<5;samepos2++){
for(int diffpos=0;diffpos<5;diffpos++){
if(diffpos==samepos||diffpos==samepos2) continue;
std::array<int,5> tempvec={same2a,same2a,same2a,same2a,same2a};
tempvec[samepos]=same2b;
tempvec[samepos2]=same2b;
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor1){
if(set_permfactor2){// full anti-symmetry
if(permcount%2==1)
permfactor5[tensor2listindex(tempvec)]=-1.;
} else { // diff is sigma0
if( ((samepos2-samepos-1)%3>0)
&& (abs(abs(samepos2-diffpos)-abs(samepos-diffpos))%3>0) )
permfactor5[tensor2listindex(tempvec)]=-1.;
}
} else { // same2a is sigma0
if((diffpos>samepos) && (diffpos<samepos2))
permfactor5[tensor2listindex(tempvec)]=-1.;
}
permcount++;
}
}
}
}
// generate all unique perms of vectors {a,a,b,b,c}, return in perms
// there must be a sigma zero if we have 4 different matrices in string
// bool is true if sigma0 is the repeated matrix
void perms2111(int same2, int diff1,int diff2,int diff3,
std::vector<std::array<int,5>> & perms
){
bool twozero(false);
if(same2==0)
twozero=true;
else if (diff1!=0){
std::cerr<<"One of first or second argurments must be a zero"<<std::endl;
return;
} else if(diff2==0|| diff3==0){
std::cerr<<"Only first and second arguments may be a zero."<<std::endl;
return;
}
int permcount = 0;
for(int diffpos1=0;diffpos1<5;diffpos1++){
for(int diffpos2=0;diffpos2<5;diffpos2++){
if(diffpos2==diffpos1) continue;
for(int diffpos3=0;diffpos3<5;diffpos3++){
if(diffpos3==diffpos2||diffpos3==diffpos1) continue;
std::array<int,5> tempvec={same2,same2,same2,same2,same2};
tempvec[diffpos1]=diff1;
tempvec[diffpos2]=diff2;
tempvec[diffpos3]=diff3;
perms.push_back(tempvec);
if(twozero){// don't care about exact positions of singles, just order
if(diffpos2>diffpos3 && diffpos3>diffpos1)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos1>diffpos2 && diffpos2>diffpos3)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos3>diffpos1 && diffpos1>diffpos2)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
else
permfactor5[tensor2listindex(tempvec)]=COM(0,1);// evwn
} else {
if(permcount%2==1)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor2listindex(tempvec)]=COM(0,1);
}
permcount++;
}
}
}
}
void perms21(int same, int diff, std::vector<std::array<int,3>> & perms){
bool set_permfactor(true);
if(same==0||diff==0)
set_permfactor=false;
for(int diffpos=0; diffpos<3;diffpos++){
std::array<int,3> tempvec={same,same,same};
tempvec[diffpos]=diff;
perms.push_back(tempvec);
if(set_permfactor && diffpos==1)
permfactor3[tensor2listindex(tempvec)]=-1.;
}
}
void perms111(int diff1, int diff2, int diff3,
std::vector<std::array<int,3>> & perms
){
bool sig_zero(false);
if(diff1==0)
sig_zero=true;
else if(diff2==0||diff3==0){
std::cerr<<"Only first argument may be a zero."<<std::endl;
return;
}
int permcount=0;
for(int pos1=0;pos1<3;pos1++){
for(int pos2=pos1+1;pos2<3;pos2++){
std::array<int,3> tempvec={diff1,diff1,diff1};
tempvec[pos1]=diff2;
tempvec[pos2]=diff3;
perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
} else if(permcount%2==1){
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
} else {
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
}
tempvec[pos1]=diff3;
tempvec[pos2]=diff2;
perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
} else if(permcount%2==1){
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
} else {
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
}
permcount++;
}
}
}
void SpinorO(CLHEP::HepLorentzVector p, Helicity hel, COM *sp){
// sp is pointer to COM sp[2]
COM pplus = p.e() +p.z();
COM pminus = p.e() -p.z();
COM sqpp= sqrt(pplus);
COM sqpm= sqrt(pminus);
COM pperp = p.x() + COM(0,1)*p.y();
if(hel == helicity::plus){
sp[0] = sqpp;
sp[1] = sqpm*pperp/abs(pperp);
} else {
sp[0] = sqpm*conj(pperp)/abs(pperp);
sp[1] = -sqpp;
}
}
void SpinorIp(COM sqpp, Helicity hel, COM *sp){
// if hel=+
if(hel == helicity::plus){
sp[0] = sqpp;
sp[1] = 0.;
} else {
sp[0] = 0.;
sp[1] = -sqpp;
}
}
void SpinorIm(COM sqpm, Helicity hel, COM *sp){
// if hel=+
if(hel == helicity::plus){
sp[0] = 0;
sp[1] = -sqpm;
} else {
sp[0] = -sqpm;
sp[1] = 0.;
}
}
void Spinor(CLHEP::HepLorentzVector p, Helicity hel, COM *sp){
COM pplus = p.e() +p.z();
COM pminus = p.e() -p.z();
// If incoming along +ve z
if (pminus==0.){
COM sqpp= sqrt(pplus);
SpinorIp(sqpp,hel,sp);
}
// If incoming along -ve z
else if(pplus==0.){
COM sqpm= sqrt(pminus);
SpinorIm(sqpm,hel,sp);
}
// Outgoing
else {
SpinorO(p,hel,sp);
}
}
} // anonymous namespace
Tensor<2> metric(){
static const Tensor<2> g = [](){
Tensor<2> g(0.);
g(0,0) = 1.;
g(1,1) = -1.;
g(2,2) = -1.;
g(3,3) = -1.;
return g;
}();
return g;
}
// <1|mu|2>
- Tensor<1> TCurrent(
+ Tensor<1> current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
){
COM sp1[2];
COM sp2[2];
Tensor<1> newT(0.);
CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
Spinor(p1new, h, sp1);
Spinor(p2new, h, sp2);
for(int i=0;i<2;i++){
for(int j=0; j<2; j++){
newT+=(Sigma(i,j,h)*sp2[j])*conj(sp1[i]);
}
}
return newT;
}
// <1|mu|2>
- Tensor<1> TCurrent(CLHEP::HepLorentzVector p1, Helicity h1,
+ Tensor<1> current(CLHEP::HepLorentzVector p1, Helicity h1,
CLHEP::HepLorentzVector p2, Helicity h2
){
COM sp1[2];
COM sp2[2];
Tensor<1> newT(0.);
CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
if(h1!=h2){
return newT;
}
Spinor(p1new, h1, sp1);
Spinor(p2new, h2, sp2);
for(int i=0;i<2;i++){
for(int j=0; j<2; j++){
newT+=(Sigma(i,j,h2)*sp2[j])*conj(sp1[i]);
}
}
return newT;
}
// <1|mu nu sigma|2>
- Tensor<3> T3Current(CLHEP::HepLorentzVector p1, Helicity h1,
+ Tensor<3> rank3_current(CLHEP::HepLorentzVector p1, Helicity h1,
CLHEP::HepLorentzVector p2, Helicity h2
){
COM sp1[2];
COM sp2[2];
Tensor<3> newT(0.);
CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
if(h1!=h2){
return newT;
}
Spinor(p1new, h1, sp1);
Spinor(p2new, h2, sp2);
COM current[4];
for(int i=0; i<2;i++){
for(int j=0; j<2;j++){
current[0]+=conj(sp1[i])*sigma0[i][j]*sp2[j];
current[1]+=conj(sp1[i])*sigma1[i][j]*sp2[j];
current[2]+=conj(sp1[i])*sigma2[i][j]*sp2[j];
current[3]+=conj(sp1[i])*sigma3[i][j]*sp2[j];
}
}
for( std::size_t itensor=0; itensor<newT.size(); itensor++ ){
double hfact = double( helfactor3[h2][itensor] );
newT.components[itensor] = current[sigma_index3[itensor]] * hfact
* permfactor3[itensor];
}
return newT;
}
// <1|mu nu sigma tau rho|2>
- Tensor<5> T5Current(CLHEP::HepLorentzVector p1, Helicity h1,
+ Tensor<5> rank5_current(CLHEP::HepLorentzVector p1, Helicity h1,
CLHEP::HepLorentzVector p2, Helicity h2
){
COM sp1[2];
COM sp2[2];
Tensor<5> newT(0.);
CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
if(h1!=h2){
return newT;
}
Spinor(p1new, h1, sp1);
Spinor(p2new, h2, sp2);
COM current[4];
for(int i=0; i<2;i++){
for(int j=0; j<2;j++){
current[0]+=conj(sp1[i])*sigma0[i][j]*sp2[j];
current[1]+=conj(sp1[i])*sigma1[i][j]*sp2[j];
current[2]+=conj(sp1[i])*sigma2[i][j]*sp2[j];
current[3]+=conj(sp1[i])*sigma3[i][j]*sp2[j];
}
}
for(std::size_t itensor=0;itensor<newT.size();itensor++){
double hfact = double(helfactor5[h2][itensor]);
newT.components[itensor] = current[sigma_index5[itensor]] * hfact
* permfactor5[itensor];
}
return newT;
}
Tensor<1> Construct1Tensor(CCurrent j){
Tensor<1> newT;
newT.components={j.c0,j.c1,j.c2,j.c3};
return newT;
}
Tensor<1> Construct1Tensor(CLHEP::HepLorentzVector p){
Tensor<1> newT;
newT.components={p.e(),p.x(),p.y(),p.z()};
return newT;
}
Tensor<1> eps(CLHEP::HepLorentzVector k, CLHEP::HepLorentzVector ref, Helicity pol){
Tensor<1> polvec;
COM spk[2];
COM spr[2];
COM denom;
CLHEP::HepLorentzVector knew{ k.e()<0?-k:k };
Spinor(knew, pol, spk);
Spinor(ref, flip(pol), spr);
denom = pow(-1.,pol)*sqrt(2)*(conj(spr[0])*spk[0] + conj(spr[1])*spk[1]);
- polvec = TCurrent(ref,knew,flip(pol))/denom;
+ polvec = current(ref,knew,flip(pol))/denom;
return polvec;
}
// slow function! - but only need to evaluate once.
bool init_sigma_index(){
// initialize permfactor(3) to list of ones (change to minus one for each odd
// permutation and multiply by i for all permutations in perms2111, perms311,
// perms111)
permfactor5.resize(1024,1.);
permfactor3.resize(64,1.);
// first set sigma_index (5)
std::vector<std::array<int,5>> sigma0indices;
std::vector<std::array<int,5>> sigma1indices;
std::vector<std::array<int,5>> sigma2indices;
std::vector<std::array<int,5>> sigma3indices;
// need to generate all possible permuations of {i,j,k,l,m}
// where each index can be {0,1,2,3,4}
// 1024 possibilities
// perms with 5 same
sigma0indices.push_back({0,0,0,0,0});
sigma1indices.push_back({1,1,1,1,1});
sigma2indices.push_back({2,2,2,2,2});
sigma3indices.push_back({3,3,3,3,3});
// perms with 4 same
perms41(1,0,sigma0indices);
perms41(2,0,sigma0indices);
perms41(3,0,sigma0indices);
perms41(0,1,sigma1indices);
perms41(2,1,sigma1indices);
perms41(3,1,sigma1indices);
perms41(0,2,sigma2indices);
perms41(1,2,sigma2indices);
perms41(3,2,sigma2indices);
perms41(0,3,sigma3indices);
perms41(1,3,sigma3indices);
perms41(2,3,sigma3indices);
// perms with 3 same, 2 same
perms32(0,1,sigma0indices);
perms32(0,2,sigma0indices);
perms32(0,3,sigma0indices);
perms32(1,0,sigma1indices);
perms32(1,2,sigma1indices);
perms32(1,3,sigma1indices);
perms32(2,0,sigma2indices);
perms32(2,1,sigma2indices);
perms32(2,3,sigma2indices);
perms32(3,0,sigma3indices);
perms32(3,1,sigma3indices);
perms32(3,2,sigma3indices);
// perms with 3 same, 2 different
perms311(1,2,3,sigma0indices);
perms311(2,3,1,sigma0indices);
perms311(3,1,2,sigma0indices);
perms311(0,2,3,sigma1indices);
perms311(2,3,0,sigma1indices);
perms311(3,2,0,sigma1indices);
perms311(0,3,1,sigma2indices);
perms311(1,3,0,sigma2indices);
perms311(3,1,0,sigma2indices);
perms311(0,1,2,sigma3indices);
perms311(1,2,0,sigma3indices);
perms311(2,1,0,sigma3indices);
perms221(1,2,0,sigma0indices);
perms221(1,3,0,sigma0indices);
perms221(2,3,0,sigma0indices);
perms221(0,2,1,sigma1indices);
perms221(0,3,1,sigma1indices);
perms221(2,3,1,sigma1indices);
perms221(0,1,2,sigma2indices);
perms221(0,3,2,sigma2indices);
perms221(1,3,2,sigma2indices);
perms221(0,1,3,sigma3indices);
perms221(0,2,3,sigma3indices);
perms221(1,2,3,sigma3indices);
perms2111(0,1,2,3,sigma0indices);
perms2111(1,0,2,3,sigma1indices);
perms2111(2,0,3,1,sigma2indices);
perms2111(3,0,1,2,sigma3indices);
if(sigma0indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma0indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma1indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma1indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma2indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma2indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
} else if(sigma3indices.size()!=256){
std::cerr<<"sigma_index not set: ";
std::cerr<<"sigma3indices has "<< sigma0indices.size() << " components" << std::endl;
return false;
}
for(int i=0;i<256;i++){
// map each unique set of tensor indices to its position in a list
int index0 = tensor2listindex(sigma0indices.at(i));
int index1 = tensor2listindex(sigma1indices.at(i));
int index2 = tensor2listindex(sigma2indices.at(i));
int index3 = tensor2listindex(sigma3indices.at(i));
sigma_index5[index0]=0;
sigma_index5[index1]=1;
sigma_index5[index2]=2;
sigma_index5[index3]=3;
short int sign[4]={1,-1,-1,-1};
// plus->true->1
helfactor5[1][index0] = sign[sigma0indices.at(i)[1]]
* sign[sigma0indices.at(i)[3]];
helfactor5[1][index1] = sign[sigma1indices.at(i)[1]]
* sign[sigma1indices.at(i)[3]];
helfactor5[1][index2] = sign[sigma2indices.at(i)[1]]
* sign[sigma2indices.at(i)[3]];
helfactor5[1][index3] = sign[sigma3indices.at(i)[1]]
* sign[sigma3indices.at(i)[3]];
// minus->false->0
helfactor5[0][index0] = sign[sigma0indices.at(i)[0]]
* sign[sigma0indices.at(i)[2]]
* sign[sigma0indices.at(i)[4]];
helfactor5[0][index1] = sign[sigma1indices.at(i)[0]]
* sign[sigma1indices.at(i)[2]]
* sign[sigma1indices.at(i)[4]];
helfactor5[0][index2] = sign[sigma2indices.at(i)[0]]
* sign[sigma2indices.at(i)[2]]
* sign[sigma2indices.at(i)[4]];
helfactor5[0][index3] = sign[sigma3indices.at(i)[0]]
* sign[sigma3indices.at(i)[2]]
* sign[sigma3indices.at(i)[4]];
}
// now set sigma_index3
std::vector<std::array<int,3>> sigma0indices3;
std::vector<std::array<int,3>> sigma1indices3;
std::vector<std::array<int,3>> sigma2indices3;
std::vector<std::array<int,3>> sigma3indices3;
// perms with 3 same
sigma0indices3.push_back({0,0,0});
sigma1indices3.push_back({1,1,1});
sigma2indices3.push_back({2,2,2});
sigma3indices3.push_back({3,3,3});
// 2 same
perms21(1,0,sigma0indices3);
perms21(2,0,sigma0indices3);
perms21(3,0,sigma0indices3);
perms21(0,1,sigma1indices3);
perms21(2,1,sigma1indices3);
perms21(3,1,sigma1indices3);
perms21(0,2,sigma2indices3);
perms21(1,2,sigma2indices3);
perms21(3,2,sigma2indices3);
perms21(0,3,sigma3indices3);
perms21(1,3,sigma3indices3);
perms21(2,3,sigma3indices3);
// none same
perms111(1,2,3,sigma0indices3);
perms111(0,2,3,sigma1indices3);
perms111(0,3,1,sigma2indices3);
perms111(0,1,2,sigma3indices3);
if(sigma0indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma0indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma1indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma1indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma2indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma2indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
} else if(sigma3indices3.size()!=16){
std::cerr<<"sigma_index3 not set: ";
std::cerr<<"sigma3indices3 has "<< sigma0indices3.size() << " components" << std::endl;
return false;
}
for(int i=0;i<16;i++){
int index0 = tensor2listindex(sigma0indices3.at(i));
int index1 = tensor2listindex(sigma1indices3.at(i));
int index2 = tensor2listindex(sigma2indices3.at(i));
int index3 = tensor2listindex(sigma3indices3.at(i));
sigma_index3[index0]=0;
sigma_index3[index1]=1;
sigma_index3[index2]=2;
sigma_index3[index3]=3;
short int sign[4]={1,-1,-1,-1};
// plus->true->1
helfactor3[1][index0] = sign[sigma0indices3.at(i)[1]];
helfactor3[1][index1] = sign[sigma1indices3.at(i)[1]];
helfactor3[1][index2] = sign[sigma2indices3.at(i)[1]];
helfactor3[1][index3] = sign[sigma3indices3.at(i)[1]];
// minus->false->0
helfactor3[0][index0] = sign[sigma0indices3.at(i)[0]]
* sign[sigma0indices3.at(i)[2]];
helfactor3[0][index1] = sign[sigma1indices3.at(i)[0]]
* sign[sigma1indices3.at(i)[2]];
helfactor3[0][index2] = sign[sigma2indices3.at(i)[0]]
* sign[sigma2indices3.at(i)[2]];
helfactor3[0][index3] = sign[sigma3indices3.at(i)[0]]
* sign[sigma3indices3.at(i)[2]];
}
return true;
} // end init_sigma_index
}
diff --git a/src/Wjets.cc b/src/Wjets.cc
index 55f6ebc..9d84c8f 100644
--- a/src/Wjets.cc
+++ b/src/Wjets.cc
@@ -1,1315 +1,1314 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include "HEJ/currents.hh"
#include "HEJ/utility.hh"
#include "HEJ/Tensor.hh"
#include "HEJ/Constants.hh"
#include <array>
#include <iostream>
using HEJ::Tensor;
using HEJ::init_sigma_index;
using HEJ::metric;
-using HEJ::TCurrent;
-using HEJ::T3Current;
-using HEJ::T5Current;
+using HEJ::rank3_current;
+using HEJ::rank5_current;
using HEJ::eps;
using HEJ::Construct1Tensor;
using HEJ::Helicity;
namespace helicity = HEJ::helicity;
namespace { // Helper Functions
// FKL W Helper Functions
double WProp (const HLV & plbar, const HLV & 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 (
HLV pout, Helicity helout,
HLV plbar, Helicity hellbar,
HLV pl, Helicity hell,
HLV pin, Helicity 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 && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
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,pl,helout);
COM prod1665=temp2.dot(t65);
temp3 = joi(plbar,helin,pin,helin);
COM prod5465=temp3.dot(t65);
temp2=joo(pout,helout,plbar,helout);
temp3=joi(pl,hell,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 (
HLV pout, Helicity helout,
HLV plbar, Helicity hellbar,
HLV pl, Helicity hell,
HLV pin, Helicity 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 && hellbar==hell) {
HLV qa=pout+plbar+pl;
HLV qb=pin-plbar-pl;
double ta(qa.m2()),tb(qb.m2());
CCurrent temp2,temp3,temp5;
CCurrent t65 = joo(pl,hell,plbar,hellbar);
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(plbar,helout,pout,helout); // temp2 is <5|alpha|1>
COM prod5165=temp2.dot(t65);
temp3 = jio(pin,helin,pl,helin); // temp3 is <4|alpha|6>
COM prod4665=temp3.dot(t65);
temp2=joo(pl,helout,pout,helout); // temp2 is now <6|mu|1>
temp3=jio(pin,helin,plbar,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> jW(HLV pin, HLV pout, HLV plbar, HLV pl, bool aqline){
// Build the external quark line W Emmision
- Tensor<1> ABCurr = TCurrent(pl, plbar, helicity::minus);
+ Tensor<1> ABCurr = HEJ::current(pl, plbar, helicity::minus);
Tensor<1> Tp4W = Construct1Tensor((pout+pl+plbar));//p4+pw
Tensor<1> TpbW = Construct1Tensor((pin-pl-plbar));//pb-pw
Tensor<3> J4bBlank;
if (aqline){
- J4bBlank = T3Current(pin,helicity::minus,pout,helicity::minus);
+ J4bBlank = rank3_current(pin,helicity::minus,pout,helicity::minus);
}
else{
- J4bBlank = T3Current(pout,helicity::minus,pin,helicity::minus);
+ J4bBlank = rank3_current(pout,helicity::minus,pin,helicity::minus);
}
double t4AB = (pout+pl+plbar).m2();
double tbAB = (pin-pl-plbar).m2();
Tensor<2> J4b1 = (J4bBlank.contract(Tp4W,2))/t4AB;
Tensor<2> J4b2 = (J4bBlank.contract(TpbW,2))/tbAB;
Tensor<2> T4bmMom(0.);
if (aqline){
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(mu, nu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,-1);
}
}
}
else{
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
T4bmMom(nu,mu) = (J4b1(nu,mu) + J4b2(mu,nu))*COM(0,1);
}
}
}
Tensor<1> T4bm = T4bmMom.contract(ABCurr,1);
return T4bm;
}
// Relevant W+Jets Unordered Contribution Helper Functions
double jM2Wuno(HLV pg, HLV p1,HLV plbar, HLV pl, HLV pa, Helicity h1,
HLV p2, HLV pb, Helicity h2, Helicity pol
){
//@TODO Simplify the below (less Tensor class?)
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 pW = pl+plbar;
HLV q1g=pa-pW-p1-pg;
HLV q1 = pa-p1-pW;
HLV 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> epsg = eps(pg,p2,pol);
- Tensor<1> epsW = TCurrent(pl,plbar,helicity::minus);
- Tensor<1> j2b = TCurrent(p2,pb,h2);
+ Tensor<1> epsW = HEJ::current(pl,plbar,helicity::minus);
+ Tensor<1> j2b = HEJ::current(p2,pb,h2);
Tensor<1> Tq1q2 = Construct1Tensor((q1+q2)/taW1 + (pb/pb.dot(pg)
+p2/p2.dot(pg)) * tb2/(2*tb2g));
Tensor<1> Tq1g = Construct1Tensor((-pg-q1))/taW1;
Tensor<1> Tq2g = Construct1Tensor((pg-q2));
Tensor<1> TqaW = Construct1Tensor((pa-pW));//pa-pw
Tensor<1> Tqag = Construct1Tensor((pa-pg));
Tensor<1> TqaWg = Construct1Tensor((pa-pg-pW));
Tensor<1> Tp1g = Construct1Tensor((p1+pg));
Tensor<1> Tp1W = Construct1Tensor((p1+pW));//p1+pw
Tensor<1> Tp1gW = Construct1Tensor((p1+pg+pW));//p1+pw+pg
Tensor<2> g=metric();
- Tensor<3> J31a = T3Current(p1, h1, pa, h1);
+ Tensor<3> J31a = rank3_current(p1, h1, pa, h1);
Tensor<2> J2_qaW =J31a.contract(TqaW/taW, 2);
Tensor<2> J2_p1W =J31a.contract(Tp1W/s1W, 2);
Tensor<3> L1a = outer(Tq1q2, J2_qaW);
Tensor<3> L1b = outer(Tq1q2, J2_p1W);
Tensor<3> L2a = outer(Tq1g,J2_qaW);
Tensor<3> L2b = outer(Tq1g, J2_p1W);
Tensor<3> L3 = outer(g, J2_qaW.contract(Tq2g,1)+J2_p1W.contract(Tq2g,2))/taW1;
Tensor<3> L(0.);
- Tensor<5> J51a = T5Current(p1, h1, pa, h1);
+ Tensor<5> J51a = rank5_current(p1, h1, pa, h1);
Tensor<4> J_qaW = J51a.contract(TqaW,4);
Tensor<4> J_qag = J51a.contract(Tqag,4);
Tensor<4> J_p1gW = J51a.contract(Tp1gW,4);
Tensor<3> U1a = J_qaW.contract(Tp1g,2);
Tensor<3> U1b = J_p1gW.contract(Tp1g,2);
Tensor<3> U1c = J_p1gW.contract(Tp1W,2);
Tensor<3> U1(0.);
Tensor<3> U2a = J_qaW.contract(TqaWg,2);
Tensor<3> U2b = J_qag.contract(TqaWg,2);
Tensor<3> U2c = J_qag.contract(Tp1W,2);
Tensor<3> U2(0.);
for(int nu=0; nu<4;nu++){
for(int mu=0;mu<4;mu++){
for(int rho=0;rho<4;rho++){
L(nu, mu, rho) = L1a(nu,mu,rho) + L1b(nu,rho,mu)
+ L2a(mu,nu,rho) + L2b(mu,rho,nu) + L3(mu,nu,rho);
U1(nu, mu, rho) = U1a(nu, mu, rho) / (s1g*taW)
+ U1b(nu,rho,mu) / (s1g*s1gW) + U1c(rho,nu,mu) / (s1W*s1gW);
U2(nu,mu,rho) = U2a(mu,nu,rho) / (taWg*taW)
+ U2b(mu,rho,nu) / (taWg*tag) + U2c(rho,mu,nu) / (s1W*tag);
}
}
}
COM X = ((((U1-L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
COM Y = ((((U2+L).contract(epsW,3)).contract(j2b,2)).contract(epsg,1));
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> gtqqxW(HLV pq,HLV pqbar,HLV pl,HLV plbar){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
// Define llx current.
- Tensor<1> ABCur = TCurrent(pl, plbar, helicity::minus);
+ Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//blank 3 Gamma Current
- Tensor<3> JV23 = T3Current(pq,helicity::minus,pqbar,helicity::minus);
+ Tensor<3> JV23 = rank3_current(pq,helicity::minus,pqbar,helicity::minus);
// Components of g->qqW before W Contraction
Tensor<2> JV1 = JV23.contract((Tpq + TAB),2)/(s2AB);
Tensor<2> 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> JVCur = (JV1.contract(ABCur,1) - JV2.contract(ABCur,2))*COM(0.,-1.);
return JVCur;
}
// Helper Functions Calculate the Crossed Contribution
Tensor <2> MCrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar, HLV pl,
HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCross?
// Useful propagator factors
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV 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> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
// Define llx current.
- Tensor<1> ABCur = TCurrent(pl, plbar,helicity::minus);
+ Tensor<1> ABCur = HEJ::current(pl, plbar,helicity::minus);
//Blank 5 gamma Current
- Tensor<5> J523 = T5Current(pq,helicity::minus,pqbar,helicity::minus);
+ Tensor<5> J523 = rank5_current(pq,helicity::minus,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_q3q = J523.contract((Tq3+Tpq),2);
Tensor<4> J_2AB = J523.contract((Tpq+TAB),2);
// Components of Crossed Vertex Contribution
Tensor<3> Xcro1 = J_q3q.contract((Tpqbar + TAB),3);
Tensor<3> Xcro2 = J_q3q.contract((Tq1-Tpqbar),3);
Tensor<3> Xcro3 = J_2AB.contract((Tq1-Tpqbar),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xcro1Cont = Xcro1.contract(ABCur,3)/(tcro1*s3AB);
Tensor<2> Xcro2Cont = Xcro2.contract(ABCur,2)/(tcro1*tcro2);
Tensor<2> Xcro3Cont = Xcro3.contract(ABCur,1)/(s2AB*tcro2);
//Initialise the Crossed Vertex Object
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = -(-Xcro1Cont(nu,mu)+Xcro2Cont(nu,mu)+Xcro3Cont(nu,mu));
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncrossW(HLV pa, HLV, HLV, HLV, HLV pq, HLV pqbar,
HLV pl, HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncross?
double s2AB=(pl+plbar+pq).m2();
double s3AB=(pl+plbar+pqbar).m2();
HLV 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> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
// Define llx current.
- Tensor<1> ABCur = TCurrent(pl, plbar, helicity::minus);
+ Tensor<1> ABCur = HEJ::current(pl, plbar, helicity::minus);
//Blank 5 gamma Current
- Tensor<5> J523 = T5Current(pq,helicity::minus,pqbar,helicity::minus);
+ Tensor<5> J523 = rank5_current(pq,helicity::minus,pqbar,helicity::minus);
// 4 gamma currents (with 1 contraction already).
Tensor<4> J_2AB = J523.contract((Tpq+TAB),2);
Tensor<4> J_q1q = J523.contract((Tq1-Tpq),2);
// 2 Contractions taken care of.
Tensor<3> Xunc1 = J_2AB.contract((Tq3+Tpqbar),3);
Tensor<3> Xunc2 = J_q1q.contract((Tq3+Tpqbar),3);
Tensor<3> Xunc3 = J_q1q.contract((Tpqbar+TAB),3);
// Term Denominators Taken Care of at this stage
Tensor<2> Xunc1Cont = Xunc1.contract(ABCur,1)/(s2AB*tunc2);
Tensor<2> Xunc2Cont = Xunc2.contract(ABCur,2)/(tunc1*tunc2);
Tensor<2> Xunc3Cont = Xunc3.contract(ABCur,3)/(tunc1*s3AB);
//Initialise the Uncrossed Vertex Object
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -(- Xunc1Cont(mu,nu)+Xunc2Cont(mu,nu) +Xunc3Cont(mu,nu));
}
}
return Xunc;
}
// Helper Functions Calculate the g->qqxW (Eikonal) Contributions
Tensor <2> MSymW(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
HLV pl,HLV plbar, std::vector<HLV> partons, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MSym?
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();
HLV 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> Tp1 = Construct1Tensor(p1);
Tensor<1> Tp4 = Construct1Tensor(p4);
Tensor<1> Tpa = Construct1Tensor(pa);
Tensor<1> Tpb = Construct1Tensor(pb);
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> TAB = Construct1Tensor(pl+plbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
Tensor<2> g=metric();
// g->qqW Current (Factors of sqrt2 dealt with in this function.)
Tensor<1> JV = gtqqxW(pq,pqbar,pl,plbar);
// 1a gluon emisson Contribution
Tensor<3> X1a = outer(g, Tp1*(t1/(s12+s13+s1A+s1B))
+ Tpa*(t1/(sa2+sa3+saA+saB)) );
Tensor<2> X1aCont = X1a.contract(JV,3);
//4b gluon emission Contribution
Tensor<3> X4b = outer(g, Tp4*(t3/(s42+s43+s4A+s4B))
+ Tpb*(t3/(sb2+sb3+sbA+sbB)) );
Tensor<2> X4bCont = X4b.contract(JV,3);
//Set up each term of 3G diagram.
Tensor<3> X3g1 = outer(Tq1+Tpq+Tpqbar+TAB, g);
Tensor<3> X3g2 = outer(Tq3-Tpq-Tpqbar-TAB, g);
Tensor<3> X3g3 = outer(Tq1+Tq3, g);
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(JV,3);
Tensor<2> X3g2Cont = X3g2.contract(JV,2);
Tensor<2> 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>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu,nu) = (X3g1Cont(nu,mu) + X3g2Cont(mu,nu) - X3g3Cont(nu,mu))
+ (X1aCont(mu,nu) - X4bCont(mu,nu));
}
}
return Xsym/s23AB;
}
Tensor <2> MCross(HLV pa, HLV pq, HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MCrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2=(q1-pqbar).m2();
Tensor<1> Tq1 = Construct1Tensor(q1-pqbar);
//Blank 3 gamma Current
- Tensor<3> J323 = T3Current(pq,hq,pqbar,hq);
+ Tensor<3> J323 = rank3_current(pq,hq,pqbar,hq);
// 2 gamma current (with 1 contraction already).
Tensor<2> XCroCont = J323.contract((Tq1),2)/(t2);
//Initialise the Crossed Vertex
Tensor<2> Xcro(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xcro(mu,nu) = XCroCont(nu,mu);
}
}
return Xcro;
}
// Helper Functions Calculate the Uncrossed Contribution
Tensor <2> MUncross(HLV pa, HLV pq,HLV pqbar, std::vector<HLV> partons,
Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MUncrossW?
HLV q1;
q1=pa;
for(int i=0;i<nabove+1;i++){
q1-=partons.at(i);
}
double t2 = (q1-pq).m2();
Tensor<1> Tq1 = Construct1Tensor(q1-pq);
//Blank 3 gamma Current
- Tensor<3> J323 = T3Current(pq,hq,pqbar,hq);
+ Tensor<3> J323 = rank3_current(pq,hq,pqbar,hq);
// 2 gamma currents (with 1 contraction already).
Tensor<2> XUncCont = J323.contract((Tq1),2)/t2;
//Initialise the Uncrossed Vertex
Tensor<2> Xunc(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xunc(mu,nu) = -XUncCont(mu,nu);
}
}
return Xunc;
}
// Helper Functions Calculate the Eikonal Contributions
Tensor <2> MSym(HLV pa, HLV p1, HLV pb, HLV p4, HLV pq, HLV pqbar,
std::vector<HLV> partons, Helicity hq, int nabove
){
//@TODO Simplify the calculation below Maybe combine with MsymW?
HLV 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> Tp1 = Construct1Tensor(p1);
Tensor<1> Tp4 = Construct1Tensor(p4);
Tensor<1> Tpa = Construct1Tensor(pa);
Tensor<1> Tpb = Construct1Tensor(pb);
Tensor<1> Tpq = Construct1Tensor(pq);
Tensor<1> Tpqbar = Construct1Tensor(pqbar);
Tensor<1> Tq1 = Construct1Tensor(q1);
Tensor<1> Tq3 = Construct1Tensor(q3);
Tensor<2> g=metric();
- Tensor<1> qqxCur = TCurrent(pq, pqbar, hq);
+ Tensor<1> qqxCur = HEJ::current(pq, pqbar, hq);
// // 1a gluon emisson Contribution
Tensor<3> X1a = outer(g, Tp1*(t1/(s12+s13))+Tpa*(t1/(sa2+sa3)));
Tensor<2> X1aCont = X1a.contract(qqxCur,3);
// //4b gluon emission Contribution
Tensor<3> X4b = outer(g, Tp4*(t3/(s42+s43)) + Tpb*(t3/(sb2+sb3)));
Tensor<2> X4bCont = X4b.contract(qqxCur,3);
// New Formulation Corresponding to New Analytics
Tensor<3> X3g1 = outer(Tq1+Tpq+Tpqbar, g);
Tensor<3> X3g2 = outer(Tq3-Tpq-Tpqbar, g);
Tensor<3> X3g3 = outer(Tq1+Tq3, g);
// Note the contraction of indices changes term by term
Tensor<2> X3g1Cont = X3g1.contract(qqxCur,3);
Tensor<2> X3g2Cont = X3g2.contract(qqxCur,2);
Tensor<2> X3g3Cont = X3g3.contract(qqxCur,1);
Tensor<2>Xsym(0.);
for(int mu=0; mu<4;mu++){
for(int nu=0;nu<4;nu++){
Xsym(mu, nu) = COM(0,1) * ( (X3g1Cont(nu,mu) + X3g2Cont(mu,nu)
- X3g3Cont(nu,mu)) + (X1aCont(mu,nu) - X4bCont(mu,nu)) );
}
}
return Xsym/s23;
}
} // Anonymous Namespace helper functions
//! W+Jets FKL Contributions
/**
* @brief W+Jets FKL Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_\mu.
* Handles all possible incoming states.
*/
double jW_j( HLV p1out, HLV plbar, HLV pl, HLV p1in, HLV p2out, HLV p2in,
bool aqlineb, bool aqlinef
){
CCurrent mj1m,mj2p,mj2m;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out);
if(aqlineb) mj1m=jWbar(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
else mj1m=jW(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
if(aqlinef){
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
} else{
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
}
COM Mmp=mj1m.dot(mj2p);
COM Mmm=mj1m.dot(mj2m);
// sum of spinor strings ||^2
double a2Mmp=abs2(Mmp);
double a2Mmm=abs2(Mmm);
double WPropfact = WProp(plbar, pl);
// 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 ME_W_qQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false);
}
double ME_W_qQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, true);
}
double ME_W_qbarQ (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false);
}
double ME_W_qbarQbar (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, true);
}
double ME_W_qg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, false, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_W_qbarg (HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out, HLV p2in){
return jW_j(p1out, plbar, pl, p1in, p2out, p2in, true, false)*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param p1out Outgoing Particle 1. (W emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (W emission)
* @param p2out Outgoing Particle 2 (Quark, unordered emission this side.)
* @param p2in Incoming Particle 2 (Quark, unordered emission this side.)
* @param pg Unordered Gluon momenta
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W ^\mu j_{uno}_\mu. Ie, unordered with W emission opposite side.
* Handles all possible incoming states.
*/
double jW_juno(HLV p1out, HLV plbar, HLV pl,HLV p1in, HLV p2out,
HLV p2in, HLV pg, bool aqlineb, bool aqlinef){
CCurrent mj1m,mj2p,mj2m, jgbm,jgbp,j2gm,j2gp;
HLV q1=p1in-p1out-plbar-pl;
HLV q2=-(p2in-p2out-pg);
HLV q3=-(p2in-p2out);
if(aqlineb) mj1m=jWbar(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
else mj1m=jW(p1out,helicity::minus,plbar,helicity::minus,pl,helicity::minus,p1in,helicity::minus);
//@TODO Is aqlinef necessary? Gives same results.
if(aqlinef){
mj2p=jio(p2in,true,p2out,true);
mj2m=jio(p2in,false,p2out,false);
j2gp=joo(pg,true,p2out,true);
j2gm=joo(pg,false,p2out,false);
jgbp=jio(p2in,true,pg,true);
jgbm=jio(p2in,false,pg,false);
} else{
mj2p=joi(p2out,true,p2in,true);
mj2m=joi(p2out,false,p2in,false);
j2gp=joo(p2out,true,pg,true);
j2gm=joo(p2out,false,pg,false);
jgbp=joi(pg,true,p2in,true);
jgbm=joi(pg,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 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
ampsq*=WProp(plbar, pl);
return ampsq/((16)*(q2.m2()*q1.m2()));
}
double ME_W_unob_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, false);
}
double ME_W_unob_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, false, true);
}
double ME_W_unob_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, false);
}
double ME_W_unob_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p2out, plbar, pl, p2in, p1out, p1in, pg, true, true);
}
double ME_W_unof_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, false);
}
double ME_W_unof_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, false);
}
double ME_W_unof_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, false, true);
}
double ME_W_unof_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jW_juno(p1out, plbar, pl, p1in, p2out, p2in, pg, true, true);
}
/**
* @brief W+Jets Unordered Contributions, function to handle all incoming types.
* @param pg Unordered Gluon momenta
* @param p1out Outgoing Particle 1. (Quark - W and Uno emission)
* @param plbar Outgoing election momenta
* @param pl Outgoing neutrino momenta
* @param p1in Incoming Particle 1. (Quark - W and Uno emission)
* @param p2out Outgoing Particle 2
* @param p2in Incoming Particle 2
* @param aqlineb Bool. Is Backwards quark line an anti-quark line?
*
* Calculates j_W_{uno} ^\mu j_\mu. Ie, unordered with W emission same side.
* Handles all possible incoming states. Note this handles both forward and back-
* -ward Wuno emission. For forward, ensure p1out is the uno and W emission parton.
* @TODO: Include separate wrapper functions for forward and backward to clean up
* ME_W_unof_current in `MatrixElement.cc`.
*/
double jWuno_j(HLV pg, HLV p1out, HLV plbar, HLV pl, HLV p1in,
HLV p2out, HLV p2in, bool aqlineb
){
//Calculate different Helicity choices
const Helicity h = aqlineb?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::plus);
double ME2mpm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::plus,helicity::minus);
double ME2mmp = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::plus);
double ME2mmm = jM2Wuno(pg, p1out,plbar,pl,p1in,h,p2out,p2in,helicity::minus,helicity::minus);
return ME2mpp + ME2mpm + ME2mmp + ME2mmm;
}
double ME_Wuno_qQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false);
}
double ME_Wuno_qbarQ(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qbarQbar(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true);
}
double ME_Wuno_qg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, false)*K_g(p2out, p2in)/HEJ::C_F;
}
double ME_Wuno_qbarg(HLV p1out, HLV p1in, HLV p2out, HLV p2in,
HLV pg, HLV plbar, HLV pl
){
return jWuno_j(pg, p1out, plbar, pl, p1in, p2out, p2in, true)*K_g(p2out, p2in)/HEJ::C_F;
}
/**
* @brief W+Jets Extremal qqx Contributions, function to handle all incoming types.
* @param pgin Incoming gluon which will split into qqx.
* @param pqout Quark of extremal qqx outgoing (W-Emission).
* @param plbar Outgoing anti-lepton momenta
* @param pl Outgoing lepton momenta
* @param pqbarout Anti-quark of extremal qqx pair. (W-Emission)
* @param pout Outgoing Particle 2 (end of FKL chain)
* @param p2in Incoming Particle 2
* @param aqlinef Bool. Is Forwards quark line an anti-quark line?
*
* Calculates j_W_{qqx} ^\mu j_\mu. Ie, Ex-QQX with W emission same side.
* Handles all possible incoming states. Calculated via crossing symmetry from jWuno_j.
*/
double jWqqx_j(HLV pgin, HLV pqout,HLV plbar,HLV pl,
HLV pqbarout, HLV p2out, HLV p2in, bool aqlinef){
//Calculate Different Helicity Configurations.
const Helicity h = aqlinef?helicity::plus:helicity::minus;
double ME2mpp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::plus,helicity::plus);
double ME2mpm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::plus,helicity::minus);
double ME2mmp = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::plus);
double ME2mmm = jM2Wuno(-pgin, pqout,plbar,pl,-pqbarout,h,p2out,p2in,helicity::minus,helicity::minus);
//Helicity sum
double ME2 = ME2mpp + ME2mpm + ME2mmp + ME2mmm;
//Correct colour averaging after crossing:
ME2*=(3.0/8.0);
return ME2;
}
double ME_WExqqx_qbarqQ(HLV pgin, HLV pqout,HLV plbar,HLV pl,
HLV pqbarout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false);
}
double ME_WExqqx_qqbarQ(HLV pgin, HLV pqbarout,HLV plbar,HLV pl,
HLV pqout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true);
}
double ME_WExqqx_qbarqg(HLV pgin, HLV pqout,HLV plbar,HLV pl,
HLV pqbarout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqout, plbar, pl, pqbarout, p2out, p2in, false)*K_g(p2out,p2in)/HEJ::C_F;
}
double ME_WExqqx_qqbarg(HLV pgin, HLV pqbarout, HLV plbar, HLV pl,
HLV pqout, HLV p2out, HLV p2in){
return jWqqx_j(pgin, pqbarout, plbar, pl, pqout, p2out, p2in, true)*K_g(p2out,p2in)/HEJ::C_F;
}
namespace {
//Function to calculate Term 1 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm1(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below. (Less Tensor class use?)
double t1 = (p3-pb)*(p3-pb);
Tensor<1> Tp3 = Construct1Tensor((p3));//p3
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
- Tensor<3> qqCurBlank = T3Current(p2,hel2,p3,hel2);
+ Tensor<3> qqCurBlank = rank3_current(p2,hel2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(Tp3-Tpb,2);
Tensor<1> 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> qggm2(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double t1 = (p2-pb)*(p2-pb);
Tensor<1> Tp2 = Construct1Tensor((p2));//p2
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb,refmom, helg);
- Tensor<3> qqCurBlank = T3Current(p2,hel2,p3,hel2);
+ Tensor<3> qqCurBlank = rank3_current(p2,hel2,p3,hel2);
Tensor<2> qqCur = qqCurBlank.contract(Tp2-Tpb,2);
Tensor<1> gqqCur = qqCur.contract(epsg,1)/t1;
return gqqCur;
}
//Function to calculate Term 3 in Equation 3.23 in James Cockburn's Thesis.
Tensor<1> qggm3(HLV pb, HLV p2, HLV p3, Helicity hel2, Helicity helg, HLV refmom){
//@TODO Simplify the calculation below (Less Tensor class use?)
double s23 = (p2+p3)*(p2+p3);
Tensor<1> Tp2 = Construct1Tensor((p2));//p2
Tensor<1> Tp3 = Construct1Tensor((p3));//p3
Tensor<1> Tpb = Construct1Tensor((pb));//pb
// Gauge choice in polarisation tensor. (see JC's Thesis)
Tensor<1> epsg = eps(pb, refmom, helg);
Tensor<2> g=metric();
Tensor<3> qqCurBlank1 = outer(Tp2+Tp3, g)/s23;
Tensor<3> qqCurBlank2 = outer(Tpb, g)/s23;
- Tensor<1> Cur23 = TCurrent(p2, p3,hel2);
+ Tensor<1> Cur23 = HEJ::current(p2, p3,hel2);
Tensor<2> qqCur1 = qqCurBlank1.contract(Cur23,3);
Tensor<2> qqCur2 = qqCurBlank2.contract(Cur23,3);
Tensor<2> qqCur3 = qqCurBlank2.contract(Cur23,1);
Tensor<1> gqqCur = (qqCur1.contract(epsg,1)
- qqCur2.contract(epsg,2)
+ qqCur3.contract(epsg,1))*2*COM(0,1);
return gqqCur;
}
}
// no wqq emission
double ME_W_Exqqx_QQq(HLV pa, HLV pb, HLV p1, HLV p2,
HLV p3,HLV plbar, HLV 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> TMmmm1 = qggm1(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMmmm2 = qggm2(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMmmm3 = qggm3(pb,p2,p3,helicity::minus,helicity::minus, pa);
Tensor<1> TMpmm1 = qggm1(pb,p2,p3,helicity::minus,helicity::plus, pa);
Tensor<1> TMpmm2 = qggm2(pb,p2,p3,helicity::minus,helicity::plus, pa);
Tensor<1> TMpmm3 = qggm3(pb,p2,p3,helicity::minus,helicity::plus, pa);
// Build the external quark line W Emmision
Tensor<1> cur1a = jW(pa,p1,plbar,pl, aqlinepa);
//Contract with the qqxCurrent.
COM Mmmm1 = TMmmm1.contract(cur1a,1);
COM Mmmm2 = TMmmm2.contract(cur1a,1);
COM Mmmm3 = TMmmm3.contract(cur1a,1);
COM Mpmm1 = TMpmm1.contract(cur1a,1);
COM Mpmm2 = TMpmm2.contract(cur1a,1);
COM Mpmm3 = TMpmm3.contract(cur1a,1);
//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 ME_WCenqqx_qq(HLV pa, HLV pb,HLV pl, HLV 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> T1am, T4bm, T1ap, T4bp;
if(!(aqlinepa)){
- T1ap = TCurrent(p1, pa, helicity::plus);
- T1am = TCurrent(p1, pa, helicity::minus);}
+ T1ap = HEJ::current(p1, pa, helicity::plus);
+ T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
- T1ap = TCurrent(pa, p1, helicity::plus);
- T1am = TCurrent(pa, p1, helicity::minus);}
+ T1ap = HEJ::current(pa, p1, helicity::plus);
+ T1am = HEJ::current(pa, p1, helicity::minus);}
if(!(aqlinepb)){
- T4bp = TCurrent(p4, pb, helicity::plus);
- T4bm = TCurrent(p4, pb, helicity::minus);}
+ T4bp = HEJ::current(p4, pb, helicity::plus);
+ T4bm = HEJ::current(p4, pb, helicity::minus);}
else if(aqlinepb){
- T4bp = TCurrent(pb, p4, helicity::plus);
- T4bm = TCurrent(pb, p4, helicity::minus);}
+ T4bp = HEJ::current(pb, p4, helicity::plus);
+ T4bm = HEJ::current(pb, p4, helicity::minus);}
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc = MUncrossW(pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> Xcro = MCrossW( pa, p1, pb, p4, pq, pqbar, pl, plbar, partons, nabove);
Tensor<2> 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));
COM M_mmCro = (((Xcro).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc).contract(T1am,1)).contract(T4bp,1));
COM M_mpCro = (((Xcro).contract(T1am,1)).contract(T4bp,1));
COM M_mpSym = (((Xsym).contract(T1am,1)).contract(T4bp,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc).contract(T1ap,1)).contract(T4bp,1));
COM M_ppCro = (((Xcro).contract(T1ap,1)).contract(T4bp,1));
COM M_ppSym = (((Xsym).contract(T1ap,1)).contract(T4bp,1));
//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.));
HLV 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 ME_W_Cenqqx_qq(HLV pa, HLV pb,HLV pl,HLV plbar, std::vector<HLV> 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> T1am(0.), T1ap(0.);
if(!(aqlinepa)){
- T1ap = TCurrent(p1, pa, helicity::plus);
- T1am = TCurrent(p1, pa, helicity::minus);}
+ T1ap = HEJ::current(p1, pa, helicity::plus);
+ T1am = HEJ::current(p1, pa, helicity::minus);}
else if(aqlinepa){
- T1ap = TCurrent(pa, p1, helicity::plus);
- T1am = TCurrent(pa, p1, helicity::minus);}
+ T1ap = HEJ::current(pa, p1, helicity::plus);
+ T1am = HEJ::current(pa, p1, helicity::minus);}
Tensor <1> T4bm = jW(pb, p4, plbar, pl, aqlinepb);
// Calculate the 3 separate contributions to the effective vertex
Tensor<2> Xunc_m = MUncross(pa, pq, pqbar,partons, helicity::minus, nabove);
Tensor<2> Xcro_m = MCross( pa, pq, pqbar,partons, helicity::minus, nabove);
Tensor<2> Xsym_m = MSym( pa, p1, pb, p4, pq, pqbar, partons, helicity::minus, nabove);
Tensor<2> Xunc_p = MUncross(pa, pq, pqbar,partons, helicity::plus, nabove);
Tensor<2> Xcro_p = MCross( pa, pq, pqbar,partons, helicity::plus, nabove);
Tensor<2> Xsym_p = MSym( pa, p1, pb, p4, pq, pqbar, partons, helicity::plus, nabove);
// (- - hel choice)
COM M_mmUnc = (((Xunc_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmCro = (((Xcro_m).contract(T1am,1)).contract(T4bm,1));
COM M_mmSym = (((Xsym_m).contract(T1am,1)).contract(T4bm,1));
// (- + hel choice)
COM M_mpUnc = (((Xunc_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpCro = (((Xcro_p).contract(T1am,1)).contract(T4bm,1));
COM M_mpSym = (((Xsym_p).contract(T1am,1)).contract(T4bm,1));
// (+ - hel choice)
COM M_pmUnc = (((Xunc_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmCro = (((Xcro_m).contract(T1ap,1)).contract(T4bm,1));
COM M_pmSym = (((Xsym_m).contract(T1ap,1)).contract(T4bm,1));
// (+ + hel choice)
COM M_ppUnc = (((Xunc_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppCro = (((Xcro_p).contract(T1ap,1)).contract(T4bm,1));
COM M_ppSym = (((Xsym_p).contract(T1ap,1)).contract(T4bm,1));
//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.));
HLV 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

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:02 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805744
Default Alt Text
(80 KB)

Event Timeline