No OneTemporary
Actions

Size

32 KB

Subscribers

None

View Options

	diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
	index a3fd0e8..d2aa6eb 100644
	--- a/include/HEJ/Tensor.hh
	+++ b/include/HEJ/Tensor.hh
	@@ -1,214 +1,218 @@
	/** \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)?basesqrtsqrt: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> 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 p2 Momentum of Particle 2
	* @param h Helicity of the Spinors
	* @returns Tensor T^mu^nu^rho = <p1\|mu nu rho\|p2>
	*
	* @note in/out configuration considered in calculation
	*/
	Tensor<3> rank3_current(
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2,
	Helicity h
	);
	/**
	* \brief Calculates current <p1\|mu nu rho tau sigma\|p2>
	* @param p1 Momentum of Particle 1
	* @param p2 Momentum of Particle 2
	* @param h Helicity of the Spinors
	* @returns Tensor T^mu^nu^rho = <p1\|mu nu rho tau sigma\|p2>
	*
	* @note in/out configuration considered in calculation
	*/
	Tensor<5> rank5_current(
	CLHEP::HepLorentzVector const & p1,
	CLHEP::HepLorentzVector const & p2,
	Helicity h
	);

	/**
	* \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);
	+ Tensor<1> eps(
	+ CLHEP::HepLorentzVector const & k,
	+ CLHEP::HepLorentzVector const & 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 704a175..e416b25 100644
	--- a/src/Tensor.cc
	+++ b/src/Tensor.cc
	@@ -1,738 +1,743 @@
	/**
	* \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];

	// 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 const & p) {
	return COM{p.x(), p.y()};
	}

	COM perp_hat(CLHEP::HepLorentzVector const & p) {
	const COM p_perp = perp(p);
	if(p_perp == COM{0., 0.}) return -1.;
	return p_perp/std::abs(p_perp);
	}

	+ COM angle(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj) {
	+ return
	+ + std::sqrt(COM{pi.minus()pj.plus()})perp_hat(pi)
	+ - std::sqrt(COM{pi.plus()pj.minus()})perp_hat(pj);
	+ }
	+
	+ COM square(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj) {
	+ return std::conj(angle(pi, pj));
	+ }
	+
	} // 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_hat = perp_hat(p);
	const COM q_perp_hat = perp_hat(q);
	std::array<std::array<COM, 2>, 2> sqrt_pq;
	sqrt_pq[minus][minus] = std::sqrt(COM{p.minus()q.minus()})p_perp_hat*conj(q_perp_hat);
	sqrt_pq[minus][plus] = std::sqrt(COM{p.minus()q.plus()})p_perp_hat;
	sqrt_pq[plus][minus] = std::sqrt(COM{p.plus()q.minus()})conj(q_perp_hat);
	sqrt_pq[plus][plus] = std::sqrt(COM{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 };

	auto j = HEJ::current(p1new, p2new, h);
	if(h == helicity::minus) {
	j *= -1;
	j(0) *= -1;
	}

	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 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 };

	auto j = HEJ::current(p1new, p2new, h);
	if(h == helicity::minus) {
	j *= -1;
	j(0) *= -1;
	}

	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 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;
	+ // polarisation vector according to eq. (55) in developer manual
	+ Tensor<1> eps(
	+ CLHEP::HepLorentzVector const & pg,
	+ CLHEP::HepLorentzVector const & pr,
	+ Helicity pol
	+ ){
	+ if(pol == helicity::plus) {
	+ return current(pr, pg, pol)/(std::sqrt(2)*square(pg, pr));
	+ }
	+ return current(pr, pg, pol)/(std::sqrt(2)*angle(pg, pr));
	}

	// 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
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 1:51 PM (15 h, 38 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023027
Default Alt Text: (32 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions