include/Pythia8/PhaseSpace.h

0001 // PhaseSpace.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Header file for phase space generators in kinematics selection.
0007 // PhaseSpace: base class for phase space generators.
0008 // Base class for derived classes> 2 ->1 , 2 -> 2, 2 -> 2 elastic/diffractive,
0009 // 2 -> 2 nondiffractive, 2 -> 3, Les Houches.
0010
0011 #ifndef Pythia8_PhaseSpace_H
0012 #define Pythia8_PhaseSpace_H
0013
0014 #include "Pythia8/Basics.h"
0015 #include "Pythia8/BeamParticle.h"
0016 #include "Pythia8/Info.h"
0017 #include "Pythia8/LesHouches.h"
0018 #include "Pythia8/MultipartonInteractions.h"
0019 #include "Pythia8/ParticleData.h"
0020 #include "Pythia8/PartonDistributions.h"
0021 #include "Pythia8/PhysicsBase.h"
0022 #include "Pythia8/PythiaStdlib.h"
0023 #include "Pythia8/SigmaProcess.h"
0024 #include "Pythia8/SigmaTotal.h"
0025 #include "Pythia8/Settings.h"
0026 #include "Pythia8/StandardModel.h"
0027 #include "Pythia8/UserHooks.h"
0028 #include "Pythia8/GammaKinematics.h"
0029
0030 namespace Pythia8 {
0031
0032 //==========================================================================
0033
0034 // Forward reference to the UserHooks class.
0035 class UserHooks;
0036
0037 //==========================================================================
0038
0039 // PhaseSpace is a base class for  phase space generators
0040 // used in the selection of hard-process kinematics.
0041
0042 class PhaseSpace : public PhysicsBase {
0043
0044 public:
0045
0046   // Destructor.
0047   virtual ~PhaseSpace() {}
0048
0049   // Perform simple initialization and store pointers.
0050   void init(bool isFirst, SigmaProcessPtr sigmaProcessPtrIn);
0051
0052   // Switch to new beam particle identities; for similar hadrons only.
0053   void updateBeamIDs() { idAold = idA; idBold = idB; idA = beamAPtr->id();
0054     idB = beamBPtr->id(); mA = beamAPtr->m(); mB = beamBPtr->m();
0055     sigmaProcessPtr->updateBeamIDs();}
0056
0057   // Update the CM energy of the event.
0058   void newECM(double eCMin) {eCM = eCMin; s = eCM * eCM;}
0059
0060   // Store or replace Les Houches pointer.
0061   void setLHAPtr(LHAupPtr lhaUpPtrIn) {lhaUpPtr = lhaUpPtrIn;}
0062
0063   // A pure virtual method, wherein an optimization procedure
0064   // is used to determine how phase space should be sampled.
0065   virtual bool setupSampling() = 0;
0066
0067   // A pure virtual method, wherein a trial event kinematics
0068   // is to be selected in the derived class.
0069   virtual bool trialKin(bool inEvent = true, bool repeatSame = false) = 0;
0070
0071   // A pure virtual method, wherein the accepted event kinematics
0072   // is to be constructed in the derived class.
0073   virtual bool finalKin() = 0;
0074
0075   // Allow for nonisotropic decays when ME's available.
0076   void   decayKinematics( Event& process);
0077
0078   // Give back current or maximum cross section, or set latter.
0079   double sigmaNow() const {return sigmaNw;}
0080   double sigmaMax() const {return sigmaMx;}
0081   double biasSelectionWeight()  const {return biasWt;}
0082   bool   newSigmaMax() const {return newSigmaMx;}
0083   void   setSigmaMax(double sigmaMaxIn) {sigmaMx = sigmaMaxIn;}
0084
0085   // New value for switched beam identity or energy (for SoftQCD processes).
0086   double sigmaMaxSwitch() {sigmaNw = sigmaProcessPtr->sigmaHatWrap();
0087     sigmaMx = sigmaNw; return sigmaMx;}
0088
0089   // For Les Houches with negative event weight needs
0090   virtual double sigmaSumSigned() const {return sigmaMx;}
0091
0092   // Give back constructed four-vectors and known masses.
0093   Vec4   p(int i)   const {return pH[i];}
0094   double m(int i)   const {return mH[i];}
0095
0096   // Reset the four-momentum.
0097   void   setP(int i, Vec4 pNew) { pH[i] = pNew; }
0098
0099   // Give back other event properties.
0100   double ecm()      const {return eCM;}
0101   double x1()       const {return x1H;}
0102   double x2()       const {return x2H;}
0103   double sHat()     const {return sH;}
0104   double tHat()     const {return tH;}
0105   double uHat()     const {return uH;}
0106   double pTHat()    const {return pTH;}
0107   double thetaHat() const {return theta;}
0108   double phiHat()   const {return phi;}
0109   double runBW3()   const {return runBW3H;}
0110   double runBW4()   const {return runBW4H;}
0111   double runBW5()   const {return runBW5H;}
0112
0113   // Inform whether beam particles are resolved in partons or scatter directly.
0114   virtual bool isResolved() const {return true;}
0115
0116   // Functions to rescale momenta and cross section for new sHat
0117   // Currently implemented only for PhaseSpace2to2tauyz class.
0118   virtual void rescaleSigma( double){}
0119   virtual void rescaleMomenta( double){}
0120
0121   // Calculate the weight for over-estimated cross section.
0122   virtual double weightGammaPDFApprox(){ return 1.;}
0123
0124   // Set the GammaKinematics pointer needed for soft photoproduction.
0125   virtual void setGammaKinPtr( GammaKinematics* gammaKinPtrIn) {
0126     gammaKinPtr = gammaKinPtrIn; }
0127
0128 protected:
0129
0130   // Constructor.
0131   PhaseSpace() : sigmaProcessPtr(), lhaUpPtr(), gammaKinPtr(),
0132     useBreitWigners(), doEnergySpread(), showSearch(), showViolation(),
0133     increaseMaximum(), hasQ2Min(), gmZmodeGlobal(), mHatGlobalMin(),
0134     mHatGlobalMax(), pTHatGlobalMin(), pTHatGlobalMax(), Q2GlobalMin(),
0135     pTHatMinDiverge(), minWidthBreitWigners(), minWidthNarrowBW(), idA(),
0136     idB(), idAold(), idBold(), idAgm(), idBgm(), mA(), mB(), eCM(), s(),
0137     sigmaMxGm(), hasLeptonBeamA(), hasLeptonBeamB(), hasOneLeptonBeam(),
0138     hasTwoLeptonBeams(), hasPointGammaA(), hasPointGammaB(),
0139     hasOnePointParticle(), hasTwoPointParticles(), hasGamma(), hasVMD(),
0140     newSigmaMx(), canModifySigma(), canBiasSelection(), canBias2Sel(),
0141     gmZmode(), bias2SelPow(), bias2SelRef(), wtBW(), sigmaNw(),
0142     sigmaMx(), sigmaPos(), sigmaNeg(), biasWt(), mHatMin(), mHatMax(),
0143     sHatMin(), sHatMax(), pTHatMin(), pTHatMax(), pT2HatMin(), pT2HatMax(),
0144     x1H(), x2H(), m3(), m4(), m5(), s3(), s4(), s5(), mHat(), sH(), tH(), uH(),
0145     pAbs(), p2Abs(), pTH(), theta(), phi(), betaZ(), mH(), idResA(), idResB(),
0146     mResA(), mResB(), GammaResA(), GammaResB(), tauResA(), tauResB(),
0147     widResA(), widResB(), sameResMass(), useMirrorWeight(), hasNegZ(),
0148     hasPosZ(), tau(), y(), z(), tauMin(), tauMax(), yMax(), zMin(), zMax(),
0149     ratio34(), unity34(), zNeg(), zPos(), wtTau(), wtY(), wtZ(), wt3Body(),
0150     runBW3H(), runBW4H(), runBW5H(), intTau0(), intTau1(), intTau2(),
0151     intTau3(), intTau4(), intTau5(), intTau6(), intY0(), intY12(), intY34(),
0152     intY56(), mTchan1(), sTchan1(), mTchan2(), sTchan2(), frac3Flat(),
0153     frac3Pow1(), frac3Pow2(), zNegMin(), zNegMax(), zPosMin(), zPosMax(),
0154     nTau(), nY(), nZ(), tauCoef(), yCoef(), zCoef(), tauCoefSum(), yCoefSum(),
0155     zCoefSum(), useBW(), useNarrowBW(), idMass(), mPeak(), sPeak(), mWidth(),
0156     mMin(), mMax(), mw(), wmRat(), mLower(), mUpper(), sLower(), sUpper(),
0157     fracFlatS(), fracFlatM(), fracInv(), fracInv2(), atanLower(), atanUpper(),
0158     intBW(), intFlatS(), intFlatM(), intInv(), intInv2() {}
0159
0160   // Constants: could only be changed in the code itself.
0161   static const int    NMAXTRY, NTRY3BODY;
0162   static const double SAFETYMARGIN, TINY, EVENFRAC, SAMESIGMA, MRESMINABS,
0163                       WIDTHMARGIN, SAMEMASS, MASSMARGIN, EXTRABWWTMAX,
0164                       THRESHOLDSIZE, THRESHOLDSTEP, YRANGEMARGIN, LEPTONXMIN,
0165                       LEPTONXMAX, LEPTONXLOGMIN, LEPTONXLOGMAX, LEPTONTAUMIN,
0166                       SHATMINZ, PT2RATMINZ, WTCORRECTION[11];
0167
0168   // Pointer to cross section.
0169   SigmaProcessPtr sigmaProcessPtr;
0170
0171   // Pointer to LHAup for generating external events.
0172   LHAupPtr      lhaUpPtr;
0173
0174   // Pointer to object that samples photon kinematics from leptons.
0175   GammaKinematics* gammaKinPtr;
0176
0177   // Initialization data, normally only set once.
0178   bool   useBreitWigners, doEnergySpread, showSearch, showViolation,
0179          increaseMaximum, hasQ2Min;
0180   int    gmZmodeGlobal;
0181   double mHatGlobalMin, mHatGlobalMax, pTHatGlobalMin, pTHatGlobalMax,
0182          Q2GlobalMin, pTHatMinDiverge, minWidthBreitWigners, minWidthNarrowBW;
0183
0184   // Information on incoming beams.
0185   int    idA, idB, idAold, idBold, idAgm, idBgm;
0186   double mA, mB, eCM, s, sigmaMxGm;
0187   bool   hasLeptonBeamA, hasLeptonBeamB, hasOneLeptonBeam, hasTwoLeptonBeams,
0188          hasPointGammaA, hasPointGammaB, hasOnePointParticle,
0189          hasTwoPointParticles, hasGamma, hasVMD;
0190
0191   // Cross section information.
0192   bool   newSigmaMx, canModifySigma, canBiasSelection, canBias2Sel;
0193   int    gmZmode;
0194   double bias2SelPow, bias2SelRef, wtBW, sigmaNw, sigmaMx, sigmaPos,
0195          sigmaNeg, biasWt;
0196
0197   // Process-specific kinematics properties, almost always available.
0198   double mHatMin, mHatMax, sHatMin, sHatMax, pTHatMin, pTHatMax,
0199          pT2HatMin, pT2HatMax;
0200
0201   // Event-specific kinematics properties, almost always available.
0202   double x1H, x2H, m3, m4, m5, s3, s4, s5, mHat, sH, tH, uH, pAbs, p2Abs,
0203          pTH, theta, phi, betaZ;
0204   Vec4   pH[12];
0205   double mH[12];
0206
0207   // Reselect decay products momenta isotropically in phase space.
0208   void decayKinematicsStep( Event& process, int iRes);
0209
0210   // Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
0211
0212   // Determine how phase space should be sampled.
0213   void setup3Body();
0214   bool setupSampling123(bool is2, bool is3);
0215
0216   // Select a trial kinematics phase space point.
0217   bool trialKin123(bool is2, bool is3, bool inEvent = true);
0218
0219   // Presence and properties of any s-channel resonances.
0220   int    idResA, idResB;
0221   double mResA, mResB, GammaResA, GammaResB, tauResA, tauResB, widResA,
0222          widResB;
0223   bool   sameResMass;
0224
0225   // Kinematics properties specific to 2 -> 1/2/3.
0226   bool   useMirrorWeight, hasNegZ, hasPosZ;
0227   double tau, y, z, tauMin, tauMax, yMax, zMin, zMax, ratio34, unity34,
0228          zNeg, zPos, wtTau, wtY, wtZ, wt3Body, runBW3H, runBW4H, runBW5H,
0229          intTau0, intTau1, intTau2, intTau3, intTau4, intTau5, intTau6,
0230          intY0, intY12, intY34, intY56, mTchan1, sTchan1, mTchan2, sTchan2,
0231          frac3Flat, frac3Pow1, frac3Pow2, zNegMin, zNegMax, zPosMin, zPosMax;
0232   Vec4   p3cm, p4cm, p5cm;
0233
0234   // Coefficients for optimized selection in 2 -> 1/2/3.
0235   int    nTau, nY, nZ;
0236   double tauCoef[8], yCoef[8], zCoef[8], tauCoefSum[8], yCoefSum[8],
0237          zCoefSum[8];
0238
0239   // Calculate kinematical limits for 2 -> 1/2/3.
0240   bool limitTau(bool is2, bool is3);
0241   bool limitY();
0242   bool limitZ();
0243
0244   // Select kinematical variable between defined limits for 2 -> 1/2/3.
0245   void selectTau(int iTau, double tauVal, bool is2);
0246   void selectY(int iY, double yVal);
0247   void selectZ(int iZ, double zVal);
0248   bool select3Body();
0249
0250   // Solve equation system for better phase space coefficients in 2 -> 1/2/3.
0251   void solveSys( int n, int bin[8], double vec[8], double mat[8][8],
0252     double coef[8]);
0253
0254   // Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
0255   bool   useBW[6], useNarrowBW[6];
0256   int    idMass[6];
0257   double mPeak[6], sPeak[6], mWidth[6], mMin[6], mMax[6], mw[6], wmRat[6],
0258          mLower[6], mUpper[6], sLower[6], sUpper[6], fracFlatS[6],
0259          fracFlatM[6], fracInv[6], fracInv2[6], atanLower[6], atanUpper[6],
0260          intBW[6], intFlatS[6], intFlatM[6], intInv[6], intInv2[6];
0261
0262   // Setup mass selection for one resonance at a time. Split in two parts.
0263   void   setupMass1(int iM);
0264   void   setupMass2(int iM, double distToThresh);
0265
0266   // Do mass selection and find the associated weight.
0267   void   trialMass(int iM);
0268   double weightMass(int iM);
0269
0270   // Standard methods to find t range of a 2 -> 2 process
0271   // and to check whether a given t value is in that range.
0272   pair<double,double> tRange( double sIn, double s1In, double s2In,
0273     double s3In,  double s4In) {
0274     double lambda12 = pow2( sIn - s1In - s2In) - 4. * s1In * s2In;
0275     double lambda34 = pow2( sIn - s3In - s4In) - 4. * s3In * s4In;
0276     if (lambda12 < 0. || lambda34 < 0.) return make_pair( 0., 0.);
0277     double tLow = -0.5 * (sIn - (s1In + s2In + s3In + s4In) + (s1In - s2In)
0278       * (s3In - s4In) / sIn + sqrtpos(lambda12 *  lambda34) / sIn);
0279     double tUpp = ( (s3In - s1In) * (s4In - s2In) + (s1In + s4In - s2In - s3In)
0280       * (s1In * s4In - s2In * s3In) / sIn ) / tLow;
0281     return make_pair( tLow, tUpp); }
0282   bool tInRange( double tIn, double sIn, double s1In, double s2In,
0283     double s3In,  double s4In) {
0284     pair<double, double> tRng = tRange( sIn, s1In, s2In, s3In, s4In);
0285     return (tIn > tRng.first && tIn < tRng.second); }
0286
0287   // The error function erf(x) should normally be in your math library,
0288   // but if not uncomment this simple parametrization by Sergei Winitzki.
0289   //double erf(double x) { double x2 = x * x; double kx2 = 0.147 * x2;
0290   //  double tmp = sqrt(1. - exp(-x2 * (4./M_PI + kx2) / (1. + kx2)));
0291   //  return ((x >= 0.) ? tmp : -tmp); }
0292
0293 };
0294
0295 //==========================================================================
0296
0297 // A derived class with 2 -> 1 kinematics set up in tau, y.
0298
0299 class PhaseSpace2to1tauy : public PhaseSpace {
0300
0301 public:
0302
0303   // Constructor.
0304   PhaseSpace2to1tauy() {}
0305
0306   // Optimize subsequent kinematics selection.
0307   virtual bool setupSampling() {if (!setupMass()) return false;
0308     return setupSampling123(false, false);}
0309
0310   // Construct the trial kinematics.
0311   virtual bool trialKin(bool inEvent = true, bool = false) {wtBW = 1.;
0312     return trialKin123(false, false, inEvent);}
0313
0314   // Construct the final event kinematics.
0315   virtual bool finalKin();
0316
0317 private:
0318
0319   // Set up allowed mass range.
0320   bool setupMass();
0321
0322 };
0323
0324 //==========================================================================
0325
0326 // A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
0327
0328 class PhaseSpace2to2tauyz : public PhaseSpace {
0329
0330 public:
0331
0332   // Constructor.
0333   PhaseSpace2to2tauyz() {}
0334
0335   // Optimize subsequent kinematics selection.
0336   virtual bool setupSampling() {if (!setupMasses()) return false;
0337     return setupSampling123(true, false);}
0338
0339   // Construct the trial kinematics.
0340   virtual bool trialKin(bool inEvent = true, bool = false) {
0341     if (!trialMasses()) return false;
0342     return trialKin123(true, false, inEvent);}
0343
0344   // Construct the final event kinematics.
0345   virtual bool finalKin();
0346
0347   // Rescales the momenta of incoming and outgoing partons according to
0348   // new sHat.
0349   virtual void rescaleMomenta ( double sHatNew);
0350
0351   // Recalculates cross section with rescaled sHat.
0352   virtual void rescaleSigma ( double sHatNew);
0353
0354   // Calculate the weight for over-estimated cross section.
0355   virtual double weightGammaPDFApprox();
0356
0357 private:
0358
0359   // Set up for fixed or Breit-Wigner mass selection.
0360   bool setupMasses();
0361
0362   // Select fixed or Breit-Wigner-distributed masses.
0363   bool trialMasses();
0364
0365   // Pick off-shell initialization masses when on-shell not allowed.
0366   bool constrainedM3M4();
0367   bool constrainedM3();
0368   bool constrainedM4();
0369
0370 };
0371
0372 //==========================================================================
0373
0374 // A derived class with 2 -> 2 kinematics set up for elastic scattering.
0375
0376 class PhaseSpace2to2elastic : public PhaseSpace {
0377
0378 public:
0379
0380   // Constructor.
0381   PhaseSpace2to2elastic() : isOneExp(), useCoulomb(), s1(), s2(), alphaEM0(),
0382     lambda12S(), lambda12(), lambda34(), tempA(), tempB(), tempC(),
0383     tLow(), tUpp(), bSlope1(), bSlope2(), sigRef1(), sigRef2(),
0384     sigRef(), sigNorm1(), sigNorm2(), sigNorm3(), sigNormSum(), rel2() {}
0385
0386   // Construct the trial or final event kinematics.
0387   virtual bool setupSampling();
0388   virtual bool trialKin(bool inEvent = true, bool = false);
0389   virtual bool finalKin();
0390
0391   // Are beam particles resolved in partons or scatter directly?
0392   virtual bool isResolved() const {return false;}
0393
0394 private:
0395
0396   // Constants: could only be changed in the code itself.
0397   static const int    NTRY;
0398   static const double BNARROW, BWIDE, WIDEFRAC, TOFFSET;
0399
0400   // Kinematics properties specific to 2 -> 2 elastic.
0401   bool   isOneExp, useCoulomb;
0402   double s1, s2, alphaEM0, lambda12S, lambda12, lambda34, tempA, tempB, tempC,
0403          tLow, tUpp, bSlope1, bSlope2, sigRef1, sigRef2, sigRef,
0404          sigNorm1, sigNorm2, sigNorm3, sigNormSum, rel2;
0405
0406 };
0407
0408 //==========================================================================
0409
0410 // A derived class with 2 -> 2 kinematics set up for diffractive scattering.
0411
0412 class PhaseSpace2to2diffractive : public PhaseSpace {
0413
0414 public:
0415
0416   // Constructor.
0417   PhaseSpace2to2diffractive(bool isDiffAin = false, bool isDiffBin = false)
0418     : isDiffA(isDiffAin), isDiffB(isDiffBin), splitxit(), m3ElDiff(),
0419     m4ElDiff(), s1(), s2(), xiMin(), xiMax(), xiNow(), sigNow(), sigMax(),
0420     sigMaxNow(), lambda12(), lambda34(), bNow(), tempA(), tempB(), tempC(),
0421     tLow(), tUpp(), tWeight(), fWid1(), fWid2(), fWid3(), fWid4(), fbWid1(),
0422     fbWid2(), fbWid3(), fbWid4(), fbWid1234() {isSD = !isDiffA || !isDiffB;}
0423
0424   // Construct the trial or final event kinematics.
0425   virtual bool setupSampling();
0426   virtual bool trialKin(bool inEvent = true, bool = false);
0427   virtual bool finalKin();
0428
0429   // Are beam particles resolved in partons or scatter directly?
0430   virtual bool isResolved() const {return false;}
0431
0432 private:
0433
0434   // Constants: could only be changed in the code itself.
0435   static const int    NTRY;
0436   static const double BWID1, BWID2, BWID3, BWID4, FWID1SD, FWID2SD, FWID3SD,
0437                       FWID4SD, FWID1DD, FWID2DD, FWID3DD, FWID4DD, MAXFUDGESD,
0438                       MAXFUDGEDD, MAXFUDGET, DIFFMASSMARGIN, SPROTON;
0439
0440   // Initialization data.
0441   bool   isDiffA, isDiffB, isSD, splitxit;
0442
0443   // Pion mass, cached for efficiency.
0444   double mPi;
0445
0446   // Kinematics properties specific to 2 -> 2 diffraction.
0447   double m3ElDiff, m4ElDiff, s1, s2, xiMin, xiMax, xiNow, sigNow, sigMax,
0448          sigMaxNow, lambda12, lambda34, bNow, tempA, tempB, tempC,
0449          tLow, tUpp, tWeight, fWid1, fWid2, fWid3, fWid4, fbWid1, fbWid2,
0450          fbWid3, fbWid4, fbWid1234;
0451
0452 };
0453
0454 //==========================================================================
0455
0456 // A derived class with 2 -> 3 kinematics set up for central diffractive
0457 // scattering.
0458
0459 class PhaseSpace2to3diffractive : public PhaseSpace {
0460
0461 public:
0462
0463   // Constructor.
0464  PhaseSpace2to3diffractive() : PhaseSpace(), splitxit(), s1(), s2(), m5min(),
0465     s5min(), sigNow(), sigMax(), sigMaxNow(), xiMin(), xi1(), xi2(), fWid1(),
0466     fWid2(), fWid3(), fbWid1(), fbWid2(), fbWid3(), fbWid123() {}
0467
0468   // Construct the trial or final event kinematics.
0469   virtual bool setupSampling();
0470   virtual bool trialKin(bool inEvent = true, bool = false);
0471   virtual bool finalKin();
0472
0473   // Are beam particles resolved in partons or scatter directly?
0474   virtual bool isResolved() const {return false;}
0475
0476  private:
0477
0478   // Constants: could only be changed in the code itself.
0479   static const int    NTRY;
0480   static const double BWID1, BWID2, BWID3, BWID4, FWID1, FWID2, FWID3,
0481                       MAXFUDGECD, MAXFUDGET, DIFFMASSMARGIN;
0482
0483   // Initialization data.
0484   bool   splitxit;
0485
0486   // Local variables to calculate DPE kinematics.
0487   double s1, s2, m5min, s5min, sigNow, sigMax, sigMaxNow, xiMin, xi1, xi2,
0488          fWid1, fWid2, fWid3, fbWid1, fbWid2, fbWid3, fbWid123;
0489   Vec4   p1, p2, p3, p4, p5;
0490
0491 };
0492
0493 //==========================================================================
0494
0495 class PhaseSpace2to2nondiffractive : public PhaseSpace {
0496
0497 // A derived class for nondiffractive events. Hardly does anything, since
0498 // the real action is taken care of by the MultipartonInteractions class.
0499
0500 public:
0501
0502   // Constructor.
0503   PhaseSpace2to2nondiffractive() {}
0504
0505   // Construct the trial or final event kinematics.
0506   virtual bool setupSampling();
0507   virtual bool trialKin( bool , bool = false);
0508   virtual bool finalKin() { if (hasGamma) gammaKinPtr->finalize();
0509     return true;}
0510
0511 private:
0512
0513 };
0514
0515 //==========================================================================
0516
0517 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
0518 // tau, y, pT2_4, pT2_5, phi_4, phi_5 and y_3 (partial cylindrical symmetry).
0519
0520 class PhaseSpace2to3tauycyl : public PhaseSpace {
0521
0522 public:
0523
0524   // Constructor.
0525   PhaseSpace2to3tauycyl() {}
0526
0527   // Optimize subsequent kinematics selection.
0528   virtual bool setupSampling() {if (!setupMasses()) return false;
0529     setup3Body(); return setupSampling123(false, true);}
0530
0531   // Construct the trial kinematics.
0532   virtual bool trialKin(bool inEvent = true, bool = false) {
0533     if (!trialMasses()) return false;
0534     return trialKin123(false, true, inEvent);}
0535
0536   // Construct the final event kinematics.
0537   virtual bool finalKin();
0538
0539 private:
0540
0541   // Constants: could only be changed in the code itself.
0542   static const int    NITERNR;
0543
0544   // Set up for fixed or Breit-Wigner mass selection.
0545   bool setupMasses();
0546
0547   // Select fixed or Breit-Wigner-distributed masses.
0548   bool trialMasses();
0549
0550 };
0551
0552 //==========================================================================
0553
0554 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
0555 // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
0556 // Intended specifically for (essentially massless) 2 -> 3 QCD processes.
0557
0558 class PhaseSpace2to3yyycyl : public PhaseSpace {
0559
0560 public:
0561
0562   // Constructor.
0563   PhaseSpace2to3yyycyl() : pTHat3Min(), pTHat3Max(), pTHat5Min(), pTHat5Max(),
0564     RsepMin(), R2sepMin(), hasBaryonBeams(), pT3Min(), pT3Max(), pT5Min(),
0565     pT5Max(), y3Max(), y4Max(), y5Max(), pT3(), pT4(), pT5(), phi3(), phi4(),
0566     phi5(), y3(), y4(), y5(), dphi() {}
0567
0568   // Optimize subsequent kinematics selection.
0569   virtual bool setupSampling();
0570
0571   // Construct the trial kinematics.
0572   virtual bool trialKin(bool inEvent = true, bool = false);
0573
0574   // Construct the final event kinematics.
0575   virtual bool finalKin();
0576
0577 private:
0578
0579   // Phase space cuts specifically for 2 -> 3 QCD processes.
0580   double pTHat3Min, pTHat3Max, pTHat5Min, pTHat5Max, RsepMin, R2sepMin;
0581   bool   hasBaryonBeams;
0582
0583   // Event kinematics choices.
0584   double pT3Min, pT3Max, pT5Min, pT5Max, y3Max, y4Max, y5Max,
0585          pT3, pT4, pT5, phi3, phi4, phi5, y3, y4, y5, dphi;
0586   Vec4   pInSum;
0587
0588 };
0589
0590 //==========================================================================
0591
0592 // A derived class for Les Houches events.
0593
0594 class PhaseSpaceLHA : public PhaseSpace {
0595
0596 public:
0597
0598   // Constructor.
0599   PhaseSpaceLHA() : strategy(), stratAbs(), nProc(), idProcSave(0),
0600     xMaxAbsSum(), xSecSgnSum(), sigmaSgn() {}
0601
0602   // Find maximal cross section for comparison with internal processes.
0603   virtual bool setupSampling();
0604
0605   // Construct the next process, by interface to Les Houches class.
0606   virtual bool trialKin( bool , bool repeatSame = false);
0607
0608   // Set scale, alpha_s and alpha_em if not done.
0609   virtual bool finalKin() {sigmaProcessPtr->setScale(); return true;}
0610
0611   // For Les Houches with negative event weight needs
0612   virtual double sigmaSumSigned() const {return sigmaSgn;}
0613
0614 private:
0615
0616   // Constants.
0617   static const double CONVERTPB2MB;
0618
0619   // Local properties.
0620   int    strategy, stratAbs, nProc, idProcSave;
0621   double xMaxAbsSum, xSecSgnSum, sigmaSgn;
0622   vector<int>    idProc;
0623   vector<double> xMaxAbsProc;
0624
0625 };
0626
0627 //==========================================================================
0628
0629 // Rambo flat phase-space generator.
0630
0631 // This is an implementation of the Rambo phase-space generator as
0632 // presented in A New Monte Carlo Treatment Of Multiparticle Phase
0633 // Space At High-Energies, R. Kleiss, W.J. Stirling, S.D. Ellis, CPC40
0634 // (1986) 359.
0635
0636 class Rambo {
0637
0638  public:
0639
0640   // Deafult constructor.
0641   Rambo() { rndmPtr=nullptr; isInitPtr=false;}
0642
0643   // Initializing constructor.
0644   Rambo(Rndm* rndmPtrIn) { initPtr(rndmPtrIn); }
0645
0646   // Destructor.
0647   virtual ~Rambo() {}
0648
0649   // Initialize pointers.
0650   void initPtr(Rndm* rndmPtrIn) {rndmPtr = rndmPtrIn; isInitPtr = true;}
0651
0652   // Rambo phase space generator. Generates nOut uniformly distributed
0653   // massless 4-vectors with sqrt(s) = eCM. Output in pOut.
0654   double genPoint(double eCM,int nOut,vector<Vec4>& pOut);
0655
0656   // Massive generalisation, weights NOT 1 anymore - literal implementation
0657   // of original RAMBO paper by Ellis, Kleiss and Stirling. Number of particles
0658   // determined from size of mIn vector.
0659   double genPoint(double eCM,vector<double> mIn,vector<Vec4>& pOut);
0660
0661  private:
0662
0663   // Is initialized.
0664   bool isInitPtr;
0665
0666   // Pointer to the random number generator.
0667   Rndm*  rndmPtr;
0668
0669 };
0670
0671 //==========================================================================
0672
0673 } // end namespace Pythia8
0674
0675 #endif // Pythia8_PhaseSpace_H