Page MenuHomeHEPForge
Files F8309617

No OneTemporary
Actions

View Options

diff --git a/include/HEJ/detail/Tensor_impl.hh b/include/HEJ/detail/Tensor_impl.hh
index 00c9758..5dd3f55 100644
--- a/include/HEJ/detail/Tensor_impl.hh
+++ b/include/HEJ/detail/Tensor_impl.hh
@@ -1,250 +1,203 @@
/** \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 <vector>
#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 K>
Tensor<K+1> outer(CLHEP::HepLorentzVector const & p1, Tensor<K> const & T2) {
- Tensor<K+1> product;
- auto target = begin(product.components);
- //Need to rearrange ordering from HLV standard (x,y,z,t) -> HEJ Tensor (t,x,y,z)
- std::array<double,4> v1({p1[3], p1[0], p1[1], p1[2]});
- for(auto && t1: v1) {
- std::transform(
- begin(T2.components), end(T2.components),
- target,
- [t1](auto && t2) { return t1*t2; }
- );
- std::advance(target, T2.components.size());
- }
- return product;
+ return outer(Construct1Tensor(p1), T2);
}
template<unsigned K>
Tensor<K+1> outer(Tensor<K> const & T1, CLHEP::HepLorentzVector const & p2) {
- Tensor<K+1> product;
- auto target = begin(product.components);
- //Need to rearrange ordering from HLV standard (x,y,z,t) -> HEJ Tensor (t,x,y,z)
- std::array<double,4> v2({p2[3], p2[0], p2[1], p2[2]});
- for(auto && t1: T1.components) {
- std::transform(
- begin(v2), end(v2), target,
- [t1](auto && t2) { return t1*t2; }
- );
- std::advance(target, v2.size());
- }
- return product;
+ return outer(T1, Construct1Tensor(p2));
}
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 const & 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;
+ ){
+ return contract(Construct1Tensor(p2), k);
}
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, 3:57 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4016353
Default Alt Text
(7 KB)

Event Timeline