diff --git a/src/Tensor.cc b/src/Tensor.cc
index 4972e69..5aad2f4 100644
--- a/src/Tensor.cc
+++ b/src/Tensor.cc
@@ -1,766 +1,741 @@
/**
* \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.} };
// map from a list of tensor lorentz indices onto a single index
// in d dimensional spacetime
// TODO: basically the same as detail::accumulate_idx
template<std::size_t N>
std::size_t tensor_to_list_idx(std::array<int, N> const & indices) {
std::size_t list_idx = 0;
for(std::size_t i = 1, factor = 1; i <= N; ++i) {
list_idx += factor*indices[N-i];
factor *= detail::d;
}
assert(list_idx < (1u << 2*N));
return list_idx;
}
// 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[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd perm
// if this perm is odd, swapping diff1<->diff2 automatically even
set_permfactor2=false;
} else {
permfactor5[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor_to_list_idx(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[tensor_to_list_idx(tempvec)]=-1.;
} else { // diff is sigma0
if( ((samepos2-samepos-1)%3>0)
&& (abs(abs(samepos2-diffpos)-abs(samepos-diffpos))%3>0) )
permfactor5[tensor_to_list_idx(tempvec)]=-1.;
}
} else { // same2a is sigma0
if((diffpos>samepos) && (diffpos<samepos2))
permfactor5[tensor_to_list_idx(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[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos1>diffpos2 && diffpos2>diffpos3)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else if(diffpos3>diffpos1 && diffpos1>diffpos2)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);// odd
else
permfactor5[tensor_to_list_idx(tempvec)]=COM(0,1);// evwn
} else {
if(permcount%2==1)
permfactor5[tensor_to_list_idx(tempvec)]=-1.*COM(0,1);
else
permfactor5[tensor_to_list_idx(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[tensor_to_list_idx(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[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
} else if(permcount%2==1){
permfactor3[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd
} else {
permfactor3[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
}
tempvec[pos1]=diff3;
tempvec[pos2]=diff2;
perms.push_back(tempvec);
if(sig_zero){
permfactor3[tensor_to_list_idx(tempvec)]=-1.*COM(0,1); // odd
} else if(permcount%2==1){
permfactor3[tensor_to_list_idx(tempvec)]=1.*COM(0,1); // even
} else {
permfactor3[tensor_to_list_idx(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);
}
}
COM perp(CLHEP::HepLorentzVector p) {
return COM{p.x(), p.y()};
}
} // 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>
// see eqs. (48), (49) in developer manual
Tensor<1> current(
CLHEP::HepLorentzVector const & p,
CLHEP::HepLorentzVector const & q,
Helicity h
){
using namespace helicity;
const COM p_perp = perp(p);
const COM q_perp = perp(q);
const COM p_perp_hat = (p_perp != COM{0., 0.})?p_perp/std::abs(p_perp):-1.;
const COM q_perp_hat = (q_perp != COM{0., 0.})?q_perp/std::abs(q_perp):-1.;
std::array<std::array<COM, 2>, 2> sqrt_pq;
sqrt_pq[minus][minus] = std::sqrt(p.minus()*q.minus())*p_perp_hat*conj(q_perp_hat);
sqrt_pq[minus][plus] = std::sqrt(p.minus()*q.plus())*p_perp_hat;
sqrt_pq[plus][minus] = std::sqrt(p.plus()*q.minus())*conj(q_perp_hat);
sqrt_pq[plus][plus] = std::sqrt(p.plus()*q.plus());
Tensor<1> j;
j(0) = sqrt_pq[plus][plus] + sqrt_pq[minus][minus];
j(1) = sqrt_pq[plus][minus] + sqrt_pq[minus][plus];
j(2) = COM{0., 1.}*(sqrt_pq[plus][minus] - sqrt_pq[minus][plus]);
j(3) = sqrt_pq[plus][plus] - sqrt_pq[minus][minus];
return (h == minus)?j:j.complex_conj();
}
// <1|mu nu sigma|2>
+ // TODO: how does this work?
Tensor<3> rank3_current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
){
+ const CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
+ const CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
- 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 };
-
- Spinor(p1new, h, sp1);
- Spinor(p2new, h, 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];
- }
+ auto j = HEJ::current(p1new, p2new, h);
+ if(h == helicity::minus) {
+ j *= -1;
+ j(0) *= -1;
}
- for( std::size_t itensor=0; itensor<newT.size(); itensor++ ){
- double hfact = double( helfactor3[h][itensor] );
- newT.components[itensor] = current[sigma_index3[itensor]] * hfact
- * permfactor3[itensor];
+ Tensor<3> j3;
+ for(std::size_t tensor_idx=0; tensor_idx < j3.size(); ++tensor_idx){
+ const double hfact = helfactor3[h][tensor_idx];
+ j3.components[tensor_idx] = j(sigma_index3[tensor_idx]) * hfact
+ * permfactor3[tensor_idx];
}
- return newT;
+ return j3;
}
// <1|mu nu sigma tau rho|2>
Tensor<5> rank5_current(
CLHEP::HepLorentzVector const & p1,
CLHEP::HepLorentzVector const & p2,
Helicity h
){
+ const CLHEP::HepLorentzVector p1new{ p1.e()<0?-p1:p1 };
+ const CLHEP::HepLorentzVector p2new{ p2.e()<0?-p2:p2 };
- 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 };
-
- Spinor(p1new, h, sp1);
- Spinor(p2new, h, 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];
- }
+ auto j = HEJ::current(p1new, p2new, h);
+ if(h == helicity::minus) {
+ j *= -1;
+ j(0) *= -1;
}
- for(std::size_t itensor=0;itensor<newT.size();itensor++){
- double hfact = double(helfactor5[h][itensor]);
- newT.components[itensor] = current[sigma_index5[itensor]] * hfact
- * permfactor5[itensor];
+ Tensor<5> j5;
+ for(std::size_t tensor_idx=0; tensor_idx < j5.size(); ++tensor_idx){
+ const double hfact = helfactor5[h][tensor_idx];
+ j5.components[tensor_idx] = j(sigma_index5[tensor_idx]) * hfact
+ * permfactor5[tensor_idx];
}
- return newT;
+ return j5;
}
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 = 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 = tensor_to_list_idx(sigma0indices.at(i));
int index1 = tensor_to_list_idx(sigma1indices.at(i));
int index2 = tensor_to_list_idx(sigma2indices.at(i));
int index3 = tensor_to_list_idx(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 = tensor_to_list_idx(sigma0indices3.at(i));
int index1 = tensor_to_list_idx(sigma1indices3.at(i));
int index2 = tensor_to_list_idx(sigma2indices3.at(i));
int index3 = tensor_to_list_idx(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
}