Page MenuHomeHEPForge
Files F10664337

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
index 8ee8b4e..3472bd3 100644
--- a/include/HEJ/Tensor.hh
+++ b/include/HEJ/Tensor.hh
@@ -1,264 +1,269 @@
/** \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;
}
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);
//! Outer product of tensor with HLV
template<unsigned K>
friend Tensor<K+1> outer(Tensor<K> const & T1, CLHEP::HepLorentzVector const & p2);
//! Outer product of HLV with tensor
template<unsigned K>
friend Tensor<K+1> outer(CLHEP::HepLorentzVector const & p1, Tensor<K> 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) const;
/**
* \brief T^(mu1...mk..mN)p2_(muk) contract kth index, where k member of [1,N]
* @param p2 HepLorentzVector to contract with.
* @param k int to contract HLV p2 with from original Tensor.
* @returns T1.contract(p2,1) -> T1^(mu,nu,rho...)p2_mu
*/
Tensor<N-1> contract(CLHEP::HepLorentzVector const & p2, int k) const;
/**
* \brief "Dot" product, contracting the last open index
* @param p2 Vector (e.g. Rank 1 Tensor) to contract with
* @returns Same as T1.contract(T2, N)
*/
template<typename T>
Tensor<N-1> dot(T const & T2) const {
return contract(T2, N);
}
//! Complex conjugate tensor
Tensor<N> & complex_conj();
std::array<COM, detail::power(d, N)> components;
private:
};
template<unsigned N>
Tensor<N> operator*(COM const & x, Tensor<N> t);
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;
+ /**
+ * \brief Flip helicity
+ * @param h Helicity to flip
+ * @returns plus iff the input helicity was minus and vice versa
+ */
inline
constexpr Helicity flip(Helicity h) {
return (h == helicity::plus)?helicity::minus:helicity::plus;
}
//! Absolute square of a vector
inline
double abs2(Tensor<1> const & v) {
double res = std::norm(v(0));
for(size_t i = 1; i < detail::d; ++i) {
res -= std::norm(v(i));
}
return res;
}
//! Real part of a.dot(b.complex_conj())
// @TODO needs better name
inline
double vre(const Tensor<1> & a, const Tensor<1> & b) {
COM res = a(0)*conj(b(0));
for(size_t i = 1; i < detail::d; ++i) {
res -= a(i)*conj(b(i));
}
return real(res);
}
/**
* \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
);
//! "angle" product angle(pi, pj) = <i j>
COM angle(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj);
//! "square" product square(pi, pj) = [i j]
COM square(CLHEP::HepLorentzVector const & pi, CLHEP::HepLorentzVector const & pj);
/**
* \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 HLV class
* @param p Current in HLV format
* @returns Current in Tensor Format
*/
Tensor<1> to_tensor(CLHEP::HepLorentzVector const & 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 const & k,
CLHEP::HepLorentzVector const & ref,
Helicity pol
);
//! Initialises Tensor values by iterating over permutations of gamma matrices.
void init_sigma_index();
}
// implementation of template functions
#include "HEJ/detail/Tensor_impl.hh"

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:36 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887776
Default Alt Text
(8 KB)

Event Timeline