Page MenuHomeHEPForge
Files F8310408

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
index 55e0ce8..b02ae30 100644
--- a/include/HEJ/Tensor.hh
+++ b/include/HEJ/Tensor.hh
@@ -1,203 +1,201 @@
/** \file
* \brief Tensor Template Class declaration.
*
* This file contains the declaration of the Tensor Template class. This
* is used to calculate some of the more complex currents within the
* W+Jets implementation particularly.
*/
#pragma once
#include <array>
///@TODO remove function implementation from header
///@TODO put in some namespace
template <unsigned int N, unsigned int D>
class Tensor{
public:
//Constructor
Tensor();
Tensor(COM x);
//Destructor
virtual ~Tensor();
int rank(){
return N;
}
int dim(){
return D;
}
int len(){
return size;
}
COM at(int i){
return components[i];
}
COM at(int i, int j) {
return components[D*i +j];
}
COM at(int i, int j, int k) {
return components[D*(D*i + j)+ k];
}
COM at(int i,int j, int k,int l) {
return components[D*(D*(D*i +j) + k) + l];
}
COM at(int i,int j, int k,int l, int m){
return components[D*(D*(D*(D*i + j) + k) + l) + m];
}
bool isSet(){
if(components.size()==0)
return false;
else
return true;
}
void Fill(COM x){
components=x;
}
//Set component indexed by i,j,k,l,m
void Set(int i,COM x){
components[i] = x;
}
void Set(int i, int j, COM x) {
components[D*i +j] = x;
}
void Set(int i, int j, int k, COM x) {
components[D*(D*i + j)+ k] = x;
}
void Set(int i,int j, int k,int l,COM x) {
components[D*(D*(D*i +j) + k) + l] = x;
}
void Set(int i,int j, int k,int l, int m, COM x){
components[D*(D*(D*(D*i + j) + k) + l) + m] = x;
}
Tensor<N,D> operator*(const double x){
Tensor<N,D> newT;
newT.components=components*COM(x,0);
return newT;
}
Tensor<N,D> operator*(const COM x){
Tensor<N,D> newT;
newT.components=components*x;
return newT;
}
Tensor<N,D> operator/(const double x){
Tensor<N,D> newT;
newT.components=components/COM(x,0);
return newT;
}
Tensor<N,D> operator/(const COM x){
Tensor<N,D> newT;
newT.components=components/x;
return newT;
}
Tensor<N,D> operator+(const Tensor<N,D> T2){
Tensor<N,D> newT;
newT.components=components+T2.components;
return newT;
}
Tensor<N,D> operator-(const Tensor<N,D> T2){
Tensor<N,D> newT;
newT.components=components-T2.components;
return newT;
}
void operator+=(const Tensor<N,D> T2){
components = components+T2.components;
}
void operator-=(const Tensor<N,D> T2){
components=components-T2.components;
}
Tensor<N+1,D> rightprod(const Tensor<1,D> T2){
Tensor<N+1,D> newT;
for(int i=0; i<size;i++){
for(unsigned int j=0;j<D;j++){
newT.components[i*D+j]=components[i]*T2.components[j];
}
}
return newT;
}
Tensor<N+1,D> leftprod(const Tensor<1,D> T2){
Tensor<N+1,D> newT;
for(unsigned int j=0;j<D;j++){
for(int i=0; i<size;i++){
newT.components[j*size+i]=components[i]*T2.components[j];
}
}
return newT;
}
//T^(mu1...mk..mN)T2_(muk) contract kth index, where k member of [1,N]
Tensor<N-1,D> contract(const Tensor<1,D> T2, int k){
Tensor<N-1,D> newT;
for(int j=0; j<newT.len(); j++){
COM temp;
int itemp = pow(D,(N-k));
for (unsigned int i=0; i<D; i++){
int index = D*itemp*floor(j/itemp) + itemp*i +j%(itemp);
temp+=components[index]*T2.components[i]*sign(i);
}
newT.components[j]=temp;
}
return newT;
}
std::valarray<COM> components;
private:
int size;
COM sign(unsigned int i){
if(i==0)
return 1.;
else
return -1.;
}
};
template <unsigned int N, unsigned int D> Tensor<N,D>::Tensor()
{
size = pow(D,N);
components.resize(size);
}
template <unsigned int N, unsigned int D> Tensor<N,D>::Tensor(COM x) {
size = pow(D,N);
components.resize(size, x);
}
template <unsigned int N, unsigned int D> Tensor<N,D>::~Tensor() {}
-
// Tensor Functions:
// Tensor<1,4> Sigma(int i, int j, bool hel);
// Tensor<2,4> Metric();
// int tensor2listindex(std::array<int,5> indexlist);
// int tensor2listindex(std::array<int,3> indexlist);
// void perms41(int same4, int diff, std::vector<std::array<int,5>> * perms);
// void perms32(int same3, int diff, std::vector<std::array<int,5>> * perms);
// void perms311(int same3, int diff1, int diff2, std::vector<std::array<int,5>> * perms);
// void perms221(int same2a, int same2b, int diff, std::vector<std::array<int,5>> * perms);
// void perms2111(int same2, int diff1,int diff2,int diff3, std::vector<std::array<int,5>> * perms);
// void perms21(int same, int diff, std::vector<std::array<int,3>> * perms);
// void perms111(int diff1, int diff2, int diff3, std::vector<std::array<int,3>> * perms);
-
Tensor<2,4> Metric();
Tensor<1,4> TCurrent(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2);
Tensor<3,4> T3Current(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2);
Tensor<5,4> T5Current(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2);
Tensor<1,4> Construct1Tensor(CCurrent j);
Tensor<1,4> Construct1Tensor(CLHEP::HepLorentzVector p);
Tensor<1,4> eps(CLHEP::HepLorentzVector k, CLHEP::HepLorentzVector ref, bool pol);
bool init_sigma_index();
diff --git a/src/Tensor.cc b/src/Tensor.cc
index 40032f9..0fba3bb 100644
--- a/src/Tensor.cc
+++ b/src/Tensor.cc
@@ -1,824 +1,790 @@
#include "HEJ/currents.hh"
#include "HEJ/Tensor.hh"
#include <array>
#include <iostream>
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,4> Sigma(int i, int j, bool hel){
Tensor<1,4> newT;
if(hel){
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];
+ } 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){
+ 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);
+ 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){
+ 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);
+ 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){
+ 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);
+ perms.push_back(tempvec);
- if(!same0&&!diff0){
+ 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{
+ } 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{
+ } 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);
+ 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
+ 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);
+ 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
+ } 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))
+ } 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
+ 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);
+ perms.push_back(tempvec);
if(twozero){// don't care about exact positions of singles, just order
- if(diffpos2>diffpos3&&diffpos3>diffpos1)
+ if(diffpos2>diffpos3 && diffpos3>diffpos1)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
- else if(diffpos1>diffpos2&&diffpos2>diffpos3)
+ else if(diffpos1>diffpos2 && diffpos2>diffpos3)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
- else if(diffpos3>diffpos1&&diffpos1>diffpos2)
+ else if(diffpos3>diffpos1 && diffpos1>diffpos2)
permfactor5[tensor2listindex(tempvec)]=-1.*COM(0,1);// odd
else
permfactor5[tensor2listindex(tempvec)]=COM(0,1);// evwn
- }else{
+ } 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){
+ 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)
+ 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
+ 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);
+ perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
- }else if(permcount%2==1){
+ } else if(permcount%2==1){
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
- }else{
+ } else {
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
}
tempvec[pos1]=diff3;
tempvec[pos2]=diff2;
- perms->push_back(tempvec);
+ perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
- }else if(permcount%2==1){
+ } else if(permcount%2==1){
permfactor3[tensor2listindex(tempvec)]=1.*COM(0,1); // even
- }else{
+ } else {
permfactor3[tensor2listindex(tempvec)]=-1.*COM(0,1); // odd
}
permcount++;
}
-
}
-
}
-
void SpinorO(CLHEP::HepLorentzVector p, bool 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=+
if(hel){
sp[0] = sqpp;
sp[1] = sqpm*pperp/abs(pperp);
- }else{
+ } else {
sp[0] = sqpm*conj(pperp)/abs(pperp);
sp[1] = -sqpp;
}
}
void SpinorIp(COM sqpp, bool hel, COM *sp){
// if hel=+
if(hel){
sp[0] = sqpp;
sp[1] = 0.;
- }else{
+ } else {
sp[0] = 0.;
sp[1] = -sqpp;
}
}
void SpinorIm(COM sqpm, bool hel, COM *sp){
// if hel=+
if(hel){
sp[0] = 0;
sp[1] = -sqpm;
- }else{
+ } else {
sp[0] = -sqpm;
sp[1] = 0.;
}
}
-
void Spinor(CLHEP::HepLorentzVector p, bool 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{
+ else {
SpinorO(p,hel,sp);
}
}
} // anonymous namespace
Tensor<2,4> Metric(){
Tensor<2,4> g(0.);
g.Set(0,0, 1.);
g.Set(1,1, -1.);
g.Set(2,2, -1.);
g.Set(3,3, -1.);
return g;
}
// <1|mu|2>
Tensor<1,4> TCurrent(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2
){
COM sp1[2];
COM sp2[2];
Tensor<1,4> newT(0.);
- CLHEP::HepLorentzVector p1new=p1;
- CLHEP::HepLorentzVector p2new=p2;
-
- if(p1.e()<0){
- p1new=-p1;
- }
-
- if(p2.e()<0){
- p2new=-p2;
- }
+ 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,4> T3Current(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2
){
COM sp1[2];
COM sp2[2];
Tensor<3,4> newT(0.);
- CLHEP::HepLorentzVector p1new=p1;
- CLHEP::HepLorentzVector p2new=p2;
-
- if(p1.e()<0){
- p1new=-p1;
- }
-
- if(p2.e()<0){
- p2new=-p2;
- }
+ 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( int itensor=0; itensor<newT.len(); 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,4> T5Current(CLHEP::HepLorentzVector p1, bool h1,
CLHEP::HepLorentzVector p2, bool h2
){
COM sp1[2];
COM sp2[2];
Tensor<5,4> newT(0.);
- CLHEP::HepLorentzVector p1new=p1;
- CLHEP::HepLorentzVector p2new=p2;
-
- if(p1.e()<0){
- p1new=-p1;
- }
-
- if(p2.e()<0){
- p2new=-p2;
- }
+ 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(int itensor=0;itensor<newT.len();itensor++){
double hfact = double(helfactor5[h2][itensor]);
newT.components[itensor] = current[sigma_index5[itensor]] * hfact
* permfactor5[itensor];
}
return newT;
}
Tensor<1,4> Construct1Tensor(CCurrent j){
Tensor<1,4> newT;
newT.components={j.c0,j.c1,j.c2,j.c3};
return newT;
}
Tensor<1,4> Construct1Tensor(CLHEP::HepLorentzVector p){
Tensor<1,4> newT;
newT.components={p.e(),p.x(),p.y(),p.z()};
return newT;
}
Tensor<1,4> eps(CLHEP::HepLorentzVector k, CLHEP::HepLorentzVector ref, bool pol){
Tensor<1,4> polvec;
COM spk[2];
COM spr[2];
COM denom;
- CLHEP::HepLorentzVector knew=k;
-
- if(k.e()<0){
- knew=-k;
- }
+ CLHEP::HepLorentzVector knew{ k.e()<0?-k:k };
Spinor(knew, pol, spk);
Spinor(ref, !pol, spr);
denom = pow(-1.,pol)*sqrt(2)*(conj(spr[0])*spk[0] + conj(spr[1])*spk[1]);
polvec = TCurrent(ref,!pol,knew,!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);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:24 PM (8 h, 1 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023788
Default Alt Text
(35 KB)

Event Timeline