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