Page MenuHomeHEPForge
Files F8309980

No OneTemporary
Actions

View Options

diff --git a/Config/Complex.h b/Config/Complex.h
--- a/Config/Complex.h
+++ b/Config/Complex.h
@@ -1,82 +1,28 @@
// -*- C++ -*-
//
// Complex.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Complex_H
#define ThePEG_Complex_H
//
// This is file wraps the standard complex header and makes some
// convenient typedefs in the ThePEG namespace.
//
#include <complex>
namespace ThePEG {
using std::complex;
-/** ThePEG code should use Complex for all complex scalars */
-typedef std::complex<double> Complex;
-
-/** @cond TRAITSPECIALIZATIONS */
-
-template <typename T, typename U>
-struct BinaryOpTraits;
-
-
-template <typename T, typename U>
-struct BinaryOpTraits<complex<T>, U> {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef complex<typename BinaryOpTraits<T,U>::MulT> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef complex<typename BinaryOpTraits<T,U>::DivT> DivT;
-
-};
-
-template <typename T, typename U>
-struct BinaryOpTraits<T, complex<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef complex<typename BinaryOpTraits<T,U>::MulT> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef complex<typename BinaryOpTraits<T,U>::DivT> DivT;
-
-};
-
-template <typename T, typename U>
-struct BinaryOpTraits<complex<T>, complex<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef complex<typename BinaryOpTraits<T,U>::MulT> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef complex<typename BinaryOpTraits<T,U>::DivT> DivT;
-
-};
-
-template <typename T>
-struct BinaryOpTraits<complex<T>, complex<T> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef complex<typename BinaryOpTraits<T,T>::MulT> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef complex<typename BinaryOpTraits<T,T>::DivT> DivT;
-
- /** @endcond */
-
-};
-
-
+/// ThePEG code should use Complex for all complex scalars
+using Complex = std::complex<double>;
}
#endif /* ThePEG_Complex_H */
diff --git a/Config/PhysicalQty.h b/Config/PhysicalQty.h
--- a/Config/PhysicalQty.h
+++ b/Config/PhysicalQty.h
@@ -1,337 +1,321 @@
// -*- C++ -*-
//
// PhysicalQty.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Physical_Qty_H
#define Physical_Qty_H
#include "TemplateTools.h"
#include <sstream>
+#include <ratio>
+#include <type_traits>
/** @file
*
* The PhysicalQty class allows compile-time checking of dimensional
* correctness. Mathematical operations that are inconsistent are
* flagged as type errors.
*
* Do not use the classes directly in ThePEG, use the wrappers defined
* in Units.h or Phys_Qty.h instead.
*/
namespace ThePEG {
/// Helper class to construct zero unitful quantities.
struct ZeroUnit {
/** Automatic conversion to double. */
constexpr operator double() const { return 0.0; }
};
/// ZERO can be used as zero for any unitful quantity.
constexpr ZeroUnit ZERO = ZeroUnit();
-/// Helper classes to extend or shorten fractions
-//@{
-/**
- * Template to help with fractional powers of dimensions
- */
-template <int M, int II>
-struct QtyHelper
-{
- /// The numerator, indicating failure.
- static constexpr int I = -999999;
-};
-
-/**
- * Template to help with fractional powers of dimensions
- */
-template <int II>
-struct QtyHelper<0,II>
-{
- /// The new numerator.
- static constexpr int I = II;
-};
-
-/**
- * Template to help with fractional powers of dimensions
- */
-template <int II, int DI, int DI2>
-struct QtyInt
-{
- /// The new numerator.
- static constexpr int I = QtyHelper<(DI2*II)%DI,(DI2*II)/DI>::I;
-};
-//@}
-
/**
* This template class allows the compiler to check calculations with
* physical quantities for dimensional correctness. A quantity can be
* composed of arbitrary fractional powers of length L, energy E and
* charge Q. Commonly used quantities should be typedef'ed (see Units.h).
*
* Some member functions can break dimensional consistency if the user
* is not careful; these are marked explicitly.
*
* Do not use this class directly in ThePEG, use the pre-defined quantities
* from Units.h or the wrapper in Phys_Qty.h instead.
*/
-template<int L, int E, int Q, int DL = 1, int DE = 1, int DQ = 1>
-class Qty
+
+// only specialization is with std::ratio below
+template <typename L, typename E, typename T>
+class Qty;
+
+template <typename T, typename U>
+struct qty_equal {
+ static constexpr bool value = false;
+};
+
+template<typename L1, typename L2, typename E1, typename E2, typename Q1, typename Q2>
+struct qty_equal<Qty<L1,E1,Q1>, Qty<L2,E2,Q2>> {
+ static constexpr bool value
+ = std::ratio_equal<L1,L2>::value
+ && std::ratio_equal<E1,E2>::value
+ && std::ratio_equal<Q1,Q2>::value;
+};
+
+template <typename T>
+struct is_qty {
+ static constexpr bool value = qty_equal<T,T>::value;
+};
+
+template <typename ResultT, typename T, typename U = T>
+using enable_if_same_qty = typename std::enable_if<qty_equal<T,U>::value, ResultT>::type;
+
+template<int L, int E, int Q, int DL, int DE, int DQ>
+class Qty<std::ratio<L,DL>, std::ratio<E,DE>, std::ratio<Q,DQ>>
{
private:
/// Constructor from raw values. Breaks consistency.
constexpr Qty(double val) : rawValue_(val) {}
public:
/// The name of the class for persistent IO
static std::string className() {
std::ostringstream os;
os << "Qty<"
<< L << ','
<< E << ','
<< Q << ','
<< DL << ','
<< DE << ','
<< DQ << '>';
return os.str();
}
-
+ /// General power type
+ template <int Num, int Den>
+ using Power = Qty<typename std::ratio<Num*L,Den*DL>::type,
+ typename std::ratio<Num*E,Den*DE>::type,
+ typename std::ratio<Num*Q,Den*DQ>::type>;
+ /// Our type
+ using Type = Power<1,1>;
/// The squared type.
- typedef Qty<2*L,2*E,2*Q,DL,DE,DQ> Squared;
+ using Squared = Power<2,1>;
+ /// The inverse type.
+ using Inverse = Power<-1,1>;
+ /// The sqrt type.
+ using Sqrt = Power<1,2>;
/// Basic unit of this quantity.
- static constexpr Qty<L,E,Q,DL,DE,DQ> baseunit()
+ static constexpr Type baseunit()
{
- return Qty<L,E,Q,DL,DE,DQ>(1.0);
+ return Type(1.0);
}
/// Default constructor to 0.
constexpr Qty() : rawValue_(0.0) {}
/// Default constructor to 0.
constexpr Qty(ZeroUnit) : rawValue_(0.0) {}
/// Constructor from a compatible quantity
- template <int DL2, int DE2, int DQ2>
- constexpr
- Qty(const Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> & q,
- double factor = 1.0)
+ template <typename U>
+ constexpr Qty(const U & q, double factor = 1.0,
+ enable_if_same_qty<void,Type,U> * = nullptr)
: rawValue_(q.rawValue() * factor) {}
/// Access to the raw value. Breaks consistency.
constexpr double rawValue() const { return rawValue_; }
/// Assignment multiplication by dimensionless number.
- Qty<L,E,Q,DL,DE,DQ> & operator*=(double x) { rawValue_ *= x; return *this; }
+ Type & operator*=(double x) { rawValue_ *= x; return *this; }
/// Assignment division by dimensionless number.
- Qty<L,E,Q,DL,DE,DQ> & operator/=(double x) { rawValue_ /= x; return *this; }
+ Type & operator/=(double x) { rawValue_ /= x; return *this; }
/// Assignment addition with compatible quantity.
- template <int DL2, int DE2, int DQ2>
- Qty<L,E,Q,DL,DE,DQ> &
- operator+=(const Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> x)
+ Type & operator+=(const Type & x)
{
rawValue_ += x.rawValue();
return *this;
}
/// Assignment subtraction with compatible quantity.
- template <int DL2, int DE2, int DQ2>
- Qty<L,E,Q,DL,DE,DQ> &
- operator-=(const Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> x)
+ Type & operator-=(const Type & x)
{
rawValue_ -= x.rawValue();
return *this;
}
private:
/// The raw value in units of Qty::baseunit().
double rawValue_;
};
/// Specialization of Qty for <0,0,0> with conversions to double.
-template<int DL, int DE, int DQ>
-class Qty<0,0,0,DL,DE,DQ>
+template<>
+class Qty<std::ratio<0>,std::ratio<0>,std::ratio<0>>
{
public:
+ /// Our type
+ using Type = Qty<std::ratio<0>,std::ratio<0>,std::ratio<0>>;
+ /// General power type
+ template <int Num, int Den>
+ using Power = Type;
/// The squared type.
- typedef double Squared;
+ using Squared = Type;
+ /// The inverse type.
+ using Inverse = Type;
+ /// The sqrt type.
+ using Sqrt = Type;
+
/// Basic unit of this quantity.
- static constexpr double baseunit() {
+ static constexpr Type baseunit() {
return 1.0;
}
/// Default constructor to 0.
constexpr Qty(ZeroUnit) : rawValue_(0.0) {}
/// Default constructor from a double.
constexpr Qty(double x = 0.0, double factor=1.0)
: rawValue_(x * factor) {}
/// Constructor from a compatible quantity
- template <int DL2, int DE2, int DQ2>
- constexpr
- Qty(const Qty<0,0,0,DL2,DE2,DQ2> & q, double factor=1.0)
+ template <typename U>
+ constexpr Qty(const U & q, double factor=1.0,
+ enable_if_same_qty<void,Type,U> * = nullptr)
: rawValue_(q.rawValue() * factor) {}
/// Access to the raw value.
constexpr double rawValue() const { return rawValue_; }
/// Cast to double.
constexpr operator double() const { return rawValue_; }
/// Assignment multiplication by dimensionless number.
- Qty<0,0,0,DL,DE,DQ> & operator*=(double x) { rawValue_ *= x; return *this; }
+ Type & operator*=(double x) { rawValue_ *= x; return *this; }
/// Assignment division by dimensionless number.
- Qty<0,0,0,DL,DE,DQ> & operator/=(double x) { rawValue_ /= x; return *this; }
+ Type & operator/=(double x) { rawValue_ /= x; return *this; }
/// Assignment addition with compatible quantity.
- template <int DL2, int DE2, int DQ2>
- Qty<0,0,0,DL,DE,DQ> & operator+=(const Qty<0,0,0,DL2,DE2,DQ2> x) {
+ Type & operator+=(const Type & x) {
rawValue_ += x.rawValue();
return *this;
}
/// Assignment subtraction with compatible quantity.
- template <int DL2, int DE2, int DQ2>
- Qty<0,0,0,DL,DE,DQ> & operator-=(const Qty<0,0,0,DL2,DE2,DQ2> x) {
+ Type & operator-=(const Type & x) {
rawValue_ -= x.rawValue();
return *this;
}
/// Assignment addition with double.
- Qty<0,0,0,DL,DE,DQ> & operator+=(double x) {
+ Type & operator+=(double x) {
rawValue_ += x;
return *this;
}
/// Assignment subtraction with double.
- Qty<0,0,0,DL,DE,DQ> & operator-=(double x) {
+ Type & operator-=(double x) {
rawValue_ -= x;
return *this;
}
private:
/// The raw value.
double rawValue_;
};
+using QtyDouble = Qty<std::ratio<0>,std::ratio<0>,std::ratio<0>>;
+
/// @name Result types for binary operations.
//@{
/**
* BinaryOpTraits should be specialized with typdefs called MulT and
* DivT which gives the type resulting when multiplying and dividing
* the template argument types respectively.
*/
template <typename T, typename U>
struct BinaryOpTraits;
/** @cond TRAITSPECIALIZATIONS */
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-struct BinaryOpTraits<Qty<L1,E1,Q1,DL1,DE1,DQ1>,
- Qty<L2,E2,Q2,DL2,DE2,DQ2> > {
+template<typename L1, typename L2,
+ typename E1, typename E2,
+ typename Q1, typename Q2>
+struct BinaryOpTraits<Qty<L1,E1,Q1>, Qty<L2,E2,Q2>> {
/** The type resulting from multiplication of the template type with
itself. */
- typedef Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> MulT;
+ typedef Qty<std::ratio_add<L1,L2>,
+ std::ratio_add<E1,E2>,
+ std::ratio_add<Q1,Q2>> MulT;
/** The type resulting from division of one template type with
another. */
- typedef Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> DivT;
-};
-
-
-template<int L1, int E1, int Q1, int DL1, int DE1, int DQ1>
-struct BinaryOpTraits<Qty<L1,E1,Q1,DL1,DE1,DQ1>,
- Qty<L1,E1,Q1,DL1,DE1,DQ1> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef Qty<2*L1,2*E1,2*Q1,
- DL1,DE1,DQ1> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef double DivT;
+ typedef Qty<std::ratio_subtract<L1,L2>,
+ std::ratio_subtract<E1,E2>,
+ std::ratio_subtract<Q1,Q2>> DivT;
};
/**
* Multiplication template
*/
-template<int L1, int E1, int Q1, int DL1, int DE1, int DQ1>
-struct BinaryOpTraits<double,
- Qty<L1,E1,Q1,DL1,DE1,DQ1> > {
+template<typename L, typename E, typename Q>
+struct BinaryOpTraits<double, Qty<L,E,Q>> {
/** The type resulting from multiplication of the template type */
- typedef Qty<L1,E1,Q1,
- DL1,DE1,DQ1> MulT;
+ typedef Qty<L,E,Q> MulT;
/** The type resulting from division of the template type */
- typedef Qty<-L1,-E1,-Q1,
- DL1,DE1,DQ1> DivT;
+ typedef typename BinaryOpTraits<QtyDouble, Qty<L,E,Q>>::DivT DivT;
};
/**
* Multiplication template
*/
-template<int L1, int E1, int Q1, int DL1, int DE1, int DQ1>
-struct BinaryOpTraits<Qty<L1,E1,Q1,DL1,DE1,DQ1>,
- double> {
+template<typename L, typename E, typename Q>
+struct BinaryOpTraits<Qty<L,E,Q>, double> {
/** The type resulting from multiplication of the template type */
- typedef Qty<L1,E1,Q1,
- DL1,DE1,DQ1> MulT;
+ typedef Qty<L,E,Q> MulT;
/** The type resulting from division of the template type */
- typedef Qty<L1,E1,Q1,
- DL1,DE1,DQ1> DivT;
+ typedef Qty<L,E,Q> DivT;
};
//@}
/// @name Type traits for alternative code generation.
//@{
/** Type traits for alternative code generation*/
-template <int L, int E, int Q, int DL, int DE, int DQ>
-struct TypeTraits<Qty<L,E,Q,DL,DE,DQ> >
+template <typename L, typename E, typename Q>
+struct TypeTraits<Qty<L,E,Q>>
{
/** Enum for dimensions*/
enum { hasDimension = true };
/// Type switch set to dimensioned type.
typedef DimensionT DimType;
/// Base unit
- static constexpr Qty<L,E,Q,DL,DE,DQ> baseunit()
- { return Qty<L,E,Q,DL,DE,DQ>::baseunit(); }
+ static constexpr Qty<L,E,Q> baseunit()
+ { return Qty<L,E,Q>::baseunit(); }
};
/** Type traits for alternative code generation*/
-template <int DL, int DE, int DQ>
-struct TypeTraits<Qty<0,0,0,DL,DE,DQ> >
+template <>
+struct TypeTraits<QtyDouble>
{
/** Enum for dimensions*/
enum { hasDimension = false };
/// Type switch set to standard type.
typedef StandardT DimType;
/// Base unit
- static constexpr double baseunit() { return 1.0; }
+ static constexpr QtyDouble baseunit() { return 1.0; }
};
//@}
/** @endcond */
}
#endif
diff --git a/Config/PhysicalQtyComplex.h b/Config/PhysicalQtyComplex.h
--- a/Config/PhysicalQtyComplex.h
+++ b/Config/PhysicalQtyComplex.h
@@ -1,330 +1,284 @@
// -*- C++ -*-
//
// PhysicalQtyComplex.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Physical_Qty_Complex_H
#define Physical_Qty_Complex_H
+#include "PhysicalQty.h"
+#include "PhysicalQtyOps.h"
#include <complex>
/** @file PhysicalQtyComplex.h
* Overloads for operations on complex physical quantities.
*/
namespace std {
/**
* Template specialization for std::complex<Qty<0,0,0> >
* with conversions to complex<double>
*/
- template<int DL, int DE, int DQ>
- class complex<ThePEG::Qty<0,0,0,DL,DE,DQ> >
+ template<>
+ class complex<ThePEG::QtyDouble>
{
public:
/// Default constructor
constexpr complex(double r=0.0, double i=0.0)
: rawValue_(r,i) {}
/// Constructor from complex<double>
constexpr complex(complex<double> C)
: rawValue_(C) {}
/**
* The internal representation of the dimensionful quantity.
* Using this will break dimension-consistency.
*/
constexpr complex<double> rawValue() const { return rawValue_; }
/// Real part
constexpr double real() const { return rawValue_.real(); }
/// Imaginary part
constexpr double imag() const { return rawValue_.imag(); }
/// Cast to complex<double>
constexpr operator complex<double>() const {
return rawValue_;
}
/// Addition-assignment
- complex<ThePEG::Qty<0,0,0,DL,DE,DQ> > &
- operator+=(const complex<ThePEG::Qty<0,0,0,DL,DE,DQ> > x) {
+ complex<ThePEG::QtyDouble> &
+ operator+=(const complex<ThePEG::QtyDouble> x) {
rawValue_ += x.rawValue();
return *this;
}
/// Subtraction-assignment
- complex<ThePEG::Qty<0,0,0,DL,DE,DQ> > &
- operator-=(const complex<ThePEG::Qty<0,0,0,DL,DE,DQ> > x) {
+ complex<ThePEG::QtyDouble> &
+ operator-=(const complex<ThePEG::QtyDouble> x) {
rawValue_ -= x.rawValue();
return *this;
}
private:
/// Internal value of the dimensioned quantity
complex<double> rawValue_;
};
}
// =========================================
namespace ThePEG {
/// @name Overloads for mathematical operations
//@{
// complex qty = complex qty * complex qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator*(std::complex<Qty<L1,E1,Q1,DL1,DE1,DQ1> > q1,
- std::complex<Qty<L2,E2,Q2,DL2,DE2,DQ2> > q2) {
- typedef std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > RetT;
- return RetT(q1.real()*q2.real() - q1.imag()*q2.imag(),
- q1.real()*q2.imag() + q1.imag()*q2.real());
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline constexpr auto
+operator*(std::complex<Qty<L1,E1,Q1>> q1,
+ std::complex<Qty<L2,E2,Q2>> q2)
+-> std::complex<decltype(q1.real()*q2.real())>
+{
+ return {q1.real()*q2.real() - q1.imag()*q2.imag(),
+ q1.real()*q2.imag() + q1.imag()*q2.real()};
}
// complex qty = complex qty * complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> >
-operator*(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- typedef std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> > RetT;
- return RetT(q1.real()*q2.real() - q1.imag()*q2.imag(),
- q1.real()*q2.imag() + q1.imag()*q2.real());
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<typename Qty<L,E,Q>::Squared>
+operator*(std::complex<Qty<L,E,Q>> q1,
+ std::complex<Qty<L,E,Q>> q2)
+{
+ return {q1.real()*q2.real() - q1.imag()*q2.imag(),
+ q1.real()*q2.imag() + q1.imag()*q2.real()};
}
// complex qty = complex double - complex qty
-template<int DL, int DE, int DQ>
inline constexpr std::complex<double>
-operator-(std::complex<double> q1,
- std::complex<Qty<0,0,0,DL,DE,DQ> > q2) {
- typedef std::complex<double> RetT;
- return RetT(q1.real()-q2.real(),q1.imag()-q2.imag());
+operator-(std::complex<double> q1, std::complex<QtyDouble> q2) {
+ return {q1.real()-q2.real(), q1.imag()-q2.imag()};
}
// complex qty = complex double * complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator*(std::complex<double> q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- typedef std::complex<Qty<L,E,Q,DL,DE,DQ> > RetT;
- return RetT(q1.real()*q2.real() - q1.imag()*q2.imag(),
- q1.real()*q2.imag() + q1.imag()*q2.real());
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>>
+operator*(std::complex<double> q1, std::complex<Qty<L,E,Q>> q2) {
+ return {q1.real()*q2.real() - q1.imag()*q2.imag(),
+ q1.real()*q2.imag() + q1.imag()*q2.real()};
}
// complex qty = complex double / complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline std::complex<Qty<-L,-E,-Q,DL,DE,DQ> >
-operator/(std::complex<double> q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- typedef std::complex<Qty<-L,-E,-Q,DL,DE,DQ> > RetT;
- std::complex<Qty<L,E,Q,DL,DE,DQ> > tmp = q1*conj(q2);
- Qty<2*L,2*E,2*Q,DL,DE,DQ> norm = (q2*conj(q2)).real();
- return RetT(tmp.real()/norm,tmp.imag()/norm);
+template<typename L, typename E, typename Q>
+inline std::complex<typename Qty<L,E,Q>::Inverse>
+operator/(std::complex<double> q1, std::complex<Qty<L,E,Q>> q2) {
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex qty = complex double / qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<-L,-E,-Q,DL,DE,DQ> >
-operator/(std::complex<double> q1,
- Qty<L,E,Q,DL,DE,DQ> q2) {
- typedef std::complex<Qty<-L,-E,-Q,DL,DE,DQ> > RetT;
- return RetT(q1.real()/q2,q1.imag()/q2);
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<typename Qty<L,E,Q>::Inverse>
+operator/(std::complex<double> q1, Qty<L,E,Q> q2) {
+ return {q1.real()/q2, q1.imag()/q2};
}
// complex qty = complex qty / complex double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator/(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- std::complex<double> q2) {
- std::complex<Qty<L,E,Q,DL,DE,DQ> > tmp = q1*conj(q2);
- double norm = (q2*conj(q2)).real();
- return std::complex<Qty<L,E,Q,DL,DE,DQ> >(tmp.real()/norm,tmp.imag()/norm);
+template<typename L, typename E, typename Q>
+inline std::complex<Qty<L,E,Q>>
+operator/(std::complex<Qty<L,E,Q>> q1, std::complex<double> q2) {
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex qty = qty / complex double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator/(Qty<L,E,Q,DL,DE,DQ> q1,
- std::complex<double> q2) {
- std::complex<Qty<L,E,Q,DL,DE,DQ> > tmp = q1*conj(q2);
- double norm = (q2*conj(q2)).real();
- return std::complex<Qty<L,E,Q,DL,DE,DQ> >(tmp.real()/norm,tmp.imag()/norm);
+template<typename L, typename E, typename Q>
+inline std::complex<Qty<L,E,Q>>
+operator/(Qty<L,E,Q> q1, std::complex<double> q2) {
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex double = complex qty / complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
+template<typename L, typename E, typename Q>
inline std::complex<double>
-operator/(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> > tmp = q1*conj(q2);
- Qty<2*L,2*E,2*Q,DL,DE,DQ> norm = (q2*conj(q2)).real();
- return std::complex<double>(tmp.real()/norm,tmp.imag()/norm);
+operator/(std::complex<Qty<L,E,Q>> q1,
+ std::complex<Qty<L,E,Q>> q2) {
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex double = qty / complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
+template<typename L, typename E, typename Q>
inline std::complex<double>
-operator/(Qty<L,E,Q,DL,DE,DQ> q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> > tmp = q1*conj(q2);
- Qty<2*L,2*E,2*Q,DL,DE,DQ> norm = (q2*conj(q2)).real();
- return std::complex<double>(tmp.real()/norm,tmp.imag()/norm);
+operator/(Qty<L,E,Q> q1, std::complex<Qty<L,E,Q>> q2) {
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex double = complex qty / qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
+template<typename L, typename E, typename Q>
inline constexpr std::complex<double>
-operator/(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- Qty<L,E,Q,DL,DE,DQ> q2) {
- return std::complex<double>(q1.real()/q2,q1.imag()/q2);
+operator/(std::complex<Qty<L,E,Q>> q1, Qty<L,E,Q> q2) {
+ return {q1.real()/q2, q1.imag()/q2};
}
// complex qty = complex qty / complex qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator/(std::complex<Qty<L1,E1,Q1,DL1,DE1,DQ1> > q1,
- std::complex<Qty<L2,E2,Q2,DL2,DE2,DQ2> > q2) {
- typedef std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > RetT;
- std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > tmp = q1*conj(q2);
- Qty<2*L2,2*E2,2*Q2,DL2,DE2,DQ2> norm = (q2*conj(q2)).real();
- return RetT(tmp.real()/norm,tmp.imag()/norm);
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline auto
+operator/(std::complex<Qty<L1,E1,Q1>> q1,
+ std::complex<Qty<L2,E2,Q2>> q2)
+-> std::complex<decltype(q1.real()/q2.real())>
+{
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex qty = qty / complex qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator/(Qty<L1,E1,Q1,DL1,DE1,DQ1> q1,
- std::complex<Qty<L2,E2,Q2,DL2,DE2,DQ2> > q2) {
- typedef std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > RetT;
- std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > tmp = q1*conj(q2);
- Qty<2*L2,2*E2,2*Q2,DL2,DE2,DQ2> norm = (q2*conj(q2)).real();
- return RetT(tmp.real()/norm,tmp.imag()/norm);
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline auto
+operator/(Qty<L1,E1,Q1> q1,
+ std::complex<Qty<L2,E2,Q2>> q2)
+-> std::complex<decltype(q1/q2.real())>
+{
+ auto tmp = q1*conj(q2);
+ auto norm = (q2*conj(q2)).real();
+ return {tmp.real()/norm, tmp.imag()/norm};
}
// complex qty = complex qty / qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator/(std::complex<Qty<L1,E1,Q1,DL1,DE1,DQ1> > q1,
- Qty<L2,E2,Q2,DL2,DE2,DQ2> q2) {
- typedef std::complex<Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > RetT;
- return RetT(q1.real()/q2,q1.imag()/q2);
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline constexpr auto
+operator/(std::complex<Qty<L1,E1,Q1>> q1, Qty<L2,E2,Q2> q2)
+-> std::complex<decltype(q1.real()/q2)>
+{
+ return {q1.real()/q2, q1.imag()/q2};
}
// complex qty = complex qty * complex double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator*(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- std::complex<double> q2) {
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>>
+operator*(std::complex<Qty<L,E,Q>> q1, std::complex<double> q2) {
return q2 * q1;
}
// complex qty = qty * complex qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator*(Qty<L1,E1,Q1,DL1,DE1,DQ1> q1,
- std::complex<Qty<L2,E2,Q2,DL2,DE2,DQ2> > q2) {
- typedef std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> > RetT;
- return RetT(q1*q2.real(), q1*q2.imag());
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline constexpr auto
+operator*(Qty<L1,E1,Q1> q1, std::complex<Qty<L2,E2,Q2>> q2)
+-> std::complex<decltype(q1*q2.real())>
+{
+ return {q1*q2.real(), q1*q2.imag()};
}
// complex qty = qty * complex qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> >
-operator*(Qty<L,E,Q,DL,DE,DQ> q1,
- std::complex<Qty<L,E,Q,DL,DE,DQ> > q2) {
- typedef std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> > RetT;
- return RetT(q1*q2.real(), q1*q2.imag());
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<typename Qty<L,E,Q>::Squared>
+operator*(Qty<L,E,Q> q1, std::complex<Qty<L,E,Q>> q2) {
+ return {q1*q2.real(), q1*q2.imag()};
}
// complex qty = qty * complex double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator*(Qty<L,E,Q,DL,DE,DQ> q1,
- std::complex<double> q2) {
- typedef std::complex<Qty<L,E,Q,DL,DE,DQ> > RetT;
- return RetT(q1*q2.real(), q1*q2.imag());
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>>
+operator*(Qty<L,E,Q> q1, std::complex<double> q2) {
+ return {q1*q2.real(), q1*q2.imag()};
}
// complex qty = complex double * qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<L,E,Q,DL,DE,DQ> >
-operator*(std::complex<double> q1,
- Qty<L,E,Q,DL,DE,DQ> q2) {
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>>
+operator*(std::complex<double> q1, Qty<L,E,Q> q2) {
return q2 * q1;
}
// complex qty = complex qty * qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr std::complex<Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,
- DL1*DL2,DE1*DE2,DQ1*DQ2> >
-operator*(std::complex<Qty<L1,E1,Q1,DL1,DE1,DQ1> > q1,
- Qty<L2,E2,Q2,DL2,DE2,DQ2> q2) {
+template<typename L1, typename E1, typename Q1,
+ typename L2, typename E2, typename Q2>
+inline constexpr auto
+operator*(std::complex<Qty<L1,E1,Q1>> q1, Qty<L2,E2,Q2> q2)
+-> decltype(q2*q1)
+{
return q2 * q1;
}
// complex qty = complex qty * qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr std::complex<Qty<2*L,2*E,2*Q,DL,DE,DQ> >
-operator*(std::complex<Qty<L,E,Q,DL,DE,DQ> > q1,
- Qty<L,E,Q,DL,DE,DQ> q2) {
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<typename Qty<L,E,Q>::Squared>
+operator*(std::complex<Qty<L,E,Q>> q1, Qty<L,E,Q> q2) {
return q2 * q1;
}
-// // complex qty *= complex double
-// template<int L, int E, int Q, int DL, int DE, int DQ>
-// inline std::complex<Qty<L,E,Q,DL,DE,DQ> > &
-// operator*=(std::complex<Qty<L,E,Q,DL,DE,DQ> > & q1,
-// std::complex<double> q2) {
-// q1 = q1 * q2;
-// return q1;
-// }
-
// complex qty *= double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline std::complex<Qty<L,E,Q,DL,DE,DQ> > &
-operator*=(std::complex<Qty<L,E,Q,DL,DE,DQ> > & q1,
- double q2) {
- q1 = q1 * q2;
- return q1;
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>> &
+operator*=(std::complex<Qty<L,E,Q>> & q1, double q2) {
+ return (q1 = q1 * q2);
}
-// // complex qty /= complex double
-// template<int L, int E, int Q, int DL, int DE, int DQ>
-// inline std::complex<Qty<L,E,Q,DL,DE,DQ> > &
-// operator/=(std::complex<Qty<L,E,Q,DL,DE,DQ> > & q1,
-// std::complex<double> q2) {
-// q1 = q1 / q2;
-// return q1;
-// }
-
// complex qty /= double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline std::complex<Qty<L,E,Q,DL,DE,DQ> > &
-operator/=(std::complex<Qty<L,E,Q,DL,DE,DQ> > & q1,
- double q2) {
- q1 = q1 / q2;
- return q1;
+template<typename L, typename E, typename Q>
+inline constexpr std::complex<Qty<L,E,Q>> &
+operator/=(std::complex<Qty<L,E,Q>> & q1, double q2) {
+ return (q1 = q1 / q2);
}
//@}
}
#endif
diff --git a/Config/PhysicalQtyOps.h b/Config/PhysicalQtyOps.h
--- a/Config/PhysicalQtyOps.h
+++ b/Config/PhysicalQtyOps.h
@@ -1,276 +1,226 @@
// -*- C++ -*-
//
// PhysicalQtyOps.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Physical_Qty_Ops_H
#define Physical_Qty_Ops_H
+#include "PhysicalQty.h"
#include <cmath>
/** @file PhysicalQtyOps.h
* Overloads for mathematical operations on physical quantities.
*/
namespace ThePEG {
/// @name Overloads for mathematical operations on physical quantities.
//@{
// qty = qty * qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2>
-operator*(Qty<L1,E1,Q1,DL1,DE1,DQ1> q1, Qty<L2,E2,Q2,DL2,DE2,DQ2> q2) {
- typedef
- Qty<L1*DL2+L2*DL1,E1*DE2+E2*DE1,Q1*DQ2+Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2> RetT;
+template<typename T, typename U>
+inline constexpr typename BinaryOpTraits<T,U>::MulT
+operator*(T q1, U q2) {
+ typedef typename BinaryOpTraits<T,U>::MulT RetT;
return RetT{RetT::baseunit(), q1.rawValue()*q2.rawValue()};
}
-
// qty = qty / qty
-template<int L1, int L2, int E1, int E2, int Q1, int Q2,
- int DL1, int DL2, int DE1, int DE2, int DQ1, int DQ2>
-inline constexpr Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2>
-operator/(Qty<L1,E1,Q1,DL1,DE1,DQ1> q1, Qty<L2,E2,Q2,DL2,DE2,DQ2> q2) {
- typedef
- Qty<L1*DL2-L2*DL1,E1*DE2-E2*DE1,Q1*DQ2-Q2*DQ1,DL1*DL2,DE1*DE2,DQ1*DQ2> RetT;
+template<typename T, typename U>
+inline constexpr typename BinaryOpTraits<T,U>::DivT
+operator/(T q1, U q2) {
+ typedef typename BinaryOpTraits<T,U>::DivT RetT;
return RetT{RetT::baseunit(), q1.rawValue()/q2.rawValue()};
}
// qty = qty + qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline Qty<L,E,Q,DL,DE,DQ>
-operator+(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
- Qty<L,E,Q,DL,DE,DQ> q = q1;
- q += q2;
- return q;
+template<typename T, typename U>
+inline enable_if_same_qty<T,T,U>
+operator+(T q1, U q2) {
+ q1 += q2;
+ return q1;
}
// qty = qty - qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline Qty<L,E,Q,DL,DE,DQ>
-operator-(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
- Qty<L,E,Q,DL,DE,DQ> q = q1;
- q -= q2;
- return q;
+template<typename T, typename U>
+inline enable_if_same_qty<T,T,U>
+operator-(T q1, U q2) {
+ q1 -= q2;
+ return q1;
}
// qty == qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator==(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator==(T q1, U q2) {
return q1.rawValue()==q2.rawValue();
}
// qty != qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator!=(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator!=(T q1, U q2) {
return q1.rawValue()!=q2.rawValue();
}
// qty < qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator<(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator<(T q1, U q2) {
return q1.rawValue()<q2.rawValue();
}
// qty <= qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator<=(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator<=(T q1, U q2) {
return q1.rawValue()<=q2.rawValue();
}
// qty > qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator>(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator>(T q1, U q2) {
return q1.rawValue()>q2.rawValue();
}
// qty >= qty
-template<int L, int E, int Q, int DL, int DE, int DQ, int DL2, int DE2, int DQ2>
-inline constexpr bool
-operator>=(Qty<L,E,Q,DL,DE,DQ> q1,
- Qty<QtyInt<L,DL,DL2>::I,
- QtyInt<E,DE,DE2>::I,
- QtyInt<Q,DQ,DQ2>::I,
- DL2,DE2,DQ2> q2) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<bool,T,U>
+operator>=(T q1, U q2) {
return q1.rawValue()>=q2.rawValue();
}
// comparisons with ZERO
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator==(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator==(T q1, ZeroUnit) {
return q1.rawValue() == 0.0;
}
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator!=(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator!=(T q1, ZeroUnit) {
return q1.rawValue() != 0.0;
}
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator<(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator<(T q1, ZeroUnit) {
return q1.rawValue() < 0.0;
}
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator>(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator>(T q1, ZeroUnit) {
return q1.rawValue() > 0.0;
}
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator<=(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator<=(T q1, ZeroUnit) {
return q1.rawValue() <= 0.0;
}
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr bool
-operator>=(Qty<L,E,Q,DL,DE,DQ> q1, ZeroUnit) {
+template<typename T>
+inline constexpr enable_if_same_qty<bool, T>
+operator>=(T q1, ZeroUnit) {
return q1.rawValue() >= 0.0;
}
// qty = qty * double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<L,E,Q,DL,DE,DQ>
-operator*(Qty<L,E,Q,DL,DE,DQ> q,double x) {
- return Qty<L,E,Q,DL,DE,DQ>{q,x};
+template<typename T>
+inline constexpr enable_if_same_qty<T, T>
+operator*(T q,double x) {
+ return T{q,x};
}
// qty = double * qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<L,E,Q,DL,DE,DQ>
-operator*(double x,Qty<L,E,Q,DL,DE,DQ> q) {
- return Qty<L,E,Q,DL,DE,DQ>{q,x};
+template<typename T>
+inline constexpr enable_if_same_qty<T, T>
+operator*(double x,T q) {
+ return T{q,x};
}
// qty = qty / double
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<L,E,Q,DL,DE,DQ>
-operator/(Qty<L,E,Q,DL,DE,DQ> q,double x) {
- return Qty<L,E,Q,DL,DE,DQ>{q, 1./x};
+template<typename T>
+inline constexpr enable_if_same_qty<T, T>
+operator/(T q,double x) {
+ return T{q, 1./x};
}
// qty = double / qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<-L,-E,-Q,DL,DE,DQ>
-operator/(double x, Qty<L,E,Q,DL,DE,DQ> q) {
- typedef Qty<-L,-E,-Q,DL,DE,DQ> RetT;
+template<typename T>
+inline constexpr enable_if_same_qty<typename T::Inverse, T>
+operator/(double x, T q) {
+ typedef typename T::Inverse RetT;
return RetT{RetT::baseunit(), x/q.rawValue()};
}
// qty = -qty
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<L,E,Q,DL,DE,DQ>
-operator-(Qty<L,E,Q,DL,DE,DQ> q) {
- typedef Qty<L,E,Q,DL,DE,DQ> RetT;
- return RetT{RetT::baseunit(), -q.rawValue()};
-}
-
-// qty = sqr(qty)
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<2*L,2*E,2*Q,DL,DE,DQ>
-sqr(Qty<L,E,Q,DL,DE,DQ> q) {
- return q*q;
+template<typename T>
+inline constexpr enable_if_same_qty<T, T>
+operator-(T q) {
+ typedef T RetT;
+ return RetT{q, -1.0};
}
// qty = sqrt(qty) // std::sqrt is not constexpr
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline Qty<L,E,Q,DL*2,DE*2,DQ*2>
-sqrt(Qty<L,E,Q,DL,DE,DQ> q) {
- typedef Qty<L,E,Q,DL*2,DE*2,DQ*2> RetT;
- return RetT(std::sqrt(q.rawValue())*RetT::baseunit());
+template<typename T>
+inline enable_if_same_qty<typename T::Sqrt, T>
+sqrt(T q) {
+ typedef typename T::Sqrt RetT;
+ return RetT{RetT::baseunit(), std::sqrt(q.rawValue())};
}
// double = atan2(y,x)
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr double
-atan2(Qty<L,E,Q,DL,DE,DQ> y, Qty<L,E,Q,DL,DE,DQ> x) {
+template<typename T, typename U>
+inline constexpr enable_if_same_qty<double,T,U>
+atan2(T y, U x) {
return std::atan2(y.rawValue(), x.rawValue());
}
// qty = abs(qty)
-template<int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<L,E,Q,DL,DE,DQ>
-abs(Qty<L,E,Q,DL,DE,DQ> q) {
- typedef Qty<L,E,Q,DL,DE,DQ> RetT;
- return RetT(std::abs(q.rawValue())*RetT::baseunit());
+template<typename T>
+inline constexpr enable_if_same_qty<T, T>
+abs(T q) {
+ return T{T::baseunit(), std::abs(q.rawValue())};
}
// qty = pow<P,R>(qty)
-template<int P, int R, int L, int E, int Q, int DL, int DE, int DQ>
-inline constexpr Qty<P*L,P*E,P*Q,R*DL,R*DE,R*DQ>
-pow(Qty<L,E,Q,DL,DE,DQ> q) {
- typedef Qty<P*L,P*E,P*Q,R*DL,R*DE,R*DQ> RetT;
- return RetT(std::pow(q.rawValue(),double(P)/double(R))*RetT::baseunit());
+template<int Num, int Den, typename T>
+inline constexpr enable_if_same_qty<typename T::template Power<Num,Den>, T>
+pow(T q) {
+ typedef typename T::template Power<Num,Den> RetT;
+ return RetT{RetT::baseunit(), std::pow(q.rawValue(),double(Num)/double(Den))};
}
// max for T,U types
template<typename T, typename U>
-inline
-T max(const T & t, const U & u) {
+inline T max(const T & t, const U & u) {
const T & utmp = u;
return std::max(t, utmp);
}
// ZeroUnit in front should take U type
template<typename U>
-inline
-U max(const ZeroUnit & t, const U & u) {
+inline U max(const ZeroUnit & t, const U & u) {
const U & ttmp = t;
return std::max(ttmp, u);
}
// min for T,U types
template<typename T, typename U>
-inline
-T min(const T & t, const U & u) {
+inline T min(const T & t, const U & u) {
const T & utmp = u;
return std::min(t, utmp);
}
// ZeroUnit in front should take U type
template<typename U>
-inline
-U min(const ZeroUnit & t, const U & u) {
+inline U min(const ZeroUnit & t, const U & u) {
const U & ttmp = t;
return std::min(ttmp, u);
}
//@}
}
#endif
diff --git a/Config/TemplateTools.h b/Config/TemplateTools.h
--- a/Config/TemplateTools.h
+++ b/Config/TemplateTools.h
@@ -1,97 +1,64 @@
// -*- C++ -*-
//
// TemplateTools.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Template_Tools_H
#define Template_Tools_H
#include <type_traits>
/**
* @file TemplateTools.h
* Useful template machinery. Based on Alexandrescu, "Modern C++ Design".
*/
namespace ThePEG {
/// Conversion between integers and types.
template <int v>
struct Int2Type
{
enum { value = v };
};
/// Dummy type for ambiguous function signatures.
struct DummyType {};
-/// Result type calculations for binary operators.
-template <typename T, typename U>
-struct BinaryOpTraits;
-
-/** @cond TRAITSPECIALIZATIONS */
-
-template <>
-struct BinaryOpTraits<double,double> {
- /** The type resulting from multiplication of the template types. */
- typedef double MulT;
- /** The type resulting from division of the first template type by
- the second. */
- typedef double DivT;
-};
-
-template <>
-struct BinaryOpTraits<long double, long double> {
- /** The type resulting from multiplication of the template types. */
- typedef long double MulT;
- /** The type resulting from division of the first template type by
- the second. */
- typedef long double DivT;
-};
-
-template <>
-struct BinaryOpTraits<int,int> {
- /** The type resulting from multiplication of the template types. */
- typedef int MulT;
- /** The type resulting from division of the first template type by
- the second. */
- typedef int DivT;
-};
-
/** @endcond */
/// Selection mechanism for type-dependent implementations.
enum ImplSelector { Dimensioned, Standard };
/// Typedef for dimensioned types.
typedef Int2Type<Dimensioned> DimensionT;
/// Typedef for non-dimensioned types.
typedef Int2Type<Standard> StandardT;
/// Type traits for built-in types
template <typename T>
struct TypeTraits
{
/// Boolean flag. Is true for physical quantities.
enum { hasDimension = false };
/// Implementation selector
typedef StandardT DimType;
/// Base unit for arithmetic types
// construction with extra U type is necessary to make
// enable_if work before concepts are supported properly
template <typename U = T>
static constexpr typename
std::enable_if< (std::is_arithmetic<U>::value &&
std::is_same<U, T>::value), U>::type
baseunit()
{ return static_cast<U>(1); }
};
}
#endif
diff --git a/Config/ThePEG.h b/Config/ThePEG.h
--- a/Config/ThePEG.h
+++ b/Config/ThePEG.h
@@ -1,149 +1,148 @@
// -*- C++ -*-
//
// ThePEG.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_H
#define ThePEG_H
/** \file ThePEG.h
* This is the main config header file for ThePEG. Do not make
* changes in this file. If you need to modify anything, edit a copy
* of the file which can be included instead of this file using the
* macro <code>ThePEG_ALTERNATE_CONFIG</code>.
*/
#ifndef ThePEG_ALTERNATE_CONFIG
#include "ThePEG/Pointer/Ptr.h"
#include "ThePEG/Pointer/PtrTraits.h"
#include "ThePEG/Pointer/RCPtr.h"
#include "ThePEG/Utilities/Rebinder.fh"
#include "ThePEG/Utilities/Interval.fh"
#include "ThePEG/Utilities/ClassDescription.fh"
#include "ThePEG/Interface/InterfaceBase.fh"
#include "ThePEG/Persistency/PersistentOStream.fh"
#include "ThePEG/Persistency/PersistentIStream.fh"
#include "TemplateTools.h"
#include "Complex.h"
#include "Unitsystem.h"
#include "Constants.h"
#include "std.h"
/**
* This is the main namespace within which all identifiers in ThePEG
* are declared. External packages based on ThePEG should not
* introduce identifiers in the ThePEG namespace, but in a separate
* namespace which need not be nested within the ThePEG namespace.
*/
namespace ThePEG {
// Introduce some identifiers in the ThePEG namespace/
using namespace ThePEG::Pointer;
using ThePEG::Pointer::Ptr;
using namespace ThePEG::Units;
/**
* Define the base class from which all (polymorphic) classes in
* ThePEG are derived.
*/
struct Base: public ReferenceCounted {
/** The virtual destructor */
virtual ~Base() {}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init() {}
/**
* Print out debugging information for this object on std::cerr. To
* be called from within a debugger. Simply calls the virtual
* debugme() function.
*/
void debug() const;
/**
* Print out debugging information for this object on std::cerr. To
* be called from within a debugger via the debug() function.
*/
virtual void debugme() const;
};
/**
* Define the base class from which all persistent classes in
* ThePEG are derived.
*/
typedef Base PersistentBase;
/**
* TraitsType is an empty, non-polymorphic, base class. It is used as
* a base class of all traits classes in ThePEG in order to group them
* together in the documentation. It currently serves no other
* purpose.
*/
struct TraitsType {};
/**
* A standard exception class to be used for vetoing a whole event.
*/
struct Veto {
/** the default constructor. */
Veto();
};
/**
* A standard exception class to be used to temporarily stop the
* generation of an event.
*/
struct Stop {};
/**
* The square function should really have been included in the
* standard C++ library.
*/
template <typename T>
-inline constexpr
-typename BinaryOpTraits<T,T>::MulT
-sqr(const T& x) {
+inline constexpr
+auto sqr(const T& x) -> decltype(x*x) {
return x*x;
}
// Debugging in ThePEG may be swithced off completely by this
// compilation swithc, eliminating possible overhead in error
// checking.
#ifndef ThePEG_NO_DEBUG
/** Macro for accessing debug functions to enable compile-time disabling. */
#define ThePEG_DEBUG_LEVEL Debug::level
/** Macro for accessing debug functions to enable compile-time disabling. */
#define ThePEG_DEBUG_ITEM(item) Debug::debugItem(item)
#else
/** Macro for accessing debug functions to enable compile-time disabling. */
#define ThePEG_DEBUG_LEVEL 0
/** Macro for accessing debug functions to enable compile-time disabling. */
#define ThePEG_DEBUG_ITEM(item) false
#endif
}
#include "Pointers.h"
#include "Containers.h"
#else
#include ThePEG_ALTERNATIVE_CONFIG
#endif
#endif /* ThePEG_H */
diff --git a/Config/Unitsystem.h b/Config/Unitsystem.h
--- a/Config/Unitsystem.h
+++ b/Config/Unitsystem.h
@@ -1,264 +1,268 @@
// -*- C++ -*-
//
// Unitsystem.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad, David Grellscheid
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Units_H
#define ThePEG_Units_H
#include "ThePEG/Vectors/Lorentz5Vector.fh"
#include "ThePEG/Vectors/LorentzVector.fh"
#include "ThePEG/Vectors/ThreeVector.fh"
#include "ThePEG/Vectors/Transverse.fh"
#include "PhysicalQty.h"
#include "PhysicalQtyOps.h"
#include "PhysicalQtyComplex.h"
namespace ThePEG {
/**
* The Units namespace contains the declaration of a number of classes
* for variables with dimension. Currently they are all typedefs of
* double, but in the future the SIUnits package will be used.
*
* The file Utilities/UnitIO.h defines helper-classes and helper
* functions to read and write variables with dimensions. As an
* example, to read and write an energy variable <code>e</code> in
* units of GeV, use: <code>os << ounit(e, GeV)</code> and <code>is >>
* iunit(e, GeV)</code>
*/
namespace Units {
+/// adapter for the old style of naming quantities
+template<int L, int E, int Q, int DL=1, int DE=1, int DQ=1>
+using Qty = ThePEG::Qty<std::ratio<L,DL>, std::ratio<E,DE>, std::ratio<Q,DQ>>;
+
/** Energy. */
typedef Qty<0,1,0> Energy;
/** Mass has the same unit as Energy <=> c == 1. */
typedef Energy Mass;
/** Length. */
typedef Qty<1,0,0> Length;
/** Time has the same unit as Length. <=> c == 1. */
typedef Length Time;
/** Inverse Length. */
typedef Qty<-1,0,0> InvLength;
/** Velocities are dimensionless fractions of c. */
typedef double Velocity;
/** Charge. */
typedef Qty<0,0,1> Charge;
/** Angular momentum. */
typedef Qty<1,1,0> AngularMomentum;
/** Tension. */
typedef Qty<-1,1,0> Tension;
/** Area will be assumed to be Length\f$^2\f$. */
typedef Qty<2,0,0> Area;
/** Inverse Area. */
typedef Qty<-2,0,0> InvArea;
/** Cross section is an area. */
typedef Area CrossSection;
/**
* @name Higher powers of energy.
* Even higher powers can be created with similar typedefs.
*/
//@{
typedef Qty<0, 2, 0> Energy2;
typedef Qty<0, 3, 0> Energy3;
typedef Qty<0, 4, 0> Energy4;
typedef Qty<0, 5, 0> Energy5;
typedef Qty<0, 6, 0> Energy6;
typedef Qty<0, 7, 0> Energy7;
typedef Qty<0, 8, 0> Energy8;
typedef Qty<0, 9, 0> Energy9;
typedef Qty<0,10, 0> Energy10;
typedef Qty<0,11, 0> Energy11;
typedef Qty<0,12, 0> Energy12;
typedef Qty<0, 1,0, 1,2,1> SqrtEnergy;
typedef Qty<0,-1,0, 1,2,1> InvSqrtEnergy;
typedef Qty<0, -1, 0> InvEnergy;
typedef Qty<0, -2, 0> InvEnergy2;
typedef Qty<0, -3, 0> InvEnergy3;
typedef Qty<0, -4, 0> InvEnergy4;
typedef Qty<0, -5, 0> InvEnergy5;
typedef Qty<0, -6, 0> InvEnergy6;
typedef Qty<0, -7, 0> InvEnergy7;
typedef Qty<0, -8, 0> InvEnergy8;
typedef Qty<0, -9, 0> InvEnergy9;
typedef Qty<0,-10, 0> InvEnergy10;
typedef Qty<0,-11, 0> InvEnergy11;
typedef Qty<0,-12, 0> InvEnergy12;
//@}
/** CrossSection*Energy2. */
typedef Qty<2,2,0> Energy2XSec;
/** CrossSection/Energy2. */
typedef Qty<2,-2,0> DiffXSec;
/** CrossSection/Energy4. */
typedef Qty<2,-4,0> Diff2XSec;
/** CrossSection/Energy6 */
typedef Qty<2,-6,0> Diff3XSec;
/** Scale is the same as a squared energy. */
typedef Energy2 Scale;
/** A point in three-dimensional euclidean space. */
typedef ThreeVector<Length> Point;
/** A distance in three-dimensional euclidean space. */
typedef ThreeVector<Length> Distance;
/** A direction in three-dimensional euclidean space. */
typedef ThreeVector<double> Axis;
/** A momentum in three-dimensional euclidean space. */
typedef ThreeVector<Energy> Momentum3;
/** A three-dimensional boost vector. */
typedef ThreeVector<double> Boost;
/** A distance in four-dimensional space-time. */
typedef LorentzVector<Length> LorentzDistance;
/** A distance in four-dimensional space-time with an explicit
* invariant time component. */
typedef Lorentz5Vector<Length> Lorentz5Distance;
/** A point in four-dimensional space-time. */
typedef LorentzVector<Length> LorentzPoint;
/** A momentum in four-dimensional space-time. */
typedef LorentzVector<Energy> LorentzMomentum;
/** A momentum in four-dimensional space-time with an explicit
* invariant mass component. */
typedef Lorentz5Vector<Energy> Lorentz5Momentum;
/** Transverse components of a momentum. */
typedef Transverse<Energy> TransverseMomentum;
/// @name Pre-defined basic units.
//@{
constexpr Length operator "" _mm( long double x ) {
return Length{Length::baseunit(), static_cast<double>(x)};
}
constexpr Length operator "" _mm( unsigned long long x ) {
return Length{Length::baseunit(), static_cast<double>(x)};
}
constexpr Length meter = 1.0e+3_mm;
constexpr Length millimeter = 1_mm;
constexpr Length mm = 1_mm;
constexpr Length centimeter = 10_mm;
constexpr Length micrometer = 1.0e-3_mm;
constexpr Length nanometer = 1.0e-6_mm;
constexpr Length picometer = 1.0e-9_mm;
constexpr Length femtometer = 1.0e-12_mm;
constexpr Energy operator "" _MeV( long double x ) {
return Energy{Energy::baseunit(), static_cast<double>(x)};
}
constexpr Energy operator "" _MeV( unsigned long long x ) {
return Energy{Energy::baseunit(), static_cast<double>(x)};
}
constexpr Energy operator "" _GeV( long double x ) {
return Energy{1000_MeV, static_cast<double>(x)};
}
constexpr Energy operator "" _GeV( unsigned long long x ) {
return Energy{1000_MeV, static_cast<double>(x)};
}
constexpr Energy keV = 1.0e-3_MeV;
constexpr Energy MeV = 1_MeV;
constexpr Energy GeV = 1_GeV;
constexpr Energy TeV = 1.0e+6_MeV;
constexpr Energy2 operator "" _MeV2( long double x ) {
return Energy2{Energy2::baseunit(), static_cast<double>(x)};
}
constexpr Energy2 operator "" _MeV2( unsigned long long x ) {
return Energy2{Energy2::baseunit(), static_cast<double>(x)};
}
constexpr Energy2 operator "" _GeV2( long double x ) {
return Energy2{1.0e+6_MeV2, static_cast<double>(x)};
}
constexpr Energy2 operator "" _GeV2( unsigned long long x ) {
return Energy2{1.0e+6_MeV2, static_cast<double>(x)};
}
constexpr Energy2 MeV2 = 1_MeV2;
constexpr Energy2 GeV2 = 1_GeV2;
constexpr InvEnergy InvGeV = 1/GeV;
constexpr Area operator "" _pb( long double x ) {
return Area{1.0e-34 * Area::baseunit(), static_cast<double>(x)};
}
constexpr Area operator "" _pb( unsigned long long x ) {
return Area{1.0e-34 * Area::baseunit(), static_cast<double>(x)};
}
constexpr Area femtobarn = 1.0e-03_pb;
constexpr Area picobarn = 1_pb;
constexpr Area nanobarn = 1.0e+03_pb;
constexpr Area microbarn = 1.0e+06_pb;
constexpr Area millibarn = 1.0e+09_pb;
constexpr Area barn = 1.0e+12_pb;
constexpr Charge eplus = Charge::baseunit();
//@}
/// Planck's constant times c (PDG 2006 value 197.326968(17) MeV fm)
constexpr Qty<1,1,0> hbarc = 197.326968e-15 * MeV * meter;
/// Planck's constant (PDG 2006 value 197.326968(17) MeV fm)
constexpr Qty<1,1,0> hbar_Planck = hbarc / 1.0; // c is one
}
/**
* Use symbols from this namespace to make forced breaks of unit
* consistency explicit.
*/
namespace UnitRemoval {
/// @name Helper units to make breaks of unit consistency explicit.
//@{
constexpr Units::Energy E = Units::Energy::baseunit();
constexpr Units::Energy2 E2 = E*E;
constexpr Units::Energy3 E3 = E*E2;
constexpr Units::Energy4 E4 = E2*E2;
constexpr Units::InvEnergy InvE = 1.0/E;
constexpr Units::InvEnergy2 InvE2 = 1.0/E2;
constexpr Units::InvEnergy3 InvE3 = 1.0/E3;
constexpr Units::InvEnergy4 InvE4 = 1.0/E4;
constexpr Units::SqrtEnergy SqrtE = Units::SqrtEnergy::baseunit();
constexpr Units::InvSqrtEnergy InvSqrtE = Units::InvSqrtEnergy::baseunit();
//@}
}
}
#endif /* ThePEG_Units_H */
diff --git a/Helicity/LorentzRSSpinor.h b/Helicity/LorentzRSSpinor.h
--- a/Helicity/LorentzRSSpinor.h
+++ b/Helicity/LorentzRSSpinor.h
@@ -1,409 +1,413 @@
// -*- C++ -*-
//
// LorentzRSSpinor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzRSSpinor_H
#define ThePEG_LorentzRSSpinor_H
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzRSSpinor.fh"
#include "LorentzRSSpinorBar.h"
#include "LorentzSpinorBar.h"
#include "LorentzSpinor.h"
#include "LorentzPolarizationVector.h"
namespace ThePEG{
namespace Helicity{
/**
* The LorentzRSSpinor class is designed to store a Rarita-Schwinger
* spinor for a spin-3/2 particle. In addition to storing the
* components of the spinor information is stored on the type of
* spinor, for example u or v type.
*
* At the moment only one choice of the Dirac matrix representation
* is supported. For high-energy calculations the choice made by the
* HELAS collaboration is more efficient for numerical
* calculations. In this representation
*
* \f[
* \gamma_{i=1,2,3}=\left(\begin{array}{cc}
* 0 & \sigma_i \\
* -\sigma_i & 0
* \end{array}\right)
* \quad
* \gamma_0=\left(\begin{array}{cc}
* 0 & 1 \\
* 1 & 0
* \end{array}\right)
* \quad
* \gamma_5=\left(\begin{array}{cc}
* -1 & 0 \\
* 0 & 1
* \end{array}\right)
* \f]
*
* The type of the spinor is also stored using the SpinorType
* enumeration. There are three types supported SpinorType::u,
* SpinorType::v, SpinorType::unknown. This information is intended
* mainly for use in the case of Majorana particles where matrix
* elements can be calculated with either u or v type spinors and
* knowledge of which was used will be needed in order to give the
* correct correlations. The SpinorType::unknowne is intended for
* cases where either the spinor for an off-shell line in a matrix
* element calculation or the information is genuinely unknown.
*
* The LorentzRSSpinorBar class is also provided to store the barred
* spinor.
*
*
* @see HelicityDefinitions
* @see LorentzRSSpinorBar
*
* \author Peter Richardson
*
*/
template<typename Value>
class LorentzRSSpinor {
public:
/** @name Standard constructors. */
//@{
/**
* Default zero constructor, optionally specifying \a t, the type.
*/
LorentzRSSpinor(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
/**
* Constructor with complex numbers specifying the components,
* optionally specifying \a t, the type.
*/
LorentzRSSpinor(complex<Value> a1, complex<Value> b1,
complex<Value> c1, complex<Value> d1,
complex<Value> a2, complex<Value> b2,
complex<Value> c2, complex<Value> d2,
complex<Value> a3, complex<Value> b3,
complex<Value> c3, complex<Value> d3,
complex<Value> a4, complex<Value> b4,
complex<Value> c4, complex<Value> d4,
SpinorType t=SpinorType::unknown)
: _type(t), _spin{{ {{a1,b1,c1,d1}},
{{a2,b2,c2,d2}},
{{a3,b3,c3,d3}},
{{a4,b4,c4,d4}}
}} {}
//@}
/** @name Access the components. */
//@{
/**
* Subscript operator to return spinor components
*/
complex<Value> operator()(int i, int j) const {
assert( i >= 0 && i <= 3 && j>=0 && j<=3);
return _spin[i][j];
}
/**
* Set components by index
*/
complex<Value> & operator () (int i, int j) {
assert( i >= 0 && i <= 3 && j>=0 && j<=3);
return _spin[i][j];
}
/**
* Get first spinor component for the x vector
*/
complex<Value> xs1() const {return _spin[0][0];}
/**
* Get second spinor component for the x vector
*/
complex<Value> xs2() const {return _spin[0][1];}
/**
* Get third spinor component for the x vector
*/
complex<Value> xs3() const {return _spin[0][2];}
/**
* Get fourth spinor component for the x vector
*/
complex<Value> xs4() const {return _spin[0][3];}
/**
* Get first spinor component for the y vector
*/
complex<Value> ys1() const {return _spin[1][0];}
/**
* Get second spinor component for the y vector
*/
complex<Value> ys2() const {return _spin[1][1];}
/**
* Get third spinor component for the y vector
*/
complex<Value> ys3() const {return _spin[1][2];}
/**
* Get fourth spinor component for the y vector
*/
complex<Value> ys4() const {return _spin[1][3];}
/**
* Get first spinor component for the z vector
*/
complex<Value> zs1() const {return _spin[2][0];}
/**
* Get second spinor component for the z vector
*/
complex<Value> zs2() const {return _spin[2][1];}
/**
* Get third spinor component for the z vector
*/
complex<Value> zs3() const {return _spin[2][2];}
/**
* Get fourth spinor component for the z vector
*/
complex<Value> zs4() const {return _spin[2][3];}
/**
* Get first spinor component for the t vector
*/
complex<Value> ts1() const {return _spin[3][0];}
/**
* Get second spinor component for the t vector
*/
complex<Value> ts2() const {return _spin[3][1];}
/**
* Get third spinor component for the t vector
*/
complex<Value> ts3() const {return _spin[3][2];}
/**
* Get fourth spinor component for the t vector
*/
complex<Value> ts4() const {return _spin[3][3];}
/**
* Set first spinor component for the x vector
*/
void setXS1(complex<Value> in) {_spin[0][0]=in;}
/**
* Set second spinor component for the x vector
*/
void setXS2(complex<Value> in) {_spin[0][1]=in;}
/**
* Set third spinor component for the x vector
*/
void setXS3(complex<Value> in) {_spin[0][2]=in;}
/**
* Set fourth spinor component for the x vector
*/
void setXS4(complex<Value> in) {_spin[0][3]=in;}
/**
* Set first spinor component for the y vector
*/
void setYS1(complex<Value> in) {_spin[1][0]=in;}
/**
* Set second spinor component for the y vector
*/
void setYS2(complex<Value> in) {_spin[1][1]=in;}
/**
* Set third spinor component for the y vector
*/
void setYS3(complex<Value> in) {_spin[1][2]=in;}
/**
* Set fourth spinor component for the y vector
*/
void setYS4(complex<Value> in) {_spin[1][3]=in;}
/**
* Set first spinor component for the z vector
*/
void setZS1(complex<Value> in) {_spin[2][0]=in;}
/**
* Set second spinor component for the z vector
*/
void setZS2(complex<Value> in) {_spin[2][1]=in;}
/**
* Set third spinor component for the z vector
*/
void setZS3(complex<Value> in) {_spin[2][2]=in;}
/**
* Set fourth spinor component for the z vector
*/
void setZS4(complex<Value> in) {_spin[2][3]=in;}
/**
* Set first spinor component for the t vector
*/
void setTS1(complex<Value> in) {_spin[3][0]=in;}
/**
* Set second spinor component for the t vector
*/
void setTS2(complex<Value> in) {_spin[3][1]=in;}
/**
* Set third spinor component for the t vector
*/
void setTS3(complex<Value> in) {_spin[3][2]=in;}
/**
* Set fourth spinor component for the t vector
*/
void setTS4(complex<Value> in) {_spin[3][3]=in;}
//@}
/** @name Arithmetic operators. */
//@{
/**
* dot product with a polarization vector
*/
LorentzSpinor<Value> dot(const LorentzPolarizationVector & vec) const {
LorentzSpinor<Value> output(_type);
complex<Value> temp;
unsigned int ix;
for(ix=0;ix<4;++ix) {
temp = _spin[3][ix]*vec.t();
temp -= _spin[0][ix]*vec.x();
temp -= _spin[1][ix]*vec.y();
temp -= _spin[2][ix]*vec.z();
output[ix]=temp;
}
return output;
}
/**
* dot product with a 4-vector
*/
LorentzSpinor<Value> dot(const LorentzMomentum & invec) const {
LorentzSpinor<Value> output(_type);
complex<Value> temp;
LorentzVector<double> vec = UnitRemoval::InvE * invec;
unsigned int ix;
for(ix=0;ix<4;++ix) {
temp = - ( _spin[0][ix]*vec.x() + _spin[1][ix]*vec.y()+
_spin[2][ix]*vec.z() ) + _spin[3][ix]*vec.t();
output[ix]=temp;
}
return output;
}
//@}
/** @name Transformations. */
//@{
/**
* return the barred spinor
*/
LorentzRSSpinorBar<Value> bar() const;
/**
* Standard Lorentz boost specifying the components of the beta vector.
*/
LorentzRSSpinor & boost(double,double,double);
/**
* Standard Lorentz boost specifying the beta vector.
*/
LorentzRSSpinor & boost(const Boost &);
/**
* General transform
*/
LorentzRSSpinor & transform(const LorentzRotation &);
//@}
/** @name Functions related to type. */
//@{
/**
* Return the type of the spinor.
*/
SpinorType Type() const {return _type;}
//@}
/**
* Scalar product \f$\bar{f}^\alpha(c_LP_L+c_RP_R)f_\alpha\f$for general couplings
* @param fbar The barred spinor
* @param left The left-handed coupling, \f$c_L\f$.
* @param right The right-handed coupling, \f$c_R\f$.
*/
template <typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- generalScalar(LorentzRSSpinorBar<ValueB>& fbar, Complex left, Complex right) {
- complex<typename BinaryOpTraits<Value,ValueB>::MulT> output;
+ auto generalScalar(LorentzRSSpinorBar<ValueB>& fbar,
+ Complex left, Complex right)
+ -> decltype(left*fbar(3,0)*ts1())
+ {
+ decltype(left*fbar(3,0)*ts1()) output;
unsigned int iz;
output =
left*(fbar(3,0)*_spin[3][0]+fbar(3,1)*_spin[3][1])
+right*(fbar(3,2)*_spin[3][2]+fbar(3,3)*_spin[3][3]);
for(iz=0;iz<3;++iz) {
output -=
left*(fbar(iz,0)*_spin[iz][0]+fbar(iz,1)*_spin[iz][1])
+right*(fbar(iz,2)*_spin[iz][2]+fbar(iz,3)*_spin[iz][3]);
}
return output;
}
/**
* Current \f$\bar{f}(c_LP_L+c_RP_R)f^\alpha\f$ for general couplings.
* @param fbar The barred spinor
* @param left The left-handed coupling, \f$c_L\f$.
* @param right The right-handed coupling, \f$c_R\f$.
*/
template <typename ValueB>
- LorentzVector<complex<typename BinaryOpTraits<Value,ValueB>::MulT> >
- generalCurrent(LorentzSpinorBar<ValueB>& fbar, Complex left, Complex right) {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto generalCurrent(LorentzSpinorBar<ValueB>& fbar,
+ Complex left, Complex right)
+ -> LorentzVector<decltype(left*fbar.s1()*ts1())>
+ {
+ typedef decltype(left*fbar.s1()*ts1()) ResultT;
ResultT output[4];
for(size_t iz=0;iz<4;++iz)
output[iz]= left*(fbar.s1()*_spin[iz][0]+fbar.s2()*_spin[iz][1])
+right*(fbar.s3()*_spin[iz][2]+fbar.s4()*_spin[iz][3]);
return LorentzVector<ResultT>(output[0],output[1],output[2],output[3]);
}
private:
/**
* Type of spinor
*/
SpinorType _type;
/**
* Storage of the components.
*/
std::array<std::array<complex<Value>,4>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzRSSpinor.tcc"
#endif
#endif
diff --git a/Helicity/LorentzRSSpinorBar.h b/Helicity/LorentzRSSpinorBar.h
--- a/Helicity/LorentzRSSpinorBar.h
+++ b/Helicity/LorentzRSSpinorBar.h
@@ -1,346 +1,347 @@
// -*- C++ -*-
//
// LorentzRSSpinorBar.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzRSSpinorBar_H
#define ThePEG_LorentzRSSpinorBar_H
// This is the declaration of the LorentzRSSpinorBar class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzRSSpinor.fh"
#include "LorentzRSSpinorBar.fh"
#include "LorentzSpinorBar.h"
#include "LorentzSpinor.h"
#include "LorentzPolarizationVector.h"
namespace ThePEG {
namespace Helicity {
/**
* The <code>LorentzRSSpinorBar</code> class implements the storage of a
* barred Lorentz Rarita-Schwinger Spinor for a spin-3/2 particle.
* The design is based on that of the
* <code>LorentzRSSpinor</code> class and the details of the implemented
* are discussed in more detail in the header file for that class.
*
* @see LorentzRSSpinor
*
* \author Peter Richardson
*
*/
template <typename Value>
class LorentzRSSpinorBar {
public:
/** @name Standard constructors. */
//@{
/**
* Default zero constructor, optionally specifying \a t, the type.
*/
LorentzRSSpinorBar(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
/**
* Constructor with complex numbers specifying the components,
* optionally specifying \a t, the type.
*/
LorentzRSSpinorBar(complex<Value> a1,complex<Value> b1,
complex<Value> c1,complex<Value> d1,
complex<Value> a2,complex<Value> b2,
complex<Value> c2,complex<Value> d2,
complex<Value> a3,complex<Value> b3,
complex<Value> c3,complex<Value> d3,
complex<Value> a4,complex<Value> b4,
complex<Value> c4,complex<Value> d4,
SpinorType t=SpinorType::unknown)
: _type(t), _spin{{ {{a1,b1,c1,d1}},
{{a2,b2,c2,d2}},
{{a3,b3,c3,d3}},
{{a4,b4,c4,d4}}
}} {}
//@}
/** @name Access the components. */
//@{
/**
* Subscript operator to return spinor components
*/
complex<Value> operator()(int i, int j) const {
assert( i >= 0 && i <= 3 && j>=0 && j<=3 );
return _spin[i][j];
}
/**
* Set components by index
*/
complex<Value> & operator () (int i, int j) {
assert( i >= 0 && i <= 3 && j>=0 && j<=3 );
return _spin[i][j];
}
/**
* Get first spinor component for the x vector
*/
complex<Value> xs1() const {return _spin[0][0];}
/**
* Get second spinor component for the x vector
*/
complex<Value> xs2() const {return _spin[0][1];}
/**
* Get third spinor component for the x vector
*/
complex<Value> xs3() const {return _spin[0][2];}
/**
* Get fourth spinor component for the x vector
*/
complex<Value> xs4() const {return _spin[0][3];}
/**
* Get first spinor component for the y vector
*/
complex<Value> ys1() const {return _spin[1][0];}
/**
* Get second spinor component for the y vector
*/
complex<Value> ys2() const {return _spin[1][1];}
/**
* Get third spinor component for the y vector
*/
complex<Value> ys3() const {return _spin[1][2];}
/**
* Get fourth spinor component for the y vector
*/
complex<Value> ys4() const {return _spin[1][3];}
/**
* Get first spinor component for the z vector
*/
complex<Value> zs1() const {return _spin[2][0];}
/**
* Get second spinor component for the z vector
*/
complex<Value> zs2() const {return _spin[2][1];}
/**
* Get third spinor component for the z vector
*/
complex<Value> zs3() const {return _spin[2][2];}
/**
* Get fourth spinor component for the z vector
*/
complex<Value> zs4() const {return _spin[2][3];}
/**
* Get first spinor component for the t vector
*/
complex<Value> ts1() const {return _spin[3][0];}
/**
* Get second spinor component for the t vector
*/
complex<Value> ts2() const {return _spin[3][1];}
/**
* Get third spinor component for the t vector
*/
complex<Value> ts3() const {return _spin[3][2];}
/**
* Get fourth spinor component for the t vector
*/
complex<Value> ts4() const {return _spin[3][3];}
/**
* Set first spinor component for the x vector
*/
void setXS1(complex<Value> in) {_spin[0][0]=in;}
/**
* Set second spinor component for the x vector
*/
void setXS2(complex<Value> in) {_spin[0][1]=in;}
/**
* Set third spinor component for the x vector
*/
void setXS3(complex<Value> in) {_spin[0][2]=in;}
/**
* Set fourth spinor component for the x vector
*/
void setXS4(complex<Value> in) {_spin[0][3]=in;}
/**
* Set first spinor component for the y vector
*/
void setYS1(complex<Value> in) {_spin[1][0]=in;}
/**
* Set second spinor component for the y vector
*/
void setYS2(complex<Value> in) {_spin[1][1]=in;}
/**
* Set third spinor component for the y vector
*/
void setYS3(complex<Value> in) {_spin[1][2]=in;}
/**
* Set fourth spinor component for the y vector
*/
void setYS4(complex<Value> in) {_spin[1][3]=in;}
/**
* Set first spinor component for the z vector
*/
void setZS1(complex<Value> in) {_spin[2][0]=in;}
/**
* Set second spinor component for the z vector
*/
void setZS2(complex<Value> in) {_spin[2][1]=in;}
/**
* Set third spinor component for the z vector
*/
void setZS3(complex<Value> in) {_spin[2][2]=in;}
/**
* Set fourth spinor component for the z vector
*/
void setZS4(complex<Value> in) {_spin[2][3]=in;}
/**
* Set first spinor component for the t vector
*/
void setTS1(complex<Value> in) {_spin[3][0]=in;}
/**
* Set second spinor component for the t vector
*/
void setTS2(complex<Value> in) {_spin[3][1]=in;}
/**
* Set third spinor component for the t vector
*/
void setTS3(complex<Value> in) {_spin[3][2]=in;}
/**
* Set fourth spinor component for the t vector
*/
void setTS4(complex<Value> in ) {_spin[3][3]=in;}
//@}
/** @name Arithmetic operators. */
//@{
/**
* dot product with a polarization vector
*/
LorentzSpinorBar<Value> dot(const LorentzPolarizationVector & vec) const {
LorentzSpinorBar<Value> output(_type);
for(unsigned int ix=0;ix<4;++ix) {
output[ix]=_spin[3][ix]*vec.t()-_spin[0][ix]*vec.x()
-_spin[1][ix]*vec.y()-_spin[2][ix]*vec.z();
}
return output;
}
/**
* dot product with a 4-momentum
*/
LorentzSpinorBar<Value> dot(const LorentzMomentum & invec) const {
LorentzSpinorBar<Value> output(_type);
LorentzVector<double> vec = UnitRemoval::InvE * invec;
unsigned int ix;
for(ix=0;ix<4;++ix) {
output[ix]=_spin[3][ix]*vec.t()-_spin[0][ix]*vec.x()
-_spin[1][ix]*vec.y()-_spin[2][ix]*vec.z();
}
return output;
}
//@}
/** @name Transformations. */
//@{
/**
* return the barred spinor
*/
LorentzRSSpinor<Value> bar() const;
/**
* Standard Lorentz boost specifying the components of the beta vector.
*/
LorentzRSSpinorBar & boost(double,double,double);
/**
* Standard Lorentz boost specifying the beta vector.
*/
LorentzRSSpinorBar & boost(const Boost &);
/**
* General transform
*/
LorentzRSSpinorBar & transform(const LorentzRotation &);
//@}
/** @name Functions related to type. */
//@{
/**
* Return the type of the spinor.
*/
SpinorType Type() const {return _type;}
//@}
/**
* Current \f$\bar{f}^\alpha(c_LP_L+c_RP_R)f\f$ for general couplings.
* @param f The unbarred spinor
* @param left The left-handed coupling, \f$c_L\f$.
* @param right The right-handed coupling, \f$c_R\f$.
*/
template <typename ValueB>
- LorentzVector<complex<
- typename BinaryOpTraits<Value,ValueB>::MulT> >
- generalCurrent(LorentzSpinor<ValueB>& f, Complex left, Complex right) {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto generalCurrent(LorentzSpinor<ValueB>& f,
+ Complex left, Complex right)
+ -> LorentzVector<decltype(left*ts1()*f.s1())>
+ {
+ typedef decltype(left*ts1()*f.s1()) ResultT;
ResultT output[4];
unsigned int iz;
for(iz=0;iz<4;++iz){
output[iz]= left*(_spin[iz][0]*f.s1()+_spin[iz][1]*f.s2())
+right*(_spin[iz][2]*f.s3()+_spin[iz][3]*f.s4());
}
return LorentzVector<ResultT>(output[0],output[1],
output[2],output[3]);
}
private:
/**
* Type of spinor.
*/
SpinorType _type;
/**
* Storage of the components.
*/
std::array<std::array<complex<Value>,4>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzRSSpinorBar.tcc"
#endif
#endif
diff --git a/Helicity/LorentzSpinor.h b/Helicity/LorentzSpinor.h
--- a/Helicity/LorentzSpinor.h
+++ b/Helicity/LorentzSpinor.h
@@ -1,494 +1,507 @@
// -*- C++ -*-
//
// LorentzSpinor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzSpinor_H
#define ThePEG_LorentzSpinor_H
// This is the declaration of the LorentzSpinor class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzSpinor.fh"
#include "LorentzSpinorBar.h"
#include "LorentzPolarizationVector.h"
#include "LorentzTensor.h"
#include <array>
namespace ThePEG{
namespace Helicity{
/**
* The LorentzSpinor class is designed to store a spinor. In addition
* to storing the components of the spinor, information is stored on
* the representation of the type of spinor, for example u or v type.
*
* At the moment only one choice of the Dirac matrix representation
* is supported. For high-energy calculations the choice made by the
* HELAS collaboration is more efficient for numerical
* calculations. In this representation
*
* \f[
* \gamma_{i=1,2,3}=\left(\begin{array}{cc}
* 0 & \sigma_i \\
* -\sigma_i & 0
* \end{array}\right)
* \quad
* \gamma_0=\left(\begin{array}{cc}
* 0 & 1 \\
* 1 & 0
* \end{array}\right)
* \quad
* \gamma_5=\left(\begin{array}{cc}
* -1 & 0 \\
* 0 & 1
* \end{array}\right)
* \f]
*
* The type of the spinor is also stored using the SpinorType
* enumeration. There are three types supported SpinorType::u,
* SpinorType::v, SpinorType::unknown. This information is intended
* mainly for use in the case of Majorana particles where matrix
* elements can be calculated with either u or v type spinors and
* knowledge of which was used will be needed in order to give the
* correct correlations. The SpinorType::unknowne is intended for
* cases where either the spinor for an off-shell line in a matrix
* element calculation or the information is genuinely unknown.
*
* The LorentzSpinorBar class is also provided to store the barred
* spinor.
*
* @see LorentzSpinorBar
*
* @author Peter Richardson
*
*/
template<typename Value>
class LorentzSpinor {
public:
/** @name Standard constructors. */
//@{
/**
* Default zero constructor, optionally specifying \a t, the type.
*/
LorentzSpinor(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
/**
* Constructor with complex numbers specifying the components,
* optionally specifying \a s, the type.
*/
LorentzSpinor(complex<Value> a,complex<Value> b,
complex<Value> c,complex<Value> d,
SpinorType s = SpinorType::unknown) : _type(s), _spin{{a,b,c,d}} {}
//@}
/** @name Access the components. */
//@{
/**
* Subscript operator to return spinor components
*/
complex<Value> operator[](int i) const {
assert( i >= 0 && i <= 3 );
return _spin[i];
}
/**
* Subscript operator to return spinor components
*/
complex<Value> operator()(int i) const {
assert( i >= 0 && i <= 3 );
return _spin[i];
}
/**
* Set components by index.
*/
complex<Value> & operator()(int i) {
assert( i >= 0 && i <= 3 );
return _spin[i];
}
/**
* Set components by index.
*/
complex<Value> & operator[](int i) {
assert( i >= 0 && i <= 3 );
return _spin[i];
}
/**
* Get first component.
*/
complex<Value> s1() const {return _spin[0];}
/**
* Get second component.
*/
complex<Value> s2() const {return _spin[1];}
/**
* Get third component.
*/
complex<Value> s3() const {return _spin[2];}
/**
* Get fourth component.
*/
complex<Value> s4() const {return _spin[3];}
/**
* Set first component.
*/
void setS1(complex<Value> in) {_spin[0]=in;}
/**
* Set second component.
*/
void setS2(complex<Value> in) {_spin[1]=in;}
/**
* Set third component.
*/
void setS3(complex<Value> in) {_spin[2]=in;}
/**
* Set fourth component.
*/
void setS4(complex<Value> in) {_spin[3]=in;}
//@}
/** @name Transformations. */
//@{
/**
* Return the barred spinor
*/
LorentzSpinorBar<Value> bar() const;
/**
* Return the conjugated spinor \f$u_c=C\bar{u}^T\f$. This operation
* transforms u-spinors to v-spinors and vice-versa and is required when
* dealing with majorana particles.
*/
LorentzSpinor conjugate() const;
/**
* Standard Lorentz boost specifying the components of the beta vector.
*/
LorentzSpinor & boost(double,double,double);
/**
* Standard Lorentz boost specifying the beta vector.
*/
LorentzSpinor & boost(const Boost &);
/**
* General Lorentz transformation
*/
LorentzSpinor & transform(const SpinHalfLorentzRotation & );
/**
* General Lorentz transformation
*/
LorentzSpinor & transform(const LorentzRotation & r) {
transform(r.half());
return *this;
}
//@}
/** @name Functions related to type. */
//@{
/**
* Return the type of the spinor.
*/
SpinorType Type() const {return _type;}
//@}
/**
* @name Functions to apply the projection operator
*/
//@{
/**
* Apply \f$p\!\!\!\!\!\not\,\,\,+m\f$
*/
template<typename ValueB>
- LorentzSpinor<typename BinaryOpTraits<Value,ValueB>::MulT>
- projectionOperator(const LorentzVector<ValueB> & p, const ValueB & m) const {
- typedef typename BinaryOpTraits<Value,ValueB>::MulT ResultT;
- LorentzSpinor<ResultT> spin;
+ auto projectionOperator(const LorentzVector<ValueB> & p,
+ const ValueB & m) const
+ -> LorentzSpinor<decltype(m*s1())>
+ {
+ LorentzSpinor<decltype(m*s1())> spin;
static const Complex ii(0.,1.);
complex<ValueB> p0pp3=p.t()+p.z();
complex<ValueB> p0mp3=p.t()-p.z();
complex<ValueB> p1pp2=p.x()+ii*p.y();
complex<ValueB> p1mp2=p.x()-ii*p.y();
spin.setS1(m*s1()+p0mp3*s3()-p1mp2*s4());
spin.setS2(m*s2()+p0pp3*s4()-p1pp2*s3());
spin.setS3(m*s3()+p0pp3*s1()+p1mp2*s2());
spin.setS4(m*s4()+p0mp3*s2()+p1pp2*s1());
return spin;
}
/**
* Apply \f$g^LP_L+g^RP_R\f$
*/
LorentzSpinor
helicityProjectionOperator(const Complex & gL, const Complex & gR) const {
LorentzSpinor spin;
spin.setS1(gL*s1());
spin.setS2(gL*s2());
spin.setS3(gR*s3());
spin.setS4(gR*s4());
return spin;
}
//@}
/** @name Functions to calculate certain currents. */
//@{
/**
* Apply \f$p\!\!\!\!\!\not\f$
*/
template<typename ValueB>
- LorentzSpinor<typename BinaryOpTraits<Value,ValueB>::MulT>
- slash(const LorentzVector<ValueB> & p) const {
- typedef typename BinaryOpTraits<Value,ValueB>::MulT ResultT;
- LorentzSpinor<ResultT> spin;
+ auto slash(const LorentzVector<ValueB> & p) const
+ -> LorentzSpinor<decltype(p.t()*s3())>
+ {
+ LorentzSpinor<decltype(p.t()*s3())> spin;
static const Complex ii(0.,1.);
complex<ValueB> p0pp3=p.t()+p.z();
complex<ValueB> p0mp3=p.t()-p.z();
complex<ValueB> p1pp2=p.x()+ii*p.y();
complex<ValueB> p1mp2=p.x()-ii*p.y();
spin.setS1(p0mp3*s3()-p1mp2*s4());
spin.setS2(p0pp3*s4()-p1pp2*s3());
spin.setS3(p0pp3*s1()+p1mp2*s2());
spin.setS4(p0mp3*s2()+p1pp2*s1());
return spin;
}
/**
* Apply \f$p\!\!\!\!\!\not\f$
*/
template<typename ValueB>
- LorentzSpinor<typename BinaryOpTraits<Value,ValueB>::MulT>
- slash(const LorentzVector<complex<ValueB> > & p) const {
- typedef typename BinaryOpTraits<Value,ValueB>::MulT ResultT;
- LorentzSpinor<ResultT> spin;
+ auto slash(const LorentzVector<complex<ValueB> > & p) const
+ -> LorentzSpinor<decltype(p.t()*s3())>
+ {
+ LorentzSpinor<decltype(p.t()*s3())> spin;
static const Complex ii(0.,1.);
complex<ValueB> p0pp3=p.t()+p.z();
complex<ValueB> p0mp3=p.t()-p.z();
complex<ValueB> p1pp2=p.x()+ii*p.y();
complex<ValueB> p1mp2=p.x()-ii*p.y();
spin.setS1(p0mp3*s3()-p1mp2*s4());
spin.setS2(p0pp3*s4()-p1pp2*s3());
spin.setS3(p0pp3*s1()+p1mp2*s2());
spin.setS4(p0mp3*s2()+p1pp2*s1());
return spin;
}
/**
* Calculate the left-handed current \f$\bar{f}\gamma^\mu P_Lf\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- LorentzVector<complex<typename BinaryOpTraits<Value,ValueB>::MulT> >
- leftCurrent(const LorentzSpinorBar<ValueB>& fb) const {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto leftCurrent(const LorentzSpinorBar<ValueB>& fb) const
+ -> LorentzVector<decltype(fb.s3()*s2())>
+ {
+ typedef decltype(fb.s3()*s2()) ResultT;
LorentzVector<ResultT> vec;
Complex ii(0.,1.);
ResultT p1(fb.s3()*s2()),p2(fb.s4()*s1());
vec.setX( -(p1+p2) );
vec.setY( ii*(p1-p2) );
p1 = fb.s3()*s1();p2 = fb.s4()*s2();
vec.setZ( -(p1-p2) );
vec.setT( (p1+p2) );
return vec;
}
/**
* Calculate the right-handed current \f$\bar{f}\gamma^\mu P_Rf\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- LorentzVector<complex<typename BinaryOpTraits<Value,ValueB>::MulT> >
- rightCurrent(const LorentzSpinorBar<ValueB>& fb) const {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto rightCurrent(const LorentzSpinorBar<ValueB>& fb) const
+ -> LorentzVector<decltype(fb.s1()*s4())>
+ {
+ typedef decltype(fb.s1()*s4()) ResultT;
LorentzVector<ResultT> vec;
Complex ii(0.,1.);
ResultT p1(fb.s1()*s4()),p2(fb.s2()*s3());
vec.setX( (p1+p2));
vec.setY( -ii*(p1-p2));
p1 = fb.s1()*s3();p2 = fb.s2()*s4();
vec.setZ( (p1-p2));
vec.setT( (p1+p2));
return vec;
}
/**
* Calculate the vector current \f$\bar{f}\gamma^\mu f\f$
* @param fb The barred spinor.
*/
template<typename ValueB>
- LorentzVector<complex<typename BinaryOpTraits<Value,ValueB>::MulT> >
- vectorCurrent(const LorentzSpinorBar<ValueB>& fb) const {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto vectorCurrent(const LorentzSpinorBar<ValueB>& fb) const
+ -> LorentzVector<decltype(fb.s1()*s4())>
+ {
+ typedef decltype(fb.s1()*s4()) ResultT;
LorentzVector<ResultT> vec;
Complex ii(0.,1.);
ResultT s1s4(fb.s1()*s4()),s2s3(fb.s2()*s3()),
s3s2(fb.s3()*s2()),s4s1(fb.s4()*s1()),
s1s3(fb.s1()*s3()),s2s4(fb.s2()*s4()),
s3s1(fb.s3()*s1()),s4s2(fb.s4()*s2());
vec.setX( s1s4+s2s3-s3s2-s4s1 );
vec.setY( -ii*(s1s4-s2s3-s3s2+s4s1));
vec.setZ( s1s3-s2s4-s3s1+s4s2 );
vec.setT( s1s3+s2s4+s3s1+s4s2);
return vec;
}
/**
* Calculate a general current with arbitary left and right couplings,
* i.e. \f$\bar{f}\gamma^\mu(c_lP_L+c_RP_R)f\f$
* @param fb The barred spinor.
* @param left The left coupling, \f$c_L\f$.
* @param right The right coupling, \f$c_R\f$.
*/
template<typename ValueB>
- LorentzVector<complex<typename BinaryOpTraits<Value,ValueB>::MulT> >
- generalCurrent(const LorentzSpinorBar<ValueB>& fb,
- Complex left, Complex right) const {
- typedef complex<typename BinaryOpTraits<Value,ValueB>::MulT> ResultT;
+ auto generalCurrent(const LorentzSpinorBar<ValueB>& fb,
+ Complex left, Complex right) const
+ -> LorentzVector<decltype(fb.s3()*s2())>
+ {
+ typedef decltype(fb.s3()*s2()) ResultT;
LorentzVector<ResultT> vec;
Complex ii(0.,1.);
ResultT p1(fb.s3()*s2()),p2(fb.s4()*s1());
vec.setX( -left*(p1+p2));
vec.setY( ii*left*(p1-p2));
p1 = fb.s3()*s1();p2 = fb.s4()*s2();
vec.setZ( -left*(p1-p2));
vec.setT( left*(p1+p2));
p1=fb.s1()*s4();p2=fb.s2()*s3();
vec.setX(vec.x()+right*(p1+p2));
vec.setY(vec.y()-ii*right*(p1-p2));
p1 = fb.s1()*s3();p2 = fb.s2()*s4();
vec.setZ(vec.z()+right*(p1-p2));
vec.setT(vec.t()+right*(p1+p2));
return vec;
}
//@}
/** @name Functions to calculate certain scalars. */
//@{
/**
* Calculate the left-handed scalar \f$\bar{f}P_Lf\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- leftScalar(const LorentzSpinorBar<ValueB>& fb) const {
+ auto leftScalar(const LorentzSpinorBar<ValueB>& fb) const
+ -> decltype(fb.s1()*s1())
+ {
return fb.s1()*s1()+fb.s2()*s2();
}
/**
* Calculate the right-handed scalar \f$\bar{f}P_Rf\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- rightScalar(const LorentzSpinorBar<ValueB>& fb) const {
+ auto rightScalar(const LorentzSpinorBar<ValueB>& fb) const
+ -> decltype(fb.s3()*s3())
+ {
return fb.s3()*s3()+fb.s4()*s4();
}
/**
* Calculate the scalar \f$\bar{f}f\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- scalar(const LorentzSpinorBar<ValueB>& fb) const {
- return fb.s1()*s1()+fb.s2()*s2()+fb.s3()*s3()+fb.s4()*s4();
+ auto scalar(const LorentzSpinorBar<ValueB>& fb) const
+ -> decltype(fb.s1()*s1())
+ {
+ return fb.s1()*s1()+fb.s2()*s2()
+ +fb.s3()*s3()+fb.s4()*s4();
}
/**
* Calculate the pseudoscalar \f$\bar{f}\gamma_5f\f$.
* @param fb The barred spinor.
*/
template<typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- pseudoScalar(const LorentzSpinorBar<ValueB>& fb) const {
- return -fb.s1()*s1()-fb.s2()*s2()+fb.s3()*s3()+fb.s4()*s4();
+ auto pseudoScalar(const LorentzSpinorBar<ValueB>& fb) const
+ -> decltype(fb.s1()*s1())
+ {
+ return -fb.s1()*s1()-fb.s2()*s2()
+ +fb.s3()*s3()+fb.s4()*s4();
}
/**
* Calculate a general scalar product with arbitary left and right couplings,
* i.e. \f$\bar{f}c_lP_L+c_RP_Rf\f$
* @param fb The barred spinor.
* @param left The left coupling, \f$c_L\f$.
* @param right The right coupling, \f$c_R\f$.
*/
template<typename ValueB>
- complex<typename BinaryOpTraits<Value,ValueB>::MulT>
- generalScalar(const LorentzSpinorBar<ValueB>& fb,
- Complex left, Complex right) const {
- return left*(fb.s1()*s1()+fb.s2()*s2())
+ auto generalScalar(const LorentzSpinorBar<ValueB>& fb,
+ Complex left, Complex right) const
+ -> decltype(left*fb.s1()*s1())
+ {
+ return left*(fb.s1()*s1()+fb.s2()*s2())
+ right*(fb.s3()*s3()+fb.s4()*s4());
}
//@}
/**
* Calculate the product with \f$\sigma^{\mu\nu}\f$, i.e.
* \f$\bar{f}\sigma^{\mu\nu}f\f$
*/
template<typename ValueB>
- LorentzTensor<typename BinaryOpTraits<Value,ValueB>::MulT>
- sigma(const LorentzSpinorBar<ValueB>& fb) const {
- typedef typename BinaryOpTraits<Value,ValueB>::MulT ResultT;
+ auto sigma(const LorentzSpinorBar<ValueB>& fb) const
+ -> LorentzTensor<decltype(fb.s1()*s1())>
+ {
+ typedef decltype(fb.s1()*s1()) ResultT;
LorentzTensor<ResultT> output;
- complex<ResultT> s11(fb.s1()*s1()),s22(fb.s2()*s2()),
+ ResultT s11(fb.s1()*s1()),s22(fb.s2()*s2()),
s33(fb.s3()*s3()),s44(fb.s4()*s4()),
s12(fb.s1()*s2()),s21(fb.s2()*s1()),
s34(fb.s3()*s4()),s43(fb.s4()*s3());
Complex ii(0.,1.);
- complex<ResultT> zero;
+ ResultT zero;
zero = ZERO;
output.setTT( zero );
output.setTX(-ii*( s12+s21-s34-s43));
output.setTY( -s12+s21+s34-s43 );
output.setTZ(-ii*( s11-s22-s33+s44));
output.setXT( -output.tx() );
output.setXX( zero );
output.setXY( s11-s22+s33-s44 );
output.setXZ(-ii*(-s12+s21-s34+s43));
output.setYT( -output.ty() );
output.setYX( -output.xy() );
output.setYY( zero );
output.setYZ( s12+s21+s34+s43 );
output.setZT( -output.tz() );
output.setZX( -output.xz() );
output.setZY( -output.yz() );
output.setZZ( zero );
return output;
}
private:
/**
* Type of spinor
*/
SpinorType _type;
/**
* Storage of the components.
*/
std::array<complex<Value>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzSpinor.tcc"
#endif
#endif
diff --git a/Helicity/LorentzSpinorBar.h b/Helicity/LorentzSpinorBar.h
--- a/Helicity/LorentzSpinorBar.h
+++ b/Helicity/LorentzSpinorBar.h
@@ -1,233 +1,235 @@
// -*- C++ -*-
//
// LorentzSpinorBar.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzSpinorBar_H
#define ThePEG_LorentzSpinorBar_H
// This is the declaration of the LorentzSpinorBar class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "HelicityDefinitions.h"
#include "LorentzSpinor.fh"
#include "LorentzSpinorBar.fh"
namespace ThePEG {
namespace Helicity {
/**
* The LorentzSpinorBar class implements the storage of a barred
* LorentzSpinor. The design is based on that of the LorentzSpinor
* class where the details of the implemented are discussed in more
* detail.
*
* @see LorentzSpinor
*
* @author Peter Richardson
*/
template<typename Value>
class LorentzSpinorBar {
public:
/** @name Standard constructors. */
//@{
/**
* Default zero constructor, optionally specifying \a t, the type
*/
LorentzSpinorBar(SpinorType t = SpinorType::unknown) : _type(t), _spin() {}
/**
* Constructor with complex numbers specifying the components,
* optionally specifying \a t, the type
*/
LorentzSpinorBar(complex<Value> a, complex<Value> b,
complex<Value> c, complex<Value> d,
SpinorType t = SpinorType::unknown)
: _type(t), _spin{{a,b,c,d}} {}
//@}
/** @name Access the components. */
//@{
/**
* Subscript operator to return spinor components
*/
complex<Value> operator[](int i) const {
assert( i>= 0 && i <= 3 );
return _spin[i];
}
/**
* Subscript operator to return spinor components
*/
complex<Value> operator()(int i) const {
assert( i>= 0 && i <= 3 );
return _spin[i];
}
/**
* Set components by index.
*/
complex<Value> & operator()(int i) {
assert( i>= 0 && i <= 3 );
return _spin[i];
}
/**
* Set components by index.
*/
complex<Value> & operator[](int i) {
assert( i>= 0 && i <= 3 );
return _spin[i];
}
/**
* Get first component.
*/
complex<Value> s1() const {return _spin[0];}
/**
* Get second component.
*/
complex<Value> s2() const {return _spin[1];}
/**
* Get third component.
*/
complex<Value> s3() const {return _spin[2];}
/**
* Get fourth component.
*/
complex<Value> s4() const {return _spin[3];}
/**
* Set first component.
*/
void setS1(complex<Value> in) {_spin[0]=in;}
/**
* Set second component.
*/
void setS2(complex<Value> in) {_spin[1]=in;}
/**
* Set third component.
*/
void setS3(complex<Value> in) {_spin[2]=in;}
/**
* Set fourth component.
*/
void setS4(complex<Value> in) {_spin[3]=in;}
//@}
/** @name Transformations. */
//@{
/**
* Return the barred spinor
*/
LorentzSpinor<Value> bar() const;
/**
* Return the conjugated spinor \f$u_c=C\bar{u}^T\f$. This operation
* transforms u-spinors to v-spinors and vice-versa and is required when
* dealing with majorana particles.
*/
LorentzSpinorBar conjugate() const;
/**
* Standard Lorentz boost specifying the components of the beta vector.
*/
LorentzSpinorBar & boost(double,double,double);
/**
* Standard Lorentz boost specifying the beta vector.
*/
LorentzSpinorBar & boost(const Boost &);
/**
* General Lorentz transformation
*/
LorentzSpinorBar & transform(const SpinHalfLorentzRotation &) ;
/**
* General Lorentz transformation
*/
LorentzSpinorBar & transform(const LorentzRotation & r) {
transform(r.half());
return *this;
}
//@}
/** @name Functions related to type. */
//@{
/**
* Return the type of the spinor.
*/
SpinorType Type() const {return _type;}
//@}
/**
* @name Functions to apply the projection operator
*/
//@{
/**
* Apply \f$p\!\!\!\!\!\not\,\,\,+m\f$
*/
template<typename ValueB>
- LorentzSpinorBar<typename BinaryOpTraits<Value,ValueB>::MulT>
- projectionOperator(const LorentzVector<ValueB> & p, const ValueB & m) const {
- typedef typename BinaryOpTraits<Value,ValueB>::MulT ResultT;
+ auto projectionOperator(const LorentzVector<ValueB> & p,
+ const ValueB & m) const
+ -> LorentzSpinorBar<decltype(m*s1())>
+ {
+ typedef decltype(m*s1()) ResultT;
LorentzSpinorBar<ResultT> spin;
static const Complex ii(0.,1.);
complex<ValueB> p0pp3=p.t()+p.z();
complex<ValueB> p0mp3=p.t()-p.z();
complex<ValueB> p1pp2=p.x()+ii*p.y();
complex<ValueB> p1mp2=p.x()-ii*p.y();
spin.setS1(m*s1()+p0pp3*s3()+p1pp2*s4());
spin.setS2(m*s2()+p0mp3*s4()+p1mp2*s3());
spin.setS3(m*s3()+p0mp3*s1()-p1pp2*s2());
spin.setS4(m*s4()+p0pp3*s2()-p1mp2*s1());
return spin;
}
/**
* Apply \f$g^LP_L+g^RP_R\f$
*/
LorentzSpinorBar
helicityProjectionOperator(const Complex & gL, const Complex & gR) const {
LorentzSpinorBar spin;
spin.setS1(gL*s1());
spin.setS2(gL*s2());
spin.setS3(gR*s3());
spin.setS4(gR*s4());
return spin;
}
//@}
private:
/**
* Type of spinor
*/
SpinorType _type;
/**
* Storage of the components.
*/
std::array<complex<Value>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzSpinorBar.tcc"
#endif
#endif
diff --git a/Helicity/LorentzTensor.h b/Helicity/LorentzTensor.h
--- a/Helicity/LorentzTensor.h
+++ b/Helicity/LorentzTensor.h
@@ -1,505 +1,518 @@
// -*- C++ -*-
//
// LorentzTensor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzTensor_H
#define ThePEG_LorentzTensor_H
// This is the declaration of the LorentzTensor class.
-
+#include "ThePEG/Config/PhysicalQtyComplex.h"
#include "ThePEG/Config/ThePEG.h"
#include "LorentzPolarizationVector.h"
namespace ThePEG {
namespace Helicity {
// compiler magic needs these pre-declarations to make friend templates work
template<typename Value> class LorentzTensor;
/**
* The LorentzTensor class is designed to implement the storage of a
* complex tensor to be used to representation the wavefunction of a
* spin-2 particle.
*
* At the moment it only implements the storage of the tensor
* components but it is envisaged that it will be extended to include
* boost methods etc.
*
* @author Peter Richardson
*
*/
template<typename Value>
class LorentzTensor {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default zero constructor.
*/
LorentzTensor() = default;
/**
* Constructor specifyign all components.
*/
LorentzTensor(complex<Value> xx, complex<Value> xy,
- complex<Value> xz, complex<Value> xt,
- complex<Value> yx, complex<Value> yy,
- complex<Value> yz, complex<Value> yt,
- complex<Value> zx, complex<Value> zy,
- complex<Value> zz, complex<Value> zt,
- complex<Value> tx, complex<Value> ty,
- complex<Value> tz, complex<Value> tt)
+ complex<Value> xz, complex<Value> xt,
+ complex<Value> yx, complex<Value> yy,
+ complex<Value> yz, complex<Value> yt,
+ complex<Value> zx, complex<Value> zy,
+ complex<Value> zz, complex<Value> zt,
+ complex<Value> tx, complex<Value> ty,
+ complex<Value> tz, complex<Value> tt)
: _tensor{{ {{xx,xy,xz,xt}},
{{yx,yy,yz,yt}},
{{zx,zy,zz,zt}},
{{tx,ty,tz,tt}} }} {}
/**
* Constructor in terms of two polarization vectors.
*/
LorentzTensor(const LorentzPolarizationVector & p,
- const LorentzPolarizationVector & q) {
+ const LorentzPolarizationVector & q) {
setXX(p.x() * q.x()); setYX(p.y() * q.x());
setZX(p.z() * q.x()); setTX(p.t() * q.x());
setXY(p.x() * q.y()); setYY(p.y() * q.y());
setZY(p.z() * q.y()); setTY(p.t() * q.y());
setXZ(p.x() * q.z()); setYZ(p.y() * q.z());
setZZ(p.z() * q.z()); setTZ(p.t() * q.z());
setXT(p.x() * q.t()); setYT(p.y() * q.t());
setZT(p.z() * q.t()); setTT(p.t() * q.t());
}
//@}
/** @name Access individual components. */
//@{
/**
* Get x,x component.
*/
complex<Value> xx() const {return _tensor[0][0];}
/**
* Get y,x component.
*/
complex<Value> yx() const {return _tensor[1][0];}
/**
* Get z,x component.
*/
complex<Value> zx() const {return _tensor[2][0];}
/**
* Get t,x component.
*/
complex<Value> tx() const {return _tensor[3][0];}
/**
* Get x,y component.
*/
complex<Value> xy() const {return _tensor[0][1];}
/**
* Get y,y component.
*/
complex<Value> yy() const {return _tensor[1][1];}
/**
* Get z,y component.
*/
complex<Value> zy() const {return _tensor[2][1];}
/**
* Get t,y component.
*/
complex<Value> ty() const {return _tensor[3][1];}
/**
* Get x,z component.
*/
complex<Value> xz() const {return _tensor[0][2];}
/**
* Get y,z component.
*/
complex<Value> yz() const {return _tensor[1][2];}
/**
* Get z,z component.
*/
complex<Value> zz() const {return _tensor[2][2];}
/**
* Get t,z component.
*/
complex<Value> tz() const {return _tensor[3][2];}
/**
* Get x,t component.
*/
complex<Value> xt() const {return _tensor[0][3];}
/**
* Get y,t component.
*/
complex<Value> yt() const {return _tensor[1][3];}
/**
* Get z,t component.
*/
complex<Value> zt() const {return _tensor[2][3];}
/**
* Get t,t component.
*/
complex<Value> tt() const {return _tensor[3][3];}
/**
* Set x,x component.
*/
void setXX(complex<Value> a) {_tensor[0][0]=a;}
/**
* Set y,x component.
*/
void setYX(complex<Value> a) {_tensor[1][0]=a;}
/**
* Set z,x component.
*/
void setZX(complex<Value> a) {_tensor[2][0]=a;}
/**
* Set t,x component.
*/
void setTX(complex<Value> a) {_tensor[3][0]=a;}
/**
* Set x,y component.
*/
void setXY(complex<Value> a) {_tensor[0][1]=a;}
/**
* Set y,y component.
*/
void setYY(complex<Value> a) {_tensor[1][1]=a;}
/**
* Set z,y component.
*/
void setZY(complex<Value> a) {_tensor[2][1]=a;}
/**
* Set t,y component.
*/
void setTY(complex<Value> a) {_tensor[3][1]=a;}
/**
* Set x,z component.
*/
void setXZ(complex<Value> a) {_tensor[0][2]=a;}
/**
* Set y,z component.
*/
void setYZ(complex<Value> a) {_tensor[1][2]=a;}
/**
* Set z,z component.
*/
void setZZ(complex<Value> a) {_tensor[2][2]=a;}
/**
* Set t,z component.
*/
void setTZ(complex<Value> a) {_tensor[3][2]=a;}
/**
* Set x,t component.
*/
void setXT(complex<Value> a) {_tensor[0][3]=a;}
/**
* Set y,t component.
*/
void setYT(complex<Value> a) {_tensor[1][3]=a;}
/**
* Set z,t component.
*/
void setZT(complex<Value> a) {_tensor[2][3]=a;}
/**
* Set t,t component.
*/
void setTT(complex<Value> a) {_tensor[3][3]=a;}
/**
* Get components by indices.
*/
complex<Value> operator () (int i, int j) const {
assert( i>=0 && i<=3 && j>=0 && j<=3);
return _tensor[i][j];
}
/**
* Set components by indices.
*/
complex<Value> & operator () (int i, int j) {
assert( i>=0 && i<=3 && j>=0 && j<=3);
return _tensor[i][j];
}
//@}
/** @name Transformations. */
//@{
/**
* Standard Lorentz boost specifying the components of the beta vector.
*/
LorentzTensor & boost(double,double,double);
/**
* Standard Lorentz boost specifying the beta vector.
*/
LorentzTensor<Value> & boost(const Boost & b) {
return boost(b.x(), b.y(), b.z());
}
/**
* General Lorentz transformation
*/
LorentzTensor & transform(const SpinOneLorentzRotation & r){
unsigned int ix,iy,ixa,iya;
LorentzTensor<Value> output;
complex<Value> temp;
for(ix=0;ix<4;++ix) {
for(iy=0;iy<4;++iy) {
- temp=complex<Value>();
- for(ixa=0;ixa<4;++ixa) {
- for(iya=0;iya<4;++iya)
- temp+=r(ix,ixa)*r(iy,iya)*(*this)(ixa,iya);
- }
- output(ix,iy)=temp;
+ temp=complex<Value>();
+ for(ixa=0;ixa<4;++ixa) {
+ for(iya=0;iya<4;++iya)
+ temp+=r(ix,ixa)*r(iy,iya)*(*this)(ixa,iya);
+ }
+ output(ix,iy)=temp;
}
}
*this=output;
return *this;
}
/**
* Return the complex conjugate.
*/
LorentzTensor<Value> conjugate() {
return LorentzTensor<Value>(conj(xx()), conj(xy()), conj(xz()), conj(xt()),
- conj(yx()), conj(yy()), conj(yz()), conj(yt()),
- conj(zx()), conj(zy()), conj(zz()), conj(zt()),
- conj(tx()), conj(ty()), conj(tz()), conj(tt()));
+ conj(yx()), conj(yy()), conj(yz()), conj(yt()),
+ conj(zx()), conj(zy()), conj(zz()), conj(zt()),
+ conj(tx()), conj(ty()), conj(tz()), conj(tt()));
}
//@}
/** @name Arithmetic operators. */
//@{
/**
* Scaling with a complex number
*/
LorentzTensor<Value> operator*=(Complex a) {
for(int ix=0;ix<4;++ix)
for(int iy=0;iy<4;++iy) _tensor[ix][iy]*=a;
return *this;
}
/**
* Scalar product with other tensor
*/
template <typename T, typename U>
- friend complex<typename BinaryOpTraits<T,U>::MulT>
- operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u);
+ friend auto
+ operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u)
+ -> decltype(t.xx()*u.xx())
+ ;
/**
* Addition.
*/
LorentzTensor<Value> operator+(const LorentzTensor<Value> & in) const {
return LorentzTensor<Value>(xx()+in.xx(),xy()+in.xy(),xz()+in.xz(),xt()+in.xt(),
- yx()+in.yx(),yy()+in.yy(),yz()+in.yz(),yt()+in.yt(),
- zx()+in.zx(),zy()+in.zy(),zz()+in.zz(),zt()+in.zt(),
- tx()+in.tx(),ty()+in.ty(),tz()+in.tz(),tt()+in.tt());
+ yx()+in.yx(),yy()+in.yy(),yz()+in.yz(),yt()+in.yt(),
+ zx()+in.zx(),zy()+in.zy(),zz()+in.zz(),zt()+in.zt(),
+ tx()+in.tx(),ty()+in.ty(),tz()+in.tz(),tt()+in.tt());
}
/**
* Subtraction.
*/
LorentzTensor<Value> operator-(const LorentzTensor<Value> & in) const {
return LorentzTensor<Value>(xx()-in.xx(),xy()-in.xy(),xz()-in.xz(),xt()-in.xt(),
- yx()-in.yx(),yy()-in.yy(),yz()-in.yz(),yt()-in.yt(),
- zx()-in.zx(),zy()-in.zy(),zz()-in.zz(),zt()-in.zt(),
- tx()-in.tx(),ty()-in.ty(),tz()-in.tz(),tt()-in.tt());
+ yx()-in.yx(),yy()-in.yy(),yz()-in.yz(),yt()-in.yt(),
+ zx()-in.zx(),zy()-in.zy(),zz()-in.zz(),zt()-in.zt(),
+ tx()-in.tx(),ty()-in.ty(),tz()-in.tz(),tt()-in.tt());
}
/**
* Trace
*/
complex<Value> trace() const {
return _tensor[3][3]-_tensor[0][0]-_tensor[1][1]-_tensor[2][2];
}
//@}
/**
* Various dot products
*/
//@{
/**
* First index dot product with polarization vector
*/
LorentzVector<complex<Value> > preDot (const LorentzPolarizationVector & vec) const {
LorentzVector<complex<Value> > output;
output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
- vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
+ vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
- vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
+ vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
- vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
+ vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
- vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
+ vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
return output;
}
/**
* Second index dot product with polarization vector
*/
LorentzVector<complex<Value> > postDot(const LorentzPolarizationVector & vec) const {
LorentzVector<complex<Value> > output;
output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
- vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
+ vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
- vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
+ vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
- vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
+ vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
- vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
+ vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
return output;
}
/**
* First index dot product with momentum
- */
- LorentzVector<complex<typename BinaryOpTraits<Value,Energy>::MulT> >
- preDot (const Lorentz5Momentum & vec) const {
- LorentzVector<complex<typename BinaryOpTraits<Value,Energy>::MulT> > output;
+ */
+ auto preDot (const Lorentz5Momentum & vec) const
+ -> LorentzVector<decltype(vec.x()*xx())>
+ {
+ LorentzVector<decltype(vec.x()*xx())> output;
output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][0]-
- vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
+ vec.y()*_tensor[1][0]-vec.z()*_tensor[2][0]);
output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][1]-
- vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
+ vec.y()*_tensor[1][1]-vec.z()*_tensor[2][1]);
output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][2]-
- vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
+ vec.y()*_tensor[1][2]-vec.z()*_tensor[2][2]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][3]-
- vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
+ vec.y()*_tensor[1][3]-vec.z()*_tensor[2][3]);
return output;
}
/**
* Second index dot product with momentum
- */
- LorentzVector<complex<typename BinaryOpTraits<Value,Energy>::MulT> >
- postDot(const Lorentz5Momentum & vec) const {
- LorentzVector<complex<typename BinaryOpTraits<Value,Energy>::MulT> > output;
+ */
+ auto postDot(const Lorentz5Momentum & vec) const
+ -> LorentzVector<decltype(vec.x()*xx())>
+ {
+ LorentzVector<decltype(vec.x()*xx())> output;
output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
- vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
+ vec.y()*_tensor[0][1]-vec.z()*_tensor[0][2]);
output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
- vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
+ vec.y()*_tensor[1][1]-vec.z()*_tensor[1][2]);
output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
- vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
+ vec.y()*_tensor[2][1]-vec.z()*_tensor[2][2]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
- vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
+ vec.y()*_tensor[3][1]-vec.z()*_tensor[3][2]);
return output;
}
//@}
private:
/**
* The components.
*/
std::array<std::array<complex<Value>,4>,4> _tensor;
};
/**
* Multiplication by a complex number.
*/
template<typename T, typename U>
-inline LorentzTensor<typename BinaryOpTraits<T,U>::MulT>
-operator*(complex<U> a, const LorentzTensor<T> & t) {
- return LorentzTensor<typename BinaryOpTraits<T,U>::MulT>
- (a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
+inline auto
+operator*(complex<U> a, const LorentzTensor<T> & t)
+-> LorentzTensor<decltype(a.real()*t.xx().real())>
+{
+ return
+ {a*t.xx(), a*t.xy(), a*t.xz(), a*t.xt(),
a*t.yx(), a*t.yy(), a*t.yz(), a*t.yt(),
a*t.zx(), a*t.zy(), a*t.zz(), a*t.zt(),
- a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt());
+ a*t.tx(), a*t.ty(), a*t.tz(), a*t.tt()};
}
/**
* Multiply a LorentzVector by a LorentzTensor.
*/
template<typename T, typename U>
-inline LorentzVector<typename BinaryOpTraits<complex<T>,U>::MulT>
-operator*(const LorentzVector<U> & invec,
- const LorentzTensor<T> & inten) {
- LorentzVector<typename BinaryOpTraits<complex<T>,U>::MulT> outvec;
- outvec.setX(invec.t()*inten(3,0)-invec.x()*inten(0,0)
- -invec.y()*inten(1,0)-invec.z()*inten(2,0));
- outvec.setY(invec.t()*inten(3,1)-invec.x()*inten(0,1)
- -invec.y()*inten(1,1)-invec.z()*inten(2,1));
- outvec.setZ(invec.t()*inten(3,2)-invec.x()*inten(0,2)
- -invec.y()*inten(1,2)-invec.z()*inten(2,2));
- outvec.setT(invec.t()*inten(3,3)-invec.x()*inten(0,3)
- -invec.y()*inten(1,3)-invec.z()*inten(2,3));
+inline auto
+operator*(const LorentzVector<U> & v,
+ const LorentzTensor<T> & t)
+-> LorentzVector<decltype(v.t()*t(3,0))>
+{
+ LorentzVector<decltype(v.t()*t(3,0))> outvec;
+ outvec.setX( v.t()*t(3,0)-v.x()*t(0,0)
+ -v.y()*t(1,0)-v.z()*t(2,0));
+ outvec.setY( v.t()*t(3,1)-v.x()*t(0,1)
+ -v.y()*t(1,1)-v.z()*t(2,1));
+ outvec.setZ( v.t()*t(3,2)-v.x()*t(0,2)
+ -v.y()*t(1,2)-v.z()*t(2,2));
+ outvec.setT( v.t()*t(3,3)-v.x()*t(0,3)
+ -v.y()*t(1,3)-v.z()*t(2,3));
return outvec;
}
/**
* Multiply a LorentzTensor by a LorentzVector.
*/
template<typename T, typename U>
-inline LorentzVector<typename BinaryOpTraits<complex<T>,U>::MulT>
-operator*(const LorentzTensor<T> & inten, const LorentzVector<U> & invec){
- LorentzVector<typename BinaryOpTraits<complex<T>,U>::MulT> outvec;
- outvec.setX(invec.t()*inten(0,3)-invec.x()*inten(0,0)
- -invec.y()*inten(0,1)-invec.z()*inten(0,2));
- outvec.setY(invec.t()*inten(1,3)-invec.x()*inten(1,0)
- -invec.y()*inten(1,1)-invec.z()*inten(1,2));
- outvec.setZ(invec.t()*inten(2,3)-invec.x()*inten(2,0)
- -invec.y()*inten(2,1)-invec.z()*inten(2,2));
- outvec.setT(invec.t()*inten(3,3)-invec.x()*inten(3,0)
- -invec.y()*inten(3,1)-invec.z()*inten(3,2));
+inline auto
+operator*(const LorentzTensor<T> & t, const LorentzVector<U> & v)
+-> LorentzVector<decltype(v.t()*t(0,3))>
+{
+ LorentzVector<decltype(v.t()*t(0,3))> outvec;
+ outvec.setX( v.t()*t(0,3)-v.x()*t(0,0)
+ -v.y()*t(0,1)-v.z()*t(0,2));
+ outvec.setY( v.t()*t(1,3)-v.x()*t(1,0)
+ -v.y()*t(1,1)-v.z()*t(1,2));
+ outvec.setZ( v.t()*t(2,3)-v.x()*t(2,0)
+ -v.y()*t(2,1)-v.z()*t(2,2));
+ outvec.setT( v.t()*t(3,3)-v.x()*t(3,0)
+ -v.y()*t(3,1)-v.z()*t(3,2));
return outvec;
}
/**
* Multiply a LorentzTensor by a LorentzTensor
*/
template <typename T, typename U>
-inline complex<typename BinaryOpTraits<T,U>::MulT>
-operator*(const LorentzTensor<T> & t, const LorentzTensor<U> & u) {
- typedef complex<typename BinaryOpTraits<T,U>::MulT> RetT;
+inline auto
+operator*(const LorentzTensor<T> & t,
+ const LorentzTensor<U> & u)
+-> decltype(t.xx()*u.xx())
+{
+ using RetT = decltype(t.xx()*u.xx());
RetT output=RetT(),temp;
for(unsigned int ix=0;ix<4;++ix) {
temp = t._tensor[ix][3]*u._tensor[ix][3];
for(unsigned int iy=0;iy<3;++iy) {
temp+= t._tensor[ix][iy]*u._tensor[ix][iy];
}
if(ix<3) output-=temp;
else output+=temp;
}
return output;
}
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzTensor.tcc"
#endif
#endif
diff --git a/Helicity/epsilon.h b/Helicity/epsilon.h
--- a/Helicity/epsilon.h
+++ b/Helicity/epsilon.h
@@ -1,69 +1,66 @@
// -*- C++ -*-
//
// epsilon.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_epsilon_H
#define ThePEG_epsilon_H
//
// This is the declaration of the epsilon class.
#include "ThePEG/Vectors/LorentzVector.h"
namespace ThePEG {
namespace Helicity {
/** \ingroup Helicity
* \author Peter Richardson
*
* This class is designed to combine 5-momenta and polarization
* vectors together with the result being the product with the
* eps function. The class is purely static and contains no data.
*
* @see LorentzPolarizationVector
* @see Lorentz5Vector
*/
/**
* Return the product
* \f$\epsilon^{\mu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}v_{3\gamma}\f$.
* @param a The first vector \f$v_{1\alpha}\f$.
* @param b The second vector \f$v_{2\alpha}\f$.
* @param c The third vector \f$v_{3\alpha}\f$.
* @return The product
* \f$\epsilon^{\mu\alpha\beta\gamma}v_{1\alpha}v_{2\beta}v_{3\gamma}\f$.
*/
template <typename A, typename B, typename C>
- inline LorentzVector<typename BinaryOpTraits<
- typename BinaryOpTraits<A,B>::MulT, C>::MulT>
- epsilon(const LorentzVector<A> & a,
- const LorentzVector<B> & b,
- const LorentzVector<C> & c) {
- typedef typename BinaryOpTraits<A,B>::MulT ABType;
- typedef typename BinaryOpTraits<ABType, C>::MulT ABCType;
- typedef LorentzVector<ABCType> ResultType;
-
- ABType diffxy = a.x() * b.y() - a.y() * b.x();
- ABType diffxz = a.x() * b.z() - a.z() * b.x();
- ABType diffxt = a.x() * b.t() - a.t() * b.x();
- ABType diffyz = a.y() * b.z() - a.z() * b.y();
- ABType diffyt = a.y() * b.t() - a.t() * b.y();
- ABType diffzt = a.z() * b.t() - a.t() * b.z();
-
+ auto epsilon(const LorentzVector<A> & a,
+ const LorentzVector<B> & b,
+ const LorentzVector<C> & c)
+ -> LorentzVector<decltype(a.x()*b.y()*c.z())>
+ {
+ auto diffxy = a.x() * b.y() - a.y() * b.x();
+ auto diffxz = a.x() * b.z() - a.z() * b.x();
+ auto diffxt = a.x() * b.t() - a.t() * b.x();
+ auto diffyz = a.y() * b.z() - a.z() * b.y();
+ auto diffyt = a.y() * b.t() - a.t() * b.y();
+ auto diffzt = a.z() * b.t() - a.t() * b.z();
+
+ using ResultType = LorentzVector<decltype(a.x()*b.x()*c.x())>;
ResultType result;
result.setX( c.z() * diffyt - c.t() * diffyz - c.y() * diffzt);
result.setY( c.t() * diffxz - c.z() * diffxt + c.x() * diffzt);
result.setZ(-c.t() * diffxy + c.y() * diffxt - c.x() * diffyt);
result.setT(-c.z() * diffxy + c.y() * diffxz - c.x() * diffyz);
return result;
}
}
}
#endif /* ThePEG_epsilon_H */
diff --git a/Vectors/Lorentz5Vector.h b/Vectors/Lorentz5Vector.h
--- a/Vectors/Lorentz5Vector.h
+++ b/Vectors/Lorentz5Vector.h
@@ -1,395 +1,325 @@
// -*- C++ -*-
//
// Lorentz5Vector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Lorentz5Vector_H
#define ThePEG_Lorentz5Vector_H
// This is the declaration of the Lorentz5vector class.
#include "LorentzVector.h"
#include "Lorentz5Vector.fh"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/Utilities/Direction.h"
#include "ThePEG/Utilities/UnitIO.h"
#include "LorentzRotation.h"
namespace ThePEG {
template <typename Value>
/**
* The Lorentz5Vector inherits from the
* <code>LorentzVector</code> class. It is templated on the
* type of the member variables. The <code>Lorentz5Vector</code> class is a
* <code>LorentzVector</code> with an extra member for the invariant
* length/mass of the vector. Note that an object of the
* <code>Lorentz5Vector</code> class may be internally inconsistent in
* that the invariant length/mass of the <code>LorentzVector</code>
* class need not be the same as the member variable representing the
* invariant length/mass. The degree of inconsistency can be accessed
* with the <code>massError()</code>, <code>energyError()</code> and
* <code>rhoError()</code> methods and an object can be made consistent
* using the <code>rescaleMass()</code>, <code>rescaleEnergy()</code> or
* <code>rescaleRho()</code> methods.
*
* @see Math
*
*/
class Lorentz5Vector: public LorentzVector<Value> {
public:
/** Template argument typedef. */
- typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
+ using Value2 = decltype(sqr(std::declval<Value>()));
public:
/// Component access.
//@{
Value x() const { return LorentzVector<Value>::x(); }
Value y() const { return LorentzVector<Value>::y(); }
Value z() const { return LorentzVector<Value>::z(); }
Value t() const { return LorentzVector<Value>::t(); }
//@}
public:
/** @name Constructors and destructor. */
//@{
/**
* Constructor giving the null vector.
*/
Lorentz5Vector() : mm() {}
/**
* Constructor giving the invariant length.
*/
Lorentz5Vector(Value m)
: LorentzVector<Value>(Value(), Value(), Value(), m), mm(m) {}
/**
* Constructor giving the components x, y, z, t. The invariant
* length is set to LorentzVector::mag().
*/
Lorentz5Vector(Value x, Value y, Value z, Value t = Value())
: LorentzVector<Value>(x, y, z, t) { rescaleMass(); }
/**
* Constructor giving the components x, y, z, t and invariant length.
* May result in an inconsistent Lorentz5Vector.
*/
Lorentz5Vector(Value x, Value y, Value z, Value t, Value tau)
: LorentzVector<Value>(x, y, z, t), mm(tau) {}
/**
* Constructor giving a 3-Vector and a time component. The invariant
* length is set to LorentzVector::mag().
*/
Lorentz5Vector(const ThreeVector<Value> & p, Value e)
: LorentzVector<Value>(p, e) { rescaleMass(); }
/**
* Constructor giving an invariant length and a 3-Vector
* component. The time component is set to the corresponding value.
*/
Lorentz5Vector(Value m, const ThreeVector<Value> & p)
: LorentzVector<Value>(p, sqrt(p.mag2() + m*m)), mm(m) {}
/**
* Constructor giving a 3-Vector, a time component and an invariant
* length. May result in an inconsistent Lorentz5Vector.
*/
Lorentz5Vector(const ThreeVector<Value> & p, Value t, Value tau)
: LorentzVector<Value>(p, t), mm(tau) {}
/**
* Constructor giving a LorentzVector and an invariant length.
* May result in an inconsistent Lorentz5Vector.
*/
Lorentz5Vector(const LorentzVector<Value> & p, Value m)
: LorentzVector<Value>(p), mm(m) {}
/**
* Copy from HepLorentzVector constructor. The invariant
* length is set to LorentzVector::mag().
*/
Lorentz5Vector(const LorentzVector<Value> & p)
: LorentzVector<Value>(p) { rescaleMass(); }
/**
* Construct from value type U convertible to Value.
*/
template<class U>
Lorentz5Vector(const Lorentz5Vector<U> & p)
: LorentzVector<Value>(p), mm(p.m) {}
//@}
/** @name Assignment and set functions. */
//@{
/**
* Set invariant length/mass.
*/
void setTau(Value a) { mm = a; }
/**
* Set invariant length/mass.
*/
void setMass(Value a) { mm = a; }
/**
* Assignment. The invariant length is kept fixed. May result in an
* inconsistent Lorentz5Vector.
*/
Lorentz5Vector & operator=(const LorentzVector<Value> & q) {
LorentzVector<Value>::operator=(q);
return *this;
}
//@}
/** @name Rescale functions to make consistent. */
//@{
/**
* Rescale energy, so that the invariant length/mass of the
* LorentzVector agrees with the current one.
*/
void rescaleEnergy() {
LorentzVector<Value>::setT(sqrt(LorentzVector<Value>::vect().mag2() + mass2()));
}
/**
* Rescale spatial component, so that the invariant length/mass of
* the LorentzVector agrees with the current one.
*/
void rescaleRho() {
LorentzVector<Value>::setRho(sqrt(t()*t() - mass2()));
}
/**
* Set the invariant length/mass member, so that it agrees with the
* invariant length/mass of the LorentzVector.
*/
void rescaleMass() {
mm = LorentzVector<Value>::m();
}
//@}
/** @name Check consistency. */
//@{
/**
* Return the relative inconsistency in the mass component.
*/
double massError() const {
return sqrt(abs(Math::relativeError(mass2(),
LorentzVector<Value>::m2())));
}
/**
* Return the relative inconsistency in the energy component.
*/
double energyError() const {
return sqrt(abs(Math::relativeError(t()*t(), mass2()
+ LorentzVector<Value>::vect().mag2())));
}
/**
* Return the relative inconsistency in the spatial components.
*/
double rhoError() const {
return sqrt(abs(Math::relativeError(LorentzVector<Value>::vect().mag2(),
t()*t() - mass2())));
}
//@}
/** @name Access components. */
//@{
/**
* Mass/invariant length component squared. m2() gives
* the same calculated from the LorentzVector
*/
Value2 mass2() const { return mm > Value() ? mm*mm: -mm*mm; }
/**
* Mass/invariant length component squared. m2() gives
* the same calculated from the LorentzVector
*/
Value2 tau2() const { return mass2(); }
/**
* Mass/invariant length component. m() gives the same
* calculated from the LorentzVector
*/
Value mass() const { return mm; }
/**
* Mass/invariant length component. m() gives the same
* calculated from the LorentzVector
*/
Value tau() const { return mass(); }
/**
* Return the positive negative light-cone components (depending on
* the value of Direction<0>.
*/
Value dirPlus() const {
return Direction<0>::pos() ?
LorentzVector<Value>::plus()
:
LorentzVector<Value>::minus();
}
/**
* Return the positive negative light-cone components (depending on
* the value of Direction<0>.
*/
Value dirMinus() const {
return Direction<0>::neg() ?
LorentzVector<Value>::plus()
:
LorentzVector<Value>::minus();
}
//@}
/**
* Perform a Lorentz transformation
*/
Lorentz5Vector & transform(const LorentzRotation & r)
{
LorentzVector<Value>::transform(r.one());
return *this;
}
private:
/** The invariant mass/length member. */
Value mm;
};
/** Output a Lorentz5Vector to a stream. */
template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const Lorentz5Vector<T> & p, UT & u) {
os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
<< ounit(p.e(), u) << ounit(p.mass(), u);
}
/** Input a Lorentz5Vector from a stream. */
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, Lorentz5Vector<T> & p, UT & u) {
T x, y, z, e, mass;
is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u)
>> iunit(mass, u);
p = Lorentz5Vector<T>(x, y, z, e, mass);
}
-template <typename T, typename U>
-struct BinaryOpTraits;
-
-/**
- * Template for multiplication by scalar
- */
-template <typename T>
-struct BinaryOpTraits<Lorentz5Vector<T>, double> {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef Lorentz5Vector<T> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef Lorentz5Vector<T> DivT;
-};
-
-/**
- * Template for multiplication by scalar
- */
-template <typename U>
-struct BinaryOpTraits<double, Lorentz5Vector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef Lorentz5Vector<U> MulT;
-};
-
-/**
- * Template for multiplication for complex and Lorentz5Vector
- */
-template <typename T, typename U>
-struct BinaryOpTraits<Lorentz5Vector<T>, std::complex<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef Lorentz5Vector<std::complex<typename BinaryOpTraits<T,U>::MulT> > MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef Lorentz5Vector<std::complex<typename BinaryOpTraits<T,U>::DivT> > DivT;
-};
-
-/**
- * Template for multiplication by scalar
- */
-template <typename T, typename U>
-struct BinaryOpTraits<std::complex<T>, Lorentz5Vector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef Lorentz5Vector<std::complex<typename BinaryOpTraits<T,U>::MulT> > MulT;
-};
-
-/**
- * Template for multiplication of two Lorentz5Vectors
- */
-template <typename T, typename U>
-struct BinaryOpTraits<Lorentz5Vector<T>, Lorentz5Vector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef typename BinaryOpTraits<T,U>::MulT MulT;
-};
-
-/**
- * Template for multiplication for LorentzVector and Lorentz5Vector
- */
-template <typename T, typename U>
-struct BinaryOpTraits<LorentzVector<T>, Lorentz5Vector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef typename BinaryOpTraits<T,U>::MulT MulT;
-};
-
-/**
- * Template for multiplication for LorentzVector and Lorentz5Vector
- */
-template <typename T, typename U>
-struct BinaryOpTraits<Lorentz5Vector<T>, LorentzVector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef typename BinaryOpTraits<T,U>::MulT MulT;
-};
/// @name Dot product overloads.
//@{
template <typename ValueA, typename ValueB>
-inline typename BinaryOpTraits<ValueA,ValueB>::MulT
-operator*(const Lorentz5Vector<ValueA> & a, const Lorentz5Vector<ValueB> & b) {
+inline auto
+operator*(const Lorentz5Vector<ValueA> & a, const Lorentz5Vector<ValueB> & b)
+-> decltype(a.dot(b))
+{
return a.dot(b);
}
template <typename ValueA, typename ValueB>
-inline typename BinaryOpTraits<ValueA,ValueB>::MulT
-operator*(const LorentzVector<ValueA> & a, const Lorentz5Vector<ValueB> & b) {
+inline auto
+operator*(const LorentzVector<ValueA> & a, const Lorentz5Vector<ValueB> & b)
+-> decltype(a.dot(b))
+{
return a.dot(b);
}
template <typename ValueA, typename ValueB>
-inline typename BinaryOpTraits<ValueA,ValueB>::MulT
-operator*(const Lorentz5Vector<ValueA> & a, const LorentzVector<ValueB> & b) {
+inline auto
+operator*(const Lorentz5Vector<ValueA> & a, const LorentzVector<ValueB> & b)
+-> decltype(a.dot(b))
+{
return a.dot(b);
}
template <typename Value>
-inline typename BinaryOpTraits<Value,Value>::MulT
-operator*(const Lorentz5Vector<Value> & a, const Lorentz5Vector<Value> & b) {
+inline auto
+operator*(const Lorentz5Vector<Value> & a, const Lorentz5Vector<Value> & b)
+-> decltype(a.dot(b))
+{
return a.dot(b);
}
//@}
}
#endif /* ThePEG_Particle_H */
diff --git a/Vectors/LorentzVector.h b/Vectors/LorentzVector.h
--- a/Vectors/LorentzVector.h
+++ b/Vectors/LorentzVector.h
@@ -1,783 +1,749 @@
// -*- C++ -*-
//
// LorentzVector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_LorentzVector_H
#define ThePEG_LorentzVector_H
/**
* @file LorentzVector.h contains the LorentzVector class. Lorentz
* vectors can be created with any unit type as template parameter.
* All basic mathematical operations are supported, as well as a
* subset of the CLHEP LorentzVector functionality.
*/
#include "LorentzVector.fh"
#include "ThePEG/Utilities/Direction.h"
#include "ThePEG/Utilities/UnitIO.h"
#include "LorentzRotation.h"
#include "ThreeVector.h"
/// Debug helper function
#ifdef NDEBUG
#define ERROR_IF(condition,message) if (false) {}
#else
#define ERROR_IF(condition,message) \
if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror)
#endif
namespace ThePEG {
template <typename Value> class LorentzVector;
/**
* A 4-component Lorentz vector. It can be created with any unit type
* as template parameter. All basic mathematical operations are
* supported, as well as a subset of the CLHEP LorentzVector
* functionality.
*/
template <typename Value> class LorentzVector
{
private:
/// Value squared
- typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
+ using Value2 = decltype(sqr(std::declval<Value>()));
public:
/** @name Constructors. */
//@{
LorentzVector()
: theX(), theY(), theZ(), theT() {}
LorentzVector(Value x, Value y, Value z, Value t)
: theX(x), theY(y), theZ(z), theT(t) {}
LorentzVector(const ThreeVector<Value> & v, Value t)
: theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {}
template<typename U>
LorentzVector(const LorentzVector<U> & v)
: theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {}
//@}
/// Assignment operator
template <typename ValueB>
LorentzVector<Value> & operator=(const LorentzVector<ValueB> & b) {
setX(b.x());
setY(b.y());
setZ(b.z());
setT(b.t());
return *this;
}
public:
/// @name Component access methods.
//@{
Value x() const { return theX; }
Value y() const { return theY; }
Value z() const { return theZ; }
Value t() const { return theT; }
Value e() const { return t(); }
//@}
/// @name Component set methods.
//@{
void setX(Value x) { theX = x; }
void setY(Value y) { theY = y; }
void setZ(Value z) { theZ = z; }
void setT(Value t) { theT = t; }
void setE(Value e) { setT(e); }
//@}
public:
/// Access to the 3-component part.
ThreeVector<Value> vect() const {
return ThreeVector<Value>(x(),y(),z());
}
/// Cast to the 3-component part.
operator ThreeVector<Value>() const { return vect(); }
/// Set the 3-component part.
void setVect(const ThreeVector<Value> & p) {
theX = p.x();
theY = p.y();
theZ = p.z();
}
public:
/// The complex conjugate vector.
LorentzVector<Value> conjugate() const
{
return LorentzVector<Value>(conj(x()),conj(y()),conj(z()),conj(t()));
}
/// Squared magnitude \f$x^\mu\,x_\mu=t^2 - \vec{x}^2\f$.
Value2 m2() const
{
return (t()-z())*(t()+z()) - sqr(x()) - sqr(y());
}
/// Squared magnitude with another vector
Value2 m2(const LorentzVector<Value> & a) const {
Value tt(a.t()+t()),zz(a.z()+z());
return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y());
}
/// Magnitude (signed) \f$\pm\sqrt{|t^2 - \vec{x}^2|}\f$.
Value m() const
{
Value2 tmp = m2();
return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
}
/// Transverse mass squared \f$t^2-z^2\f$.
Value2 mt2() const { return (t()-z())*(t()+z()); }
/// Transverse mass (signed) \f$\pm\sqrt{|t^2 - z^2|}\f$.
Value mt() const
{
Value2 tmp = mt2();
return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp));
}
/// Squared transverse component of the spatial vector \f$x^2+y^2\f$.
Value2 perp2() const { return sqr(x()) + sqr(y()); }
/// Transverse component of the spatial vector \f$\pm\sqrt{x^2 + y^2}\f$.
Value perp() const { return sqrt(perp2()); }
/**
* Squared transverse component of the spatial vector with respect to the
* given axis.
*/
template <typename U>
Value2 perp2(const ThreeVector<U> & p) const
{
return vect().perp2(p);
}
/**
* Transverse component of the spatial vector with respect to the
* given axis.
*/
template <typename U>
Value perp(const ThreeVector<U> & p) const
{
return vect().perp(p);
}
/// Transverse energy squared.
Value2 et2() const
{
Value2 pt2 = vect().perp2();
return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z());
}
/// Transverse energy (signed).
Value et() const
{
Value2 etet = et2();
return e() < Value() ? -sqrt(etet) : sqrt(etet);
}
/// Transverse energy squared with respect to the given axis.
Value2 et2(const ThreeVector<double> & v) const
{
Value2 pt2 = vect().perp2(v);
Value pv = vect().dot(v.unit());
return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv);
}
/// Transverse energy with respect to the given axis (signed).
Value et(const ThreeVector<double> & v) const
{
Value2 etet = et2(v);
return e() < Value() ? -sqrt(etet) : sqrt(etet);
}
/// @name Spherical coordinates for the spatial part.
//@{
/// Radius squared.
Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
/// Radius.
Value rho() const { return sqrt(rho2()); }
/// Set new radius.
void setRho(Value newRho)
{
Value oldRho = rho();
if (oldRho == Value())
return;
double factor = newRho / oldRho;
setX(x()*factor);
setY(y()*factor);
setZ(z()*factor);
}
/// Polar angle.
double theta() const
{
assert(!(x() == Value() && y() == Value() && z() == Value()));
return atan2(perp(),z());
}
/// Cosine of the polar angle.
double cosTheta() const
{
Value ptot = rho();
assert( ptot > Value() );
return z() / ptot;
}
/// Azimuthal angle.
double phi() const {
return atan2(y(),x()) ;
}
//@}
/// Pseudorapidity of spatial part.
double eta() const {
Value m = rho();
if ( m == Value() ) return 0.0;
Value pt = max(Constants::epsilon*m, perp());
double rap = log((m + abs(z()))/pt);
return z() > ZERO? rap: -rap;
}
/// Spatial angle with another vector.
double angle(const LorentzVector<Value> & w) const
{
return vect().angle(w.vect());
}
/// Rapidity \f$\frac{1}{2}\ln\frac{t+z}{t-z} \f$
double rapidity() const {
if ( z() == ZERO ) return 0.0;
ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2()));
double rap = log((t() + abs(z()))/pt);
return z() > ZERO? rap: -rap;
}
/// Rapidity with respect to another vector
double rapidity(const Axis & ref) const {
double r = ref.mag2();
ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity");
Value vdotu = vect().dot(ref)/sqrt(r);
if ( vdotu == ZERO ) return 0.0;
ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector");
Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2()));
double rap = log((t() + abs(z()))/pt);
return z() > ZERO? rap: -rap;
}
/**
* Boost from reference frame into this vector's rest
* frame: \f$\frac{\vec{x}}{t}\f$.
*/
Boost boostVector() const {
if (t() == Value()) {
if (rho2() == Value2())
return Boost();
else
ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result");
}
// result will make analytic sense but is physically meaningless
ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector");
return vect() * (1./t());
}
/**
* Boost from reference frame into this vector's rest
* frame: \f$-\frac{\vec{x}}{t}\f$.
*/
Boost findBoostToCM() const
{
return -boostVector();
}
/// Returns the positive light-cone component \f$t + z\f$.
Value plus() const { return t() + z(); }
/// Returns the negative light-cone component \f$t - z\f$.
Value minus() const { return t() - z(); }
/// Are two vectors nearby, using Euclidean measure \f$t^2 + |\vec{x}|^2\f$?
bool isNear(const LorentzVector<Value> & w, double epsilon) const
{
Value2 limit = abs(vect().dot(w.vect()));
limit += 0.25 * sqr( t() + w.t() );
limit *= sqr(epsilon);
Value2 delta = (vect() - w.vect()).mag2();
delta += sqr( t() - w.t() );
return (delta <= limit);
}
/// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$.
LorentzVector<Value> & transform(const SpinOneLorentzRotation & m)
{
return *this = m.operator*(*this);
}
/// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$.
LorentzVector<Value> & operator*=(const SpinOneLorentzRotation & m)
{
return transform(m);
}
/// Dot product with metric \f$(+,-,-,-)\f$
template <typename U>
- typename BinaryOpTraits<Value,U>::MulT
- dot(const LorentzVector<U> & a) const
+ auto dot(const LorentzVector<U> & a) const -> decltype(t() * a.t())
{
return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() );
}
public:
/**
* Apply boost.
*
* @param bx Component x of the boost.
* @param by Component y of the boost.
* @param bz Component z of the boost.
* @param gamma Optional gamma parameter for higher numerical
* accuracy. The user has to ensure consistency. If not given, it
* will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$.
*
*/
LorentzVector<Value> &
boost(double bx, double by, double bz, double gamma=-1.)
{
const double b2 = bx*bx + by*by + bz*bz;
if ( b2 == 0.0 ) return *this;
if ( gamma < 0.0 ) {
gamma = 1.0 / sqrt(1.0 - b2);
}
const Value bp = bx*x() + by*y() + bz*z();
const double gamma2 = (gamma - 1.0)/b2;
setX(x() + gamma2*bp*bx + gamma*bx*t());
setY(y() + gamma2*bp*by + gamma*by*t());
setZ(z() + gamma2*bp*bz + gamma*bz*t());
setT(gamma*(t() + bp));
return *this;
}
/**
* Apply boost.
*
* @param b Three-vector giving the boost.
*
* @param gamma Optional gamma parameter for higher numerical
* accuracy. The user has to ensure consistency. If not given, it
* will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$.
*
*/
LorentzVector<Value> & boost(Boost b, double gamma=-1.) {
return boost(b.x(), b.y(), b.z(),gamma);
}
/**
* Apply rotation around the x-axis.
*
* @param phi Angle in radians.
*/
LorentzVector<Value> & rotateX (double phi) {
double sinphi = sin(phi);
double cosphi = cos(phi);
Value ty = y() * cosphi - z() * sinphi;
theZ = z() * cosphi + y() * sinphi;
theY = ty;
return *this;
}
/**
* Apply rotation around the y-axis.
*
* @param phi Angle in radians.
*/
LorentzVector<Value> & rotateY (double phi) {
double sinphi = sin(phi);
double cosphi = cos(phi);
Value tz = z() * cosphi - x() * sinphi;
theX = x() * cosphi + z() * sinphi;
theZ = tz;
return *this;
}
/**
* Apply rotation around the z-axis.
*
* @param phi Angle in radians.
*/
LorentzVector<Value> & rotateZ (double phi) {
double sinphi = sin(phi);
double cosphi = cos(phi);
Value tx = x() * cosphi - y() * sinphi;
theY = y() * cosphi + x() * sinphi;
theX = tx;
return *this;
}
/**
* Rotate the reference frame to a new z-axis.
*/
LorentzVector<Value> & rotateUz (const Axis & axis) {
Axis ax = axis.unit();
double u1 = ax.x();
double u2 = ax.y();
double u3 = ax.z();
double up = u1*u1 + u2*u2;
if (up>0) {
up = sqrt(up);
Value px = x(), py = y(), pz = z();
setX( (u1*u3*px - u2*py)/up + u1*pz );
setY( (u2*u3*px + u1*py)/up + u2*pz );
setZ( -up*px + u3*pz );
}
else if (u3 < 0.) {
setX(-x());
setZ(-z());
}
return *this;
}
/**
* Apply a rotation.
* @param angle Rotation angle in radians.
* @param axis Rotation axis.
*/
template <typename U>
LorentzVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
if (angle == 0.0)
return *this;
const U ll = axis.mag();
assert( ll > U() );
const double sa = sin(angle), ca = cos(angle);
const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
const Value xx = x(), yy = y(), zz = z();
setX((ca+(1-ca)*dx*dx) * xx
+((1-ca)*dx*dy-sa*dz) * yy
+((1-ca)*dx*dz+sa*dy) * zz
);
setY(((1-ca)*dy*dx+sa*dz) * xx
+(ca+(1-ca)*dy*dy) * yy
+((1-ca)*dy*dz-sa*dx) * zz
);
setZ(((1-ca)*dz*dx-sa*dy) * xx
+((1-ca)*dz*dy+sa*dx) * yy
+(ca+(1-ca)*dz*dz) * zz
);
return *this;
}
public:
/// @name Mathematical assignment operators.
//@{
template <typename ValueB>
LorentzVector<Value> & operator+=(const LorentzVector<ValueB> & a) {
theX += a.x();
theY += a.y();
theZ += a.z();
theT += a.t();
return *this;
}
template <typename ValueB>
LorentzVector<Value> & operator-=(const LorentzVector<ValueB> & a) {
theX -= a.x();
theY -= a.y();
theZ -= a.z();
theT -= a.t();
return *this;
}
LorentzVector<Value> & operator*=(double a) {
theX *= a;
theY *= a;
theZ *= a;
theT *= a;
return *this;
}
LorentzVector<Value> & operator/=(double a) {
theX /= a;
theY /= a;
theZ /= a;
theT /= a;
return *this;
}
//@}
private:
/// @name Vector components
//@{
Value theX;
Value theY;
Value theZ;
Value theT;
//@}
};
/// @name Basic mathematical operations
//@{
template <typename Value>
inline LorentzVector<double>
operator/(const LorentzVector<Value> & v, Value a) {
return LorentzVector<double>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
}
inline LorentzVector<Complex>
operator/(const LorentzVector<Complex> & v, Complex a) {
return LorentzVector<Complex>(v.x()/a, v.y()/a, v.z()/a, v.t()/a);
}
template <typename Value>
inline LorentzVector<Value> operator-(const LorentzVector<Value> & v) {
return LorentzVector<Value>(-v.x(),-v.y(),-v.z(),-v.t());
}
template <typename ValueA, typename ValueB>
inline LorentzVector<ValueA>
operator+(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
return a += b;
}
template <typename ValueA, typename ValueB>
inline LorentzVector<ValueA>
operator-(LorentzVector<ValueA> a, const LorentzVector<ValueB> & b) {
return a -= b;
}
template <typename Value>
inline LorentzVector<Value>
operator*(const LorentzVector<Value> & a, double b) {
return LorentzVector<Value>(a.x()*b, a.y()*b, a.z()*b, a.t()*b);
}
template <typename Value>
inline LorentzVector<Value>
operator*(double b, LorentzVector<Value> a) {
return a *= b;
}
template <typename ValueA, typename ValueB>
-inline
-LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(ValueB a, const LorentzVector<ValueA> & v) {
- typedef typename BinaryOpTraits<ValueB,ValueA>::MulT ResultT;
- return LorentzVector<ResultT>(a*v.x(), a*v.y(), a*v.z(), a*v.t());
+inline auto operator*(ValueB a, const LorentzVector<ValueA> & v)
+-> LorentzVector<decltype(a*v.x())>
+{
+ return {a*v.x(), a*v.y(), a*v.z(), a*v.t()};
}
template <typename ValueA, typename ValueB>
-inline
-LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(const LorentzVector<ValueA> & v, ValueB b) {
+inline auto operator*(const LorentzVector<ValueA> & v, ValueB b)
+-> LorentzVector<decltype(b*v.x())>
+{
return b*v;
}
template <typename ValueA, typename ValueB>
-inline
-LorentzVector<typename BinaryOpTraits<ValueA,ValueB>::DivT>
-operator/(const LorentzVector<ValueA> & v, ValueB b) {
- typedef typename BinaryOpTraits<ValueA,ValueB>::DivT ResultT;
- return LorentzVector<ResultT>(v.x()/b, v.y()/b, v.z()/b, v.t()/b);
+inline auto operator/(const LorentzVector<ValueA> & v, ValueB b)
+-> LorentzVector<decltype(v.x()/b)>
+{
+ return {v.x()/b, v.y()/b, v.z()/b, v.t()/b};
}
//@}
/// @name Scalar product with metric \f$(+,-,-,-)\f$
//@{
template <typename ValueA, typename ValueB>
-inline typename BinaryOpTraits<ValueA,ValueB>::MulT
-operator*(const LorentzVector<ValueA> & a, const LorentzVector<ValueB> & b) {
+inline auto
+operator*(const LorentzVector<ValueA> & a, const LorentzVector<ValueB> & b)
+-> decltype(a.dot(b))
+{
return a.dot(b);
}
-template <typename Value>
-inline typename BinaryOpTraits<Value,Value>::MulT
-operator*(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
- return a.dot(b);
-}
//@}
/// Equality
template <typename Value>
inline bool
operator==(const LorentzVector<Value> & a, const LorentzVector<Value> & b) {
return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t();
}
/// Stream output. Format \f$(x,y,z;t)\f$.
inline ostream & operator<< (ostream & os, const LorentzVector<double> & v) {
return os << "(" << v.x() << "," << v.y() << "," << v.z()
<< ";" << v.t() << ")";
}
/** Return the positive light-cone component. Or negative if the
* current Direction<0> is reversed. */
template <typename Value>
inline Value dirPlus(const LorentzVector<Value> & p) {
return Direction<0>::pos()? p.plus(): p.minus();
}
/** Return the negative light-cone component. Or positive if the
* current Direction<0> is reversed. */
template <typename Value>
inline Value dirMinus(const LorentzVector<Value> & p) {
return Direction<0>::neg()? p.plus(): p.minus();
}
/** Return the component along the positive z-axis. Or the negative
* z-axis if the current Direction<0> is reversed. */
template <typename Value>
inline Value dirZ(const LorentzVector<Value> & p) {
return Direction<0>::dir()*p.z();
}
/** Return the polar angle wrt. the positive z-axis. Or the negative
* z-axis if the current Direction<0> is reversed. */
template <typename Value>
inline double dirTheta(const LorentzVector<Value> & p) {
return Direction<0>::pos()? p.theta(): Constants::pi - p.theta();
}
/** Return the cosine of the polar angle wrt. the positive z-axis. Or
* the negative z-axis if the current Direction<0> is reversed. */
template <typename Value>
inline double dirCosTheta(const LorentzVector<Value> & p) {
return Direction<0>::pos()? p.cosTheta(): -p.cosTheta();
}
/** Get the boost vector for the LorentzVector. If the current
* Direction<0> is reversed, so is the z-component. */
template <typename Value>
inline ThreeVector<Value> dirBoostVector(const LorentzVector<Value> & p) {
ThreeVector<Value> b(p.boostVector());
if ( Direction<0>::neg() ) b.setZ(-b.z());
return b;
}
/** Create a LorentzVector giving its light-cone and transverse
* components. */
template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus, Value x, Value y) {
LorentzVector<Value> r(x, y, 0.5*(plus-minus), 0.5*(plus+minus));
return r;
}
/** Create a LorentzVector giving its light-cone components. */
template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus) {
// g++-3.3 has a problem with using Value() directly
// gcc-bug c++/3650, fixed in 3.4
static const Value zero = Value();
LorentzVector<Value> r(zero, zero,
0.5*(plus-minus), 0.5*(plus+minus));
return r;
}
}
// delayed header inclusion to break inclusion loop:
// LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec
#include "Transverse.h"
namespace ThePEG {
/** Create a LorentzVector giving its light-cone and transverse
* components. */
template <typename Value>
inline LorentzVector<Value>
lightCone(Value plus, Value minus, Transverse<Value> pt) {
LorentzVector<Value> r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus));
return r;
}
/** Create a LorentzVector giving its light-cone and transverse
* components. If the current Direction<0> is reversed, so is the
* z-component. */
template <typename Value>
inline LorentzVector<Value>
lightConeDir(Value plus, Value minus,
Value x = Value(), Value y = Value()) {
LorentzVector<Value> r(x, y, Direction<0>::dir()*0.5*(plus - minus),
0.5*(plus + minus));
return r;
}
/** Create a LorentzVector giving its light-cone and transverse
* components. If the current Direction<0> is reversed, so is the
* z-component. */
template <typename Value>
inline LorentzVector<Value>
lightConeDir(Value plus, Value minus, Transverse<Value> pt) {
LorentzVector<Value> r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus),
0.5*(plus + minus));
return r;
}
/** Output a LorentzVector with units to a stream. */
template <typename OStream, typename UnitT, typename Value>
void ounitstream(OStream & os, const LorentzVector<Value> & p, UnitT & u) {
os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u)
<< ounit(p.e(), u);
}
/** Input a LorentzVector with units from a stream. */
template <typename IStream, typename UnitT, typename Value>
void iunitstream(IStream & is, LorentzVector<Value> & p, UnitT & u) {
Value x, y, z, e;
is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u);
p = LorentzVector<Value>(x, y, z, e);
}
-/// @name Traits for binary operations
-//@{
-template <typename T, typename U>
-struct BinaryOpTraits;
-/**
- * Template for multiplication by scalar
- */
-template <typename T, typename U>
-struct BinaryOpTraits<LorentzVector<T>, U> {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef LorentzVector<typename BinaryOpTraits<T,U>::MulT> MulT;
- /** The type resulting from division of one template type with
- another. */
- typedef LorentzVector<typename BinaryOpTraits<T,U>::DivT> DivT;
-};
-
-/**
- * Template for multiplication by scalar
- */
-template <typename T, typename U>
-struct BinaryOpTraits<T, LorentzVector<U> > {
- /** The type resulting from multiplication of the template type with
- itself. */
- typedef LorentzVector<typename BinaryOpTraits<T,U>::MulT> MulT;
-};
-//@}
-
}
#undef ERROR_IF
#endif /* ThePEG_LorentzVector_H */
diff --git a/Vectors/ThreeVector.h b/Vectors/ThreeVector.h
--- a/Vectors/ThreeVector.h
+++ b/Vectors/ThreeVector.h
@@ -1,426 +1,421 @@
// -*- C++ -*-
//
// ThreeVector.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_ThreeVector_H
#define ThePEG_ThreeVector_H
/**
* @file ThreeVector.h contains the ThreeVector class. ThreeVector can be
* created with any unit type as template parameter. All basic
* mathematical operations are supported, as well as a subset of the
* CLHEP Vector3 functionality.
*/
#include "ThreeVector.fh"
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Utilities/UnitIO.h"
#include <cassert>
#include <cmath>
namespace ThePEG {
/**
* A 3-component vector. It can be created with any unit type
* as template parameter. All basic mathematical operations are
* supported, as well as a subset of the CLHEP Vector3
* functionality.
*/
template <typename Value>
class ThreeVector
-{
+{
private:
/// Value squared
- typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
- /// Value to the 4th power
- typedef typename BinaryOpTraits<Value2,Value2>::MulT Value4;
+ using Value2 = decltype(sqr(std::declval<Value>()));
public:
/** @name Constructors. */
//@{
ThreeVector()
: theX(), theY(), theZ() {}
ThreeVector(Value x, Value y, Value z)
: theX(x), theY(y), theZ(z) {}
template<typename ValueB>
ThreeVector(const ThreeVector<ValueB> & v)
: theX(v.x()), theY(v.y()), theZ(v.z()) {}
//@}
public:
/// @name Component access methods.
//@{
Value x() const { return theX; }
Value y() const { return theY; }
Value z() const { return theZ; }
//@}
/// @name Component set methods.
//@{
void setX(Value x) { theX = x; }
void setY(Value y) { theY = y; }
void setZ(Value z) { theZ = z; }
//@}
public:
/// Squared magnitude \f$x^2+y^2+z^2\f$.
Value2 mag2() const { return sqr(x()) + sqr(y()) + sqr(z()); }
/// Magnitude \f$\sqrt{x^2+y^2+z^2}\f$.
Value mag() const { return sqrt(mag2()); }
/// Squared transverse component \f$x^2+y^2\f$.
Value2 perp2() const { return sqr(x()) + sqr(y()); }
/// Transverse component \f$\sqrt{x^2+y^2}\f$.
Value perp() const { return sqrt(perp2()); }
/// Dot product.
template <typename U>
- typename BinaryOpTraits<Value,U>::MulT
- dot(const ThreeVector<U> & a) const {
+ auto dot(const ThreeVector<U> & a) const
+ -> decltype(x()*a.x())
+ {
return x()*a.x() + y()*a.y() + z()*a.z();
}
/// Squared transverse component with respect to the given axis.
template <typename U>
Value2 perp2(const ThreeVector<U> & p) const {
- typedef typename BinaryOpTraits<U,U>::MulT pSqType;
- const pSqType pMag2 = p.mag2();
- assert( pMag2 > pSqType() );
- typename BinaryOpTraits<Value,U>::MulT
- ss = this->dot(p);
+ const auto pMag2 = p.mag2();
+ assert( pMag2 > ZERO );
+ auto ss = this->dot(p);
Value2 ret = mag2() - sqr(ss)/pMag2;
- if ( ret <= Value2() )
- ret = Value2();
+ if ( ret <= ZERO )
+ ret = ZERO;
return ret;
}
/// Transverse component with respect to the given axis.
template <typename U>
Value perp(const ThreeVector<U> & p) const {
return sqrt(perp2(p));
}
/// @name Spherical coordinates.
//@{
/// Polar angle.
double theta() const {
- assert(!(x() == Value() && y() == Value() && z() == Value()));
+ assert(!(x() == ZERO && y() == ZERO && z() == ZERO));
return atan2(perp(),z());
}
/// Azimuthal angle.
double phi() const {
return atan2(y(),x());
}
/// Set the polar angle.
void setTheta(double th) {
double ma = mag();
double ph = phi();
setX(ma*sin(th)*cos(ph));
setY(ma*sin(th)*sin(ph));
setZ(ma*cos(th));
}
/// Set the azimuthal angle.
void setPhi(double ph) {
double xy = perp();
setX(xy*cos(ph));
setY(xy*sin(ph));
}
//@}
/// Parallel vector with unit length.
ThreeVector<double> unit() const {
Value2 mg2 = mag2();
- assert(mg2 > Value2());
+ assert(mg2 > ZERO);
Value mg = sqrt(mg2);
- return ThreeVector<double>(x()/mg, y()/mg, z()/mg);
+ return {x()/mg, y()/mg, z()/mg};
}
/// Orthogonal vector.
ThreeVector<Value> orthogonal() const {
Value xx = abs(x());
Value yy = abs(y());
Value zz = abs(z());
+ using TVec = ThreeVector<Value>;
if (xx < yy) {
- return xx < zz ? ThreeVector<Value>(Value(),z(),-y())
- : ThreeVector<Value>(y(),-x(),Value());
+ return xx < zz ? TVec{ZERO,z(),-y()} : TVec{y(),-x(),ZERO};
} else {
- return yy < zz ? ThreeVector<Value>(-z(),Value(),x())
- : ThreeVector<Value>(y(),-x(),Value());
+ return yy < zz ? TVec{-z(),ZERO,x()} : TVec{y(),-x(),ZERO};
}
}
/// Azimuthal angle difference, brought into the range \f$(-\pi,\pi]\f$.
template <typename U>
double deltaPhi (const ThreeVector<U> & v2) const {
double dphi = v2.phi() - phi();
if ( dphi > Constants::pi ) {
dphi -= Constants::twopi;
} else if ( dphi <= -Constants::pi ) {
dphi += Constants::twopi;
}
return dphi;
}
/**
* Apply a rotation.
* @param angle Rotation angle in radians.
* @param axis Rotation axis.
*/
template <typename U>
ThreeVector<Value> & rotate(double angle, const ThreeVector<U> & axis) {
if (angle == 0.0)
return *this;
const U ll = axis.mag();
- assert( ll > U() );
+ assert( ll > ZERO );
const double sa = sin(angle), ca = cos(angle);
const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll;
const Value xx = x(), yy = y(), zz = z();
setX((ca+(1-ca)*dx*dx) * xx
- +((1-ca)*dx*dy-sa*dz) * yy
- +((1-ca)*dx*dz+sa*dy) * zz
- );
+ +((1-ca)*dx*dy-sa*dz) * yy
+ +((1-ca)*dx*dz+sa*dy) * zz
+ );
setY(((1-ca)*dy*dx+sa*dz) * xx
- +(ca+(1-ca)*dy*dy) * yy
- +((1-ca)*dy*dz-sa*dx) * zz
- );
+ +(ca+(1-ca)*dy*dy) * yy
+ +((1-ca)*dy*dz-sa*dx) * zz
+ );
setZ(((1-ca)*dz*dx-sa*dy) * xx
- +((1-ca)*dz*dy+sa*dx) * yy
- +(ca+(1-ca)*dz*dz) * zz
- );
+ +((1-ca)*dz*dy+sa*dx) * yy
+ +(ca+(1-ca)*dz*dz) * zz
+ );
return *this;
}
/**
* Rotate the reference frame to a new z-axis.
*/
ThreeVector<Value> & rotateUz (const Axis & axis) {
Axis ax = axis.unit();
double u1 = ax.x();
double u2 = ax.y();
double u3 = ax.z();
double up = u1*u1 + u2*u2;
if (up>0) {
up = sqrt(up);
Value px = x(), py = y(), pz = z();
setX( (u1*u3*px - u2*py)/up + u1*pz );
setY( (u2*u3*px + u1*py)/up + u2*pz );
setZ( -up*px + u3*pz );
}
else if (u3 < 0.) {
setX(-x());
setZ(-z());
}
return *this;
}
/**
* Rotate from a reference frame to the z-axis.
*/
ThreeVector<Value> & rotateUzBack (const Axis & axis) {
Axis ax = axis.unit();
double u1 = ax.x();
double u2 = ax.y();
double u3 = ax.z();
double up = u1*u1 + u2*u2;
if (up>0) {
up = sqrt(up);
Value px = x(), py = y(), pz = z();
setX( ( u1*u3*px + u2*u3*py)/up - up*pz );
setY( (-u2*px + u1*py)/up );
setZ( u1*px + u2*py + u3*pz );
}
else if (u3 < 0.) {
setX(-x());
setZ(-z());
}
return *this;
}
/// Vector cross-product
template <typename U>
- ThreeVector<typename BinaryOpTraits<Value,U>::MulT>
- cross(const ThreeVector<U> & a) const {
- typedef ThreeVector<typename BinaryOpTraits<Value,U>::MulT> ResultT;
- return ResultT( y()*a.z()-z()*a.y(),
- -x()*a.z()+z()*a.x(),
- x()*a.y()-y()*a.x());
+ auto cross(const ThreeVector<U> & a) const
+ -> ThreeVector<decltype(y()*a.z())>
+ {
+ return { y()*a.z()-z()*a.y(),
+ -x()*a.z()+z()*a.x(),
+ x()*a.y()-y()*a.x() };
}
public:
/// @name Comparison operators.
//@{
bool operator==(const ThreeVector<Value> & a) const {
return (theX == a.x() && theY == a.y() && theZ == a.z());
}
bool operator!=(const ThreeVector<Value> & a) const {
return !(*this == a);
}
bool almostEqual(const ThreeVector<Value> & a, double threshold = 1e-04) const {
return ((std::abs(theX - a.x()) < threshold) && (std::abs(theY - a.y()) < threshold) && (std::abs(theZ - a.z()) < threshold));
}
bool almostUnequal(const ThreeVector<Value> & a, double threshold = 1e-04) const {
return ! this->almostEqual(a, threshold);
}
//@}
public:
/// @name Mathematical assignment operators.
//@{
ThreeVector<Value> & operator+=(const ThreeVector<Value> & a) {
theX += a.x();
theY += a.y();
theZ += a.z();
return *this;
}
ThreeVector<Value> & operator-=(const ThreeVector<Value> & a) {
theX -= a.x();
theY -= a.y();
theZ -= a.z();
return *this;
}
ThreeVector<Value> & operator*=(double a) {
theX *= a;
theY *= a;
theZ *= a;
return *this;
}
ThreeVector<Value> & operator/=(double a) {
theX /= a;
theY /= a;
theZ /= a;
return *this;
}
//@}
/// Cosine of the azimuthal angle between two vectors.
template <typename U>
double cosTheta(const ThreeVector<U> & q) const {
- typedef typename BinaryOpTraits<Value,U>::MulT
- ProdType;
- ProdType ptot = mag()*q.mag();
- assert( ptot > ProdType() );
+ auto ptot = mag()*q.mag();
+ assert( ptot > ZERO );
double arg = dot(q)/ptot;
if (arg > 1.0) arg = 1.0;
else if(arg < -1.0) arg = -1.0;
return arg;
}
/// Angle between two vectors.
template <typename U>
double angle(const ThreeVector<U> & v) const {
return acos(cosTheta(v));
}
private:
/// @name Vector components
//@{
Value theX;
Value theY;
Value theZ;
//@}
};
/// Stream output. Format \f$(x,y,z)\f$.
inline ostream &
operator<< (ostream & os, const ThreeVector<double> & v)
{
return os << '(' << v.x() << ',' << v.y() << ',' << v.z() << ')';
}
/// @name Basic mathematical operations
//@{
template <typename Value>
inline ThreeVector<Value>
operator+(ThreeVector<Value> a,
- const ThreeVector<Value> & b)
+ const ThreeVector<Value> & b)
{
return a += b;
}
template <typename Value>
inline ThreeVector<Value>
operator-(ThreeVector<Value> a,
- const ThreeVector<Value> & b)
+ const ThreeVector<Value> & b)
{
return a -= b;
}
template <typename Value>
inline ThreeVector<Value> operator-(const ThreeVector<Value> & v) {
- return ThreeVector<Value>(-v.x(),-v.y(),-v.z());
+ return {-v.x(),-v.y(),-v.z()};
}
template <typename Value>
inline ThreeVector<Value> operator*(ThreeVector<Value> v, double a) {
return v *= a;
}
template <typename Value>
inline ThreeVector<Value> operator*(double a, ThreeVector<Value> v) {
return v *= a;
}
template <typename ValueA, typename ValueB>
-inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(ValueB a, ThreeVector<ValueA> v) {
- typedef typename BinaryOpTraits<ValueA,ValueB>::MulT ResultT;
- return ThreeVector<ResultT>(a*v.x(), a*v.y(), a*v.z());
+inline auto operator*(ValueB a, ThreeVector<ValueA> v)
+-> ThreeVector<decltype(a*v.x())>
+{
+ return {a*v.x(), a*v.y(), a*v.z()};
}
template <typename ValueA, typename ValueB>
-inline ThreeVector<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(ThreeVector<ValueA> v, ValueB a) {
- return a*v;
+inline auto operator*(ThreeVector<ValueA> v, ValueB a)
+-> ThreeVector<decltype(v.x()*a)>
+{
+ return {v.x()*a, v.y()*a, v.z()*a};
}
//@}
/// Vector dot product.
template <typename ValueA, typename ValueB>
-inline typename BinaryOpTraits<ValueA,ValueB>::MulT
-operator*(const ThreeVector<ValueA> & a,
- const ThreeVector<ValueB> & b)
+inline auto operator*(const ThreeVector<ValueA> & a,
+ const ThreeVector<ValueB> & b)
+-> decltype(a.x()*b.x())
{
return a.dot(b);
}
/// A parallel vector with unit length.
template <typename Value>
ThreeVector<double> unitVector(const ThreeVector<Value> & v) {
return v.unit();
}
/** Output a ThreeVector with units to a stream. */
template <typename OStream, typename UT, typename Value>
void ounitstream(OStream & os, const ThreeVector<Value> & p, UT & u) {
os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u);
}
/** Input a ThreeVector with units from a stream. */
template <typename IStream, typename UT, typename Value>
void iunitstream(IStream & is, ThreeVector<Value> & p, UT & u) {
Value x, y, z;
is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u);
p = ThreeVector<Value>(x, y, z);
}
}
#endif /* ThePEG_ThreeVector_H */
diff --git a/Vectors/Transverse.h b/Vectors/Transverse.h
--- a/Vectors/Transverse.h
+++ b/Vectors/Transverse.h
@@ -1,262 +1,260 @@
// -*- C++ -*-
//
// Transverse.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Transverse_H
#define ThePEG_Transverse_H
// This is the declaration of the Transverse class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Vectors/Lorentz5Vector.h"
#include "Transverse.fh"
namespace ThePEG {
/**
* Transverse represents the transverse components of a
* LorentzVector. It inherits from
* <code>std::pair<Value,Value></code> and can be used
* anywhere such a pair is called for. It can also be created directly
* from a <code>ThreeVector</code>, <code>LorentzVector</code> and
* <code>Lorentz5Momentum</code>.
*
* @see Lorentz5Vector
*/
template <typename Value>
class Transverse: public pair<Value,Value> {
public:
/** Template argument typedef. */
- typedef typename BinaryOpTraits<Value,Value>::MulT Value2;
+ using Value2 = decltype(sqr(std::declval<Value>()));
+
/** Template argument typedef. */
- typedef pair<Value,Value> BasePair;
+ using BasePair = pair<Value,Value>;
public:
/** @name Constructors. */
//@{
/**
* Default constructor.
*/
Transverse() : BasePair(Value(), Value()) {}
/**
* Constructor from underlying representation.
*/
Transverse(const BasePair & p) : BasePair(p) {}
/**
* Constructor from x and y components.
*/
Transverse(Value x, Value y) : BasePair(x, y) {}
/**
* Constructor taking the transverse parts of a ThreeVector.
*/
Transverse(const ThreeVector<Value> & p) : BasePair(p.x(), p.y()) {}
/**
* Constructor taking the transverse parts of a LorentzVector.
*/
Transverse(const LorentzVector<Value> & p) : BasePair(p.x(), p.y()) {}
/**
* Constructor taking the transverse parts of a Lorentz5Vector.
*/
Transverse(const Lorentz5Vector<Value> & p) : BasePair(p.x(), p.y()) {}
//@}
/** @name Assignment operators. */
//@{
/**
* Assignment from underlying representation.
*/
const Transverse & operator=(const BasePair & p) {
BasePair::operator=(p);
return *this;
}
/**
* Assignment taking the transverse parts of a ThreeVector.
*/
const Transverse & operator=(const ThreeVector<Value> & p) {
BasePair::operator=(BasePair(p.x(), p.y()));
return *this;
}
/**
* Assignment taking the transverse parts of a LorentzVector.
*/
const Transverse & operator=(const LorentzVector<Value> & p) {
BasePair::operator=(BasePair(p.x(), p.y()));
return *this;
}
/**
* Assignment taking the transverse parts of a Lorentz5Vector.
*/
const Transverse & operator=(const Lorentz5Vector<Value> & p) {
BasePair::operator=(BasePair(p.x(), p.y()));
return *this;
}
//@}
/** @name Arithmetric operations */
//@{
/**
* Unary minus.
*/
Transverse operator-() const { return Transverse(-x(), -y()); }
/**
* Binary minus.
*/
Transverse operator-(const Transverse & pt) const {
return Transverse(x() - pt.x(), y() - pt.y());
}
/**
* Assign-subtract.
*/
Transverse & operator-=(const Transverse & pt) {
BasePair::first -= pt.x();
BasePair::second -= pt.y();
return *this;
}
/**
* Addition.
*/
Transverse operator+(const Transverse & pt) const {
return Transverse(x() + pt.x(), y() + pt.y());
}
/**
* Assign-add.
*/
Transverse & operator+=(const Transverse & pt) {
BasePair::first += pt.x();
BasePair::second += pt.y();
return *this;
}
/**
* Multiply-assign with a scalar.
*/
inline Transverse & operator*=(double a) {
BasePair::first *= a;
BasePair::second *= a;
return *this;
}
/**
* Divide-assign with a scalar.
*/
inline Transverse & operator/=(double a) {
BasePair::first /= a;
BasePair::second /= a;
return *this;
}
//@}
/** @name Access coordinates. */
//@{
/**
* The x-component.
*/
Value x() const { return BasePair::first; }
/**
* The y-component.
*/
Value y() const { return BasePair::second; }
/**
* The magnitude squared.
*/
Value2 pt2() const { return sqr(x()) + sqr(y()); }
/**
* The magnitude.
*/
Value pt() const { return sqrt(pt2()); }
/**
* The azimuth angle.
*/
double phi() const { return atan2(y(), x()); }
//@}
};
/** Output a Transverse with units to a stream. */
template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const Transverse<T> & p, UT & u) {
os << ounit(p.x(), u) << ounit(p.y(), u);
}
/** Input a Transverse with units from a stream. */
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, Transverse<T> & p, UT & u) {
T x, y;
is >> iunit(x, u) >> iunit(y, u);
p = Transverse<T>(x, y);
}
/** Multiply a Transverse with a number. */
template <typename Value>
inline Transverse<Value>
operator*(Transverse<Value> a, double b) {
return a *= b;
}
/** Multiply a number with a Transverse. */
template <typename Value>
inline Transverse<Value>
operator*(double b, Transverse<Value> a) {
return a *= b;
}
/** Multiply a (unitful) number with a Transverse. */
template <typename ValueA, typename ValueB>
-inline
-Transverse<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(ValueB a, const Transverse<ValueA> & v) {
- typedef typename BinaryOpTraits<ValueB,ValueA>::MulT ResultT;
- return Transverse<ResultT>(a*v.x(), a*v.y());
+inline auto operator*(ValueB a, const Transverse<ValueA> & v)
+-> Transverse<decltype(a*v.x())>
+{
+ return {a*v.x(), a*v.y()};
}
/** Multiply a Transverse with a (unitful) number. */
template <typename ValueA, typename ValueB>
-inline
-Transverse<typename BinaryOpTraits<ValueA,ValueB>::MulT>
-operator*(const Transverse<ValueA> & v, ValueB a) {
- typedef typename BinaryOpTraits<ValueB,ValueA>::MulT ResultT;
- return Transverse<ResultT>(a*v.x(), a*v.y());
+inline auto operator*(const Transverse<ValueA> & v, ValueB a)
+-> Transverse<decltype(a*v.x())>
+{
+ return {a*v.x(), a*v.y()};
}
/** Divide a Transverse by a number. */
template <typename Value>
inline Transverse<double>
operator/(const Transverse<Value> & v, Value a) {
- return Transverse<double>(v.x()/a, v.y()/a);
+ return {v.x()/a, v.y()/a};
}
/** Divide a Transverse by a (unitful) number. */
template <typename ValueA, typename ValueB>
-inline
-Transverse<typename BinaryOpTraits<ValueA,ValueB>::DivT>
-operator/(const Transverse<ValueA> & v, ValueB b) {
- typedef typename BinaryOpTraits<ValueA,ValueB>::DivT ResultT;
- return Transverse<ResultT>(v.x()/b, v.y()/b);
+inline auto operator/(const Transverse<ValueA> & v, ValueB b)
+-> Transverse<decltype(v.x()/b)>
+{
+ return {v.x()/b, v.y()/b};
}
}
#endif /* ThePEG_Transverse_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:02 PM (14 h, 59 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023561
Default Alt Text
(177 KB)

Event Timeline