Page MenuHomeHEPForge
Files F7879038

No OneTemporary
Actions

View Options

diff --git a/EventRecord/RhoDMatrix.h b/EventRecord/RhoDMatrix.h
--- a/EventRecord/RhoDMatrix.h
+++ b/EventRecord/RhoDMatrix.h
@@ -1,145 +1,146 @@
// -*- C++ -*-
//
// RhoDMatrix.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_RhoDMatrix_H
#define ThePEG_RhoDMatrix_H
// This is the declaration of the RhoDMatrix class.
#include "ThePEG/PDT/PDT.h"
#include "ThePEG/Helicity/HelicityDefinitions.h"
#include <cassert>
+#include <array>
namespace ThePEG {
/**
* The RhoDMatrix class is designed to implement the storage of the
* rho and D matrices which are required for the spin correlation
* algorithm. The matrix stores the spin as 2s+1.
*
* @author Peter Richardson
*
*/
class RhoDMatrix {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor with undefined spin.
*/
- RhoDMatrix() : _spin(), _ispin() {}
+ RhoDMatrix() = default;
/**
* Standard constructor giving the spin as 2s+1. The matrix starts out averaged,
* unless the second argument is false, when it is zeroed.
*/
RhoDMatrix(PDT::Spin inspin, bool average = true)
: _spin(inspin), _ispin(abs(int(inspin))) {
assert(_ispin <= MAXSPIN);
// initialize to average
- for(size_t ix=0; ix<_ispin; ++ix)
- for(size_t iy=0; iy<_ispin; ++iy)
- _matrix[ix][iy] = (average && ix==iy) ? 1./_ispin : 0.;
+ if ( average )
+ for(size_t ix=0; ix<_ispin; ++ix)
+ _matrix[ix][ix] = 1./_ispin;
}
//@}
public:
/** @name Access matrix elements. */
//@{
/**
* Return an element of the matrix.
*/
Complex operator() (size_t ix, size_t iy) const {
assert(ix < _ispin);
assert(iy < _ispin);
return _matrix[ix][iy];
}
/**
* Set an element of the matrix.
*/
Complex & operator() (size_t ix, size_t iy) {
assert(ix < _ispin);
assert(iy < _ispin);
return _matrix[ix][iy];
}
/**
* renormalise the matrix so it has unit trace
*/
void normalize() {
#ifndef NDEBUG
static const double epsa=1e-40, epsb=1e-10;
#endif
Complex norm = 0.;
for(size_t ix=0; ix<_ispin; ++ix)
norm += _matrix[ix][ix];
assert(norm.real() > epsa);
assert(norm.imag()/norm.real() < epsb);
double invnorm = 1./norm.real();
for(size_t ix=0; ix<_ispin; ++ix)
for(size_t iy=0; iy<_ispin; ++iy)
_matrix[ix][iy]*=invnorm;
}
//@}
/** @name Access the spin. */
//@{
/**
* Get the spin. The spin is returned as 2J+1 in units of hbar/2.
*/
PDT::Spin iSpin() const { return _spin; }
//@}
/**
* Output the spin density matrix for debugging purposes.
*/
friend ostream & operator<<(ostream & os, const RhoDMatrix & rd);
private:
/**
* 2s+1 for the particle.
*/
PDT::Spin _spin;
/**
* Storage of 2s+1 for speed.
*/
size_t _ispin;
/**
* Spin matrix size
*/
enum { MAXSPIN = 5 };
/**
* Storage for the matrix allowing up to spin 2 particles.
*/
// Deliberately not using vector<> to avoid calls to 'new'
// from this commonly used class.
- Complex _matrix[MAXSPIN][MAXSPIN];
+ std::array<std::array<Complex,MAXSPIN>,MAXSPIN> _matrix;
};
/** Output operator */
inline ostream & operator<<(ostream & os, const RhoDMatrix & rd) {
for (size_t ix = 0; ix < rd._ispin; ++ix) {
for (size_t iy = 0; iy < rd._ispin; ++iy)
os << rd._matrix[ix][iy] << " ";
os << '\n';
}
return os;
}
}
#endif /* ThePEG_RhoDMatrix_H */
diff --git a/Helicity/FermionSpinInfo.h b/Helicity/FermionSpinInfo.h
--- a/Helicity/FermionSpinInfo.h
+++ b/Helicity/FermionSpinInfo.h
@@ -1,211 +1,209 @@
// -*- C++ -*-
//
// FermionSpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_FermionSpinInfo_H
#define ThePEG_FermionSpinInfo_H
// This is the declaration of the FermionSpinInfo class.
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/Helicity/LorentzSpinor.h"
#include "FermionSpinInfo.fh"
+#include <array>
namespace ThePEG {
namespace Helicity {
/**
* The FermionSpinInfo class inherits from the SpinInfo class and
* implements the storage of the basis vectors for a spin-1/2
* particle. The basis states are the u-spinors for a particle and
* the v-spinors for an antiparticle. The barred spinors can be
* obtained from these.
*
* These basis states should be set by either matrixelements or
* decayers which are capable of generating spin correlation
* information.
*
* The basis states in the rest frame of the particles can then be
* accessed by decayers to produce the correct correlations.
*
* N.B. in our convention 0 is the \f$-\frac12\f$ helicity state and
* 1 is the \f$+\frac12\f$ helicity state.
*
* @author Peter Richardson
*/
class FermionSpinInfo: public SpinInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
- FermionSpinInfo() : SpinInfo(PDT::Spin1Half), _productionstates(2),
- _currentstates(2), _decaystates(2),
- _decaycalc(false) {}
+ FermionSpinInfo()
+ : SpinInfo(PDT::Spin1Half), _decaycalc(false) {}
/**
* Standard Constructor.
* @param p the production momentum.
* @param time true if the particle is time-like.
*/
FermionSpinInfo(const Lorentz5Momentum & p, bool time)
- : SpinInfo(PDT::Spin1Half, p, time),
- _productionstates(2), _currentstates(2), _decaystates(2),
- _decaycalc(false) {}
+ : SpinInfo(PDT::Spin1Half, p, time), _decaycalc(false) {}
//@}
public:
/** @name Set and get methods for the basis state. */
//@{
/**
* Set the basis state, this is production state.
* @param hel the helicity (0 or 1 as described above.)
* @param in the LorentzSpinor for the given helicity.
*/
void setBasisState(unsigned int hel,
const LorentzSpinor<SqrtEnergy> & in) const {
assert(hel<2);
_productionstates[hel] = in;
_currentstates [hel] = in;
}
/**
* Set the basis state for the decay.
* @param hel the helicity (0 or 1 as described above.)
* @param in the LorentzSpinor for the given helicity.
*/
void setDecayState(unsigned int hel,
const LorentzSpinor<SqrtEnergy> & in) const {
assert(hel<2);
_decaycalc = true;
_decaystates[hel] = in;
}
/**
* Get the basis state for the production for the given helicity, \a
* hel (which is 0 or 1 as described above.)
*/
const LorentzSpinor<SqrtEnergy> & getProductionBasisState(unsigned int hel) const {
assert(hel<2);
return _productionstates[hel];
}
/**
* Get the current basis state for the given helicity, \a
* hel (which is 0 or 1 as described above.)
*/
const LorentzSpinor<SqrtEnergy> & getCurrentBasisState(unsigned int hel) const {
assert(hel<2);
return _currentstates[hel];
}
/**
* Get the basis state for the decay for the given helicity, \a hel
* (which is 0 or 1 as described above.)
*/
const LorentzSpinor<SqrtEnergy> & getDecayBasisState(unsigned int hel) const {
assert(hel<2);
if(!_decaycalc) {
for(unsigned int ix=0;ix<2;++ix) _decaystates[ix]=_currentstates[ix];
_decaycalc=true;
}
return _decaystates[hel];
}
//@}
/**
* Perform a lorentz rotation of the spin information
*/
virtual void transform(const LorentzMomentum &,const LorentzRotation &);
public:
/**
* Standard Init function.
*/
static void Init();
/**
* Standard clone method.
*/
virtual EIPtr clone() const;
private:
/**
* Describe a concrete class without persistent data.
*/
static NoPIOClassDescription<FermionSpinInfo> initFermionSpinInfo;
/**
* Private and non-existent assignment operator.
*/
FermionSpinInfo & operator=(const FermionSpinInfo &);
private:
/**
* basis states in the frame in which the particle was produced
*/
- mutable vector<LorentzSpinor<SqrtEnergy> > _productionstates;
+ mutable std::array<LorentzSpinor<SqrtEnergy>,2> _productionstates;
/**
* basis states in the current frame of the particle
*/
- mutable vector<LorentzSpinor<SqrtEnergy> > _currentstates;
+ mutable std::array<LorentzSpinor<SqrtEnergy>,2> _currentstates;
/**
* basis states in the frame in which the particle decays
*/
- mutable vector<LorentzSpinor<SqrtEnergy> > _decaystates;
+ mutable std::array<LorentzSpinor<SqrtEnergy>,2> _decaystates;
/**
* True if the decay state has been set.
*/
mutable bool _decaycalc;
};
}
}
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the base class of
* FermionSpinInfo.
*/
template <>
struct BaseClassTrait<ThePEG::Helicity::FermionSpinInfo,1>
: public ClassTraitsType {
/** Typedef of the base class of FermionSpinInfo. */
typedef ThePEG::SpinInfo NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* FermionSpinInfo class and the shared object where it is defined.
*/
template <>
struct ClassTraits<ThePEG::Helicity::FermionSpinInfo>
: public ClassTraitsBase<ThePEG::Helicity::FermionSpinInfo> {
/**
* Return the class name.
*/
static string className() { return "ThePEG::Helicity::FermionSpinInfo"; }
};
/** @endcond */
}
#endif /* ThePEG_FermionSpinInfo_H */
diff --git a/Helicity/HelicityDefinitions.h b/Helicity/HelicityDefinitions.h
--- a/Helicity/HelicityDefinitions.h
+++ b/Helicity/HelicityDefinitions.h
@@ -1,54 +1,54 @@
// -*- C++ -*-
//
// HelicityDefinitions.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_HelicityDefinitions_H
#define THEPEG_HelicityDefinitions_H
// This is the declaration of the HelicityDefinitions class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/Utilities/Exception.h"
/** \file HelicityDefinitions.h
*
* This file contains enumerations used by LorentzSpinor and
* LorentzSpinorBar classes.
*
* @see LorentzSpinor
*
* @author Peter Richardson
*/
namespace ThePEG {
/**
* The Helicity namespace contains classes for spin representation
* classes in ThePEG.
*/
namespace Helicity {
/**
* Enumeration to specify spinor type.
*/
-enum SpinorType {
- u_spinortype, /**< u spinor. */
- v_spinortype, /**< v spinor. */
- unknown_spinortype /**< Undefined spinor type. */
+enum class SpinorType {
+ u, /**< u spinor. */
+ v, /**< v spinor. */
+ unknown /**< Undefined spinor type. */
};
/** @cond EXCEPTIONCLASSES */
/** Exception class used by Helicity classes to signal a logical error. */
class HelicityLogicalError: public Exception {};
/** Exception class used by Helicity classes to signal a inconsistencies. */
class HelicityConsistencyError: public Exception {};
/** @endcond */
}
}
#endif /* THEPEG_HelicityDefinitions_H */
diff --git a/Helicity/LorentzRSSpinor.h b/Helicity/LorentzRSSpinor.h
--- a/Helicity/LorentzRSSpinor.h
+++ b/Helicity/LorentzRSSpinor.h
@@ -1,414 +1,409 @@
// -*- C++ -*-
//
// LorentzRSSpinor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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 u_spinortype,
- * v_spinortype, unknown_spinortype. This information is intended
+ * 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 unknown_spinortypee is intended for
+ * 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 = unknown_spinortype) : _type(t) {
- for(unsigned int ix=0;ix<4;++ix)
- for(unsigned int iy=0;iy<4;++iy)
- _spin[ix][iy]=Value();
- }
+ 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=unknown_spinortype)
- : _type(t) {
- _spin[0][0]=a1;_spin[1][0]=a2;_spin[2][0]=a3;_spin[3][0]=a4;
- _spin[0][1]=b1;_spin[1][1]=b2;_spin[2][1]=b3;_spin[3][1]=b4;
- _spin[0][2]=c1;_spin[1][2]=c2;_spin[2][2]=c3;_spin[3][2]=c4;
- _spin[0][3]=d1;_spin[1][3]=d2;_spin[2][3]=d3;_spin[3][3]=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;
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;
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.
*/
- complex<Value> _spin[4][4];
+ 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,351 +1,346 @@
// -*- C++ -*-
//
// LorentzRSSpinorBar.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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 = unknown_spinortype) : _type(t) {
- for(unsigned int ix=0;ix<4;++ix)
- for(unsigned int iy=0;iy<4;++iy)
- _spin[ix][iy]=Value();
- }
+ 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=unknown_spinortype)
- : _type(t) {
- _spin[0][0]=a1;_spin[1][0]=a2;_spin[2][0]=a3;_spin[3][0]=a4;
- _spin[0][1]=b1;_spin[1][1]=b2;_spin[2][1]=b3;_spin[3][1]=b4;
- _spin[0][2]=c1;_spin[1][2]=c2;_spin[2][2]=c3;_spin[3][2]=c4;
- _spin[0][3]=d1;_spin[1][3]=d2;_spin[2][3]=d3;_spin[3][3]=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;
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.
*/
- complex<Value> _spin[4][4];
+ 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,500 +1,494 @@
// -*- C++ -*-
//
// LorentzSpinor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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 u_spinortype,
- * v_spinortype, unknown_spinortype. This information is intended
+ * 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 unknown_spinortypee is intended for
+ * 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 = unknown_spinortype) : _type(t) {
- for(unsigned int ix=0;ix<4;++ix) _spin[ix]=Value();
- }
+ 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 = unknown_spinortype) : _type(s) {
- _spin[0]=a;
- _spin[1]=b;
- _spin[2]=c;
- _spin[3]=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;
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;
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;
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;
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;
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;
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;
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 {
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 {
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();
}
/**
* 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();
}
/**
* 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())
+ 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;
LorentzTensor<ResultT> output;
complex<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;
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.
*/
- complex<Value> _spin[4];
+ std::array<complex<Value>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzSpinor.tcc"
#endif
#endif
diff --git a/Helicity/LorentzSpinor.tcc b/Helicity/LorentzSpinor.tcc
--- a/Helicity/LorentzSpinor.tcc
+++ b/Helicity/LorentzSpinor.tcc
@@ -1,112 +1,112 @@
// -*- C++ -*-
//
// LorentzSpinor.tcc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined member
// functions of the LorentzSpinor class.
//
// Author: Peter Richardson
//
#include "LorentzSpinor.h"
#include "LorentzSpinorBar.h"
using namespace ThePEG;
using namespace ThePEG::Helicity;
// return the barred spinor
template <typename Value>
LorentzSpinorBar<Value> LorentzSpinor<Value>::bar() const
{
complex<Value> output[4];
// HELAS
output[0] = conj(_spin[2]);
output[1] = conj(_spin[3]);
output[2] = conj(_spin[0]);
output[3] = conj(_spin[1]);
return LorentzSpinorBar<Value>(output[0],output[1],output[2],output[3],_type);
}
// boost the spinor
template <typename Value>
LorentzSpinor<Value> & LorentzSpinor<Value>::boost(double bx,double by,double bz)
{
// work out beta and chi
double beta=sqrt(bx*bx+by*by+bz*bz);
double chi = atanh(beta);
double sinhchi = sinh(0.5*chi)/beta, coshchi = cosh(0.5*chi);
// calculate the new spinor
complex<Value> out[4];
Complex ii(0.,1.);
Complex nxminy=bx-ii*by;
Complex nxpiny=bx+ii*by;
out[0] = coshchi*_spin[0]+sinhchi*(-bz*_spin[0]-nxminy*_spin[1]);
out[1] = coshchi*_spin[1]+sinhchi*(+bz*_spin[1]-nxpiny*_spin[0]);
out[2] = coshchi*_spin[2]+sinhchi*(+bz*_spin[2]+nxminy*_spin[3]);
out[3] = coshchi*_spin[3]+sinhchi*(-bz*_spin[3]+nxpiny*_spin[2]);
for(unsigned int ix=0;ix<4;++ix) _spin[ix]=out[ix];
return *this;
}
// boost the spinor
template <typename Value>
LorentzSpinor<Value> & LorentzSpinor<Value>::boost(const Boost & boostv)
{
double beta = boostv.mag();
double bx=boostv.x(),by=boostv.y(),bz=boostv.z();
double chi = atanh(beta);
double sinhchi = sinh(0.5*chi)/beta, coshchi = cosh(0.5*chi);
complex<Value> out[4];
Complex ii(0.,1.);
Complex nxminy=bx-ii*by;
Complex nxpiny=bx+ii*by;
out[0] = coshchi*_spin[0]+sinhchi*(-bz*_spin[0]-nxminy*_spin[1]);
out[1] = coshchi*_spin[1]+sinhchi*(+bz*_spin[1]-nxpiny*_spin[0]);
out[2] = coshchi*_spin[2]+sinhchi*(+bz*_spin[2]+nxminy*_spin[3]);
out[3] = coshchi*_spin[3]+sinhchi*(-bz*_spin[3]+nxpiny*_spin[2]);
for(unsigned int ix=0;ix<4;++ix) _spin[ix]=out[ix];
return *this;
}
// general transform
template <typename Value>
LorentzSpinor<Value> & LorentzSpinor<Value>::
transform(const SpinHalfLorentzRotation & r) {
unsigned int ix,iy;
complex<Value> out[4];
for(ix=0;ix<4;++ix) {
out[ix]=complex<Value>();
for(iy=0;iy<4;++iy) out[ix]+=r(ix,iy)*_spin[iy];
}
for(ix=0;ix<4;++ix) _spin[ix]=out[ix];
return *this;
}
// conjugation
template <typename Value>
LorentzSpinor<Value> LorentzSpinor<Value>::conjugate() const {
SpinorType new_type;
switch(_type) {
- case u_spinortype:
- new_type=v_spinortype;
+ case SpinorType::u:
+ new_type=SpinorType::v;
break;
- case v_spinortype:
- new_type=u_spinortype;
+ case SpinorType::v:
+ new_type=SpinorType::u;
break;
- case unknown_spinortype:
+ case SpinorType::unknown:
default:
- new_type=unknown_spinortype;
+ new_type=SpinorType::unknown;
break;
}
return LorentzSpinor<Value>( conj(_spin[3]),
-conj(_spin[2]),
-conj(_spin[1]),
+conj(_spin[0]),
new_type);
}
diff --git a/Helicity/LorentzSpinorBar.h b/Helicity/LorentzSpinorBar.h
--- a/Helicity/LorentzSpinorBar.h
+++ b/Helicity/LorentzSpinorBar.h
@@ -1,240 +1,233 @@
// -*- C++ -*-
//
// LorentzSpinorBar.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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 = unknown_spinortype) : _type(t) {
- for(unsigned int ix=0;ix<4;++ix) _spin[ix]=Value();
- }
+ 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 = unknown_spinortype)
- : _type(t) {
- _spin[0]=a;
- _spin[1]=b;
- _spin[2]=c;
- _spin[3]=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;
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.
*/
- complex<Value> _spin[4];
+ std::array<complex<Value>,4> _spin;
};
}
}
#ifndef ThePEG_TEMPLATES_IN_CC_FILE
#include "LorentzSpinorBar.tcc"
#endif
#endif
diff --git a/Helicity/LorentzSpinorBar.tcc b/Helicity/LorentzSpinorBar.tcc
--- a/Helicity/LorentzSpinorBar.tcc
+++ b/Helicity/LorentzSpinorBar.tcc
@@ -1,109 +1,109 @@
// -*- C++ -*-
//
// LorentzSpinorBar.tcc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined member
// functions of the LorentzSpinorBar class.
//
// Author: Peter Richardson
//
#include "LorentzSpinorBar.h"
#include "LorentzSpinor.h"
using namespace ThePEG;
using namespace ThePEG::Helicity;
// return the unbarred spinor
template <typename Value>
LorentzSpinor<Value> LorentzSpinorBar<Value>::bar() const
{
complex<Value> output[4];
// HELAS
output[0] = conj(_spin[2]);
output[1] = conj(_spin[3]);
output[2] = conj(_spin[0]);
output[3] = conj(_spin[1]);
return LorentzSpinor<Value>(output[0],output[1],output[2],output[3],_type);
}
template <typename Value>
LorentzSpinorBar<Value> & LorentzSpinorBar<Value>::boost(double bx,double by,double bz)
{
// work out beta and chi
double beta=sqrt(bx*bx+by*by+bz*bz);
double chi = atanh(beta);
double sinhchi = sinh(0.5*chi)/beta, coshchi = cosh(0.5*chi);
// calculate the new spinor
complex<Value> out[4];
Complex ii(0.,1.);
Complex nxminy=bx-ii*by;
Complex nxpiny=bx+ii*by;
out[0] = coshchi*_spin[0]+sinhchi*(-bz*_spin[0]-nxpiny*_spin[1]);
out[1] = coshchi*_spin[1]+sinhchi*(+bz*_spin[1]-nxminy*_spin[0]);
out[2] = coshchi*_spin[2]+sinhchi*(+bz*_spin[2]+nxpiny*_spin[3]);
out[3] = coshchi*_spin[3]+sinhchi*(-bz*_spin[3]+nxminy*_spin[2]);
for(unsigned int ix=0;ix<4;++ix){_spin[ix]=out[ix];}
return *this;
}
template <typename Value>
LorentzSpinorBar<Value> & LorentzSpinorBar<Value>::boost(const Boost & boostv)
{
double beta = boostv.mag();
double bx=boostv.x(),by=boostv.y(),bz=boostv.z();
double chi = atanh(beta);
double sinhchi = sinh(0.5*chi)/beta, coshchi = cosh(0.5*chi);
complex<Value> out[4];
Complex ii(0.,1.);
Complex nxminy=bx-ii*by;
Complex nxpiny=bx+ii*by;
out[0] = coshchi*_spin[0]+sinhchi*(-bz*_spin[0]-nxpiny*_spin[1]);
out[1] = coshchi*_spin[1]+sinhchi*(+bz*_spin[1]-nxminy*_spin[0]);
out[2] = coshchi*_spin[2]+sinhchi*(+bz*_spin[2]+nxpiny*_spin[3]);
out[3] = coshchi*_spin[3]+sinhchi*(-bz*_spin[3]+nxminy*_spin[2]);
for(unsigned int ix=0;ix<4;++ix){_spin[ix]=out[ix];}
return *this;
}
// general transform
template <typename Value>
LorentzSpinorBar<Value> & LorentzSpinorBar<Value>::transform(const SpinHalfLorentzRotation & r)
{
unsigned int ix,iy;
SpinHalfLorentzRotation t(r.inverse());
complex<Value> out[4];
for(ix=0;ix<4;++ix)
{
out[ix]=complex<Value>();
for(iy=0;iy<4;++iy){out[ix]+=_spin[iy]*t(iy,ix);}
}
for(ix=0;ix<4;++ix){_spin[ix]=out[ix];}
return *this;
}
// conjugation
template <typename Value>
LorentzSpinorBar<Value> LorentzSpinorBar<Value>::conjugate() const {
SpinorType new_type;
switch(_type) {
- case u_spinortype:
- new_type=v_spinortype;
+ case SpinorType::u:
+ new_type=SpinorType::v;
break;
- case v_spinortype:
- new_type=u_spinortype;
+ case SpinorType::v:
+ new_type=SpinorType::u;
break;
- case unknown_spinortype:
+ case SpinorType::unknown:
default:
- new_type=unknown_spinortype;
+ new_type=SpinorType::unknown;
break;
}
return LorentzSpinorBar<Value>(-conj(_spin[3]),+conj(_spin[2]),
+conj(_spin[1]),-conj(_spin[0]),new_type);
}
diff --git a/Helicity/LorentzTensor.h b/Helicity/LorentzTensor.h
--- a/Helicity/LorentzTensor.h
+++ b/Helicity/LorentzTensor.h
@@ -1,511 +1,505 @@
// -*- C++ -*-
//
// LorentzTensor.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 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/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() {
- for(unsigned int ix=0;ix<4;++ix)
- for(unsigned int iy=0;iy<4;++iy)
- _tensor[ix][iy]=Value();
- }
+ 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){
- _tensor[0][0]=xx;_tensor[0][1]=xy;_tensor[0][2]=xz;_tensor[0][3]=xt;
- _tensor[1][0]=yx;_tensor[1][1]=yy;_tensor[1][2]=yz;_tensor[1][3]=yt;
- _tensor[2][0]=zx;_tensor[2][1]=zy;_tensor[2][2]=zz;_tensor[2][3]=zt;
- _tensor[3][0]=tx;_tensor[3][1]=ty;_tensor[3][2]=tz;_tensor[3][3]=tt;
- }
-
+ 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) {
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;
}
}
*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()));
}
//@}
/** @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);
/**
* 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());
}
/**
* 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());
}
/**
* 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]);
output.setY(vec.t()*_tensor[3][1]-vec.x()*_tensor[0][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]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[0][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]);
output.setY(vec.t()*_tensor[1][3]-vec.x()*_tensor[1][0]-
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]);
output.setT(vec.t()*_tensor[3][3]-vec.x()*_tensor[3][0]-
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;
output.setX(vec.t()*_tensor[3][0]-vec.x()*_tensor[0][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]);
output.setZ(vec.t()*_tensor[3][2]-vec.x()*_tensor[0][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]);
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;
output.setX(vec.t()*_tensor[0][3]-vec.x()*_tensor[0][0]-
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]);
output.setZ(vec.t()*_tensor[2][3]-vec.x()*_tensor[2][0]-
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]);
return output;
}
//@}
private:
/**
* The components.
*/
- complex<Value> _tensor[4][4];
+ 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(),
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());
}
/**
* 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));
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));
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;
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/RSFermionSpinInfo.h b/Helicity/RSFermionSpinInfo.h
--- a/Helicity/RSFermionSpinInfo.h
+++ b/Helicity/RSFermionSpinInfo.h
@@ -1,207 +1,204 @@
// -*- C++ -*-
//
// RSFermionSpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_RSFermionSpinInfo_H
#define THEPEG_RSFermionSpinInfo_H
// This is the declaration of the RSFermionSpinInfo class.
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/Helicity/LorentzRSSpinor.h"
#include "RSFermionSpinInfo.fh"
+#include <array>
namespace ThePEG {
namespace Helicity {
/**
* The RSFermionSpinInfo class inherits from the SpinInfo class and
* implements the storage of the basis vector for a spin-3/2 particle.
* The basis states are the vector u spinors for a particle and the vector
* v-spinors for an antiparticle. The barred spinors can be obtained from these.
*
* These basis states should be set by either the matrixelements or decayers
* which are capable of generating spin correlation information.
*
* The basis states in the rest frame of the particles can then be accessed by
* the decayers to produce the correct correlations.
*
* N.B. in our convention 0 is the \f$-\frac32\f$ helicity state,
* 1 is the \f$-\frac12\f$ helicity state,
* 2 is the \f$+\frac12\f$ helicity state,
* 3 is the \f$+\frac32\f$ helicity state.
*
* @see SpinInfo
*
* \author Peter Richardson
*
*/
class RSFermionSpinInfo: public SpinInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
- RSFermionSpinInfo() : SpinInfo(PDT::Spin3Half), _productionstates(4),
- _decaystates(4), _currentstates(4),
- _decaycalc(false) {}
+ RSFermionSpinInfo() : SpinInfo(PDT::Spin3Half), _decaycalc(false) {}
/**
* Standard Constructor.
* @param p the production momentum.
* @param time true if the particle is time-like.
*/
RSFermionSpinInfo(const Lorentz5Momentum & p,bool time)
- : SpinInfo(PDT::Spin3Half, p, time),
- _productionstates(4), _decaystates(4), _currentstates(4),
- _decaycalc(false) {}
+ : SpinInfo(PDT::Spin3Half, p, time), _decaycalc(false) {}
//@}
public:
/** @name Set and get methods for the basis state. */
//@{
/**
* Set the basis state, this is production state.
* @param hel the helicity (0,1,2,3 as described above.)
* @param in the LorentzRSSpinor for the given helicity.
*/
void setBasisState(unsigned int hel,
const LorentzRSSpinor<SqrtEnergy> & in) const {
assert(hel<4);
_productionstates[hel] = in;
_currentstates [hel] = in;
}
/**
* Set the basis state for the decay.
* @param hel the helicity (0,1,2,3 as described above.)
* @param in the LorentzRSSpinor for the given helicity.
*/
void setDecayState(unsigned int hel,
const LorentzRSSpinor<SqrtEnergy> & in) const {
assert(hel<4);
_decaycalc = true;
_decaystates[hel] = in;
}
/**
* Get the basis state for the production for the given helicity, \a
* hel (0,1,2,3 as described above.)
*/
const LorentzRSSpinor<SqrtEnergy> & getProductionBasisState(unsigned int hel) const {
assert(hel<4);
return _productionstates[hel];
}
/**
* Get the basis state for the decay for the given helicity, \a hel
* (0,1,2,3 as described above.)
*/
const LorentzRSSpinor<SqrtEnergy> & getDecayBasisState(unsigned int hel) const {
assert(hel<4);
if(!_decaycalc) {
for(unsigned int ix=0;ix<4;++ix) _decaystates[ix]=_currentstates[ix];
_decaycalc=true;
}
return _decaystates[hel];
}
/**
* Perform a lorentz rotation of the spin information
*/
virtual void transform(const LorentzMomentum &,const LorentzRotation &);
//@}
public:
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
/**
* Standard clone method.
*/
virtual EIPtr clone() const;
private:
/**
* Describe a concrete class without persistent data.
*/
static NoPIOClassDescription<RSFermionSpinInfo> initRSFermionSpinInfo;
/**
* Private and non-existent assignment operator.
*/
RSFermionSpinInfo & operator=(const RSFermionSpinInfo &);
private:
/**
* Basis states in the frame in which the particle was produced.
*/
- mutable vector<LorentzRSSpinor<SqrtEnergy> > _productionstates;
+ mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _productionstates;
/**
* Basis states in the frame in which the particle decays.
*/
- mutable vector<LorentzRSSpinor<SqrtEnergy> > _decaystates;
+ mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _decaystates;
/**
* Basis states in the current frame of the particle
*/
- mutable vector<LorentzRSSpinor<SqrtEnergy> > _currentstates;
+ mutable std::array<LorentzRSSpinor<SqrtEnergy>,4> _currentstates;
/**
* True if the decay state has been set.
*/
mutable bool _decaycalc;
};
}
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* The following template specialization informs ThePEG about the
* base class of RSFermionSpinInfo.
*/
template <>
struct BaseClassTrait<ThePEG::Helicity::RSFermionSpinInfo,1> {
/** Typedef of the base class of RSFermionSpinInfo. */
typedef ThePEG::SpinInfo NthBase;
};
/**
* The following template specialization informs ThePEG about the
* name of this class and the shared object where it is defined.
*/
template <>
struct ClassTraits<ThePEG::Helicity::RSFermionSpinInfo>
: public ClassTraitsBase<ThePEG::Helicity::RSFermionSpinInfo> {
/**
* Return the class name.
*/
static string className() { return "ThePEG::Helicity::RSFermionSpinInfo"; }
};
/** @endcond */
}
#endif /* THEPEG_RSFermionSpinInfo_H */
diff --git a/Helicity/TensorSpinInfo.h b/Helicity/TensorSpinInfo.h
--- a/Helicity/TensorSpinInfo.h
+++ b/Helicity/TensorSpinInfo.h
@@ -1,205 +1,202 @@
// -*- C++ -*-
//
// TensorSpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_TensorSpinInfo_H
#define THEPEG_TensorSpinInfo_H
// This is the declaration of the TensorSpinInfo class.
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/Helicity/LorentzTensor.h"
#include "TensorSpinInfo.fh"
// #include "TensorSpinInfo.xh"
+#include <array>
namespace ThePEG {
namespace Helicity {
/**
* The TensorSpinInfo class is the implementation of the spin
* information for tensor particles. It inherits from the SpinInfo
* class and implements the storage of the basis tensors.
*
* These basis states should be set by either matrix elements or
* decayers which are capable of generating spin correlation
* information.
*
* The basis states in the rest frame of the particles can then be
* accessed by decayers to produce the correct correlation.
*
* N.B. in our convention 0 is the \f$-2\f$ helicity state,
* 1 is the \f$-1\f$ helicity state,
* 2 is the \f$0\f$ helicity state,
* 3 is the \f$+1\f$ helicity state and
* 4 is the \f$+2\f$ helicity state.
*
* @author Peter Richardson
*
*/
class TensorSpinInfo: public SpinInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
- TensorSpinInfo() : SpinInfo(PDT::Spin2),_productionstates(5),
- _decaystates(5), _currentstates(5),
- _decaycalc(false) {}
+ TensorSpinInfo() : SpinInfo(PDT::Spin2), _decaycalc(false) {}
/**
* Standard Constructor.
* @param p the production momentum.
* @param time true if the particle is time-like.
*/
- TensorSpinInfo(const Lorentz5Momentum & p,bool time)
- : SpinInfo(PDT::Spin2, p, time),
- _productionstates(5), _decaystates(5), _currentstates(5),
- _decaycalc(false) {}
+ TensorSpinInfo(const Lorentz5Momentum & p, bool time)
+ : SpinInfo(PDT::Spin2, p, time), _decaycalc(false) {}
//@}
public:
/** @name Access the basis states. */
//@{
/**
* Set the basis state, this is production state.
* @param hel the helicity (0,1,2,3,4 as described above.)
* @param in the LorentzTensor for the given helicity.
*/
void setBasisState(unsigned int hel, LorentzTensor<double> in) const {
assert(hel<5);
_productionstates[hel]=in;
_currentstates [hel]=in;
}
/**
* Set the basis state for the decay.
* @param hel the helicity (0,1,2,3,4 as described above.)
* @param in the LorentzTensor for the given helicity.
*/
void setDecayState(unsigned int hel, LorentzTensor<double> in) const {
assert(hel<5);
_decaycalc = true;
_decaystates[hel] = in;
}
/**
* Get the basis state for the production for the given helicity, \a
* hel (0,1,2,3,4 as described above.)
*/
const LorentzTensor<double> & getProductionBasisState(unsigned int hel) const {
assert(hel<5);
return _productionstates[hel];
}
/**
* Get the basis state for the decay for the given helicity, \a hel
* (0,1,2,3,4 as described above.)
*/
const LorentzTensor<double> & getDecayBasisState(unsigned int hel) const {
assert(hel<5);
if(!_decaycalc) {
for(unsigned int ix=0;ix<5;++ix)
_decaystates[ix]=_currentstates[ix].conjugate();
_decaycalc=true;
}
return _decaystates[hel];
}
//@}
/**
* Perform a lorentz rotation of the spin information
*/
virtual void transform(const LorentzMomentum &,const LorentzRotation &);
public:
/**
* Standard Init function.
*/
static void Init();
/**
* Standard clone method.
*/
virtual EIPtr clone() const;
private:
/**
* Describe a concrete class without persistent data.
*/
static NoPIOClassDescription<TensorSpinInfo> initTensorSpinInfo;
/**
* Private and non-existent assignment operator.
*/
TensorSpinInfo & operator=(const TensorSpinInfo &);
private:
/**
* Basis states in the frame in which the particle was produced.
*/
- mutable vector<LorentzTensor<double> > _productionstates;
+ mutable std::array<LorentzTensor<double>,5> _productionstates;
/**
* Basis states in the frame in which the particle decays.
*/
- mutable vector<LorentzTensor<double> > _decaystates;
+ mutable std::array<LorentzTensor<double>,5> _decaystates;
/**
* Basis states in the current frame of the particle
*/
- mutable vector<LorentzTensor<double> > _currentstates;
+ mutable std::array<LorentzTensor<double>,5> _currentstates;
/**
* True if the decay state has been set.
*/
mutable bool _decaycalc;
};
}
}
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the base class of
* TensorSpinInfo.
*/
template <>
struct BaseClassTrait<ThePEG::Helicity::TensorSpinInfo,1>
: public ClassTraitsType {
/** Typedef of the base class of ScalarSpinInfo. */
typedef ThePEG::SpinInfo NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* TensorSpinInfo class and the shared object where it is defined.
*/
template <>
struct ClassTraits<ThePEG::Helicity::TensorSpinInfo>
: public ClassTraitsBase<ThePEG::Helicity::TensorSpinInfo> {
/**
* Return the class name.
*/
static string className() { return "ThePEG::Helicity::TensorSpinInfo"; }
};
/** @endcond */
}
#endif /* THEPEG_TensorSpinInfo_H */
diff --git a/Helicity/VectorSpinInfo.h b/Helicity/VectorSpinInfo.h
--- a/Helicity/VectorSpinInfo.h
+++ b/Helicity/VectorSpinInfo.h
@@ -1,205 +1,201 @@
// -*- C++ -*-
//
// VectorSpinInfo.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_VectorSpinInfo_H
#define THEPEG_VectorSpinInfo_H
// This is the declaration of the VectorSpinInfo class.
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ThePEG/Helicity/LorentzPolarizationVector.h"
#include "VectorSpinInfo.fh"
namespace ThePEG {
namespace Helicity {
/**
* The VectorSpinInfo class is the implementation of the spin
* information for vector particles. It inherits from the SpinInfo
* class and implements the storage of the basis vectors.
*
* These basis states should be set by either matrixelements or
* decayers which are capable of generating spin correlation
* information.
*
* The basis states in the rest frame of the particles can then be
* accessed by decayers to produce the correct correlation.
*
* N.B. in our convention 0 is the \f$-1\f$ helicity state,
* 1 is the \f$0\f$ helicity state and
* 2 is the \f$+1\f$ helicity state.
*
* @author Peter Richardson
*
*/
class VectorSpinInfo: public SpinInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
- VectorSpinInfo() : SpinInfo(PDT::Spin1), _productionstates(3),
- _decaystates(3), _currentstates(3),
- _decaycalc(false) {}
+ VectorSpinInfo() : SpinInfo(PDT::Spin1), _decaycalc(false) {}
/**
* Standard Constructor.
* @param p the production momentum.
* @param time true if the particle is time-like.
*/
VectorSpinInfo(const Lorentz5Momentum & p, bool time)
- : SpinInfo(PDT::Spin1, p, time),
- _productionstates(3), _decaystates(3), _currentstates(3),
- _decaycalc(false) {}
+ : SpinInfo(PDT::Spin1, p, time), _decaycalc(false) {}
//@}
public:
/** @name Set and get methods for the basis state. */
//@{
/**
* Set the basis state, this is production state.
* @param hel the helicity (0,1,2 as described above.)
* @param in the LorentzPolarizationVector for the given helicity.
*/
void setBasisState(unsigned int hel,
const LorentzPolarizationVector & in) const {
assert(hel<3);
_productionstates[hel] = in;
_currentstates [hel] = in;
}
/**
* Set the basis state for the decay.
* @param hel the helicity (0,1,2 as described above.)
* @param in the LorentzPolarizationVector for the given helicity.
*/
void setDecayState(unsigned int hel,
const LorentzPolarizationVector & in) const {
assert(hel<3);
_decaycalc = true;
_decaystates[hel] = in;;
}
/**
* Get the basis state for the production for the given helicity, \a
* hel (0,1,2 as described above.)
*/
const LorentzPolarizationVector & getProductionBasisState(unsigned int hel) const {
assert(hel<3);
return _productionstates[hel];
}
/**
* Get the basis state for the decay for the given helicity, \a hel
* (0,1,2 as described above.)
*/
const LorentzPolarizationVector & getDecayBasisState(unsigned int hel) const {
assert(hel<3);
if(!_decaycalc) {
for(unsigned int ix=0;ix<3;++ix)
_decaystates[ix]=_currentstates[ix].conjugate();
_decaycalc=true;
}
// return the basis function
return _decaystates[hel];
}
//@}
/**
* Perform a Lorentz rotation of the spin information
*/
virtual void transform(const LorentzMomentum &,const LorentzRotation & );
public:
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
/**
* Standard clone method.
*/
virtual EIPtr clone() const;
private:
/**
* Describe a concrete class without persistent data.
*/
static NoPIOClassDescription<VectorSpinInfo> initVectorSpinInfo;
/**
* Private and non-existent assignment operator.
*/
VectorSpinInfo & operator=(const VectorSpinInfo &);
private:
/**
* Basis states in the frame in which the particle was produced.
*/
- mutable vector<LorentzPolarizationVector> _productionstates;
+ mutable std::array<LorentzPolarizationVector,3> _productionstates;
/**
* Basis states in the frame in which the particle decays.
*/
- mutable vector<LorentzPolarizationVector> _decaystates;
+ mutable std::array<LorentzPolarizationVector,3> _decaystates;
/**
* Basis states in the current frame of the particle
*/
- mutable vector<LorentzPolarizationVector> _currentstates;
+ mutable std::array<LorentzPolarizationVector,3> _currentstates;
/**
* True if the decay state has been set.
*/
mutable bool _decaycalc;
};
}
}
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the base class of
* VectorSpinInfo.
*/
template <>
struct BaseClassTrait<ThePEG::Helicity::VectorSpinInfo,1>
: public ClassTraitsType {
/** Typedef of the base class of VectorSpinInfo. */
typedef ThePEG::SpinInfo NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* VectorSpinInfo class and the shared object where it is defined.
*/
template <>
struct ClassTraits<ThePEG::Helicity::VectorSpinInfo>
: public ClassTraitsBase<ThePEG::Helicity::VectorSpinInfo> {
/**
* Return the class name.
*/
static string className() { return "ThePEG::Helicity::VectorSpinInfo"; }
};
/** @endcond */
}
#endif /* THEPEG_VectorSpinInfo_H */
diff --git a/Helicity/Vertex/Vector/FFVVertex.cc b/Helicity/Vertex/Vector/FFVVertex.cc
--- a/Helicity/Vertex/Vector/FFVVertex.cc
+++ b/Helicity/Vertex/Vector/FFVVertex.cc
@@ -1,473 +1,473 @@
// -*- C++ -*-
//
// FFVVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FFVVertex class.
//
#include "FFVVertex.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
using namespace Helicity;
// Definition of the static class description member
AbstractNoPIOClassDescription<FFVVertex> FFVVertex::initFFVVertex;
void FFVVertex::Init() {
static ClassDocumentation<FFVVertex> documentation
("The FFVVertex class implements the helicity amplitude"
"calculations for a fermion-fantifermion gauge boson vertex. Any "
"implementation of such a vertex should inherit from in and implement"
" the virtual setCoupling member to calculate the coupling");
}
// evalulate the full vertex
Complex FFVVertex::evaluate(Energy2 q2,
const SpinorWaveFunction & sp,
const SpinorBarWaveFunction & sbar,
const VectorWaveFunction & vec) {
// first calculate the couplings
if(kinematics()) calculateKinematics(sp.momentum(),sbar.momentum(),vec.momentum());
setCoupling(q2,sp.particle(),sbar.particle(),vec.particle());
Complex ii(0.,1.);
// useful combinations of the polarization vector components
Complex e0p3=vec.t()+vec.z();
Complex e0m3=vec.t()-vec.z();
Complex e1p2=vec.x()+ii*vec.y();
Complex e1m2=vec.x()-ii*vec.y();
Complex vertex(0.);
if(_left!=0.) {
vertex += _left*(+sbar.s3()*(sp.s1()*e0p3+sp.s2()*e1m2)
+sbar.s4()*(sp.s1()*e1p2+sp.s2()*e0m3));
}
// then the right piece (often not needed eg W vertex)
if(_right!=0.) {
vertex += _right*(+sbar.s1()*(sp.s3()*e0m3-sp.s4()*e1m2)
-sbar.s2()*(sp.s3()*e1p2-sp.s4()*e0p3));
}
vertex*=ii;
return vertex*norm();
}
// evaluate an off-shell spinor
SpinorWaveFunction FFVVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out,
const SpinorWaveFunction & sp,
const VectorWaveFunction &vec,
complex<Energy> mass,
complex<Energy> width) {
// extract the pointers to the particle data objects
tcPDPtr Psp=sp.particle();
tcPDPtr Pvec=vec.particle();
// work out the momentum of the off-shell particle
Lorentz5Momentum pout = sp.momentum()+vec.momentum();
// first calculate the couplings
if(kinematics()) calculateKinematics(sp.momentum(),pout,vec.momentum());
setCoupling(q2,Psp,out,Pvec);
Complex ii(0.,1.); // now evaluate the contribution
// polarization components
Complex e0p3 = vec.t() + vec.z();
Complex e0m3 = vec.t() - vec.z();
Complex e1p2 = vec.x()+ii*vec.y();
Complex e1m2 = vec.x()-ii*vec.y();
// overall factor
Energy2 p2 = pout.m2();
Complex fact=-normPropagator(iopt,p2,out,mass,width);
// momentum components
if(mass.real() < ZERO) mass = iopt==5 ? ZERO : out->mass();
complex<Energy> p1p2 = pout.x()+ii*pout.y();
complex<Energy> p1m2 = pout.x()-ii*pout.y();
// complex nos for for the spinor
Complex s1(0.),s2(0.),s3(0.),s4(0.);
LorentzSpinor<double> spt =sp .wave();
complex<Energy> p0p3=pout.e() + pout.z();
complex<Energy> p0m3=pout.e() - pout.z();
// left piece
if(_left!=0.) {
Complex a3=_left*fact*( spt.s1()*e0p3+spt.s2()*e1m2);
Complex a4=_left*fact*( spt.s1()*e1p2+spt.s2()*e0m3);
s1 +=UnitRemoval::InvE * (p0m3*a3-p1m2*a4);
s2 +=UnitRemoval::InvE * (-p1p2*a3+p0p3*a4);
s3 +=UnitRemoval::InvE * a3*mass;
s4 +=UnitRemoval::InvE * a4*mass;
}
// right piece
if(_right!=0.) {
Complex a1=_right*fact*( spt.s3()*e0m3-spt.s4()*e1m2);
Complex a2=_right*fact*(-spt.s3()*e1p2+spt.s4()*e0p3);
s1 +=UnitRemoval::InvE * a1*mass;
s2 +=UnitRemoval::InvE * a2*mass;
s3 +=UnitRemoval::InvE * (p0p3*a1+p1m2*a2);
s4 +=UnitRemoval::InvE * (p1p2*a1+p0m3*a2);
}
// return the wavefunction
return SpinorWaveFunction(pout,out,s1,s2,s3,s4);
}
// evaluate an off-shell SpinorBar
SpinorBarWaveFunction FFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out,
const SpinorBarWaveFunction & sbar,
const VectorWaveFunction & vec,
complex<Energy> mass,
complex<Energy> width) {
// work out the momentum of the off-shell particle
Lorentz5Momentum pout = sbar.momentum()+vec.momentum();
// first calculate the couplings
if(kinematics()) calculateKinematics(pout,sbar.momentum(),vec.momentum());
setCoupling(q2,out,sbar.particle(),vec.particle());
Complex ii(0.,1.);
// now evaluate the contribution
// polarization components
Complex e0p3=vec.t() + vec.z();
Complex e0m3=vec.t() - vec.z();
Complex e1p2=vec.x()+ii*vec.y();
Complex e1m2=vec.x()-ii*vec.y();
// overall factor
Energy2 p2 = pout.m2();
Complex fact=-normPropagator(iopt,p2,out,mass,width);
// momentum components
if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass();
complex<Energy> p1p2=pout.x()+ii*pout.y();
complex<Energy> p1m2=pout.x()-ii*pout.y();
// complex numbers for the spinor
Complex s1(0.),s2(0.),s3(0.),s4(0.);
LorentzSpinorBar<double> sbart=sbar.wave();
Energy p0p3=pout.e() + pout.z();
Energy p0m3=pout.e() - pout.z();
// left piece
if(_left!=0.) {
Complex a1 = _left*fact*( sbart.s3()*e0p3+sbart.s4()*e1p2);
Complex a2 = _left*fact*( sbart.s3()*e1m2+sbart.s4()*e0m3);
s1 += UnitRemoval::InvE*a1*mass;
s2 += UnitRemoval::InvE*a2*mass;
s3 += UnitRemoval::InvE*(-p0m3*a1+p1p2*a2);
s4 += UnitRemoval::InvE*(+p1m2*a1-p0p3*a2);
}
// right piece
if(_right!=0.) {
Complex a3 = _right*fact*( sbart.s1()*e0m3-sbart.s2()*e1p2);
Complex a4 = _right*fact*(-sbart.s1()*e1m2+sbart.s2()*e0p3);
s1 += UnitRemoval::InvE*(-p0p3*a3-p1p2*a4);
s2 += UnitRemoval::InvE*(-p1m2*a3-p0m3*a4);
s3 += UnitRemoval::InvE*a3*mass;
s4 += UnitRemoval::InvE*a4*mass;
}
return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4);
}
// off-shell vector
VectorWaveFunction FFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out,
const SpinorWaveFunction & sp,
const SpinorBarWaveFunction & sbar,
complex<Energy> mass,
complex<Energy> width) {
// work out the momentum of the off-shell particle
Lorentz5Momentum pout = sbar.momentum()+sp.momentum();
// first calculate the couplings
if(kinematics()) calculateKinematics(sp.momentum(),sbar.momentum(),pout);
setCoupling(q2,sp.particle(),sbar.particle(),out);
Complex ii(0.,1.);
// overall factor
Energy2 p2 = pout.m2();
Complex fact = normPropagator(iopt,p2,out,mass,width);
// momentum components
if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass();
complex<Energy2> mass2 = sqr(mass);
// the vector for the fermion-antifermion
Complex vec[4];
// left coupling
if(_left!=0.) {
vec[0] = -_left*(sbar.s3()*sp.s2()+sbar.s4()*sp.s1());
vec[1] = ii*_left*(sbar.s3()*sp.s2()-sbar.s4()*sp.s1());
vec[2] = -_left*(sbar.s3()*sp.s1()-sbar.s4()*sp.s2());
vec[3] = _left*(sbar.s3()*sp.s1()+sbar.s4()*sp.s2());
}
// right coupling
if(_right!=0.) {
vec[0] += +_right*(sbar.s1()*sp.s4()+sbar.s2()*sp.s3());
vec[1] += -ii*_right*(sbar.s1()*sp.s4()-sbar.s2()*sp.s3());
vec[2] += +_right*(sbar.s1()*sp.s3()-sbar.s2()*sp.s4());
vec[3] += +_right*(sbar.s1()*sp.s3()+sbar.s2()*sp.s4());
}
// massless boson
if(mass.real()==ZERO) {
for(int ix=0;ix<4;++ix){vec[ix]*=fact;}
}
// massive boson
else {
complex<InvEnergy> dot = ( pout.e() *vec[3]
-pout.x()*vec[0]
-pout.y()*vec[1]
-pout.z()*vec[2])/mass2;
vec[0]=fact*(vec[0]-dot*pout.x());
vec[1]=fact*(vec[1]-dot*pout.y());
vec[2]=fact*(vec[2]-dot*pout.z());
vec[3]=fact*(vec[3]-dot*pout.e());
}
return VectorWaveFunction(pout,out,vec[0],vec[1],vec[2],vec[3]);
}
SpinorWaveFunction FFVVertex::evaluateSmall(Energy2 q2,int iopt, tcPDPtr out,
const SpinorWaveFunction & sp,
const VectorWaveFunction & vec,
unsigned int fhel, unsigned int vhel,
double ctheta, double phi, double stheta,
bool includeEikonal,
SmallAngleDirection direction,
Energy mass, Energy) {
assert(fhel <= 1);
assert( vhel == 0 || vhel == 2 );
SpinorWaveFunction output;
// first calculate the couplings
setCoupling(q2,sp.particle(),out,vec.particle());
Complex ii(0.,1.);
if(mass < ZERO) mass = iopt==5 ? ZERO : out->mass();
Lorentz5Momentum pout = sp.momentum()+vec.momentum();
assert(sp.direction()!=intermediate);
// helicity of the boson
double lam = double(vhel)-1.;
// energy of the boson
Energy Eg = abs(vec.momentum().e());
// energy of the fermion
Energy Ef = abs(sp .momentum().e());
// energy fraction of photon
double x = Eg/Ef;
// velocity of the fermon
double beta = sqrt(1.-sqr(mass/Ef));
// dimensionless versions of the variables
double dm = mass*UnitRemoval::InvE;
double dEf = Ef*UnitRemoval::InvE;
double rE = sqrt(.5*dEf);
// calculation of propagator accurate as beta->1 and theta -> 0
Energy2 dot = 2.*Ef*Eg*(sqr(mass/Ef)/(1.+beta)*ctheta
+ sqr(stheta)/(1.+ctheta) );
Complex fact= norm()*(0.5*left()+0.5*right())*UnitRemoval::E2/dot;
// phase factor
Complex ephig = cos(phi )+ii*sin(phi );
// calculation of the spinor
Complex s1(0.),s2(0.),s3(0.),s4(0.);
// u-type spinor
- if(sp.wave().Type()==u_spinortype) {
+ if(sp.wave().Type()==SpinorType::u) {
// fermion along +z
if(direction==PostiveZDirection) {
if(fhel==0) {
s1 = +x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s2 =-rE*dEf*sqrt(1.+beta)*stheta*
(+x*(1.+lam)-(includeEikonal ? 2.*beta*lam : 0. ));
s3 = +x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta);
s4 = +rE*dm/sqrt(1.+beta)*stheta*
(-x*(1.-lam)+(includeEikonal ? 2.*beta*lam : 0. ));
}
else if(fhel==1) {
s1 = +rE*dm/sqrt(1.+beta)*stheta*
(+x*(1.+lam)+(includeEikonal ? 2.*beta*lam : 0. ));
s2 = -x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta);
s3 = -rE*dEf*sqrt(1.+beta)*stheta*
(-x*(1.-lam)-(includeEikonal ? 2.*beta*lam : 0. ));
s4 = -x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta);
}
}
// fermion along -z
else {
if(fhel==0) {
s1 = -rE*dEf*sqrt(1.+beta)*stheta*
(+x*(1.+lam)-(includeEikonal ? 2.*beta*lam : 0. ));
s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s3 = rE*dm/sqrt(1.+beta)*stheta*
(-x*(1.-lam)+(includeEikonal ? 2.*beta*lam : 0. ));
s4 = -x*rE*dm/ephig*(-1.+lam)*(1.+ctheta)/sqrt(1.+beta);
}
else if(fhel==1) {
s1 =-x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s2 =-rE*dm/sqrt(1.+beta)*stheta*
( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s3 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta);
s4 =-rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
}
}
}
// v-type spinor
- else if(sp.wave().Type()==v_spinortype) {
+ else if(sp.wave().Type()==SpinorType::v) {
// fermion along +z
if(direction==PostiveZDirection) {
if(fhel==0) {
s1 = rE*dm/sqrt(1.+beta)*stheta*
(-x*(1.+lam) + ( includeEikonal ? 2.*beta*lam : 0. ));
s2 = x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta);
s3 = rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) - ( includeEikonal ? 2.*beta*lam : 0. ));
s4 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta);
}
else if(fhel==1) {
s1 = x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s2 =-rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0.));
s3 = x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta);
s4 = rE*dm/sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
}
}
// fermion aling -z
else {
if(fhel==0) {
s1 = x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s2 = rE*dm/sqrt(1.+beta)*stheta*
( x*(1.+lam) - ( includeEikonal ? 2.*beta*lam : 0. ));
s3 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta) / (1.+ctheta);
s4 =-rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) - ( includeEikonal ? 2.*beta*lam : 0. ));
}
else if(fhel==1) {
s1 =-rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.+lam) + ( includeEikonal ? 2.*beta*lam : 0. ));
s2 =-x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s3 = rE*dm/sqrt(1.+beta)*stheta*
( x*(1.-lam) + ( includeEikonal ? 2.*beta*lam : 0. ));
s4 =-x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta);
}
}
}
s1 *= -fact;
s2 *= -fact;
s3 *= -fact;
s4 *= -fact;
return SpinorWaveFunction(pout,out,s1,s2,s3,s4);
}
SpinorBarWaveFunction FFVVertex::evaluateSmall(Energy2 q2,int iopt, tcPDPtr out,
const SpinorBarWaveFunction & sbar,
const VectorWaveFunction & vec,
unsigned int fhel, unsigned int vhel,
double ctheta, double phi, double stheta,
bool includeEikonal,
SmallAngleDirection direction,
Energy mass, Energy) {
assert(fhel <= 1);
assert( vhel == 0 || vhel == 2 );
SpinorBarWaveFunction output;
// first calculate the couplings
setCoupling(q2,out,sbar.particle(),vec.particle());
Complex ii(0.,1.);
if(mass < ZERO) mass = iopt==5 ? ZERO : out->mass();
Lorentz5Momentum pout = sbar.momentum()+vec.momentum();
assert(sbar.direction()!=intermediate);
// helicity of the boson
double lam = double(vhel)-1.;
// energies and velocities
Energy Ef = abs(sbar.momentum().e());
Energy Eg = abs(vec .momentum().e());
double x = Eg/Ef;
double beta = sqrt(1.-sqr(mass/Ef));
// calculation of propagator accurate as beta->1 and theta -> 0
Energy2 dot = 2.*Ef*Eg*(sqr(mass/Ef)/(1.+beta)*ctheta
+ sqr(stheta)/(1.+ctheta) );
Complex fact= norm()*(0.5*left()+0.5*right())*UnitRemoval::E2/dot;
// calculation of the spinor
Complex s1(0.),s2(0.),s3(0.),s4(0.);
Complex ephig = cos(phi )+ii*sin(phi );
double dm = mass*UnitRemoval::InvE;
double dEf = Ef*UnitRemoval::InvE;
double rE = sqrt(.5*dEf);
// u-type spinor
- if(sbar.wave().Type()==u_spinortype) {
+ if(sbar.wave().Type()==SpinorType::u) {
// fermion along +z
if(direction==PostiveZDirection) {
if(fhel==0) {
s1 = x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s2 = rE*dm/sqrt(1.+beta)*stheta*
(+x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s3 = -x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta);
s4 = rE*dEf*sqrt(1.+beta)*stheta*
(+x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. ));
}
else if(fhel==1) {
s1 = -rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(lam+1.)*sqr(stheta)/(1.+ctheta);
s3 =-rE*dm/sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s4 = -x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta);
s4 = rE*dm*(1.+ctheta)/ephig*x*(1.-lam)/sqrt(1.+beta);
}
}
// fermion aling -z
else {
if(fhel==0) {
s1 = rE*dm/sqrt(1.+beta)*stheta*
(+x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s2 = -x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s3 = rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s4 = -x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta);
}
else if(fhel==1) {
s1 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s2 = rE*dEf*sqrt(1.+beta)*stheta*
(+x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s3 =-x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta);
s4 = rE*dm/sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
}
}
}
// v-type spinor
- else if(sbar.wave().Type()==v_spinortype) {
+ else if(sbar.wave().Type()==SpinorType::v) {
// anti fermion along +z
if(direction==PostiveZDirection) {
if(fhel==0) {
s1 = -rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s3 = rE*dm/sqrt(1.+beta)*stheta*
( x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s4 =-x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta);
}
else if(fhel==1) {
s1 =-x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s2 =-rE*dm/sqrt(1.+beta)*stheta*
( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s3 =-x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta);
s4 = rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
}
}
// anti fermion aling -z
else {
if(fhel==0) {
s1 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta);
s2 = rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. ));
s3 = x*rE*dm/ephig*(-1.+lam)*(1.+ctheta)/sqrt(1.+beta);
s4 =-rE*dm/sqrt(1.+beta)*stheta*
(+x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. ));
}
else if(fhel==1) {
s1 =-rE*dm/sqrt(1.+beta)*stheta*
( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s2 = x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta);
s3 = rE*dEf*sqrt(1.+beta)*stheta*
( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. ));
s4 =-x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta);
}
}
}
s1 *= -fact;
s2 *= -fact;
s3 *= -fact;
s4 *= -fact;
return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4);
}
diff --git a/Helicity/Vertex/VertexBase.cc b/Helicity/Vertex/VertexBase.cc
--- a/Helicity/Vertex/VertexBase.cc
+++ b/Helicity/Vertex/VertexBase.cc
@@ -1,332 +1,312 @@
// -*- C++ -*-
//
// VertexBase.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VertexBase class.
//
#include "VertexBase.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/Rebinder.h"
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/WidthGenerator.h>
#include <iterator>
using namespace ThePEG;
using namespace ThePEG::Helicity;
VertexBase::VertexBase(VertexType::T name, bool kine)
: _npoint(1), _norm(0), _calckinematics(kine),
- _kine(5,vector<Energy2>(5)), _theName(name),
+ _kine(), _theName(name),
_ordergEM(0), _ordergS(0),
_coupopt(0), _gs(sqrt(4.*Constants::pi*0.3)),
_ee(sqrt(4.*Constants::pi/137.04)),
_sw(sqrt(0.232))
{
assert ( name != VertexType::UNDEFINED );
// Count number of lines from length of 'name'
while ( name /= 10 ) ++_npoint;
}
// setup the lists of particles
// should only be called from child class constructors
void VertexBase::addToList(long ida, long idb, long idc, long idd) {
- if ( idd == 0 ) {
- const long id[] = { ida, idb, idc };
- addToList(vector<long>(id,id+3));
- }
- else {
- const long id[] = { ida, idb, idc, idd };
- addToList(vector<long>(id,id+4));
- }
+ if ( idd == 0 ) addToList({ ida, idb, idc });
+ else addToList({ ida, idb, idc, idd });
}
void VertexBase::addToList(const vector<long> & ids) {
assert( ids.size() == _npoint );
vector<PDPtr> tmp;
int chargeSum = 0;
- for (vector<long>::const_iterator it = ids.begin();
- it != ids.end(); ++it) {
- tPDPtr p = getParticleData(*it);
+ for ( auto id : ids ) {
+ tPDPtr p = getParticleData(id);
if ( !p ) return; // needed e.g. to deal with chi_5 in MSSM
tmp.push_back(p);
chargeSum += p->iCharge();
}
assert( tmp.size() == _npoint );
if ( chargeSum != 0 ) {
cerr << "Problem with the addToList() calls in "
<< fullName() << ":\n"
<< "Vertex particles ";
copy (ids.begin(), ids.end(),
std::ostream_iterator<long>(cerr," "));
cerr << "have non-zero electric charge " << chargeSum << "/3.\n";
assert( false );
}
_particles.push_back(tmp);
}
void VertexBase::doinit() {
Interfaced::doinit();
// set up the incoming and outgoing particles
if ( !_outpart.empty() || !_inpart.empty() )
return;
- for ( unsigned int ix=0; ix<_particles.size(); ++ix ) {
- for ( vector<PDPtr>::const_iterator it = _particles[ix].begin();
- it != _particles[ix].end(); ++it ) {
- tPDPtr p = *it;
+ for ( const auto & pvec : _particles ) {
+ for ( tPDPtr p : pvec ) {
assert( p );
- assert ( p->id() == getParticleData(p->id())->id() );
tPDPtr cc = p->CC();
_inpart.insert( cc ? cc : p );
_outpart.insert(p);
}
}
// check the couplings
if(Debug::level>1&&_npoint!=2+_ordergEM+_ordergS)
generator()->log() << fullName() << " has inconsistent number of "
<< "external particles and coupling order\nQED = "
<< _ordergEM << " QCD = " << _ordergS << " for"
<< " a perturbative interaction. Either it's an"
<< " effective vertex or something is wrong.\n";
assert(_npoint<=2+_ordergEM+_ordergS);
}
void VertexBase::persistentOutput(PersistentOStream & os) const {
os << _npoint << _inpart << _outpart
<< _particles << _calckinematics
<< _coupopt << _gs << _ee << _sw;
}
void VertexBase::persistentInput(PersistentIStream & is, int) {
is >> _npoint >> _inpart >> _outpart
>> _particles >> _calckinematics
>> _coupopt >> _gs >> _ee >> _sw;
}
// Static variable needed for the type description system in ThePEG.
DescribeAbstractClass<VertexBase,Interfaced>
describeThePEGVertexBase("ThePEG::VertexBase", "");
void VertexBase::Init() {
static Switch<VertexBase,bool> interfaceCalculateKinematics
("CalculateKinematics",
"Calculate kinematic invariants at the vertices. This is"
" mainly needed for loop vertices.",
&VertexBase::_calckinematics, false, false, false);
static SwitchOption interfaceCalculateKinematicsCalculate
(interfaceCalculateKinematics,
"Calculate",
"Calculate the kinematics",
true);
static SwitchOption interfaceCalculateKinematicsNoKinematics
(interfaceCalculateKinematics,
"NoKinematics",
"Do not calculate the kinematics",
false);
static ClassDocumentation<VertexBase> documentation
("The VertexBase class is designed to be the base class"
"of all vertices.");
static Switch<VertexBase,unsigned int> interfaceCoupling
("Coupling",
"Treatment of the running couplings",
&VertexBase::_coupopt, 0, false, false);
static SwitchOption interfaceCouplingRunning
(interfaceCoupling,
"Running",
"Use the running couplings from the StandardModel object",
0);
static SwitchOption interfaceCouplingFixedSM
(interfaceCoupling,
"FixedSM",
"Use the fixed values from the StandardModel object",
1);
static SwitchOption interfaceCouplingFixedLocal
(interfaceCoupling,
"FixedLocal",
"Use the local fixed values",
2);
static Parameter<VertexBase,double> interfaceStrongCoupling
("StrongCoupling",
"The fixed value of the strong coupling to use",
&VertexBase::_gs, sqrt(4.*Constants::pi*0.3), 0.0, 10.0,
false, false, Interface::limited);
static Parameter<VertexBase,double> interfaceElectroMagneticCoupling
("ElectroMagneticCoupling",
"The fixed value of the electromagnetic coupling to use",
&VertexBase::_ee, sqrt(4.*Constants::pi/128.91), 0.0, 10.0,
false, false, Interface::limited);
static Parameter<VertexBase,double> interfaceSinThetaW
("SinThetaW",
"The fixed value of sin theta_W to use",
&VertexBase::_sw, sqrt(0.232), 0.0, 10.0,
false, false, Interface::limited);
}
// find particles with a given id
vector<long> VertexBase::search(unsigned int iloc,long idd) const {
assert( iloc < _npoint );
vector<long> out;
- for(unsigned int ix=0; ix<_particles.size(); ++ix) {
- bool found = _particles[ix][iloc]->id() == idd;
+ for(const auto & pvec : _particles ) {
+ bool found = pvec[iloc]->id() == idd;
if(found) {
- for(unsigned int iy=0;iy<_particles[ix].size();++iy) {
- out.push_back(_particles[ix][iy]->id());
+ for( tcPDPtr p : pvec ) {
+ out.push_back(p->id());
}
}
}
return out;
}
// find particles with a given id
vector<tPDPtr> VertexBase::search(unsigned int iloc,tcPDPtr idd) const {
assert( iloc < _npoint );
vector<tPDPtr> out;
- for(unsigned int ix=0; ix<_particles.size(); ++ix) {
- bool found = _particles[ix][iloc] == idd;
+ for(const auto & pvec : _particles) {
+ bool found = pvec[iloc] == idd;
if(found) {
- for(unsigned int iy=0;iy<_particles[ix].size();++iy) {
- out.push_back(_particles[ix][iy]);
- }
+ out = vector<tPDPtr>(pvec.begin(), pvec.end());
}
}
return out;
}
// check a given combination is allowed for a three point vertex
bool VertexBase::allowed(long ida, long idb, long idc, long idd) const {
assert( ( _npoint==3 && idd == 0 ) || _npoint == 4 );
vector<long> out = search(0,ida);
for ( size_t ix = 0; ix < out.size(); ix += _npoint ) {
if ( out[ix+1] == idb && out[ix+2] == idc
&& ( idd == 0 || out[ix+3] == idd ) ) {
return true;
}
}
return false;
}
// output the information
ostream & ThePEG::Helicity::operator<<(ostream & os, const VertexBase & in) {
os << "Information on Vertex" << endl;
os << "This is an " << in._npoint << " vertex\n";
os << string( in._calckinematics ?
"The kinematic invariants are calculated" :
"The kinematics invariants are not calculated" ) << "\n";
os << " Particles allowed for this Vertex\n";
for(unsigned int ix=0;ix<in._particles.size();++ix) {
for(unsigned int iy=0;iy<in._particles[ix].size();++iy) {
os << in._particles[ix][iy]->PDGName() << " ";
}
os << '\n';
}
return os;
}
// calculate the propagator for a diagram
Complex VertexBase::propagator(int iopt, Energy2 p2,tcPDPtr part,
complex<Energy> mass, complex<Energy> width) {
if(mass.real() < ZERO) mass = part->mass();
const complex<Energy2> mass2 = sqr(mass);
if(width.real() < ZERO) {
const tcWidthGeneratorPtr widthgen = part->widthGenerator();
width = widthgen && (iopt==2 || iopt==6 ) ?
widthgen->width(*part,sqrt(p2)) : part->width();
}
const Complex ii(0.,1.);
complex<Energy2> masswidth;
if(iopt==5) {
return Complex(UnitRemoval::E2/p2);
}
else if(iopt==4) {
return 1.0;
}
else if(p2 < ZERO) {
if(iopt!=7)
masswidth = ZERO;
else
masswidth = ii * mass * width;
}
else {
switch (iopt) {
case 1: case 2: case 7:
masswidth = ii * mass * width;
break;
case 3:
masswidth = ZERO;
break;
case 6:
masswidth = ii * mass2 * width / sqrt(p2);
return Complex(UnitRemoval::E2 * (mass2/p2) / (p2-mass2+masswidth));
default:
assert( false );
return -999.999;
}
}
return Complex(UnitRemoval::E2/(p2-mass2+masswidth));
}
void VertexBase::rebind(const TranslationMap & trans) {
- vector<vector<PDPtr> >::iterator cit;
- vector<PDPtr>::iterator cjt;
- for (cit = _particles.begin(); cit != _particles.end(); ++cit) {
- for (cjt = cit->begin(); cjt != cit->end(); ++cjt) {
+ for (auto cit = _particles.begin(); cit != _particles.end(); ++cit) {
+ for (auto cjt = cit->begin(); cjt != cit->end(); ++cjt) {
*cjt = trans.translate(*cjt);
}
}
- set<tPDPtr>::iterator it;
set<tPDPtr> newinpart;
- for (it = _inpart.begin(); it != _inpart.end(); ++it) {
+ for (auto it = _inpart.begin(); it != _inpart.end(); ++it) {
newinpart.insert(trans.translate(*it));
}
_inpart = newinpart;
set<tPDPtr> newoutpart;
- for (it = _outpart.begin(); it != _outpart.end(); ++it) {
+ for (auto it = _outpart.begin(); it != _outpart.end(); ++it) {
newoutpart.insert(trans.translate(*it));
}
_outpart = newoutpart;
Interfaced::rebind(trans);
}
IVector VertexBase::getReferences() {
IVector ret = Interfaced::getReferences();
- vector<vector<PDPtr> >::iterator cit;
- vector<PDPtr>::iterator cjt;
- for (cit = _particles.begin(); cit != _particles.end(); ++cit) {
- for (cjt = cit->begin(); cjt != cit->end(); ++cjt) {
+ for (auto cit = _particles.begin(); cit != _particles.end(); ++cit) {
+ for (auto cjt = cit->begin(); cjt != cit->end(); ++cjt) {
ret.push_back(*cjt);
}
}
- set<tPDPtr>::iterator it;
- for (it = _inpart.begin();
- it != _inpart.end(); ++it) {
+ for (auto it = _inpart.begin(); it != _inpart.end(); ++it) {
ret.push_back(*it);
}
- for (it = _outpart.begin();
- it != _outpart.end(); ++it) {
+ for (auto it = _outpart.begin(); it != _outpart.end(); ++it) {
ret.push_back(*it);
}
return ret;
}
diff --git a/Helicity/Vertex/VertexBase.h b/Helicity/Vertex/VertexBase.h
--- a/Helicity/Vertex/VertexBase.h
+++ b/Helicity/Vertex/VertexBase.h
@@ -1,535 +1,536 @@
// -*- C++ -*-
//
// VertexBase.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_VertexBase_H
#define ThePEG_VertexBase_H
//
// This is the declaration of the VertexBase class.
#include <ThePEG/Interface/Interfaced.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/Helicity/HelicityDefinitions.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "VertexBase.fh"
+#include <array>
+
namespace ThePEG {
namespace Helicity {
/**
* Namespace for naming of vertices. Each child class should extend this
* with its own spin configuration.
*/
namespace VertexType {
typedef unsigned T;
/**
* Undefined Enum for the Lorentz structures
*/
const T UNDEFINED = 0;
}
/** \ingroup Helicity
*
* The VertexBase class is the base class for all helicity amplitude
* vertices. In implements the storage of the particles
* which are allowed to interact at the vertex and some simple functions
* which are often needed by the classes which implement the specific
* vertices.
*
* In practice little use is made of this information and it is mainly
* included for future extensions. It can also be used at the development
* and debugging stage.
*
*/
class VertexBase : public Interfaced {
/**
* The output operator is a friend to avoid the data being public.
*/
friend ostream & operator<<(ostream &, const VertexBase &);
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Constructor for \f$n\f$-point vertices.
* @param name The type of vertex
* @param kine Whether the kinematic invariants should be calculated.
*/
VertexBase(VertexType::T name, bool kine=false);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
public:
/**
* Access to the particle information
*/
//@{
/**
* Number of different particle combinations allowed.
*/
unsigned int size() const { return _particles.size(); }
public:
/**
* Is a particle allowed as an incoming particle?
* @param p The ParticleData pointer
*/
bool isIncoming(tPDPtr p) const {
return _inpart.find(p) != _inpart.end();
}
/**
* Is a particle allowed as an outgoing particle?
* @param p The ParticleData pointer
*/
bool isOutgoing(tPDPtr p) const {
return _outpart.find(p) != _outpart.end();
}
/**
* Get the list of incoming particles.
*/
const set<tPDPtr> & incoming() const { return _inpart; }
/**
* Get the list of outgoing particles.
*/
const set<tPDPtr> & outgoing() const { return _outpart; }
/**
* Get the coupling.
*/
Complex norm() const { return _norm; }
/**
* Function to search the list.
* @param ilist Which list to search
* @param id The PDG code to look for.
*/
vector<long> search(unsigned int ilist,long id) const;
/**
* Function to search the list.
* @param ilist Which list to search
* @param id The particle to look for.
*/
vector<tPDPtr> search(unsigned int ilist,tcPDPtr id) const;
/**
* Is a given combination allowed.
* @param id1 PDG code of the first particle.
* @param id2 PDG code of the second particle.
* @param id3 PDG code of the third particle.
* @param id4 PDG code of the fourth particle.
*/
bool allowed(long id1, long id2, long id3, long id4 = 0) const;
/**
* Get name of Vertex
*/
VertexType::T getName() const { return _theName; }
/**
* Get number of lines on Vertex
*/
unsigned int getNpoint() const { return _npoint; }
/**
* Get the order in \f$g_EM\f$
*/
unsigned int orderInGem() const { return _ordergEM; }
/**
* Get the order in \f$g_s\f$
*/
unsigned int orderInGs() const { return _ordergS; }
//@}
public:
/**
* @name Calculation of the strong, electromagnetic and weak couplings
*/
//@{
/**
* Strong coupling
*/
double strongCoupling(Energy2 q2) const {
if(_coupopt==0) {
double val = 4.0*Constants::pi*generator()->standardModel()->alphaS(q2);
assert(val>=0.);
return sqrt(val);
}
else if(_coupopt==1)
return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaS());
else
return _gs;
}
/**
* Electromagnetic coupling
*/
double electroMagneticCoupling(Energy2 q2) const {
if(_coupopt==0)
return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2));
else if(_coupopt==1)
return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ());
else
return _ee;
}
/**
* Weak coupling
*/
double weakCoupling(Energy2 q2) const {
if( _coupopt == 0 )
return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMME(q2)/
generator()->standardModel()->sin2ThetaW());
else if( _coupopt == 1 )
return sqrt(4.0*Constants::pi*generator()->standardModel()->alphaEMMZ()/
generator()->standardModel()->sin2ThetaW());
else
return _ee/_sw;
}
double sin2ThetaW() const {
if( _coupopt == 0 || _coupopt == 1)
return generator()->standardModel()->sin2ThetaW();
else
return sqr(_sw);
}
//@}
public:
/**
* Set coupling methods
*/
//@{
/**
* Calculate the couplings for a three point interaction.
* This method is virtual and must be implemented in
* classes inheriting from this.
* @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
* @param part1 The ParticleData pointer for the first particle.
* @param part2 The ParticleData pointer for the second particle.
* @param part3 The ParticleData pointer for the third particle.
*/
virtual void setCoupling(Energy2 q2,tcPDPtr part1,
tcPDPtr part2,tcPDPtr part3)=0;
/**
* Calculate the couplings for a four point interaction.
* This method is virtual and must be implemented in
* classes inheriting from this.
* @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
* @param part1 The ParticleData pointer for the first particle.
* @param part2 The ParticleData pointer for the second particle.
* @param part3 The ParticleData pointer for the third particle.
* @param part4 The ParticleData pointer for the fourth particle.
*/
virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3,
tcPDPtr part4)=0;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
//@}
protected:
/**
* Members to set-up the particles
*/
//@{
/**
* Set up the lists of outer particles for the vertex.
* @param ids A vector of PDG codes for the particles.
*/
void addToList(const vector<long> & ids);
/**
* Set up the lists of outer particles for the three-/four-point vertex.
* For small vertices, this form is much easier to use.
* @param ida The PDG codes for the first set of particles.
* @param idb The PDG codes for the second set of particles.
* @param idc The PDG codes for the third set of particles.
* @param idd The PDG codes for the fourth set of particles.
*/
void addToList(long ida, long idb, long idc, long idd = 0);
//@}
protected:
/**
* Members for the amplitude calculations
*/
//@{
/**
* Set the coupling.
* @param coup The coupling.
*/
void norm(const Complex & coup) { _norm = coup; }
/**
* Calculate the propagator for a diagram.
* @param iopt The option for the Breit-Wigner shape
* @param q2 The scale
* @param part The ParticleData pointer for the off-shell particle.
* @param mass The mass if not to be taken from the ParticleData object
* @param width The width if not to be taken from the ParticleData object
*/
virtual Complex propagator(int iopt, Energy2 q2,tcPDPtr part,
complex<Energy> mass=-GeV,
complex<Energy> width=-GeV);
/**
* Calculate propagator multiplied by coupling.
* @param iopt The option for the Breit-Wigner shape
* @param q2 The scale
* @param part The ParticleData pointer for the off-shell particle.
* @param mass The mass if not to be taken from the ParticleData object
* @param width The width if not to be taken from the ParticleData object
*/
Complex normPropagator(int iopt, Energy2 q2,tcPDPtr part,
complex<Energy> mass=-GeV,
complex<Energy> width=-GeV) {
return _norm*propagator(iopt,q2,part,mass,width);
}
//@}
public:
/** @name Kinematic invariants for loop diagrams */
//@{
/**
* Whether or not to calculate the kinematics invariants
*/
bool kinematics() const { return _calckinematics; }
/**
* Set whether or not to calculate the kinematics invariants
*/
void kinematics(bool kine ) { _calckinematics=kine; }
/**
* Calculate the kinematics for a 3-point vertex
*/
void calculateKinematics(const Lorentz5Momentum & p0,
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2) {
_kine[0][0]=p0*p0;
_kine[1][1]=p1*p1;
_kine[2][2]=p2*p2;
_kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
_kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
_kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
}
/**
* Calculate the kinematics for a 4-point vertex
*/
void calculateKinematics(const Lorentz5Momentum & p0,
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2,
const Lorentz5Momentum & p3) {
_kine[0][0]=p0*p0;
_kine[1][1]=p1*p1;
_kine[2][2]=p2*p2;
_kine[3][3]=p3*p3;
_kine[0][1]=p0*p1;_kine[1][0]=_kine[0][1];
_kine[0][2]=p0*p2;_kine[2][0]=_kine[0][2];
_kine[0][3]=p0*p3;_kine[3][0]=_kine[0][3];
_kine[1][2]=p1*p2;_kine[2][1]=_kine[1][2];
_kine[1][3]=p1*p3;_kine[3][1]=_kine[1][3];
_kine[2][3]=p2*p3;_kine[3][2]=_kine[2][3];
}
/**
* Calculate the kinematics for a n-point vertex
*/
void calculateKinematics(const vector<Lorentz5Momentum> & p) {
- unsigned int ix,iy;
- for(ix=0;ix<p.size();++ix) {
- for(iy=0;iy<=ix;++ix) {
+ for(size_t ix=0;ix<p.size();++ix) {
+ for(size_t iy=0;iy<=ix;++ix) {
_kine[ix][iy]=p[ix]*p[iy];
_kine[iy][ix]=_kine[ix][iy];
}
}
}
/**
* Get one of the kinematic invariants
*/
Energy2 invariant(unsigned int ix ,unsigned int iy) const {
assert ( ix < _npoint && iy < _npoint );
return _kine[ix][iy];
}
//@}
protected:
/**
* Set the order in \f$g_EM\f$
* @param order The order of the vertex in \f$g_EM\f$
*/
void orderInGem(unsigned int order) { _ordergEM = order; }
/**
* Set the order in \f$g_s\f$
* @param order The order of the vertex in \f$g_s\f$
*/
void orderInGs(unsigned int order) { _ordergS = order; }
private:
/**
* Describe a concrete class with persistent data.
*/
static AbstractClassDescription<ThePEG::Helicity::VertexBase> initVertexBase;
/**
* Private and non-existent assignment operator.
*/
VertexBase & operator=(const VertexBase &);
private:
/**
* Storage of the particles.
*/
//@{
/**
* Particles interacting at the vertex
*/
vector<vector<PDPtr> > _particles;
/**
* Number of particles at the vertex
*/
unsigned int _npoint;
/**
* ParticleData pointers for the allowed incoming particles.
*/
set<tPDPtr> _inpart;
/**
* ParticleData pointers for the allowed outgoing particles.
*/
set<tPDPtr> _outpart;
//@}
/**
* The overall coupling.
*/
Complex _norm;
/**
* Whether or not to calculate the kinematic invariants for the vertex
*/
bool _calckinematics;
/**
* Kinematica quantities needed for loop vertices
*/
- vector<vector<Energy2> > _kine;
+ std::array<std::array<Energy2,5>,5> _kine;
/**
* Name of vertex
*/
VertexType::T _theName;
/**
* Order of vertex in \f$g_EM\f$
*/
unsigned int _ordergEM;
/**
* Order of vertex in \f$g_s\f$
*/
unsigned int _ordergS;
/**
* option for the coupling
*/
unsigned int _coupopt;
/**
* Fixed value of strong coupling to use
*/
double _gs;
/**
* Fixed value of the electromagentic coupling to use
*/
double _ee;
/**
* Fixed value of \f$\sin\theta_W\f$ to use
*/
double _sw;
};
/**
* Output the information on the vertex.
*/
ostream & operator<<(ostream &, const VertexBase &);
}
}
#endif /* ThePEG_VertexBase_H */
diff --git a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc
--- a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc
+++ b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc
@@ -1,330 +1,330 @@
// -*- C++ -*-
//
// RSSpinorBarWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the RSSpinorBarWaveFunction class.
//
#include "RSSpinorBarWaveFunction.h"
using namespace ThePEG;
using namespace ThePEG::Helicity;
// calculate the Wavefunction
void RSSpinorBarWaveFunction::calculateWaveFunction(unsigned int ihel) {
Complex ii(0.,1.);
LorentzRSSpinorBar<double> news;
- if(direction()==incoming) news=LorentzRSSpinorBar<double>(u_spinortype);
- else news=LorentzRSSpinorBar<double>(v_spinortype);
+ if(direction()==incoming) news=LorentzRSSpinorBar<double>(SpinorType::u);
+ else news=LorentzRSSpinorBar<double>(SpinorType::v);
unsigned int ix,iy;
assert(direction()!=intermediate);
assert(ihel<=3);
// only two valid helicities in massless case
assert( mass()>ZERO || (ihel == 0 || ihel == 3 ) );
// extract the momentum components
// compute the normal spinors to construct the RS spinor
Complex hel_wf[2][2];
if(direction()==incoming) {
// the + spinor
hel_wf[0][0] = 0.;
hel_wf[0][1] = 1.;
// the - spinor
hel_wf[1][0] = 1.;
hel_wf[1][1] = 0.;
}
else {
// the + spinor
hel_wf[0][0] = 1.;
hel_wf[0][1] = 0.;
// the - spinor
hel_wf[1][0] = 0.;
hel_wf[1][1] = 1.;
}
double fact = direction()==incoming ? 1. : -1.;
Energy pmm=mass(),pee=fact*e();
Energy pabs = sqrt(sqr(px())+sqr(py())+sqr(pz()));
SqrtEnergy eplusp = sqrt(pee+pabs);
SqrtEnergy eminusp = ( pmm == ZERO ) ? ZERO : pmm/eplusp;
SqrtEnergy upper[2],lower[2];
if(direction()==incoming) {
upper[0] = eminusp;
lower[0] =-eplusp ;
upper[1] =-eplusp ;
lower[1] = eminusp;
}
else {
upper[0] = eplusp ;
lower[0] = eminusp;
upper[1] = eminusp;
lower[1] = eplusp ;
}
// now construct the spinors
complex<SqrtEnergy> spinor[2][4];
for(ix=0;ix<2;++ix) {
spinor[ix][0] = upper[ix]*hel_wf[ix][0];
spinor[ix][1] = upper[ix]*hel_wf[ix][1];
spinor[ix][2] = lower[ix]*hel_wf[ix][0];
spinor[ix][3] = lower[ix]*hel_wf[ix][1];
}
// compute the polarization vectors to construct the RS spinor
Complex vec[3][4];
double ort = sqrt(0.5);
double r1 = ( pmm == ZERO ) ? 0. : double(pee /pmm);
double r2 = ( pmm == ZERO ) ? 0. : double(pabs/pmm);
if(direction()==incoming) {
vec[0][0] = ort;
vec[0][1] = ort*ii;
vec[0][2] = 0.;
vec[0][3] = 0.;
vec[1][0] = 0.;
vec[1][1] = 0.;
vec[1][2] =-r1;
vec[1][3] =-r2;
vec[2][0] =-ort;
vec[2][1] = ort*ii;
vec[2][2] = 0.;
vec[2][3] = 0.;
}
else {
vec[0][0] =-ort;
vec[0][1] = ort*ii;
vec[0][2] = 0.;
vec[0][3] = 0.;
vec[1][0] = 0.;
vec[1][1] = 0.;
vec[1][2] = r1;
vec[1][3] = r2;
vec[2][0] = ort;
vec[2][1] = ort*ii;
vec[2][2] = 0.;
vec[2][3] = 0.;
}
// now we can put the bits together to compute the RS spinor
double or3(sqrt(1./3.)),tor3(sqrt(2./3.));
if(ihel==3) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy];
}
else if(ihel==2) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*
(or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy]);
}
else if(ihel==1) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*
(or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy]);
}
else if(ihel==0) {
for(ix=0;ix<4;++ix) {
for(iy=0;iy<4;++iy) {
news(ix,iy)=UnitRemoval::InvSqrtE*(vec[2][ix]*spinor[1][iy]);
}
}
}
// spinor is currently along the z axis, rotate so in right direction
if(pabs/pmm>1e-8) {
Axis axis;
axis.setX(fact*momentum().x()/pabs);
axis.setY(fact*momentum().y()/pabs);
axis.setZ(fact*momentum().z()/pabs);
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
if(sinth>1e-8) {
LorentzRotation rot;
rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
_wf = news.transform(rot);
}
else if (axis.z()<0.) {
LorentzRotation rot;
rot.setRotateX(Constants::pi);
_wf = news.transform(rot);
}
else {
_wf = news;
}
}
else {
_wf=news;
}
}
void RSSpinorBarWaveFunction::
calculateWaveFunctions(vector<LorentzRSSpinorBar<SqrtEnergy> > & waves,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getProductionBasisState(ix).bar();
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getDecayBasisState(ix).bar();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWf();
}
}
}
void RSSpinorBarWaveFunction::
calculateWaveFunctions(vector<RSSpinorBarWaveFunction> & waves,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorBarWaveFunction(particle,
inspin->getProductionBasisState(ix).bar(),
dir);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).bar(),
dir);
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
}
}
void RSSpinorBarWaveFunction::
calculateWaveFunctions(vector<LorentzRSSpinorBar<SqrtEnergy> > & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getProductionBasisState(ix).bar();
rho = RhoDMatrix(PDT::Spin3Half);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getDecayBasisState(ix).bar();
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWf();
}
rho = RhoDMatrix(PDT::Spin3Half);
}
}
void RSSpinorBarWaveFunction::
calculateWaveFunctions(vector<RSSpinorBarWaveFunction> & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorBarWaveFunction(particle,
inspin->getProductionBasisState(ix).bar(),
dir);
rho = RhoDMatrix(PDT::Spin3Half);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).bar(),
dir);
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
rho = RhoDMatrix(PDT::Spin3Half);
}
}
void RSSpinorBarWaveFunction::
constructSpinInfo(const vector<LorentzRSSpinorBar<SqrtEnergy> > & waves,
tPPtr particle,Direction dir, bool time) {
assert(waves.size()==4);
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) inspin->setBasisState(ix,waves[ix].bar());
else inspin->setDecayState(ix,waves[ix].bar());
}
else {
RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time));
particle->spinInfo(temp);
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) temp->setBasisState(ix,waves[ix].bar());
else temp->setDecayState(ix,waves[ix].bar());
}
}
void RSSpinorBarWaveFunction::
constructSpinInfo(const vector<RSSpinorBarWaveFunction> & waves,
tPPtr part,Direction dir, bool time) {
assert(waves.size()==4);
tRSFermionSpinPtr inspin = !part->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(part->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<4;++ix)
if (dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf().bar());
else inspin->setDecayState(ix,waves[ix].dimensionedWf().bar());
}
else {
RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(part->momentum(),time));
part->spinInfo(temp);
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) temp->setBasisState(ix,waves[ix].dimensionedWf().bar());
else temp->setDecayState(ix,waves[ix].dimensionedWf().bar());
}
}
diff --git a/Helicity/WaveFunction/RSSpinorWaveFunction.cc b/Helicity/WaveFunction/RSSpinorWaveFunction.cc
--- a/Helicity/WaveFunction/RSSpinorWaveFunction.cc
+++ b/Helicity/WaveFunction/RSSpinorWaveFunction.cc
@@ -1,328 +1,328 @@
// -*- C++ -*-
//
// RSSpinorWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the RSSpinorWaveFunction class.
//
#include "RSSpinorWaveFunction.h"
using namespace ThePEG;
using namespace ThePEG::Helicity;
// calculate the Wavefunction
void RSSpinorWaveFunction::calculateWaveFunction(unsigned int ihel) {
LorentzRSSpinor<double> news;
if(direction()==incoming)
- news=LorentzRSSpinor<double>(u_spinortype);
+ news=LorentzRSSpinor<double>(SpinorType::u);
else
- news=LorentzRSSpinor<double>(v_spinortype);
+ news=LorentzRSSpinor<double>(SpinorType::v);
unsigned int ix,iy;
// check helicity and type
assert(direction()!=intermediate);
assert(ihel<=3);
// massive
// only two valid helicities in massless case
assert( mass()>ZERO || ( ihel == 0 || ihel==3 ) );
// extract the momentum components
// compute the normal spinors to construct the RS spinor
Complex hel_wf[2][2];
if(direction()==incoming) {
// the + spinor
hel_wf[0][0] = 1.;
hel_wf[0][1] = 0.;
// the - spinor
hel_wf[1][0] = 0.;
hel_wf[1][1] = 1.;
}
else {
// the + spinor
hel_wf[0][0] = 0.;
hel_wf[0][1] = 1.;
// the - spinor
hel_wf[1][0] = 1.;
hel_wf[1][1] = 0.;
}
// prefactors
double fact = direction()==incoming ? 1. : -1.;
Energy pmm=mass(),pee=fact*e();
Energy pabs = sqrt(sqr(px())+sqr(py())+sqr(pz()));
SqrtEnergy eplusp = sqrt(pee+pabs);
SqrtEnergy eminusp = ( pmm == ZERO ) ? ZERO : pmm/eplusp;
SqrtEnergy upper[2],lower[2];
if(direction()==incoming) {
upper[0] = eminusp;
lower[0] = eplusp ;
upper[1] = eplusp ;
lower[1] = eminusp;
}
else {
upper[0] =-eplusp ;
lower[0] = eminusp;
upper[1] = eminusp;
lower[1] =-eplusp ;
}
// now construct the spinors
complex<SqrtEnergy> spinor[2][4];
for(ix=0;ix<2;++ix) {
spinor[ix][0] = upper[ix]*hel_wf[ix][0];
spinor[ix][1] = upper[ix]*hel_wf[ix][1];
spinor[ix][2] = lower[ix]*hel_wf[ix][0];
spinor[ix][3] = lower[ix]*hel_wf[ix][1];
}
// compute the polarization vectors to construct the RS spinor
Complex vec[3][4],ii(0.,1.);
double ort = sqrt(0.5);
double r1 = ( pmm == ZERO ) ? 0. : double(pee /pmm);
double r2 = ( pmm == ZERO ) ? 0. : double(pabs/pmm);
if(direction()==incoming) {
vec[0][0] =-ort;
vec[0][1] =-ort*ii;
vec[0][2] = 0.;
vec[0][3] = 0.;
vec[1][0] = 0.;
vec[1][1] = 0.;
vec[1][2] = r1;
vec[1][3] = r2;
vec[2][0] = ort;
vec[2][1] =-ort*ii;
vec[2][2] = 0.;
vec[2][3] = 0.;
}
else {
vec[0][0] = ort;
vec[0][1] =-ort*ii;
vec[0][2] = 0.;
vec[0][3] = 0.;
vec[1][0] = 0.;
vec[1][1] = 0.;
vec[1][2] =-r1;
vec[1][3] =-r2;
vec[2][0] =-ort;
vec[2][1] =-ort*ii;
vec[2][2] = 0.;
vec[2][3] = 0.;
}
// now we can put the bits together to compute the RS spinor
double or3(sqrt(1./3.)),tor3(sqrt(2./3.));
if(ihel==3) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy];
}
else if(ihel==2) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*
(or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy]);
}
else if(ihel==1) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*
(or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy]);
}
else if(ihel==0) {
for(ix=0;ix<4;++ix)
for(iy=0;iy<4;++iy)
news(ix,iy)=UnitRemoval::InvSqrtE*vec[2][ix]*spinor[1][iy];
}
// spinor is currently along the z axis, rotate so in right direction
if(pabs/pmm>1e-8) {
Axis axis;
axis.setX(fact*momentum().x()/pabs);
axis.setY(fact*momentum().y()/pabs);
axis.setZ(fact*momentum().z()/pabs);
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
if(sinth>1e-8) {
LorentzRotation rot;
rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
_wf= news.transform(rot);
}
else if (axis.z()<0.) {
LorentzRotation rot;
rot.setRotateX(Constants::pi);
_wf= news.transform(rot);
}
else {
_wf = news;
}
}
else {
_wf=news;
}
}
void RSSpinorWaveFunction::
calculateWaveFunctions(vector<LorentzRSSpinor<SqrtEnergy> > & waves,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getProductionBasisState(ix);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getDecayBasisState(ix);
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWf();
}
}
}
void RSSpinorWaveFunction::
calculateWaveFunctions(vector<RSSpinorWaveFunction> & waves,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorWaveFunction(particle,
inspin->getProductionBasisState(ix),dir);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorWaveFunction(particle,
inspin->getDecayBasisState(ix),dir);
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
}
}
void RSSpinorWaveFunction::
calculateWaveFunctions(vector<LorentzRSSpinor<SqrtEnergy> > & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getProductionBasisState(ix);
rho = RhoDMatrix(PDT::Spin3Half);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = inspin->getDecayBasisState(ix);
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWf();
}
rho = RhoDMatrix(PDT::Spin3Half);
}
}
void RSSpinorWaveFunction::
calculateWaveFunctions(vector<RSSpinorWaveFunction> & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
waves.resize(4);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorWaveFunction(particle,
inspin->getProductionBasisState(ix),dir);
rho = RhoDMatrix(PDT::Spin3Half);
}
else {
inspin->decay();
for(unsigned int ix=0;ix<4;++ix)
waves[ix] = RSSpinorWaveFunction(particle,
inspin->getDecayBasisState(ix),dir);
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<4;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
rho = RhoDMatrix(PDT::Spin3Half);
}
}
void RSSpinorWaveFunction::
constructSpinInfo(const vector<LorentzRSSpinor<SqrtEnergy> > & waves,
tPPtr particle,Direction dir,bool time) {
assert(waves.size()==4);
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) inspin->setBasisState(ix,waves[ix]);
else inspin->setDecayState(ix,waves[ix]);
}
else {
RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time));
particle->spinInfo(temp);
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) temp->setBasisState(ix,waves[ix]);
else temp->setDecayState(ix,waves[ix]);
}
}
void RSSpinorWaveFunction::
constructSpinInfo(const vector<RSSpinorWaveFunction> & waves,
tPPtr particle,Direction dir,bool time) {
assert(waves.size()==4);
tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() :
dynamic_ptr_cast<tRSFermionSpinPtr>(particle->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf());
else inspin->setDecayState(ix,waves[ix].dimensionedWf());
}
else {
RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time));
particle->spinInfo(temp);
for(unsigned int ix=0;ix<4;++ix)
if(dir==outgoing) temp->setBasisState(ix,waves[ix].dimensionedWf());
else temp->setDecayState(ix,waves[ix].dimensionedWf());
}
}
diff --git a/Helicity/WaveFunction/SpinorBarWaveFunction.cc b/Helicity/WaveFunction/SpinorBarWaveFunction.cc
--- a/Helicity/WaveFunction/SpinorBarWaveFunction.cc
+++ b/Helicity/WaveFunction/SpinorBarWaveFunction.cc
@@ -1,340 +1,340 @@
// -*- C++ -*-
//
// SpinorBarWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SpinorBarWaveFunction class.
//
// Author: Peter Richardson
//
#include "SpinorBarWaveFunction.h"
#include "SpinorWaveFunction.h"
using namespace ThePEG;
using namespace Helicity;
// calculate the Wavefunction
void SpinorBarWaveFunction::calculateWaveFunction(unsigned int ihel) {
Direction dir=direction();
if(dir==intermediate) ThePEG::Helicity::HelicityConsistencyError()
<< "In SpinorBarWaveFunction::calcluateWaveFunction "
<< "particle must be incoming or outgoing not intermediate"
<< Exception::abortnow;
// check ihelicity is O.K.
if(ihel>1) ThePEG::Helicity::HelicityConsistencyError()
<< "Invalid Helicity = " << ihel << " requested for SpinorBar"
<< Exception::abortnow;
// extract the momentum components
double fact = dir==incoming ? 1. : -1.;
Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass();
// define and calculate some kinematic quantities
Energy2 ptran2 = ppx*ppx+ppy*ppy;
Energy pabs = sqrt(ptran2+ppz*ppz);
Energy ptran = sqrt(ptran2);
// first need to evalulate the 2-component helicity spinors
// this is the same regardless of which definition of the spinors
// we are using
Complex hel_wf[2];
// compute the + spinor for + helicty particles and - helicity antiparticles
if((dir==outgoing && ihel== 1) || (dir==incoming && ihel==0)) {
// no transverse momentum
if(ptran==ZERO) {
if(ppz>=ZERO) {
hel_wf[0] = 1;
hel_wf[1] = 0;
}
else {
hel_wf[0] = 0;
hel_wf[1] = 1;
}
}
else {
InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
hel_wf[0] = denominator*rtppluspz;
hel_wf[1] = denominator/rtppluspz*complex<Energy>(ppx,-ppy);
}
}
// compute the - spinor for - helicty particles and + helicity antiparticles
else {
// no transverse momentum
if(ptran==ZERO) {
if(ppz>=ZERO) {
hel_wf[0] = 0;
hel_wf[1] = 1;
}
// transverse momentum
else {
hel_wf[0] = -1;
hel_wf[1] = 0;
}
}
else {
InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
hel_wf[0] = denominator/rtppluspz*complex<Energy>(-ppx,-ppy);
hel_wf[1] = denominator*rtppluspz;
}
}
SqrtEnergy upper, lower;
SqrtEnergy eplusp = sqrt(max(pee+pabs,ZERO));
SqrtEnergy eminusp = ( pmm!=ZERO ) ? pmm/eplusp : ZERO;
// set up the coefficients for the different cases
if(dir==outgoing) {
if(ihel==1) {
upper = eplusp;
lower = eminusp;
}
else {
upper = eminusp;
lower = eplusp;
}
}
else {
if(ihel==1) {
upper = eminusp;
lower = -eplusp;
}
else {
upper =-eplusp;
lower = eminusp;
}
}
// now finally we can construct the spinors
- _wf = LorentzSpinorBar<double>((dir==incoming) ? v_spinortype : u_spinortype);
+ _wf = LorentzSpinorBar<double>((dir==incoming) ? SpinorType::v : SpinorType::u);
_wf[0] = upper*hel_wf[0]*UnitRemoval::InvSqrtE;
_wf[1] = upper*hel_wf[1]*UnitRemoval::InvSqrtE;
_wf[2] = lower*hel_wf[0]*UnitRemoval::InvSqrtE;
_wf[3] = lower*hel_wf[1]*UnitRemoval::InvSqrtE;
}
void SpinorBarWaveFunction::conjugate() {
_wf=_wf.conjugate();
}
SpinorWaveFunction SpinorBarWaveFunction::bar() {
Lorentz5Momentum p = momentum();
if(direction()==outgoing) p *= -1.;
tcPDPtr ptemp = particle();
if(direction()==incoming&&particle()->CC())
ptemp = particle()->CC();
return SpinorWaveFunction(p,ptemp,_wf.bar(),direction());
}
void SpinorBarWaveFunction::
calculateWaveFunctions(vector<LorentzSpinorBar<SqrtEnergy> > & waves,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
waves[0] = inspin->getProductionBasisState(0).bar();
waves[1] = inspin->getProductionBasisState(1).bar();
}
else {
inspin->decay();
- if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
waves[0] = inspin->getDecayBasisState(0).conjugate().bar();
waves[1] = inspin->getDecayBasisState(1).conjugate().bar();
}
else {
waves[0] = inspin->getDecayBasisState(0).bar();
waves[1] = inspin->getDecayBasisState(1).bar();
}
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWave();
}
}
}
void SpinorBarWaveFunction::
calculateWaveFunctions(vector<SpinorBarWaveFunction> & waves,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getProductionBasisState(ix).bar(),
dir);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).conjugate().bar(),dir);
}
else {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).bar(),dir);
}
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
}
}
void SpinorBarWaveFunction::
calculateWaveFunctions(vector<LorentzSpinorBar<SqrtEnergy> > & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
waves[0] = inspin->getProductionBasisState(0).bar();
waves[1] = inspin->getProductionBasisState(1).bar();
rho = RhoDMatrix(PDT::Spin1Half);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
waves[0] = inspin->getDecayBasisState(0).conjugate().bar();
waves[1] = inspin->getDecayBasisState(1).conjugate().bar();
}
else {
waves[0] = inspin->getDecayBasisState(0).bar();
waves[1] = inspin->getDecayBasisState(1).bar();
}
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWave();
}
rho = RhoDMatrix(PDT::Spin1Half);
}
}
void SpinorBarWaveFunction::
calculateWaveFunctions(vector<SpinorBarWaveFunction> & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getProductionBasisState(ix).bar(),
dir);
rho = RhoDMatrix(PDT::Spin1Half);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).conjugate().bar(),dir);
}
else {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorBarWaveFunction(particle,
inspin->getDecayBasisState(ix).bar(),dir);
}
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
rho = RhoDMatrix(PDT::Spin1Half);
}
}
void SpinorBarWaveFunction::
constructSpinInfo(const vector<LorentzSpinorBar<SqrtEnergy> > & waves,
tPPtr part,Direction dir, bool time) {
assert(waves.size()==2);
tFermionSpinPtr inspin = !part->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(part->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
inspin->setBasisState(ix,waves[ix].bar());
else
inspin->setDecayState(ix,waves[ix].bar());
}
}
else {
FermionSpinPtr temp = new_ptr(FermionSpinInfo(part->momentum(),time));
part->spinInfo(temp);
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
temp->setBasisState(ix,waves[ix].bar());
else
temp->setDecayState(ix,waves[ix].bar());
}
}
}
void SpinorBarWaveFunction::
constructSpinInfo(const vector<SpinorBarWaveFunction> & waves,
tPPtr part,Direction dir, bool time) {
assert(waves.size()==2);
tFermionSpinPtr inspin = !part->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(part->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<2;++ix)
if (dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf().bar());
else inspin->setDecayState(ix,waves[ix].dimensionedWf().bar());
}
else {
FermionSpinPtr temp = new_ptr(FermionSpinInfo(part->momentum(),time));
part->spinInfo(temp);
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
temp->setBasisState(ix,waves[ix].dimensionedWf().bar());
else
temp->setDecayState(ix,waves[ix].dimensionedWf().bar());
}
}
}
diff --git a/Helicity/WaveFunction/SpinorWaveFunction.cc b/Helicity/WaveFunction/SpinorWaveFunction.cc
--- a/Helicity/WaveFunction/SpinorWaveFunction.cc
+++ b/Helicity/WaveFunction/SpinorWaveFunction.cc
@@ -1,340 +1,340 @@
// -*- C++ -*-
//
// SpinorWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2011 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SpinorWaveFunction class.
//
// Author: Peter Richardson
//
#include "SpinorWaveFunction.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "SpinorBarWaveFunction.h"
using namespace ThePEG;
using namespace Helicity;
// calculate the Wavefunction
void SpinorWaveFunction::calculateWaveFunction(unsigned int ihel) {
// check helicity is O.K.
Direction dir = direction();
if(dir==intermediate) throw ThePEG::Helicity::HelicityConsistencyError()
<< "In SpinorWaveFunction::calcluateWaveFunction "
<< "particle must be incoming or outgoing not intermediate"
<< Exception::abortnow;
if(ihel>1) throw ThePEG::Helicity::HelicityConsistencyError()
<< "Invalid Helicity = " << ihel << " requested for Spinor"
<< Exception::abortnow;
// extract the momentum components
double fact=-1.; if(dir==incoming){fact=1.;}
Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass();
// define and calculate some kinematic quantities
Energy2 ptran2 = ppx*ppx+ppy*ppy;
Energy pabs = sqrt(ptran2+ppz*ppz);
Energy ptran = sqrt(ptran2);
// first need to evalulate the 2-component helicity spinors
// this is the same regardless of which definition of the spinors
// we are using
complex <double> hel_wf[2];
// compute the + spinor for + helicty particles and - helicity antiparticles
if((dir==incoming && ihel== 1) || (dir==outgoing && ihel==0)) {
// no transverse momentum
if(ptran==ZERO) {
if(ppz>=ZERO) {
hel_wf[0] = 1;
hel_wf[1] = 0;
}
else {
hel_wf[0] = 0;
hel_wf[1] = 1;
}
}
else {
InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
hel_wf[0] = denominator*rtppluspz;
hel_wf[1] = denominator/rtppluspz*complex<Energy>(ppx,ppy);
}
}
// compute the - spinor for - helicty particles and + helicity antiparticles
else {
// no transverse momentum
if(ptran==ZERO) {
if(ppz>=ZERO) {
hel_wf[0] = 0;
hel_wf[1] = 1;
}
// transverse momentum
else {
hel_wf[0] = -1;
hel_wf[1] = 0;
}
}
else {
InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
hel_wf[0] = denominator/rtppluspz*complex<Energy>(-ppx,ppy);
hel_wf[1] = denominator*rtppluspz;
}
}
SqrtEnergy upper,lower;
SqrtEnergy eplusp = sqrt(max(pee+pabs,ZERO));
SqrtEnergy eminusp = ( pmm != ZERO ) ? pmm/eplusp : ZERO;
// set up the coefficients for the different cases
if(dir==incoming) {
if(ihel==1) {
upper = eminusp;
lower = eplusp;
}
else {
upper = eplusp;
lower = eminusp;
}
}
else {
if(ihel==1) {
upper = -eplusp;
lower = eminusp;
}
else {
upper = eminusp;
lower =-eplusp;
}
}
// now finally we can construct the spinors
- _wf = LorentzSpinor<double>( (dir==incoming) ? u_spinortype : v_spinortype);
+ _wf = LorentzSpinor<double>( (dir==incoming) ? SpinorType::u : SpinorType::v);
_wf[0]=upper*hel_wf[0]*UnitRemoval::InvSqrtE;
_wf[1]=upper*hel_wf[1]*UnitRemoval::InvSqrtE;
_wf[2]=lower*hel_wf[0]*UnitRemoval::InvSqrtE;
_wf[3]=lower*hel_wf[1]*UnitRemoval::InvSqrtE;
}
SpinorBarWaveFunction SpinorWaveFunction::bar() {
Lorentz5Momentum p = momentum();
if(direction()==outgoing) p *= -1.;
tcPDPtr ptemp = particle();
if(direction()==incoming&&particle()->CC())
ptemp = particle()->CC();
return SpinorBarWaveFunction(p,ptemp,_wf.bar(),direction());
}
void SpinorWaveFunction::
calculateWaveFunctions(vector<LorentzSpinor<SqrtEnergy> > & waves,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
waves[0] = inspin->getProductionBasisState(0);
waves[1] = inspin->getProductionBasisState(1);
}
else {
inspin->decay();
- if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
waves[0] = inspin->getDecayBasisState(0).conjugate();
waves[1] = inspin->getDecayBasisState(1).conjugate();
}
else {
waves[0] = inspin->getDecayBasisState(0);
waves[1] = inspin->getDecayBasisState(1);
}
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWave();
}
}
}
void SpinorWaveFunction::
calculateWaveFunctions(vector<SpinorWaveFunction> & waves,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getProductionBasisState(ix),dir);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getDecayBasisState(ix).conjugate(),dir);
}
else {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getDecayBasisState(ix),dir);
}
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
}
}
void SpinorWaveFunction::
calculateWaveFunctions(vector<LorentzSpinor<SqrtEnergy> > & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
waves[0] = inspin->getProductionBasisState(0);
waves[1] = inspin->getProductionBasisState(1);
rho = RhoDMatrix(PDT::Spin1Half);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
waves[0] = inspin->getDecayBasisState(0).conjugate();
waves[1] = inspin->getDecayBasisState(1).conjugate();
}
else {
waves[0] = inspin->getDecayBasisState(0);
waves[1] = inspin->getDecayBasisState(1);
}
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave.dimensionedWave();
}
rho = RhoDMatrix(PDT::Spin1Half);
}
}
void SpinorWaveFunction::
calculateWaveFunctions(vector<SpinorWaveFunction> & waves,
RhoDMatrix & rho,
tPPtr particle,Direction dir) {
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
waves.resize(2);
// spin info object exists
if(inspin) {
if(dir==outgoing) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getProductionBasisState(ix),dir);
rho = RhoDMatrix(PDT::Spin1Half);
}
else {
inspin->decay();
- if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=u_spinortype) ||
- (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=v_spinortype)) {
+ if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) ||
+ (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getDecayBasisState(ix).conjugate(),dir);
}
else {
for(unsigned int ix=0;ix<2;++ix)
waves[ix] = SpinorWaveFunction(particle,
inspin->getDecayBasisState(ix),dir);
}
rho = inspin->rhoMatrix();
}
}
// do the calculation
else {
assert(!particle->spinInfo());
SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir);
for(unsigned int ix=0;ix<2;++ix) {
wave.reset(ix);
waves[ix] = wave;
}
rho = RhoDMatrix(PDT::Spin1Half);
}
}
void SpinorWaveFunction::
constructSpinInfo(const vector<LorentzSpinor<SqrtEnergy> > & waves,
tPPtr particle,Direction dir,bool time) {
assert(waves.size()==2);
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
inspin->setBasisState(ix,waves[ix]);
else
inspin->setDecayState(ix,waves[ix]);
}
}
else {
FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time));
particle->spinInfo(temp);
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
temp->setBasisState(ix,waves[ix]);
else
temp->setDecayState(ix,waves[ix]);
}
}
}
void SpinorWaveFunction::
constructSpinInfo(const vector<SpinorWaveFunction> & waves,
tPPtr particle,Direction dir,bool time) {
assert(waves.size()==2);
tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() :
dynamic_ptr_cast<tFermionSpinPtr>(particle->spinInfo());
if(inspin) {
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
inspin->setBasisState(ix,waves[ix].dimensionedWf());
else
inspin->setDecayState(ix,waves[ix].dimensionedWf());
}
}
else {
FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time));
particle->spinInfo(temp);
for(unsigned int ix=0;ix<2;++ix) {
if(( dir == outgoing && time) ||
( dir == incoming && !time))
temp->setBasisState(ix,waves[ix].dimensionedWf());
else
temp->setDecayState(ix,waves[ix].dimensionedWf());
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:22 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3796626
Default Alt Text
(180 KB)

Event Timeline