diff --git a/.hgignore b/.hgignore --- a/.hgignore +++ b/.hgignore @@ -1,71 +1,75 @@ Makefile$ Makefile\.in$ \.deps$ \.libs$ \.l[ao]$ \.so$ \.so\. \.o$ ~$ \.orig$ ^include/.done-all-links$ ^lib/.done-all-links$ ^src/.*\.(run|tex|out|log|rpo|spc|top|dump|dot|aux|pdf|ps|png|svg|hepmc)$ ^src/.done-all-links$ ^autom4te.cache$ ^config.thepeg$ ^config.log$ ^config.status$ ^configure$ ^Config/config.h$ ^Config/config.h.in$ ^Config/config.sub$ ^Config/depcomp$ ^Config/install-sh$ ^Config/compile$ ^Config/missing$ ^Config/stamp-h.$ ^Config/ar-lib$ ^Config/config.guess$ ^Config/test-driver$ ^Doc/fixinterfaces.pl$ ^include/ThePEG$ ^Doc/MakeDocs.in$ ^Doc/refman.conf$ ^Doc/refman.h$ ^Doc/refman-html$ ^Doc/AllInterfaces.h$ ^Doc/MoreInterfaces.h$ ^Doc/ThePEG-refman.tag$ ^java/.*\.(java|class)$ ^java/ThePEG$ ^java/ThePEG.jar$ ^java/thepeg.sh$ ^java/thepeg$ ^PDF/.done-all-links$ ^lib/ThePEGDefaults.rpo$ ^lib/Makefile.common ^lib/Makefile.dist$ ^src/runThePEG$ ^src/setupThePEG$ ^aclocal.m4$ ^libtool$ ^INSTALL$ ^src/TestDecayMode.in$ ^ThePEG.*\.tar\.(bz2|gz)$ ^src/runThePEG.bin$ ^src/setupThePEG.bin$ ^ThePEG-default.kdev4$ ^src/SimpleLEP.cmp$ ^Repository/versionstamp.inc$ ^Repository/repository_test(.log|.trs)?$ ^Repository/test-suite.log$ ^Repository/tests/.dirstamp$ # added by aHg on Wed Jan 13 13:27:00 2016 syntax: glob Config/LWH.h PDF/done-all-links include/done-all-links lib/done-all-links src/done-all-links + +# added by aHg on Tue Oct 11 14:21:42 2016 +syntax: glob +Config/ThePEG_Qty.h diff --git a/MatrixElement/ME2to2Base.cc b/MatrixElement/ME2to2Base.cc --- a/MatrixElement/ME2to2Base.cc +++ b/MatrixElement/ME2to2Base.cc @@ -1,202 +1,205 @@ // -*- C++ -*- // // ME2to2Base.cc 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, non-templated member // functions of the ME2to2Base class. // #include "ME2to2Base.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Cuts/Cuts.h" using namespace ThePEG; ME2to2Base::~ME2to2Base() {} Energy2 ME2to2Base::scale() const { switch ( scaleChoice() ) { case 1: return -tHat()*uHat()/(tHat() + uHat()); default: return tHat()*uHat()/sHat(); } } void ME2to2Base::setKinematics() { MEBase::setKinematics(); theLastTHat = (meMomenta()[0] - meMomenta()[2]).m2(); theLastUHat = (meMomenta()[1] - meMomenta()[2]).m2(); theLastPhi = meMomenta()[2].phi(); } bool ME2to2Base::generateKinematics(const double * r) { // generate the masses of the particles for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass()); } double ctmin = -1.0; double ctmax = 1.0; - Energy q = ZERO; - try { - q = SimplePhaseSpace:: - getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); - } catch ( ImpossibleKinematics ) { - return false; - } + // Energy q = ZERO; + // try { + // q = SimplePhaseSpace:: + // getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); + // } catch ( ImpossibleKinematics ) { + // return false; + // } + Energy q = SimplePhaseSpace:: + checkMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); + if ( q < ZERO ) return false; Energy e = sqrt(sHat())/2.0; Energy2 m22 = meMomenta()[2].mass2(); Energy2 m32 = meMomenta()[3].mass2(); Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22); Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32); Energy2 pq = 2.0*e*q; Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]); if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq); thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]); if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq); thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]); if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq); Energy ptmin = max(lastCuts().minKT(mePartonData()[2]), lastCuts().minKT(mePartonData()[3])); if ( ptmin > ZERO ) { double ctm = 1.0 - sqr(ptmin/q); if ( ctm <= 0.0 ) return false; ctmin = max(ctmin, -sqrt(ctm)); ctmax = min(ctmax, sqrt(ctm)); } double ymin2 = lastCuts().minYStar(mePartonData()[2]); double ymax2 = lastCuts().maxYStar(mePartonData()[2]); double ymin3 = lastCuts().minYStar(mePartonData()[3]); double ymax3 = lastCuts().maxYStar(mePartonData()[3]); double ytot = lastCuts().Y() + lastCuts().currentYHat(); if ( ymin2 + ytot > -0.9*Constants::MaxRapidity ) ctmin = max(ctmin, sqrt(sqr(q) + m22)*tanh(ymin2)/q); if ( ymax2 + ytot < 0.9*Constants::MaxRapidity ) ctmax = min(ctmax, sqrt(sqr(q) + m22)*tanh(ymax2)/q); if ( ymin3 + ytot > -0.9*Constants::MaxRapidity ) ctmax = min(ctmax, sqrt(sqr(q) + m32)*tanh(-ymin3)/q); if ( ymax3 + ytot < 0.9*Constants::MaxRapidity ) ctmin = max(ctmin, sqrt(sqr(q) + m32)*tanh(-ymax3)/q); if ( ctmin >= ctmax ) return false; double cth = getCosTheta(ctmin, ctmax, r); Energy pt = q*sqrt(1.0-sqr(cth)); theLastPhi = rnd(2.0*Constants::pi); meMomenta()[2].setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth)); meMomenta()[3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth)); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); vector out(2); out[0] = meMomenta()[2]; out[1] = meMomenta()[3]; tcPDVector tout(2); tout[0] = mePartonData()[2]; tout[1] = mePartonData()[3]; if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) ) return false; theLastTHat = pq*cth + m22 - e0e2; theLastUHat = m22 + m32 - sHat() - theLastTHat; jacobian((pq/sHat())*Constants::pi*jacobian()); return true; } double ME2to2Base::getCosTheta(double ctmin, double ctmax, const double * r) { double cth = 0.0; static const double eps = 1.0e-6; if ( 1.0 + ctmin <= eps && 1.0 - ctmax <= eps ) { jacobian(ctmax - ctmin); cth = ctmin + (*r)*jacobian(); } else if ( 1.0 + ctmin <= eps ) { cth = 1.0 - (1.0 - ctmax)*pow((1.0 - ctmin)/(1.0 - ctmax), *r); jacobian(log((1.0 - ctmin)/(1.0 - ctmax))*(1.0 - cth)); } else if ( 1.0 - ctmax <= eps ) { cth = -1.0 + (1.0 + ctmin)*pow((1.0 + ctmax)/(1.0 + ctmin), *r); jacobian(log((1.0 + ctmax)/(1.0 + ctmin))*(1.0 + cth)); } else { double zmin = 0.5*(1.0 - ctmax); double zmax = 0.5*(1.0 - ctmin); double A1 = -ctmin/(zmax*(1.0-zmax)); double A0 = -ctmax/(zmin*(1.0-zmin)); double A = *r*(A1 - A0) + A0; double z = A < 2.0? 2.0/(sqrt(sqr(A) + 4.0) + 2 - A): 0.5*(A - 2.0 + sqrt(sqr(A) + 4.0))/A; cth = 1.0 - 2.0*z; jacobian(2.0*(A1 - A0)*sqr(z)*sqr(1.0 - z)/(sqr(z) + sqr(1.0 - z))); } return cth; } CrossSection ME2to2Base::dSigHatDR() const { return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc); } void ME2to2Base::persistentOutput(PersistentOStream & os) const { os << theScaleChoice << ounit(theLastTHat, GeV2) << ounit(theLastUHat, GeV2) << theLastPhi; } void ME2to2Base::persistentInput(PersistentIStream & is, int) { is >> theScaleChoice >> iunit(theLastTHat, GeV2) >> iunit(theLastUHat, GeV2) >> theLastPhi; } AbstractClassDescription ME2to2Base::initME2to2Base; // Definition of the static class description member. Switch & ME2to2Base::interfaceScaleChoice() { static Switch dummy ("ScaleChoice", "Different options for calculating the scale of the generated " "hard sub-process.", &ME2to2Base::theScaleChoice, 0, false, false); return dummy; } void ME2to2Base::Init() { static ClassDocumentation documentation ("The ThePEG::ME2to2Base class may be used as a base class " "for all \\f$2\\rightarrow 2\\f$ matrix elements."); static SwitchOption interfaceScaleChoice0 (interfaceScaleChoice(), "that.uhat/shat", "\\f$\\hat{t}\\hat{u}/\\hat{s}\\f$", 0); static SwitchOption interfaceScaleChoice1 (interfaceScaleChoice(), "that.uhat/(that+uhat)", "\\f$-\\hat{t}\\hat{u}/(\\hat{t}+\\hat{u})\\f$", 1); } diff --git a/Utilities/SimplePhaseSpace.cc b/Utilities/SimplePhaseSpace.cc --- a/Utilities/SimplePhaseSpace.cc +++ b/Utilities/SimplePhaseSpace.cc @@ -1,92 +1,112 @@ // -*- C++ -*- // // SimplePhaseSpace.cc 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. // #include "SimplePhaseSpace.h" #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "SimplePhaseSpace.tcc" #endif using namespace ThePEG; Energy SimplePhaseSpace::getMagnitude(Energy2 s, Energy m1, Energy m2) { const Energy2 eps = 10.0*s*Constants::epsilon; if ( m1 < ZERO && sqr(m1) < eps ) m1 = ZERO; if ( m2 < ZERO && sqr(m2) < eps ) m2 = ZERO; if ( m1 >= ZERO && m2 >= ZERO ) { Energy2 aa = s - sqr(m1+m2); if ( aa < ZERO && aa > -eps ) return ZERO; if ( aa < ZERO ) throw ImpossibleKinematics(); return 0.5*sqrt(aa*(s-sqr(m1-m2))/s); } Energy2 m12 = m1 < ZERO? -sqr(m1): sqr(m1); Energy2 m22 = m2 < ZERO? -sqr(m2): sqr(m2); Energy2 r2 = 0.25*(sqr(m12) + sqr(m22 - s) -2.0*m12*(m22 + s))/s; if ( r2 < ZERO || r2 + m12 < ZERO || r2 + m22 < ZERO ) throw ImpossibleKinematics(); return sqrt(r2); } +Energy SimplePhaseSpace::checkMagnitude(Energy2 s, Energy m1, Energy m2) +{ + if ( s < ZERO ) return -1.0*GeV; + const Energy2 eps = 10.0*s*Constants::epsilon; + if ( m1 < ZERO && sqr(m1) < eps ) m1 = ZERO; + if ( m2 < ZERO && sqr(m2) < eps ) m2 = ZERO; + if ( m1 >= ZERO && m2 >= ZERO ) { + Energy2 aa = s - sqr(m1+m2); + if ( aa < ZERO && aa > -eps ) return ZERO; + if ( aa < ZERO ) return -1.0*GeV; + return 0.5*sqrt(aa*(s-sqr(m1-m2))/s); + } + Energy2 m12 = m1 < ZERO? -sqr(m1): sqr(m1); + Energy2 m22 = m2 < ZERO? -sqr(m2): sqr(m2); + Energy2 r2 = 0.25*(sqr(m12) + sqr(m22 - s) -2.0*m12*(m22 + s))/s; + if ( r2 < ZERO || r2 + m12 < ZERO || r2 + m22 < ZERO ) + return -1.0*GeV; + return sqrt(r2); +} + vector SimplePhaseSpace:: CMSn(Energy m0, const vector & m) { using Constants::pi; // Setup constants. int Np = m.size(); vector ret(Np); Energy summ = std::accumulate(m.begin(), m.end(), Energy()); if ( summ >= m0 ) throw ImpossibleKinematics(); while ( true ) { // First get an ordered list of random numbers. vector rndv(Np); rndv[0] = 1.0; rndv.back() = 0.0; for ( int i = 1; i < Np - 1; ++i ) rndv[i] = UseRandom::rnd(); std::sort(rndv.begin() + 1, rndv.end() - 1, std::greater()); // Now setup masses of subsystems. vector sm(Np); Energy tmass = m0 - summ; Energy tmp = summ; for ( int i = 0; i < Np; ++i ) { sm[i] = rndv[i]*tmass + tmp; tmp -= m[i]; } // Now the magnitude of all the momenta can be calculated. This // gives the weight. double weight = 1.0; vector p(Np); p[Np - 1] = getMagnitude(sqr(sm[Np - 2]), m[Np -2], sm[Np - 1]); for ( int i = Np - 2; i >= 0; --i ) weight *= (p[i] = getMagnitude(sqr(sm[i]), m[i], sm[i + 1]))/sm[i]; if ( weight > UseRandom::rnd() ) continue; // Now we just have to generate the angles. ret[Np - 1] = LorentzMomentum(ZERO, ZERO, ZERO, m[Np - 1]); for ( int i = Np - 2; i >= 0; --i ) { Momentum3 p3 = polar3Vector(p[i], 2.0*UseRandom::rnd() - 1.0, 2.0*pi*UseRandom::rnd()); ret[i] = LorentzMomentum(-p3, sqrt(sqr(p[i]) + sqr(m[i]))); if ( i == Np -2 ) { ret[Np - 1] = LorentzMomentum(p3, sqrt(sqr(m[Np - 1]) + p3.mag2())); } else { Boost bv = p3*(1.0/sqrt(sqr(p[i]) + sqr(sm[i + 1]))); if ( bv.mag2() >= 1.0 ) throw ImpossibleKinematics(); LorentzRotation r(bv); for ( int j = i + 1; j < Np; ++j ) ret[j]*=r.one(); } } return ret; } } diff --git a/Utilities/SimplePhaseSpace.h b/Utilities/SimplePhaseSpace.h --- a/Utilities/SimplePhaseSpace.h +++ b/Utilities/SimplePhaseSpace.h @@ -1,232 +1,244 @@ // -*- C++ -*- // // SimplePhaseSpace.h 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. // #ifndef ThePEG_SimplePhaseSpace_H #define ThePEG_SimplePhaseSpace_H #include "ThePEG/Config/ThePEG.h" #include "ThePEG/Vectors/LorentzRotation.h" #include "ThePEG/Vectors/LorentzRotation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/ParticleTraits.h" #include "ThePEG/Repository/UseRandom.h" #include "SimplePhaseSpace.xh" #include namespace ThePEG { /** * SimplePhaseSpace defines a set of static functions to be used for * distributing momenta evenly in phase space. In most cases pointers * and references to both particle and momentum objects can be used as * arguments as long as the ParticleTraits class is specialized * properly. When needed, random numbers are generated with the * generator given by the static UseRandom class. */ struct SimplePhaseSpace { /** * Set two momenta in their center of mass system. Their total * invariant mass squared is given by s, and their direction is * distributed isotropically. * @param s the total invariant mass squared. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(Energy2 s, PType & p1, PType & p2); /** * Set two momenta in their center of mass system. Their total * invariant mass squared is given by s, and their direction is * given in terms of the polar and azimuth angle of the first * momenta. * @param s the total invariant mass squared. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param cosTheta cosine of the azimuth angle of the first momentum. * @param phi azimuth angle of the first momentum. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(PType & p1, PType & p2, Energy2 s, double cosTheta, double phi); /** * Set two momenta in their center of mass system. Their total * invariant mass squared is given by s. The helper momentum p0 is * used so that afterwards \f$t=(p0-p1)^2\f$ and p1 has the azimuth * angle phi around p0. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param s the total invariant mass squared. * @param t \f$=(p0-p1)^2\f$. * @param phi azimuth angle of the first momentum around p0. * @param p0 pointer or reference to an auxiliary momentum. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(PType & p1, PType & p2, Energy2 s, Energy2 t, double phi, const PType & p0); /** * Set two momenta in their center of mass system. Their total * invariant mass squared is given by s. p1 will be along the z-axis. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param s the total invariant mass squared. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(PType & p1, PType & p2, Energy2 s); /** * Set two momenta in their center of mass system. Their total * invariant mass squared is given by s. The first will be along the * z-axis. * @param p a pair of pointers or references to the two momenta. Their * invariant masses will be preserved. * @param s the total invariant mass squared. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(const PPairType & p, Energy2 s) { CMS(*p.first, *p.second, s); } /** * Set three momenta in their center of mass system. Their total * invariant mass squared is given by s. The energy fraction of * particle p1(3) is x1(3) of the total energy and the angles of the * system is distributed isotropically. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param p3 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param s the total invariant mass squared. * @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$. * @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(PType & p1, PType & p2, PType & p3, Energy2 s, double x1, double x3); /** * Set three momenta in their center of mass system. Their total * invariant mass squared is given by s. The energy fraction of * particle p1(3) is x1(3) of the total energy. Particle p1 is * initially placed along the z-axis and particle p2 is given * azimuth angle phii. Then the system is then rotated with * theta and phi respectively. * @param p1 pointer or reference to the first momentum. Its * invariant mass will be preserved. * @param p2 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param p3 pointer or reference to the second momentum. Its * invariant mass will be preserved. * @param s the total invariant mass squared. * @param x1 the energy fraction \f$2e_1/\sqrt{s}\f$. * @param x3 the energy fraction \f$2e_3/\sqrt{s}\f$. * @param phii the azimuth angle of p2 around p1. * @param theta the polar angle of p1. * @param phi the azimuth angle of p1. * @throw ImpossibleKinematics if the sum of the invariant masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMS(PType & p1, PType & p2, PType & p3, Energy2 s, double x1, double x3, double phii = 0.0, double theta = 0.0, double phi = 0.0); /** * Calculate the absolute magnitude of the momenta of two particles * with masses m1 and m2 when put in their CMS of total invariant * mass squared s. * @param s the total invariant mass squared. * @param m1 the mass of particle 1. * @param m2 the mass of particle 2. * @throw ImpossibleKinematics if the sum of the masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ static Energy getMagnitude(Energy2 s, Energy m1, Energy m2); /** + * Calculate the absolute magnitude of the momenta of two particles + * with masses m1 and m2 when put in their CMS of total invariant + * mass squared s. + * @param s the total invariant mass squared. + * @param m1 the mass of particle 1. + * @param m2 the mass of particle 2. + * @return a negative value if the sum of the masses was + * larger than the given invariant mass (\f$\sqrt{s}\f$). + */ + static Energy checkMagnitude(Energy2 s, Energy m1, Energy m2); + + /** * Return a three-vector given the absolute momentum, cos(theta) and * phi. * @param p the magnitude of the momentum. * @param costheta the cosine of the polar angle. * @param phi the azimuth angle. */ static Momentum3 polar3Vector(Energy p, double costheta, double phi) { return Momentum3(p*sqrt(1.0 - sqr(costheta))*sin(phi), p*sqrt(1.0 - sqr(costheta))*cos(phi), p*costheta); } /** * Get a number of randomly distributed momenta. * Given a number specified invariant masses and a * total invariant mass m0, return corresponding four-momenta * randomly distributed according to phase space. * @param m0 the * total invariant mass of the resulting momenta. * @param m a vector * of invariant masses of the resulting momenta. * @return a vector * of momenta with the given masses randomly distributed. * @throw ImpossibleKinematics if the sum of the masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ static vector CMSn(Energy m0, const vector & m); /** * Set the momentum of a number of particles. Given a number of * particles and a total invariant mass m0, distribute their * four-momenta randomly according to phase space. * @param particles a container of particles or pointers to * particles. The invariant mass of these particles will not be * chaned. * @param m0 the * total invariant mass of the resulting momenta. * @throw ImpossibleKinematics if the sum of the masses was * larger than the given invariant mass (\f$\sqrt{s}\f$). */ template static void CMSn(Container & particles, Energy m0); }; } #ifndef ThePEG_TEMPLATES_IN_CC_FILE #include "SimplePhaseSpace.tcc" #endif #endif /* ThePEG_SimplePhaseSpace_H */