diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
index abc617d..7d50cdc 100644
--- a/include/HEJ/Tensor.hh
+++ b/include/HEJ/Tensor.hh
@@ -1,219 +1,206 @@
/** \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>
///@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;
}
/**
* \brief Multiply Tensor from Right: T1^(mu1...mk..mN-1)T2_(muN)
* @param T2 Tensor of Rank 1 to multiply from right with with.
* @returns T1.contract(T2,1) -> T1^(mu,nu,rho...)T2^sigma
*/
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;
}
/**
* \brief Multiply Tensor from Left: T2_(muN)T1^(mu1...mk..mN-1)
* @param T2 Tensor of Rank 1 to multiply from left with with.
* @returns T1.contract(T2,1) -> T2^sigma T1^(mu,nu,rho...)
*/
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;
}
/**
* \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,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();