Warning, /include/ThePEG/Utilities/SimplePhaseSpace.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // SimplePhaseSpace.tcc is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 //
0010 // This is the implementation of the non-inlined templated member
0011 // functions of the SimplePhaseSpace class.
0012 //
0013
0014 namespace ThePEG {
0015
0016 template <typename PType>
0017 void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s)
0018 {
0019 typedef ParticleTraits<PType> Traits;
0020 Energy m1 = Traits::mass(p1); Energy m2 = Traits::mass(p2);
0021 Energy z = getMagnitude(s, m1, m2);
0022 Energy2 m12 = m1 >= ZERO ? sqr(m1) : -sqr(m1);
0023 Energy2 m22 = m2 >= ZERO ? sqr(m2) : -sqr(m2);
0024 Energy2 c1 = (s+m12-m22);
0025 Energy2 c2 = (s-m12+m22);
0026 Traits::set5Momentum(p1, Lorentz5Momentum(ZERO, ZERO, z, (c1 > ZERO ? 1. : -1.) * 0.5*sqrt(sqr(c1)/s), m1));
0027 Traits::set5Momentum(p2, Lorentz5Momentum(ZERO, ZERO, -z, (c2 > ZERO ? 1. : -1.) * 0.5*sqrt(sqr(c2)/s), m2));
0028 }
0029
0030 template <typename PType>
0031 void SimplePhaseSpace::CMS(Energy2 s, PType & p1, PType & p2)
0032 {
0033 CMS(p1, p2, s, 2.0*UseRandom::rnd() - 1.0, Constants::twopi*UseRandom::rnd());
0034 }
0035
0036 template <typename PType>
0037 void SimplePhaseSpace::CMS(PType & p1, PType & p2, Energy2 s,
0038 double cthe, double phi)
0039 {
0040 typedef ParticleTraits<PType> Traits;
0041 Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2));
0042 double sthe = sqrt(1.0-sqr(cthe));
0043 Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe);
0044 Traits::set3Momentum(p1, p);
0045 Traits::set3Momentum(p2, -p);
0046 }
0047
0048 template <typename PType>
0049 void SimplePhaseSpace::
0050 CMS(PType & p1, PType & p2, PType & p3,
0051 Energy2 s, double x1, double x3)
0052 {
0053 CMS(p1, p2, p3, s, x1, x3,
0054 Constants::twopi*UseRandom::rnd(), acos(2.0*UseRandom::rnd() - 1.0),
0055 Constants::twopi*UseRandom::rnd());
0056 }
0057
0058 template <typename PType>
0059 void SimplePhaseSpace::
0060 CMS(PType & p1, PType & p2, PType & p3, Energy2 s,
0061 double x1, double x3, double phii, double the, double phi)
0062 {
0063 typedef ParticleTraits<PType> Traits;
0064 Energy Etot = sqrt(s);
0065 Energy m1 = Traits::mass(p1);
0066 Energy m2 = Traits::mass(p2);
0067 Energy m3 = Traits::mass(p3);
0068 Energy e1 = 0.5*x1*Etot;
0069 Energy e3 = 0.5*x3*Etot;
0070 Energy e2 = Etot - e1 - e3;
0071 if ( e1 < m1 || e2 < m2 || e3 < m3 ) throw ImpossibleKinematics();
0072 Energy r1 = sqrt(sqr(e1)-sqr(m1));
0073 Energy r2 = sqrt(sqr(e2)-sqr(m2));
0074 Energy r3 = sqrt(sqr(e3)-sqr(m3));
0075 Traits::set3Momentum(p1, Momentum3(ZERO, ZERO, r1));
0076 double cthe2 = (sqr(r3)-sqr(r2)-sqr(r1))/(2.0*r2*r1);
0077 double cthe3 = (sqr(r2)-sqr(r3)-sqr(r1))/(2.0*r3*r1);
0078 if ( abs(cthe2) > 1.0 || abs(cthe3) > 1.0 ) throw ImpossibleKinematics();
0079 double sthe2 = sqrt(1.0-sqr(cthe2));
0080 Energy px = r2*sthe2*cos(phii);
0081 Energy py = r2*sthe2*sin(phii);
0082 Traits::set3Momentum(p2, Momentum3(px, py, r2*cthe2));
0083 Traits::set3Momentum(p3, Momentum3(-px, -py, r3*cthe3));
0084 if ( the == 0.0 && phi == 0.0 ) return;
0085 LorentzRotation r;
0086 r.rotateZ(phi);
0087 r.rotateX(the);
0088 Traits::transform(p1, r);
0089 Traits::transform(p2, r);
0090 Traits::transform(p3, r);
0091 }
0092
0093 template <typename PType>
0094 void SimplePhaseSpace::
0095 CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi,
0096 const PType & p0) {
0097 typedef ParticleTraits<PType> Traits;
0098 Energy r = getMagnitude(s, Traits::mass(p1), Traits::mass(p2));
0099 Energy e = sqrt(sqr(r) + sqr(Traits::mass(p1)));
0100 Energy r0 = Traits::momentum(p0).rho();
0101 Energy e0 = Traits::momentum(p0).e();
0102 double cthe = (t + sqr(e - e0) + sqr(r) + sqr(r0))/(2.0*r*r0);
0103 if ( abs(cthe) > 1.0 ) throw ImpossibleKinematics();
0104 double sthe = sqrt(1.0-sqr(cthe));
0105 Momentum3 p(r*sthe*cos(phi), r*sthe*sin(phi), r*cthe);
0106 Traits::set3Momentum(p1, p);
0107 Traits::set3Momentum(p2, -p);
0108 if ( Traits::momentum(p0).perp2() > ZERO ) {
0109 LorentzRotation r;
0110 r.rotateX(Traits::momentum(p0).theta());
0111 r.rotateZ(Traits::momentum(p0).phi());
0112 Traits::transform(p1, r);
0113 Traits::transform(p2, r);
0114 }
0115
0116 }
0117
0118 template <typename Container>
0119 void SimplePhaseSpace::CMSn(Container & particles, Energy m0)
0120 {
0121 typedef typename Container::value_type PType;
0122 typedef typename Container::iterator Iterator;
0123 if ( particles.size() == 2 ) {
0124 Iterator it = particles.begin();
0125 PType & p1 = *it++;
0126 PType & p2 = *it;
0127 CMS(sqr(m0), p1, p2);
0128 return;
0129 }
0130 typedef ParticleTraits<PType> Traits;
0131 vector<Energy> masses(particles.size());
0132 int j = 0;
0133 for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j )
0134 masses[j] = Traits::mass(*i);
0135 vector<LorentzMomentum> p = CMSn(m0, masses);
0136 j = 0;
0137 for ( Iterator i = particles.begin();i != particles.end(); ++i, ++j )
0138 Traits::set5Momentum(*i, p[j]);
0139 }
0140
0141 }
