/usr/include/ThePEG/Utilities/SimplePhaseSpace.tcc is in libthepeg-dev 1.8.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 | // -*- C++ -*-
//
// SimplePhaseSpace.tcc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 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 templated member
// functions of the SimplePhaseSpace class.
//
namespace ThePEG {
template <typename PType>
void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s)
{
typedef ParticleTraits<PType> Traits;
Energy z = getMagnitude(s, Traits::mass(p1), Traits::mass(p2));
Traits::set3Momentum(p1, Momentum3(ZERO, ZERO, z));
Traits::set3Momentum(p2, Momentum3(ZERO, ZERO, -z));
}
template <typename PType>
void SimplePhaseSpace::CMS(Energy2 s, PType & p1, PType & p2)
{
CMS(p1, p2, s, 2.0*UseRandom::rnd() - 1.0, Constants::twopi*UseRandom::rnd());
}
template <typename PType>
void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s,
double cthe, double phi)
{
typedef ParticleTraits<PType> Traits;
Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2));
double sthe = sqrt(1.0-sqr(cthe));
Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe);
Traits::set3Momentum(p1, p);
Traits::set3Momentum(p2, -p);
}
template <typename PType>
void SimplePhaseSpace::
CMS(PType & p1, PType & p2, PType & p3,
Energy2 s, double x1, double x3)
{
CMS(p1, p2, p3, s, x1, x3,
Constants::twopi*UseRandom::rnd(), acos(2.0*UseRandom::rnd() - 1.0),
Constants::twopi*UseRandom::rnd());
}
template <typename PType>
void SimplePhaseSpace::
CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
double x1, double x3, double phii, double the, double phi)
{
typedef ParticleTraits<PType> Traits;
Energy Etot = sqrt(s);
Energy m1 = Traits::mass(p1);
Energy m2 = Traits::mass(p2);
Energy m3 = Traits::mass(p3);
Energy e1 = 0.5*x1*Etot;
Energy e3 = 0.5*x3*Etot;
Energy e2 = Etot - e1 - e3;
if ( e1 < m1 || e2 < m2 || e3 < m3 ) throw ImpossibleKinematics();
Energy r1 = sqrt(sqr(e1)-sqr(m1));
Energy r2 = sqrt(sqr(e2)-sqr(m2));
Energy r3 = sqrt(sqr(e3)-sqr(m3));
Traits::set3Momentum(p1, Momentum3(ZERO, ZERO, r1));
double cthe2 = (sqr(r3)-sqr(r2)-sqr(r1))/(2.0*r2*r1);
double cthe3 = (sqr(r2)-sqr(r3)-sqr(r1))/(2.0*r3*r1);
if ( abs(cthe2) > 1.0 || abs(cthe3) > 1.0 ) throw ImpossibleKinematics();
double sthe2 = sqrt(1.0-sqr(cthe2));
Energy px = r2*sthe2*cos(phii);
Energy py = r2*sthe2*sin(phii);
Traits::set3Momentum(p2, Momentum3(px, py, r2*cthe2));
Traits::set3Momentum(p3, Momentum3(-px, -py, r3*cthe3));
if ( the == 0.0 && phi == 0.0 ) return;
LorentzRotation r;
r.rotateZ(phi);
r.rotateX(the);
Traits::transform(p1, r);
Traits::transform(p2, r);
Traits::transform(p3, r);
}
template <typename PType>
void SimplePhaseSpace::
CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi,
const PType & p0) {
typedef ParticleTraits<PType> Traits;
Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2));
Energy e = sqrt(sqr(r) + sqr(Traits::mass(p1)));
Energy r0 = Traits::momentum(p0).rho();
Energy e0 = Traits::momentum(p0).e();
double cthe = (t + sqr(e - e0) + sqr(r) + sqr(r0))/(2.0*r*r0);
if ( abs(cthe) > 1.0 ) throw ImpossibleKinematics();
double sthe = sqrt(1.0-sqr(cthe));
Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe);
Traits::set3Momentum(p1, p);
Traits::set3Momentum(p2, -p);
if ( Traits::momentum(p0).perp2() > ZERO ) {
LorentzRotation r;
r.rotateX(Traits::momentum(p0).theta());
r.rotateZ(Traits::momentum(p0).phi());
Traits::transform(p1, r);
Traits::transform(p2, r);
}
}
template <typename Container>
void SimplePhaseSpace::CMSn(Container & particles, Energy m0)
{
typedef typename Container::value_type PType;
typedef typename Container::iterator Iterator;
if ( particles.size() == 2 ) {
Iterator it = particles.begin();
PType & p1 = *it++;
PType & p2 = *it;
CMS(sqr(m0), p1, p2);
return;
}
typedef ParticleTraits<PType> Traits;
vector<Energy> masses(particles.size());
int j = 0;
for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j )
masses[j] = Traits::mass(*i);
vector<LorentzMomentum> p = CMSn(m0, masses);
j = 0;
for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j )
Traits::set5Momentum(*i, p[j]);
}
}
|