PYTHIA  8.316
PhaseSpace.h
1 // PhaseSpace.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2025 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Header file for phase space generators in kinematics selection.
7 // PhaseSpace: base class for phase space generators.
8 // Base class for derived classes> 2 ->1 , 2 -> 2, 2 -> 2 elastic/diffractive,
9 // 2 -> 2 nondiffractive, 2 -> 3, Les Houches.
10 
11 #ifndef Pythia8_PhaseSpace_H
12 #define Pythia8_PhaseSpace_H
13 
14 #include "Pythia8/Basics.h"
15 #include "Pythia8/BeamParticle.h"
16 #include "Pythia8/Info.h"
17 #include "Pythia8/LesHouches.h"
18 #include "Pythia8/MultipartonInteractions.h"
19 #include "Pythia8/ParticleData.h"
20 #include "Pythia8/PartonDistributions.h"
21 #include "Pythia8/PhysicsBase.h"
22 #include "Pythia8/PythiaStdlib.h"
23 #include "Pythia8/SigmaProcess.h"
24 #include "Pythia8/SigmaTotal.h"
25 #include "Pythia8/Settings.h"
26 #include "Pythia8/StandardModel.h"
27 #include "Pythia8/UserHooks.h"
28 #include "Pythia8/GammaKinematics.h"
29 
30 namespace Pythia8 {
31 
32 //==========================================================================
33 
34 // Forward reference to the UserHooks class.
35 class UserHooks;
36 
37 //==========================================================================
38 
39 // PhaseSpace is a base class for phase space generators
40 // used in the selection of hard-process kinematics.
41 
42 class PhaseSpace : public PhysicsBase {
43 
44 public:
45 
46  // Destructor.
47  virtual ~PhaseSpace() {}
48 
49  // Perform simple initialization and store pointers.
50  void init(bool isFirst, SigmaProcessPtr sigmaProcessPtrIn);
51 
52  // Switch to new beam particle identities; for similar hadrons only.
53  void updateBeamIDs() { idAold = idA; idBold = idB; idA = beamAPtr->id();
54  idB = beamBPtr->id(); mA = beamAPtr->m(); mB = beamBPtr->m();
55  sigmaProcessPtr->updateBeamIDs();}
56 
57  // Update the CM energy of the event.
58  void newECM(double eCMin) {eCM = eCMin; s = eCM * eCM;}
59 
60  // Store or replace Les Houches pointer.
61  void setLHAPtr(LHAupPtr lhaUpPtrIn) {lhaUpPtr = lhaUpPtrIn;}
62 
63  // A pure virtual method, wherein an optimization procedure
64  // is used to determine how phase space should be sampled.
65  virtual bool setupSampling() = 0;
66 
67  // A pure virtual method, wherein a trial event kinematics
68  // is to be selected in the derived class.
69  virtual bool trialKin(bool inEvent = true, bool repeatSame = false) = 0;
70 
71  // A pure virtual method, wherein the accepted event kinematics
72  // is to be constructed in the derived class.
73  virtual bool finalKin() = 0;
74 
75  // Allow for nonisotropic decays when ME's available.
76  void decayKinematics( Event& process);
77 
78  // Give back current or maximum cross section, or set latter.
79  double sigmaNow() const {return sigmaNw;}
80  double sigmaMax() const {return sigmaMx;}
81  double biasSelectionWeight() const {return biasWt;}
82  bool newSigmaMax() const {return newSigmaMx;}
83  void setSigmaMax(double sigmaMaxIn) {sigmaMx = sigmaMaxIn;}
84 
85  // New value for switched beam identity or energy (for SoftQCD processes).
86  double sigmaMaxSwitch() {sigmaNw = sigmaProcessPtr->sigmaHatWrap();
87  sigmaMx = sigmaNw; return sigmaMx;}
88 
89  // For Les Houches with negative event weight needs
90  virtual double sigmaSumSigned() const {return sigmaMx;}
91 
92  // Give back constructed four-vectors and known masses.
93  Vec4 p(int i) const {return pH[i];}
94  double m(int i) const {return mH[i];}
95 
96  // Reset the four-momentum.
97  void setP(int i, Vec4 pNew) { pH[i] = pNew; }
98 
99  // Give back other event properties.
100  double ecm() const {return eCM;}
101  double x1() const {return x1H;}
102  double x2() const {return x2H;}
103  double sHat() const {return sH;}
104  double tHat() const {return tH;}
105  double uHat() const {return uH;}
106  double pTHat() const {return pTH;}
107  double thetaHat() const {return theta;}
108  double phiHat() const {return phi;}
109  double runBW3() const {return runBW3H;}
110  double runBW4() const {return runBW4H;}
111  double runBW5() const {return runBW5H;}
112 
113  // Inform whether beam particles are resolved in partons or scatter directly.
114  virtual bool isResolved() const {return true;}
115 
116  // Functions to rescale momenta and cross section for new sHat
117  // Currently implemented only for PhaseSpace2to2tauyz class.
118  virtual void rescaleSigma( double){}
119  virtual void rescaleMomenta( double){}
120 
121  // Calculate the weight for over-estimated cross section.
122  virtual double weightGammaPDFApprox(){ return 1.;}
123 
124  // Set the GammaKinematics pointer needed for soft photoproduction.
125  virtual void setGammaKinPtr( GammaKinematics* gammaKinPtrIn) {
126  gammaKinPtr = gammaKinPtrIn; }
127 
128 protected:
129 
130  // Constructor.
132  useBreitWigners(), doEnergySpread(), showSearch(),
133  showViolation(), increaseMaximum(), hasQ2Min(), gmZmodeGlobal(),
134  mHatGlobalMin(), mHatGlobalMax(), pTHatGlobalMin(),
135  pTHatGlobalMax(), Q2GlobalMin(), Q2GlobalMax(), pTHatMinDiverge(),
136  minWidthBreitWigners(), minWidthNarrowBW(), idA(), idB(),
137  idAold(), idBold(), idAgm(), idBgm(), mA(), mB(), eCM(), s(),
138  sigmaMxGm(), hasLeptonBeamA(), hasLeptonBeamB(),
139  hasOneLeptonBeam(), hasTwoLeptonBeams(), hasPointGammaA(),
140  hasPointGammaB(), hasOnePointParticle(), hasTwoPointParticles(),
141  hasGamma(), hasVMD(), newSigmaMx(), canModifySigma(),
142  canBiasSelection(), canBias2Sel(), gmZmode(), bias2SelPow(),
143  bias2SelRef(), wtBW(), sigmaNw(), sigmaMx(), sigmaPos(),
144  sigmaNeg(), biasWt(), mHatMin(), mHatMax(), sHatMin(), sHatMax(),
145  pTHatMin(), pTHatMax(), pT2HatMin(), pT2HatMax(), x1H(), x2H(),
146  m3(), m4(), m5(), s3(), s4(), s5(), mHat(), sH(), tH(), uH(),
147  pAbs(), p2Abs(), pTH(), theta(), phi(), betaZ(), mH(), idResA(),
148  idResB(), mResA(), mResB(), GammaResA(), GammaResB(), tauResA(),
149  tauResB(), widResA(), widResB(), sameResMass(), useMirrorWeight(),
150  hasNegZ(), hasPosZ(), tau(), y(), z(), tauMin(), tauMax(), yMax(),
151  zMin(), zMax(), ratio34(), unity34(), zNeg(), zPos(), wtTau(),
152  wtY(), wtZ(), wt3Body(), runBW3H(), runBW4H(), runBW5H(),
153  intTau0(), intTau1(), intTau2(), intTau3(), intTau4(), intTau5(),
154  intTau6(), intY0(), intY12(), intY34(), intY56(), mTchan1(),
155  sTchan1(), mTchan2(), sTchan2(), frac3Flat(), frac3Pow1(),
156  frac3Pow2(), zNegMin(), zNegMax(), zPosMin(), zPosMax(), nTau(),
157  nY(), nZ(), tauCoef(), yCoef(), zCoef(), tauCoefSum(), yCoefSum(),
158  zCoefSum(), useBW(), useNarrowBW(), idMass(), mPeak(), sPeak(),
159  mWidth(), mMin(), mMax(), mw(), wmRat(), mLower(), mUpper(),
160  sLower(), sUpper(), fracFlatS(), fracFlatM(), fracInv(),
161  fracInv2(), atanLower(), atanUpper(), intBW(), intFlatS(),
162  intFlatM(), intInv(), intInv2(), doTopPair(), topThresholdModel(),
163  topThresholdWidth(), eThreshold(), m3Threshold(), m4Threshold() {}
164 
165  // Constants: could only be changed in the code itself.
166  static const int NMAXTRY, NTRY3BODY;
167  static const double SAFETYMARGIN, TINY, EVENFRAC, SAMESIGMA, MRESMINABS,
170  LEPTONXMAX, LEPTONXLOGMIN, LEPTONXLOGMAX, LEPTONTAUMIN,
172 
173  // Pointer to cross section.
174  SigmaProcessPtr sigmaProcessPtr;
175 
176  // Pointer to LHAup for generating external events.
177  LHAupPtr lhaUpPtr;
178 
179  // Pointer to object that samples photon kinematics from leptons.
181 
182  // Initialization data, normally only set once.
183  bool useBreitWigners, doEnergySpread, showSearch, showViolation,
184  increaseMaximum, hasQ2Min;
185  int gmZmodeGlobal;
186  double mHatGlobalMin, mHatGlobalMax, pTHatGlobalMin, pTHatGlobalMax,
187  Q2GlobalMin, Q2GlobalMax, pTHatMinDiverge, minWidthBreitWigners,
188  minWidthNarrowBW;
189 
190  // Information on incoming beams.
191  int idA, idB, idAold, idBold, idAgm, idBgm;
192  double mA, mB, eCM, s, sigmaMxGm;
193  bool hasLeptonBeamA, hasLeptonBeamB, hasOneLeptonBeam, hasTwoLeptonBeams,
194  hasPointGammaA, hasPointGammaB, hasOnePointParticle,
195  hasTwoPointParticles, hasGamma, hasVMD;
196 
197  // Cross section information.
198  bool newSigmaMx, canModifySigma, canBiasSelection, canBias2Sel;
199  int gmZmode;
200  double bias2SelPow, bias2SelRef, wtBW, sigmaNw, sigmaMx, sigmaPos,
201  sigmaNeg, biasWt;
202 
203  // Process-specific kinematics properties, almost always available.
204  double mHatMin, mHatMax, sHatMin, sHatMax, pTHatMin, pTHatMax,
205  pT2HatMin, pT2HatMax;
206 
207  // Event-specific kinematics properties, almost always available.
208  double x1H, x2H, m3, m4, m5, s3, s4, s5, mHat, sH, tH, uH, pAbs, p2Abs,
209  pTH, theta, phi, betaZ;
210  Vec4 pH[12];
211  double mH[12];
212 
213  // Reselect decay products momenta isotropically in phase space.
214  void decayKinematicsStep( Event& process, int iRes);
215 
216  // Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
217 
218  // Determine how phase space should be sampled.
219  void setup3Body();
220  bool setupSampling123(bool is2, bool is3);
221 
222  // Select a trial kinematics phase space point.
223  bool trialKin123(bool is2, bool is3, bool inEvent = true);
224 
225  // Presence and properties of any s-channel resonances.
226  int idResA, idResB;
227  double mResA, mResB, GammaResA, GammaResB, tauResA, tauResB, widResA,
228  widResB;
229  bool sameResMass;
230 
231  // Kinematics properties specific to 2 -> 1/2/3.
232  bool useMirrorWeight, hasNegZ, hasPosZ;
233  double tau, y, z, tauMin, tauMax, yMax, zMin, zMax, ratio34, unity34,
234  zNeg, zPos, wtTau, wtY, wtZ, wt3Body, runBW3H, runBW4H, runBW5H,
235  intTau0, intTau1, intTau2, intTau3, intTau4, intTau5, intTau6,
236  intY0, intY12, intY34, intY56, mTchan1, sTchan1, mTchan2, sTchan2,
237  frac3Flat, frac3Pow1, frac3Pow2, zNegMin, zNegMax, zPosMin, zPosMax;
238  Vec4 p3cm, p4cm, p5cm;
239 
240  // Coefficients for optimized selection in 2 -> 1/2/3.
241  int nTau, nY, nZ;
242  double tauCoef[8], yCoef[8], zCoef[8], tauCoefSum[8], yCoefSum[8],
243  zCoefSum[8];
244 
245  // Calculate kinematical limits for 2 -> 1/2/3.
246  bool limitTau(bool is2, bool is3);
247  bool limitY();
248  bool limitZ();
249 
250  // Select kinematical variable between defined limits for 2 -> 1/2/3.
251  void selectTau(int iTau, double tauVal, bool is2);
252  void selectY(int iY, double yVal);
253  void selectZ(int iZ, double zVal);
254  bool select3Body();
255 
256  // Solve equation system for better phase space coefficients in 2 -> 1/2/3.
257  void solveSys( int n, int bin[8], double vec[8], double mat[8][8],
258  double coef[8]);
259 
260  // Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
261  bool useBW[6], useNarrowBW[6];
262  int idMass[6];
263  double mPeak[6], sPeak[6], mWidth[6], mMin[6], mMax[6], mw[6], wmRat[6],
264  mLower[6], mUpper[6], sLower[6], sUpper[6], fracFlatS[6],
265  fracFlatM[6], fracInv[6], fracInv2[6], atanLower[6], atanUpper[6],
266  intBW[6], intFlatS[6], intFlatM[6], intInv[6], intInv2[6];
267 
268  // Properties specific to top threshold enhancement.
269  bool doTopPair;
270  int topThresholdModel;
271  double topThresholdWidth, eThreshold, m3Threshold, m4Threshold;
272 
273  // Setup mass selection for one resonance at a time. Split in two parts.
274  void setupMass1(int iM);
275  void setupMass2(int iM, double distToThresh);
276 
277  // Do mass selection and find the associated weight.
278  void trialMass(int iM);
279  double weightMass(int iM);
280 
281  // Standard methods to find t range of a 2 -> 2 process
282  // and to check whether a given t value is in that range.
283  pair<double,double> tRange( double sIn, double s1In, double s2In,
284  double s3In, double s4In) {
285  double lambda12 = pow2( sIn - s1In - s2In) - 4. * s1In * s2In;
286  double lambda34 = pow2( sIn - s3In - s4In) - 4. * s3In * s4In;
287  if (lambda12 < 0. || lambda34 < 0.) return make_pair( 0., 0.);
288  double tLow = -0.5 * (sIn - (s1In + s2In + s3In + s4In) + (s1In - s2In)
289  * (s3In - s4In) / sIn + sqrtpos(lambda12 * lambda34) / sIn);
290  double tUpp = ( (s3In - s1In) * (s4In - s2In) + (s1In + s4In - s2In - s3In)
291  * (s1In * s4In - s2In * s3In) / sIn ) / tLow;
292  return make_pair( tLow, tUpp); }
293  bool tInRange( double tIn, double sIn, double s1In, double s2In,
294  double s3In, double s4In) {
295  pair<double, double> tRng = tRange( sIn, s1In, s2In, s3In, s4In);
296  return (tIn > tRng.first && tIn < tRng.second); }
297 
298  // The error function erf(x) should normally be in your math library,
299  // but if not uncomment this simple parametrization by Sergei Winitzki.
300  //double erf(double x) { double x2 = x * x; double kx2 = 0.147 * x2;
301  // double tmp = sqrt(1. - exp(-x2 * (4./M_PI + kx2) / (1. + kx2)));
302  // return ((x >= 0.) ? tmp : -tmp); }
303 
304 };
305 
306 //==========================================================================
307 
308 // A derived class with 2 -> 1 kinematics set up in tau, y.
309 
311 
312 public:
313 
314  // Constructor.
316 
317  // Optimize subsequent kinematics selection.
318  virtual bool setupSampling() {if (!setupMass()) return false;
319  return setupSampling123(false, false);}
320 
321  // Construct the trial kinematics.
322  virtual bool trialKin(bool inEvent = true, bool = false) {wtBW = 1.;
323  return trialKin123(false, false, inEvent);}
324 
325  // Construct the final event kinematics.
326  virtual bool finalKin();
327 
328 private:
329 
330  // Set up allowed mass range.
331  bool setupMass();
332 
333 };
334 
335 //==========================================================================
336 
337 // A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
338 
340 
341 public:
342 
343  // Constructor.
345 
346  // Optimize subsequent kinematics selection.
347  virtual bool setupSampling() {if (!setupMasses()) return false;
348  return setupSampling123(true, false);}
349 
350  // Construct the trial kinematics.
351  virtual bool trialKin(bool inEvent = true, bool = false) {
352  if (!trialMasses()) return false;
353  return trialKin123(true, false, inEvent);}
354 
355  // Construct the final event kinematics.
356  virtual bool finalKin();
357 
358  // Rescales the momenta of incoming and outgoing partons according to
359  // new sHat.
360  virtual void rescaleMomenta ( double sHatNew);
361 
362  // Recalculates cross section with rescaled sHat.
363  virtual void rescaleSigma ( double sHatNew);
364 
365  // Calculate the weight for over-estimated cross section.
366  virtual double weightGammaPDFApprox();
367 
368 private:
369 
370  // Set up for fixed or Breit-Wigner mass selection.
371  bool setupMasses();
372 
373  // Select fixed or Breit-Wigner-distributed masses.
374  bool trialMasses();
375 
376  // Pick off-shell initialization masses when on-shell not allowed.
377  bool constrainedM3M4();
378  bool constrainedM3();
379  bool constrainedM4();
380 
381 };
382 
383 //==========================================================================
384 
385 // A derived class with 2 -> 2 kinematics set up for elastic scattering.
386 
388 
389 public:
390 
391  // Constructor.
392  PhaseSpace2to2elastic() : isOneExp(), useCoulomb(), s1(), s2(), alphaEM0(),
393  lambda12S(), lambda12(), lambda34(), tempA(), tempB(), tempC(),
394  tLow(), tUpp(), bSlope1(), bSlope2(), sigRef1(), sigRef2(),
395  sigRef(), sigNorm1(), sigNorm2(), sigNorm3(), sigNormSum(), rel2() {}
396 
397  // Construct the trial or final event kinematics.
398  virtual bool setupSampling();
399  virtual bool trialKin(bool inEvent = true, bool = false);
400  virtual bool finalKin();
401 
402  // Are beam particles resolved in partons or scatter directly?
403  virtual bool isResolved() const {return false;}
404 
405 private:
406 
407  // Constants: could only be changed in the code itself.
408  static const int NTRY;
409  static const double BNARROW, BWIDE, WIDEFRAC, TOFFSET;
410 
411  // Kinematics properties specific to 2 -> 2 elastic.
412  bool isOneExp, useCoulomb;
413  double s1, s2, alphaEM0, lambda12S, lambda12, lambda34, tempA, tempB, tempC,
414  tLow, tUpp, bSlope1, bSlope2, sigRef1, sigRef2, sigRef,
415  sigNorm1, sigNorm2, sigNorm3, sigNormSum, rel2;
416 
417 };
418 
419 //==========================================================================
420 
421 // A derived class with 2 -> 2 kinematics set up for diffractive scattering.
422 
424 
425 public:
426 
427  // Constructor.
428  PhaseSpace2to2diffractive(bool isDiffAin = false, bool isDiffBin = false)
429  : isDiffA(isDiffAin), isDiffB(isDiffBin), splitxit(), m3ElDiff(),
430  m4ElDiff(), s1(), s2(), xiMin(), xiMax(), xiNow(), sigNow(), sigMax(),
431  sigMaxNow(), lambda12(), lambda34(), bNow(), tempA(), tempB(), tempC(),
432  tLow(), tUpp(), tWeight(), fWid1(), fWid2(), fWid3(), fWid4(), fbWid1(),
433  fbWid2(), fbWid3(), fbWid4(), fbWid1234() {isSD = !isDiffA || !isDiffB;}
434 
435  // Construct the trial or final event kinematics.
436  virtual bool setupSampling();
437  virtual bool trialKin(bool inEvent = true, bool = false);
438  virtual bool finalKin();
439 
440  // Are beam particles resolved in partons or scatter directly?
441  virtual bool isResolved() const {return false;}
442 
443 private:
444 
445  // Constants: could only be changed in the code itself.
446  static const int NTRY;
447  static const double BWID1, BWID2, BWID3, BWID4, FWID1SD, FWID2SD, FWID3SD,
448  FWID4SD, FWID1DD, FWID2DD, FWID3DD, FWID4DD, MAXFUDGESD,
449  MAXFUDGEDD, MAXFUDGET, DIFFMASSMARGIN, SPROTON;
450 
451  // Initialization data.
452  bool isDiffA, isDiffB, isSD, splitxit;
453 
454  // Pion mass, cached for efficiency.
455  double mPi;
456 
457  // Kinematics properties specific to 2 -> 2 diffraction.
458  double m3ElDiff, m4ElDiff, s1, s2, xiMin, xiMax, xiNow, sigNow, sigMax,
459  sigMaxNow, lambda12, lambda34, bNow, tempA, tempB, tempC,
460  tLow, tUpp, tWeight, fWid1, fWid2, fWid3, fWid4, fbWid1, fbWid2,
461  fbWid3, fbWid4, fbWid1234;
462 
463 };
464 
465 //==========================================================================
466 
467 // A derived class with 2 -> 3 kinematics set up for central diffractive
468 // scattering.
469 
471 
472 public:
473 
474  // Constructor.
475  PhaseSpace2to3diffractive() : PhaseSpace(), splitxit(), s1(), s2(), m5min(),
476  s5min(), sigNow(), sigMax(), sigMaxNow(), xiMin(), xi1(), xi2(), fWid1(),
477  fWid2(), fWid3(), fbWid1(), fbWid2(), fbWid3(), fbWid123() {}
478 
479  // Construct the trial or final event kinematics.
480  virtual bool setupSampling();
481  virtual bool trialKin(bool inEvent = true, bool = false);
482  virtual bool finalKin();
483 
484  // Are beam particles resolved in partons or scatter directly?
485  virtual bool isResolved() const {return false;}
486 
487  private:
488 
489  // Constants: could only be changed in the code itself.
490  static const int NTRY;
491  static const double BWID1, BWID2, BWID3, BWID4, FWID1, FWID2, FWID3,
492  MAXFUDGECD, MAXFUDGET, DIFFMASSMARGIN;
493 
494  // Initialization data.
495  bool splitxit;
496 
497  // Local variables to calculate DPE kinematics.
498  double s1, s2, m5min, s5min, sigNow, sigMax, sigMaxNow, xiMin, xi1, xi2,
499  fWid1, fWid2, fWid3, fbWid1, fbWid2, fbWid3, fbWid123;
500  Vec4 p1, p2, p3, p4, p5;
501 
502 };
503 
504 //==========================================================================
505 
507 
508 // A derived class for nondiffractive events. Hardly does anything, since
509 // the real action is taken care of by the MultipartonInteractions class.
510 
511 public:
512 
513  // Constructor.
515 
516  // Construct the trial or final event kinematics.
517  virtual bool setupSampling();
518  virtual bool trialKin( bool , bool = false);
519  virtual bool finalKin() { if (hasGamma) gammaKinPtr->finalize();
520  return true;}
521 
522 private:
523 
524 };
525 
526 //==========================================================================
527 
528 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
529 // tau, y, pT2_4, pT2_5, phi_4, phi_5 and y_3 (partial cylindrical symmetry).
530 
532 
533 public:
534 
535  // Constructor.
537 
538  // Optimize subsequent kinematics selection.
539  virtual bool setupSampling() {if (!setupMasses()) return false;
540  setup3Body(); return setupSampling123(false, true);}
541 
542  // Construct the trial kinematics.
543  virtual bool trialKin(bool inEvent = true, bool = false) {
544  if (!trialMasses()) return false;
545  return trialKin123(false, true, inEvent);}
546 
547  // Construct the final event kinematics.
548  virtual bool finalKin();
549 
550 private:
551 
552  // Constants: could only be changed in the code itself.
553  static const int NITERNR;
554 
555  // Set up for fixed or Breit-Wigner mass selection.
556  bool setupMasses();
557 
558  // Select fixed or Breit-Wigner-distributed masses.
559  bool trialMasses();
560 
561 };
562 
563 //==========================================================================
564 
565 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
566 // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
567 // Intended specifically for (essentially massless) 2 -> 3 QCD processes.
568 
570 
571 public:
572 
573  // Constructor.
574  PhaseSpace2to3yyycyl() : pTHat3Min(), pTHat3Max(), pTHat5Min(), pTHat5Max(),
575  RsepMin(), R2sepMin(), hasBaryonBeams(), pT3Min(), pT3Max(), pT5Min(),
576  pT5Max(), y3Max(), y4Max(), y5Max(), pT3(), pT4(), pT5(), phi3(), phi4(),
577  phi5(), y3(), y4(), y5(), dphi() {}
578 
579  // Optimize subsequent kinematics selection.
580  virtual bool setupSampling();
581 
582  // Construct the trial kinematics.
583  virtual bool trialKin(bool inEvent = true, bool = false);
584 
585  // Construct the final event kinematics.
586  virtual bool finalKin();
587 
588 private:
589 
590  // Phase space cuts specifically for 2 -> 3 QCD processes.
591  double pTHat3Min, pTHat3Max, pTHat5Min, pTHat5Max, RsepMin, R2sepMin;
592  bool hasBaryonBeams;
593 
594  // Event kinematics choices.
595  double pT3Min, pT3Max, pT5Min, pT5Max, y3Max, y4Max, y5Max,
596  pT3, pT4, pT5, phi3, phi4, phi5, y3, y4, y5, dphi;
597  Vec4 pInSum;
598 
599 };
600 
601 //==========================================================================
602 
603 // A derived class for Les Houches events.
604 
605 class PhaseSpaceLHA : public PhaseSpace {
606 
607 public:
608 
609  // Constructor.
610  PhaseSpaceLHA() : strategy(), stratAbs(), nProc(), idProcSave(0),
611  xMaxAbsSum(), xSecSgnSum(), sigmaSgn() {}
612 
613  // Find maximal cross section for comparison with internal processes.
614  virtual bool setupSampling();
615 
616  // Construct the next process, by interface to Les Houches class.
617  virtual bool trialKin( bool , bool repeatSame = false);
618 
619  // Set scale, alpha_s and alpha_em if not done.
620  virtual bool finalKin() {sigmaProcessPtr->setScale(); return true;}
621 
622  // For Les Houches with negative event weight needs
623  virtual double sigmaSumSigned() const {return sigmaSgn;}
624 
625 private:
626 
627  // Constants.
628  static const double CONVERTPB2MB;
629 
630  // Local properties.
631  int strategy, stratAbs, nProc, idProcSave;
632  double xMaxAbsSum, xSecSgnSum, sigmaSgn;
633  vector<int> idProc;
634  vector<double> xMaxAbsProc;
635 
636 };
637 
638 //==========================================================================
639 
640 // Rambo flat phase-space generator.
641 
642 // This is an implementation of the Rambo phase-space generator as
643 // presented in A New Monte Carlo Treatment Of Multiparticle Phase
644 // Space At High-Energies, R. Kleiss, W.J. Stirling, S.D. Ellis, CPC40
645 // (1986) 359.
646 
647 class Rambo {
648 
649  public:
650 
651  // Deafult constructor.
652  Rambo() { rndmPtr=nullptr; isInitPtr=false;}
653 
654  // Initializing constructor.
655  Rambo(Rndm* rndmPtrIn) { initPtr(rndmPtrIn); }
656 
657  // Destructor.
658  virtual ~Rambo() {}
659 
660  // Initialize pointers.
661  void initPtr(Rndm* rndmPtrIn) {rndmPtr = rndmPtrIn; isInitPtr = true;}
662 
663  // Rambo phase space generator. Generates nOut uniformly distributed
664  // massless 4-vectors with sqrt(s) = eCM. Output in pOut.
665  double genPoint(double eCM, int nOut, vector<Vec4>& pOut);
666 
667  // Massive generalisation, weights NOT 1 anymore - literal implementation
668  // of original RAMBO paper by Ellis, Kleiss and Stirling. Number of particles
669  // determined from size of mIn vector.
670  double genPoint(double eCM, vector<double> mIn, vector<Vec4>& pOut);
671 
672  private:
673 
674  // Is initialized.
675  bool isInitPtr;
676 
677  // Pointer to the random number generator.
678  Rndm* rndmPtr;
679 
680 };
681 
682 //==========================================================================
683 
684 } // end namespace Pythia8
685 
686 #endif // Pythia8_PhaseSpace_H
void decayKinematicsStep(Event &process, int iRes)
Reselect decay products momenta isotropically in phase space.
Definition: PhaseSpace.cc:303
A derived class with 2 -> 2 kinematics set up for diffractive scattering.
Definition: PhaseSpace.h:423
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:173
virtual double weightGammaPDFApprox()
Calculate the weight for over-estimated cross section.
Definition: PhaseSpace.h:122
Vec4 p(int i) const
Give back constructed four-vectors and known masses.
Definition: PhaseSpace.h:93
static const double SAMEMASS
Special optimization treatment when two resonances at almost same mass.
Definition: PhaseSpace.h:167
void decayKinematics(Event &process)
Allow for nonisotropic decays when ME&#39;s available.
Definition: PhaseSpace.cc:245
static const double EXTRABWWTMAX
When using Breit-Wigners in 2 -> 2 raise maximum weight estimate.
Definition: PhaseSpace.h:167
static const double SAFETYMARGIN
Maximum cross section increase, just in case true maximum not found.
Definition: PhaseSpace.h:167
virtual bool trialKin(bool inEvent=true, bool=false)
Construct the trial kinematics.
Definition: PhaseSpace.h:351
PhaseSpace2to2tauyz()
Constructor.
Definition: PhaseSpace.h:344
A derived class for Les Houches events.
Definition: PhaseSpace.h:605
bool useBW[6]
Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
Definition: PhaseSpace.h:261
double ecm() const
Give back other event properties.
Definition: PhaseSpace.h:100
bool select3Body()
Definition: PhaseSpace.cc:1575
Definition: PhysicsBase.h:26
Definition: PhaseSpace.h:506
bool doTopPair
Properties specific to top threshold enhancement.
Definition: PhaseSpace.h:269
void setupMass1(int iM)
Setup mass selection for one resonance at a time. Split in two parts.
Definition: PhaseSpace.cc:1780
void trialMass(int iM)
Do mass selection and find the associated weight.
Definition: PhaseSpace.cc:1875
bool limitTau(bool is2, bool is3)
Calculate kinematical limits for 2 -> 1/2/3.
Definition: PhaseSpace.cc:1149
The Event class holds all info on the generated event.
Definition: Event.h:408
A derived class with 2 -> 1 kinematics set up in tau, y.
Definition: PhaseSpace.h:310
PhaseSpace2to3tauycyl()
Constructor.
Definition: PhaseSpace.h:536
static const int NTRY3BODY
Number of three-body trials in phase space optimization.
Definition: PhaseSpace.h:166
A derived class with 2 -> 2 kinematics set up for elastic scattering.
Definition: PhaseSpace.h:387
static const double SHATMINZ
Safety to avoid division with unreasonably small value for z selection.
Definition: PhaseSpace.h:167
static const double LEPTONXMIN
Definition: PhaseSpace.h:167
double sigmaMaxSwitch()
New value for switched beam identity or energy (for SoftQCD processes).
Definition: PhaseSpace.h:86
virtual bool setupSampling()
Optimize subsequent kinematics selection.
Definition: PhaseSpace.h:347
bool newSigmaMx
Cross section information.
Definition: PhaseSpace.h:198
virtual bool setupSampling()
Optimize subsequent kinematics selection.
Definition: PhaseSpace.h:318
virtual bool trialKin(bool inEvent=true, bool=false)
Construct the trial kinematics.
Definition: PhaseSpace.h:322
Rndm * rndmPtr
Pointer to the random number generator.
Definition: PhysicsBase.h:97
void setupMass2(int iM, double distToThresh)
Setup mass selection for one resonance at a time - part 2.
Definition: PhaseSpace.cc:1824
virtual double sigmaSumSigned() const
For Les Houches with negative event weight needs.
Definition: PhaseSpace.h:90
bool limitZ()
Find range of allowed z = cos(theta) values.
Definition: PhaseSpace.cc:1206
virtual bool isResolved() const
Inform whether beam particles are resolved in partons or scatter directly.
Definition: PhaseSpace.h:114
void newECM(double eCMin)
Update the CM energy of the event.
Definition: PhaseSpace.h:58
static const double THRESHOLDSIZE
Size of Breit-Wigner threshold region, for mass selection biasing.
Definition: PhaseSpace.h:167
virtual bool setupSampling()
Optimize subsequent kinematics selection.
Definition: PhaseSpace.h:539
virtual void setGammaKinPtr(GammaKinematics *gammaKinPtrIn)
Set the GammaKinematics pointer needed for soft photoproduction.
Definition: PhaseSpace.h:125
GammaKinematics * gammaKinPtr
Pointer to object that samples photon kinematics from leptons.
Definition: PhaseSpace.h:180
static const double MRESMINABS
Do not allow resonance to have mass below 2 m_e.
Definition: PhaseSpace.h:167
void init(bool isFirst, SigmaProcessPtr sigmaProcessPtrIn)
Perform simple initialization and store pointers.
Definition: PhaseSpace.cc:88
virtual bool finalKin()
Definition: PhaseSpace.h:519
PhaseSpace2to3yyycyl()
Constructor.
Definition: PhaseSpace.h:574
virtual bool isResolved() const
Are beam particles resolved in partons or scatter directly?
Definition: PhaseSpace.h:403
static const double SAMESIGMA
Two cross sections with a small relative error are assumed same.
Definition: PhaseSpace.h:167
static const double WTCORRECTION[11]
Definition: PhaseSpace.h:167
Rambo(Rndm *rndmPtrIn)
Initializing constructor.
Definition: PhaseSpace.h:655
PhaseSpaceLHA()
Constructor.
Definition: PhaseSpace.h:610
Definition: Basics.h:388
void solveSys(int n, int bin[8], double vec[8], double mat[8][8], double coef[8])
Solve equation system for better phase space coefficients in 2 -> 1/2/3.
Definition: PhaseSpace.cc:1703
virtual bool finalKin()
Set scale, alpha_s and alpha_em if not done.
Definition: PhaseSpace.h:620
Definition: PhaseSpace.h:470
int idResA
Presence and properties of any s-channel resonances.
Definition: PhaseSpace.h:226
bool useMirrorWeight
Kinematics properties specific to 2 -> 1/2/3.
Definition: PhaseSpace.h:232
int idA
Information on incoming beams.
Definition: PhaseSpace.h:191
virtual bool trialKin(bool inEvent=true, bool repeatSame=false)=0
virtual ~Rambo()
Destructor.
Definition: PhaseSpace.h:658
int nTau
Coefficients for optimized selection in 2 -> 1/2/3.
Definition: PhaseSpace.h:241
virtual ~PhaseSpace()
Destructor.
Definition: PhaseSpace.h:47
double weightMass(int iM)
Definition: PhaseSpace.cc:1919
PhaseSpace2to2diffractive(bool isDiffAin=false, bool isDiffBin=false)
Constructor.
Definition: PhaseSpace.h:428
static const int NMAXTRY
Constants: could only be changed in the code itself.
Definition: PhaseSpace.h:166
bool setupSampling123(bool is2, bool is3)
Determine how phase space should be sampled.
Definition: PhaseSpace.cc:506
PhaseSpace()
Constructor.
Definition: PhaseSpace.h:131
static const double EVENFRAC
Fraction of total weight that is shared evenly between all shapes.
Definition: PhaseSpace.h:167
PhaseSpace2to3diffractive()
Constructor.
Definition: PhaseSpace.h:475
Rambo flat phase-space generator.
Definition: PhaseSpace.h:647
void selectTau(int iTau, double tauVal, bool is2)
Select kinematical variable between defined limits for 2 -> 1/2/3.
Definition: PhaseSpace.cc:1256
void setP(int i, Vec4 pNew)
Reset the four-momentum.
Definition: PhaseSpace.h:97
int id() const
Member functions for output.
Definition: BeamParticle.h:195
LHAupPtr lhaUpPtr
Pointer to LHAup for generating external events.
Definition: PhaseSpace.h:177
PhaseSpace2to2elastic()
Constructor.
Definition: PhaseSpace.h:392
void selectY(int iY, double yVal)
Select y according to a choice of shapes.
Definition: PhaseSpace.cc:1373
PhaseSpace2to2nondiffractive()
Constructor.
Definition: PhaseSpace.h:514
static const double YRANGEMARGIN
Minimal rapidity range for allowed open range (in 2 -> 3).
Definition: PhaseSpace.h:167
static const double TINY
Small number to avoid division by zero.
Definition: PhaseSpace.h:167
Class to sample the virtuality and transverse momentum of emitted photons.
Definition: GammaKinematics.h:23
virtual double sigmaSumSigned() const
For Les Houches with negative event weight needs.
Definition: PhaseSpace.h:623
pair< double, double > tRange(double sIn, double s1In, double s2In, double s3In, double s4In)
Definition: PhaseSpace.h:283
virtual void rescaleSigma(double)
Definition: PhaseSpace.h:118
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
static const double PT2RATMINZ
Regularization for small pT2min in z = cos(theta) selection.
Definition: PhaseSpace.h:167
virtual bool isResolved() const
Are beam particles resolved in partons or scatter directly?
Definition: PhaseSpace.h:441
void setup3Body()
Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
Definition: PhaseSpace.cc:482
Definition: PhaseSpace.h:531
bool limitY()
Find range of allowed y values.
Definition: PhaseSpace.cc:1183
virtual bool isResolved() const
Are beam particles resolved in partons or scatter directly?
Definition: PhaseSpace.h:485
Definition: PhaseSpace.h:42
void setLHAPtr(LHAupPtr lhaUpPtrIn)
Store or replace Les Houches pointer.
Definition: PhaseSpace.h:61
bool useBreitWigners
Initialization data, normally only set once.
Definition: PhaseSpace.h:183
static const double THRESHOLDSTEP
Step size in optimal-mass search, for mass selection biasing.
Definition: PhaseSpace.h:167
double mHatMin
Process-specific kinematics properties, almost always available.
Definition: PhaseSpace.h:204
virtual bool trialKin(bool inEvent=true, bool=false)
Construct the trial kinematics.
Definition: PhaseSpace.h:543
double x1H
Event-specific kinematics properties, almost always available.
Definition: PhaseSpace.h:208
static const double MASSMARGIN
Minimum phase space left when kinematics constraints are combined.
Definition: PhaseSpace.h:167
void updateBeamIDs()
Switch to new beam particle identities; for similar hadrons only.
Definition: PhaseSpace.h:53
bool trialKin123(bool is2, bool is3, bool inEvent=true)
Select a trial kinematics phase space point.
Definition: PhaseSpace.cc:1000
static const double WIDTHMARGIN
Do not include resonances peaked too far outside allowed mass region.
Definition: PhaseSpace.h:167
double sigmaNow() const
Give back current or maximum cross section, or set latter.
Definition: PhaseSpace.h:79
virtual bool finalKin()=0
Rambo()
Deafult constructor.
Definition: PhaseSpace.h:652
SigmaProcessPtr sigmaProcessPtr
Pointer to cross section.
Definition: PhaseSpace.h:174
Definition: Basics.h:32
PhaseSpace2to1tauy()
Constructor.
Definition: PhaseSpace.h:315
void selectZ(int iZ, double zVal)
Definition: PhaseSpace.cc:1453
double sqrtpos(const double &x)
Avoid problem with negative square root argument (from roundoff).
Definition: PythiaStdlib.h:182
void initPtr(Rndm *rndmPtrIn)
Initialize pointers.
Definition: PhaseSpace.h:661
A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
Definition: PhaseSpace.h:339
virtual bool setupSampling()=0
bool finalize()
Save the accepted values for further use.
Definition: GammaKinematics.cc:319
Definition: PhaseSpace.h:569