Page MenuHomeHEPForge
Files F8310311

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/Tensor.hh b/include/HEJ/Tensor.hh
index d379ae1..8ec5920 100644
--- a/include/HEJ/Tensor.hh
+++ b/include/HEJ/Tensor.hh
@@ -1,218 +1,226 @@
/** \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)?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);
/**
* \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);
+ /**
+ * \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 p2, 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 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"
diff --git a/include/HEJ/detail/Tensor_impl.hh b/include/HEJ/detail/Tensor_impl.hh
index 502e68b..848d04b 100644
--- a/include/HEJ/detail/Tensor_impl.hh
+++ b/include/HEJ/detail/Tensor_impl.hh
@@ -1,183 +1,214 @@
/** \file
* \brief Tensor Template Class Implementation.
*
* This file contains the implementation of the Tensor Template functions
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#pragma once
#include "HEJ/Tensor.hh"
#include <cassert>
#include <algorithm>
+#include "CLHEP/Vector/LorentzVector.h"
+
namespace HEJ {
namespace detail {
inline
constexpr std::size_t accumulate_idx(const std::size_t acc) {
return acc;
}
template<typename... Indices>
constexpr std::size_t accumulate_idx(
const std::size_t acc,
const std::size_t idx, Indices... indices
) {
return accumulate_idx(d*acc + idx, indices...);
}
}
template <unsigned int N>
template<typename... Indices>
COM & Tensor<N>::operator()(Indices... i) {
static_assert(
sizeof...(Indices) == N,
"number of indices must match tensor rank"
);
return components[detail::accumulate_idx(0, i...)];
}
template <unsigned int N>
template<typename... Indices>
COM const & Tensor<N>::operator()(Indices... i) const {
static_assert(
sizeof...(Indices) == N,
"number of indices must match tensor rank"
);
return components[detail::accumulate_idx(0, i...)];
}
template <unsigned int N>
Tensor<N>::operator COM() const{
static_assert(N==0, "Can only convert a scalar (rank 0 tensor) to a number");
assert(components.size() == 1);
return components[0];
}
template <unsigned int N>
Tensor<N>::Tensor(COM x) {
components.fill(x);
}
template<unsigned N>
Tensor<N> & Tensor<N>::operator*=(COM const & x) {
for(auto & c: components) c *= x;
return *this;
}
template<unsigned N>
Tensor<N> & Tensor<N>::operator/=(COM const & x) {
for(auto & c: components) c /= x;
return *this;
}
template<unsigned N>
Tensor<N> operator*(Tensor<N> t, COM const & x) {
return t *= x;
}
template<unsigned N>
Tensor<N> operator/(Tensor<N> t, COM const & x) {
return t /= x;
}
template <unsigned int N>
Tensor<N> & Tensor<N>::operator+=(Tensor<N> const & T2){
for(std::size_t i = 0; i < size(); ++ i) {
components[i] += T2.components[i];
}
return *this;
}
template <unsigned int N>
Tensor<N> & Tensor<N>::operator-=(Tensor<N> const & T2){
for(std::size_t i = 0; i < size(); ++ i) {
components[i] -= T2.components[i];
}
return *this;
}
template<unsigned N>
Tensor<N> operator+(Tensor<N> T1, Tensor<N> const & T2) {
return T1 += T2;
}
template<unsigned N>
Tensor<N> operator-(Tensor<N> T1, Tensor<N> const & T2) {
return T1 -= T2;
}
template<unsigned N>
bool operator==(Tensor<N> const & T1, Tensor<N> const & T2) {
return T1.components == T2.components;
}
template<unsigned N>
bool operator!=(Tensor<N> const & T1, Tensor<N> const & T2) {
return !(T1 == T2);
}
template<unsigned K, unsigned M>
Tensor<K+M> outer(Tensor<K> const & T1, Tensor<M> const & T2) {
Tensor<K+M> product;
auto target = begin(product.components);
for(auto && t1: T1.components) {
std::transform(
begin(T2.components), end(T2.components),
target,
[t1](auto && t2) { return t1*t2; }
);
std::advance(target, T2.components.size());
}
return product;
}
template<unsigned N>
template<unsigned M>
Tensor<N+M> Tensor<N>::outer(Tensor<M> const & T2) const {
return outer(*this, T2);
}
namespace detail {
inline
COM metric_sign(unsigned int i){
if(i==0)
return 1.;
else
return -1.;
}
}
template <unsigned int N>
Tensor<N-1> Tensor<N>::contract(Tensor<1> const & T2, const int k){
Tensor<N-1> result;
// distance between consecutive contracted elements in this Tensor
const std::size_t step = detail::power(d, N-k);
for(std::size_t res_idx = 0; res_idx < result.size(); ++res_idx){
COM & res = result.components[res_idx];
// Assume we compute T_{i_1,...i_N} * T2_{i_k}.
// The index into our components (T) is
// index = d^N*i_1 + d^{N-1}*i_2 + ... + i_N
//
// The corresponding index in the resulting contracted vector is
// res_idx = d^{N-1}*i_1 + ... + d^{N-k+1}*i_{k-1}
// + d^{N-k}*i_{k+1} + ... + i_N
// The second line is just (res_idx mod d^{N-k}):
const std::size_t remainder = res_idx % step;
// To get to the starting index (i_k = 0) into our components from res_idx,
// we have to multiply the first line by d and add the second line:
std::size_t our_idx = (res_idx - remainder)*d + remainder;
for (unsigned int i = 0; i < d; ++i){
res += components[our_idx]*T2.components[i]*detail::metric_sign(i);
our_idx += step;
}
}
return result;
}
template <unsigned int N>
+ Tensor<N-1> Tensor<N>::contract(CLHEP::HepLorentzVector p2, const int k){
+ Tensor<N-1> result;
+ // distance between consecutive contracted elements in this Tensor
+ const std::size_t step = detail::power(d, N-k);
+ for(std::size_t res_idx = 0; res_idx < result.size(); ++res_idx){
+ COM & res = result.components[res_idx];
+ // Assume we compute T_{i_1,...i_N} * T2_{i_k}.
+ // The index into our components (T) is
+ // index = d^N*i_1 + d^{N-1}*i_2 + ... + i_N
+ //
+ // The corresponding index in the resulting contracted vector is
+ // res_idx = d^{N-1}*i_1 + ... + d^{N-k+1}*i_{k-1}
+ // + d^{N-k}*i_{k+1} + ... + i_N
+ // The second line is just (res_idx mod d^{N-k}):
+ const std::size_t remainder = res_idx % step;
+ // To get to the starting index (i_k = 0) into our components from res_idx,
+ // we have to multiply the first line by d and add the second line:
+ std::size_t our_idx = (res_idx - remainder)*d + remainder;
+
+ //Need to rearrange ordering from HLV standard (t,x,y,z) - > HEJ Tensor (x,y,z,t)
+ for (unsigned int i = 0; i < d; ++i){
+ res += components[our_idx+i*step]*p2[(i+3)%4]*detail::metric_sign(i);
+ }
+ }
+ return result;
+ }
+
+
+ template <unsigned int N>
Tensor<N> & Tensor<N>::complex_conj() {
for(auto & c: components) {
c = std::conj(c);
}
return *this;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:04 PM (9 h, 54 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023735
Default Alt Text
(13 KB)

Event Timeline