StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PhaseSpace.h
1 // PhaseSpace.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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, SigmaProcess* sigmaProcessPtrIn);
51 
52  // Update the CM energy of the event.
53  void newECM(double eCMin) {eCM = eCMin; s = eCM * eCM;}
54 
55  // Store or replace Les Houches pointer.
56  void setLHAPtr(LHAupPtr lhaUpPtrIn) {lhaUpPtr = lhaUpPtrIn;}
57 
58  // A pure virtual method, wherein an optimization procedure
59  // is used to determine how phase space should be sampled.
60  virtual bool setupSampling() = 0;
61 
62  // A pure virtual method, wherein a trial event kinematics
63  // is to be selected in the derived class.
64  virtual bool trialKin(bool inEvent = true, bool repeatSame = false) = 0;
65 
66  // A pure virtual method, wherein the accepted event kinematics
67  // is to be constructed in the derived class.
68  virtual bool finalKin() = 0;
69 
70  // Allow for nonisotropic decays when ME's available.
71  void decayKinematics( Event& process);
72 
73  // Give back current or maximum cross section, or set latter.
74  double sigmaNow() const {return sigmaNw;}
75  double sigmaMax() const {return sigmaMx;}
76  double biasSelectionWeight() const {return biasWt;}
77  bool newSigmaMax() const {return newSigmaMx;}
78  void setSigmaMax(double sigmaMaxIn) {sigmaMx = sigmaMaxIn;}
79 
80  // For Les Houches with negative event weight needs
81  virtual double sigmaSumSigned() const {return sigmaMx;}
82 
83  // Give back constructed four-vectors and known masses.
84  Vec4 p(int i) const {return pH[i];}
85  double m(int i) const {return mH[i];}
86 
87  // Reset the four-momentum.
88  void setP(int i, Vec4 pNew) { pH[i] = pNew; }
89 
90  // Give back other event properties.
91  double ecm() const {return eCM;}
92  double x1() const {return x1H;}
93  double x2() const {return x2H;}
94  double sHat() const {return sH;}
95  double tHat() const {return tH;}
96  double uHat() const {return uH;}
97  double pTHat() const {return pTH;}
98  double thetaHat() const {return theta;}
99  double phiHat() const {return phi;}
100  double runBW3() const {return runBW3H;}
101  double runBW4() const {return runBW4H;}
102  double runBW5() const {return runBW5H;}
103 
104  // Inform whether beam particles are resolved in partons or scatter directly.
105  virtual bool isResolved() const {return true;}
106 
107  // Functions to rescale momenta and cross section for new sHat
108  // Currently implemented only for PhaseSpace2to2tauyz class.
109  virtual void rescaleSigma( double){}
110  virtual void rescaleMomenta( double){}
111 
112  // Calculate the weight for over-estimated cross section.
113  virtual double weightGammaPDFApprox(){ return 1.;}
114 
115  // Set the GammaKinematics pointer needed for soft photoproduction.
116  virtual void setGammaKinPtr( GammaKinematics* gammaKinPtrIn) {
117  gammaKinPtr = gammaKinPtrIn; }
118 
119 protected:
120 
121  // Constructor.
122  PhaseSpace() : sigmaProcessPtr(), lhaUpPtr(), gammaKinPtr(),
123  useBreitWigners(), doEnergySpread(), showSearch(), showViolation(),
124  increaseMaximum(), hasQ2Min(), gmZmodeGlobal(), mHatGlobalMin(),
125  mHatGlobalMax(), pTHatGlobalMin(), pTHatGlobalMax(), Q2GlobalMin(),
126  pTHatMinDiverge(), minWidthBreitWigners(), minWidthNarrowBW(), idA(),
127  idB(), idAgm(), idBgm(), mA(), mB(), eCM(), s(), sigmaMxGm(),
128  hasLeptonBeamA(), hasLeptonBeamB(), hasOneLeptonBeam(),
129  hasTwoLeptonBeams(), hasPointGammaA(), hasPointGammaB(),
130  hasOnePointParticle(), hasTwoPointParticles(), hasGamma(),
131  hasVMD(), newSigmaMx(),
132  canModifySigma(), canBiasSelection(), canBias2Sel(), gmZmode(),
133  bias2SelPow(), bias2SelRef(), wtBW(), sigmaNw(),
134  sigmaMx(), sigmaPos(), sigmaNeg(), biasWt(), mHatMin(), mHatMax(),
135  sHatMin(), sHatMax(), pTHatMin(), pTHatMax(), pT2HatMin(), pT2HatMax(),
136  x1H(), x2H(), m3(), m4(), m5(), s3(), s4(), s5(), mHat(), sH(), tH(), uH(),
137  pAbs(), p2Abs(), pTH(), theta(), phi(), betaZ(), mH(), idResA(), idResB(),
138  mResA(), mResB(), GammaResA(), GammaResB(), tauResA(), tauResB(),
139  widResA(), widResB(), sameResMass(), useMirrorWeight(), hasNegZ(),
140  hasPosZ(), tau(), y(), z(), tauMin(), tauMax(), yMax(), zMin(), zMax(),
141  ratio34(), unity34(), zNeg(), zPos(), wtTau(), wtY(), wtZ(), wt3Body(),
142  runBW3H(), runBW4H(), runBW5H(), intTau0(), intTau1(), intTau2(),
143  intTau3(), intTau4(), intTau5(), intTau6(), intY0(), intY12(), intY34(),
144  intY56(), mTchan1(), sTchan1(), mTchan2(), sTchan2(), frac3Flat(),
145  frac3Pow1(), frac3Pow2(), zNegMin(), zNegMax(), zPosMin(), zPosMax(),
146  nTau(), nY(), nZ(), tauCoef(), yCoef(), zCoef(), tauCoefSum(), yCoefSum(),
147  zCoefSum(), useBW(), useNarrowBW(), idMass(), mPeak(), sPeak(), mWidth(),
148  mMin(), mMax(), mw(), wmRat(), mLower(), mUpper(), sLower(), sUpper(),
149  fracFlatS(), fracFlatM(), fracInv(), fracInv2(), atanLower(), atanUpper(),
150  intBW(), intFlatS(), intFlatM(), intInv(), intInv2() {}
151 
152  // Constants: could only be changed in the code itself.
153  static const int NMAXTRY, NTRY3BODY;
154  static const double SAFETYMARGIN, TINY, EVENFRAC, SAMESIGMA, MRESMINABS,
155  WIDTHMARGIN, SAMEMASS, MASSMARGIN, EXTRABWWTMAX,
156  THRESHOLDSIZE, THRESHOLDSTEP, YRANGEMARGIN, LEPTONXMIN,
157  LEPTONXMAX, LEPTONXLOGMIN, LEPTONXLOGMAX, LEPTONTAUMIN,
158  SHATMINZ, PT2RATMINZ, WTCORRECTION[11];
159 
160  // Pointer to cross section.
161  SigmaProcess* sigmaProcessPtr;
162 
163  // Pointer to LHAup for generating external events.
164  LHAupPtr lhaUpPtr;
165 
166  // Pointer to object that samples photon kinematics from leptons.
167  GammaKinematics* gammaKinPtr;
168 
169  // Initialization data, normally only set once.
170  bool useBreitWigners, doEnergySpread, showSearch, showViolation,
171  increaseMaximum, hasQ2Min;
172  int gmZmodeGlobal;
173  double mHatGlobalMin, mHatGlobalMax, pTHatGlobalMin, pTHatGlobalMax,
174  Q2GlobalMin, pTHatMinDiverge, minWidthBreitWigners, minWidthNarrowBW;
175 
176  // Information on incoming beams.
177  int idA, idB, idAgm, idBgm;
178  double mA, mB, eCM, s, sigmaMxGm;
179  bool hasLeptonBeamA, hasLeptonBeamB, hasOneLeptonBeam, hasTwoLeptonBeams,
180  hasPointGammaA, hasPointGammaB, hasOnePointParticle,
181  hasTwoPointParticles, hasGamma, hasVMD;
182 
183  // Cross section information.
184  bool newSigmaMx, canModifySigma, canBiasSelection, canBias2Sel;
185  int gmZmode;
186  double bias2SelPow, bias2SelRef, wtBW, sigmaNw, sigmaMx, sigmaPos,
187  sigmaNeg, biasWt;
188 
189  // Process-specific kinematics properties, almost always available.
190  double mHatMin, mHatMax, sHatMin, sHatMax, pTHatMin, pTHatMax,
191  pT2HatMin, pT2HatMax;
192 
193  // Event-specific kinematics properties, almost always available.
194  double x1H, x2H, m3, m4, m5, s3, s4, s5, mHat, sH, tH, uH, pAbs, p2Abs,
195  pTH, theta, phi, betaZ;
196  Vec4 pH[12];
197  double mH[12];
198 
199  // Reselect decay products momenta isotropically in phase space.
200  void decayKinematicsStep( Event& process, int iRes);
201 
202  // Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
203 
204  // Determine how phase space should be sampled.
205  void setup3Body();
206  bool setupSampling123(bool is2, bool is3);
207 
208  // Select a trial kinematics phase space point.
209  bool trialKin123(bool is2, bool is3, bool inEvent = true);
210 
211  // Presence and properties of any s-channel resonances.
212  int idResA, idResB;
213  double mResA, mResB, GammaResA, GammaResB, tauResA, tauResB, widResA,
214  widResB;
215  bool sameResMass;
216 
217  // Kinematics properties specific to 2 -> 1/2/3.
218  bool useMirrorWeight, hasNegZ, hasPosZ;
219  double tau, y, z, tauMin, tauMax, yMax, zMin, zMax, ratio34, unity34,
220  zNeg, zPos, wtTau, wtY, wtZ, wt3Body, runBW3H, runBW4H, runBW5H,
221  intTau0, intTau1, intTau2, intTau3, intTau4, intTau5, intTau6,
222  intY0, intY12, intY34, intY56, mTchan1, sTchan1, mTchan2, sTchan2,
223  frac3Flat, frac3Pow1, frac3Pow2, zNegMin, zNegMax, zPosMin, zPosMax;
224  Vec4 p3cm, p4cm, p5cm;
225 
226  // Coefficients for optimized selection in 2 -> 1/2/3.
227  int nTau, nY, nZ;
228  double tauCoef[8], yCoef[8], zCoef[8], tauCoefSum[8], yCoefSum[8],
229  zCoefSum[8];
230 
231  // Calculate kinematical limits for 2 -> 1/2/3.
232  bool limitTau(bool is2, bool is3);
233  bool limitY();
234  bool limitZ();
235 
236  // Select kinematical variable between defined limits for 2 -> 1/2/3.
237  void selectTau(int iTau, double tauVal, bool is2);
238  void selectY(int iY, double yVal);
239  void selectZ(int iZ, double zVal);
240  bool select3Body();
241 
242  // Solve equation system for better phase space coefficients in 2 -> 1/2/3.
243  void solveSys( int n, int bin[8], double vec[8], double mat[8][8],
244  double coef[8]);
245 
246  // Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
247  bool useBW[6], useNarrowBW[6];
248  int idMass[6];
249  double mPeak[6], sPeak[6], mWidth[6], mMin[6], mMax[6], mw[6], wmRat[6],
250  mLower[6], mUpper[6], sLower[6], sUpper[6], fracFlatS[6],
251  fracFlatM[6], fracInv[6], fracInv2[6], atanLower[6], atanUpper[6],
252  intBW[6], intFlatS[6], intFlatM[6], intInv[6], intInv2[6];
253 
254  // Setup mass selection for one resonance at a time. Split in two parts.
255  void setupMass1(int iM);
256  void setupMass2(int iM, double distToThresh);
257 
258  // Do mass selection and find the associated weight.
259  void trialMass(int iM);
260  double weightMass(int iM);
261 
262  // Standard methods to find t range of a 2 -> 2 process
263  // and to check whether a given t value is in that range.
264  pair<double,double> tRange( double sIn, double s1In, double s2In,
265  double s3In, double s4In) {
266  double lambda12 = pow2( sIn - s1In - s2In) - 4. * s1In * s2In;
267  double lambda34 = pow2( sIn - s3In - s4In) - 4. * s3In * s4In;
268  if (lambda12 < 0. || lambda34 < 0.) return make_pair( 0., 0.);
269  double tLow = -0.5 * (sIn - (s1In + s2In + s3In + s4In) + (s1In - s2In)
270  * (s3In - s4In) / sIn + sqrtpos(lambda12 * lambda34) / sIn);
271  double tUpp = ( (s3In - s1In) * (s4In - s2In) + (s1In + s4In - s2In - s3In)
272  * (s1In * s4In - s2In * s3In) / sIn ) / tLow;
273  return make_pair( tLow, tUpp); }
274  bool tInRange( double tIn, double sIn, double s1In, double s2In,
275  double s3In, double s4In) {
276  pair<double, double> tRng = tRange( sIn, s1In, s2In, s3In, s4In);
277  return (tIn > tRng.first && tIn < tRng.second); }
278 
279  // The error function erf(x) should normally be in your math library,
280  // but if not uncomment this simple parametrization by Sergei Winitzki.
281  //double erf(double x) { double x2 = x * x; double kx2 = 0.147 * x2;
282  // double tmp = sqrt(1. - exp(-x2 * (4./M_PI + kx2) / (1. + kx2)));
283  // return ((x >= 0.) ? tmp : -tmp); }
284 
285 };
286 
287 //==========================================================================
288 
289 // A derived class with 2 -> 1 kinematics set up in tau, y.
290 
291 class PhaseSpace2to1tauy : public PhaseSpace {
292 
293 public:
294 
295  // Constructor.
296  PhaseSpace2to1tauy() {}
297 
298  // Optimize subsequent kinematics selection.
299  virtual bool setupSampling() {if (!setupMass()) return false;
300  return setupSampling123(false, false);}
301 
302  // Construct the trial kinematics.
303  virtual bool trialKin(bool inEvent = true, bool = false) {wtBW = 1.;
304  return trialKin123(false, false, inEvent);}
305 
306  // Construct the final event kinematics.
307  virtual bool finalKin();
308 
309 private:
310 
311  // Set up allowed mass range.
312  bool setupMass();
313 
314 };
315 
316 //==========================================================================
317 
318 // A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
319 
320 class PhaseSpace2to2tauyz : public PhaseSpace {
321 
322 public:
323 
324  // Constructor.
325  PhaseSpace2to2tauyz() {}
326 
327  // Optimize subsequent kinematics selection.
328  virtual bool setupSampling() {if (!setupMasses()) return false;
329  return setupSampling123(true, false);}
330 
331  // Construct the trial kinematics.
332  virtual bool trialKin(bool inEvent = true, bool = false) {
333  if (!trialMasses()) return false;
334  return trialKin123(true, false, inEvent);}
335 
336  // Construct the final event kinematics.
337  virtual bool finalKin();
338 
339  // Rescales the momenta of incoming and outgoing partons according to
340  // new sHat.
341  virtual void rescaleMomenta ( double sHatNew);
342 
343  // Recalculates cross section with rescaled sHat.
344  virtual void rescaleSigma ( double sHatNew);
345 
346  // Calculate the weight for over-estimated cross section.
347  virtual double weightGammaPDFApprox();
348 
349 private:
350 
351  // Set up for fixed or Breit-Wigner mass selection.
352  bool setupMasses();
353 
354  // Select fixed or Breit-Wigner-distributed masses.
355  bool trialMasses();
356 
357  // Pick off-shell initialization masses when on-shell not allowed.
358  bool constrainedM3M4();
359  bool constrainedM3();
360  bool constrainedM4();
361 
362 };
363 
364 //==========================================================================
365 
366 // A derived class with 2 -> 2 kinematics set up for elastic scattering.
367 
368 class PhaseSpace2to2elastic : public PhaseSpace {
369 
370 public:
371 
372  // Constructor.
373  PhaseSpace2to2elastic() : isOneExp(), useCoulomb(), s1(), s2(), alphaEM0(),
374  lambda12S(), lambda12(), lambda34(), tempA(), tempB(), tempC(),
375  tLow(), tUpp(), bSlope1(), bSlope2(), sigRef1(), sigRef2(),
376  sigRef(), sigNorm1(), sigNorm2(), sigNorm3(), sigNormSum(), rel2() {}
377 
378  // Construct the trial or final event kinematics.
379  virtual bool setupSampling();
380  virtual bool trialKin(bool inEvent = true, bool = false);
381  virtual bool finalKin();
382 
383  // Are beam particles resolved in partons or scatter directly?
384  virtual bool isResolved() const {return false;}
385 
386 private:
387 
388  // Constants: could only be changed in the code itself.
389  static const int NTRY;
390  static const double BNARROW, BWIDE, WIDEFRAC, TOFFSET;
391 
392  // Kinematics properties specific to 2 -> 2 elastic.
393  bool isOneExp, useCoulomb;
394  double s1, s2, alphaEM0, lambda12S, lambda12, lambda34, tempA, tempB, tempC,
395  tLow, tUpp, bSlope1, bSlope2, sigRef1, sigRef2, sigRef,
396  sigNorm1, sigNorm2, sigNorm3, sigNormSum, rel2;
397 
398 };
399 
400 //==========================================================================
401 
402 // A derived class with 2 -> 2 kinematics set up for diffractive scattering.
403 
404 class PhaseSpace2to2diffractive : public PhaseSpace {
405 
406 public:
407 
408  // Constructor.
409  PhaseSpace2to2diffractive(bool isDiffAin = false, bool isDiffBin = false)
410  : isDiffA(isDiffAin), isDiffB(isDiffBin), splitxit(), m3ElDiff(),
411  m4ElDiff(), s1(), s2(), xiMin(), xiMax(), xiNow(), sigNow(), sigMax(),
412  sigMaxNow(), lambda12(), lambda34(), bNow(), tempA(), tempB(), tempC(),
413  tLow(), tUpp(), tWeight(), fWid1(), fWid2(), fWid3(), fWid4(), fbWid1(),
414  fbWid2(), fbWid3(), fbWid4(), fbWid1234() {isSD = !isDiffA || !isDiffB;}
415 
416  // Construct the trial or final event kinematics.
417  virtual bool setupSampling();
418  virtual bool trialKin(bool inEvent = true, bool = false);
419  virtual bool finalKin();
420 
421  // Are beam particles resolved in partons or scatter directly?
422  virtual bool isResolved() const {return false;}
423 
424 private:
425 
426  // Constants: could only be changed in the code itself.
427  static const int NTRY;
428  static const double BWID1, BWID2, BWID3, BWID4, FWID1SD, FWID2SD, FWID3SD,
429  FWID4SD, FWID1DD, FWID2DD, FWID3DD, FWID4DD, MAXFUDGESD,
430  MAXFUDGEDD, MAXFUDGET, DIFFMASSMARGIN, SPROTON;
431 
432  // Initialization data.
433  bool isDiffA, isDiffB, isSD, splitxit;
434 
435  // Kinematics properties specific to 2 -> 2 diffraction.
436  double m3ElDiff, m4ElDiff, s1, s2, xiMin, xiMax, xiNow, sigNow, sigMax,
437  sigMaxNow, lambda12, lambda34, bNow, tempA, tempB, tempC,
438  tLow, tUpp, tWeight, fWid1, fWid2, fWid3, fWid4, fbWid1, fbWid2,
439  fbWid3, fbWid4, fbWid1234;
440 
441 };
442 
443 //==========================================================================
444 
445 // A derived class with 2 -> 3 kinematics set up for central diffractive
446 // scattering.
447 
448 class PhaseSpace2to3diffractive : public PhaseSpace {
449 
450 public:
451 
452  // Constructor.
453  PhaseSpace2to3diffractive() : PhaseSpace(), splitxit(), s1(), s2(), m5min(),
454  s5min(), sigNow(), sigMax(), sigMaxNow(), xiMin(), xi1(), xi2(), fWid1(),
455  fWid2(), fWid3(), fbWid1(), fbWid2(), fbWid3(), fbWid123() {}
456 
457  // Construct the trial or final event kinematics.
458  virtual bool setupSampling();
459  virtual bool trialKin(bool inEvent = true, bool = false);
460  virtual bool finalKin();
461 
462  // Are beam particles resolved in partons or scatter directly?
463  virtual bool isResolved() const {return false;}
464 
465  private:
466 
467  // Constants: could only be changed in the code itself.
468  static const int NTRY;
469  static const double BWID1, BWID2, BWID3, BWID4, FWID1, FWID2, FWID3,
470  MAXFUDGECD, MAXFUDGET, DIFFMASSMARGIN;
471 
472  // Initialization data.
473  bool splitxit;
474 
475  // Local variables to calculate DPE kinematics.
476  double s1, s2, m5min, s5min, sigNow, sigMax, sigMaxNow, xiMin, xi1, xi2,
477  fWid1, fWid2, fWid3, fbWid1, fbWid2, fbWid3, fbWid123;
478  Vec4 p1, p2, p3, p4, p5;
479 
480 };
481 
482 //==========================================================================
483 
484 class PhaseSpace2to2nondiffractive : public PhaseSpace {
485 
486 // A derived class for nondiffractive events. Hardly does anything, since
487 // the real action is taken care of by the MultipartonInteractions class.
488 
489 public:
490 
491  // Constructor.
492  PhaseSpace2to2nondiffractive() {}
493 
494  // Construct the trial or final event kinematics.
495  virtual bool setupSampling();
496  virtual bool trialKin( bool , bool = false);
497  virtual bool finalKin() { if (hasGamma) gammaKinPtr->finalize();
498  return true;}
499 
500 private:
501 
502 };
503 
504 //==========================================================================
505 
506 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
507 // tau, y, pT2_4, pT2_5, phi_4, phi_5 and y_3 (partial cylindrical symmetry).
508 
509 class PhaseSpace2to3tauycyl : public PhaseSpace {
510 
511 public:
512 
513  // Constructor.
514  PhaseSpace2to3tauycyl() {}
515 
516  // Optimize subsequent kinematics selection.
517  virtual bool setupSampling() {if (!setupMasses()) return false;
518  setup3Body(); return setupSampling123(false, true);}
519 
520  // Construct the trial kinematics.
521  virtual bool trialKin(bool inEvent = true, bool = false) {
522  if (!trialMasses()) return false;
523  return trialKin123(false, true, inEvent);}
524 
525  // Construct the final event kinematics.
526  virtual bool finalKin();
527 
528 private:
529 
530  // Constants: could only be changed in the code itself.
531  static const int NITERNR;
532 
533  // Set up for fixed or Breit-Wigner mass selection.
534  bool setupMasses();
535 
536  // Select fixed or Breit-Wigner-distributed masses.
537  bool trialMasses();
538 
539 };
540 
541 //==========================================================================
542 
543 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
544 // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
545 // Intended specifically for (essentially massless) 2 -> 3 QCD processes.
546 
547 class PhaseSpace2to3yyycyl : public PhaseSpace {
548 
549 public:
550 
551  // Constructor.
552  PhaseSpace2to3yyycyl() : pTHat3Min(), pTHat3Max(), pTHat5Min(), pTHat5Max(),
553  RsepMin(), R2sepMin(), hasBaryonBeams(), pT3Min(), pT3Max(), pT5Min(),
554  pT5Max(), y3Max(), y4Max(), y5Max(), pT3(), pT4(), pT5(), phi3(), phi4(),
555  phi5(), y3(), y4(), y5(), dphi() {}
556 
557  // Optimize subsequent kinematics selection.
558  virtual bool setupSampling();
559 
560  // Construct the trial kinematics.
561  virtual bool trialKin(bool inEvent = true, bool = false);
562 
563  // Construct the final event kinematics.
564  virtual bool finalKin();
565 
566 private:
567 
568  // Phase space cuts specifically for 2 -> 3 QCD processes.
569  double pTHat3Min, pTHat3Max, pTHat5Min, pTHat5Max, RsepMin, R2sepMin;
570  bool hasBaryonBeams;
571 
572  // Event kinematics choices.
573  double pT3Min, pT3Max, pT5Min, pT5Max, y3Max, y4Max, y5Max,
574  pT3, pT4, pT5, phi3, phi4, phi5, y3, y4, y5, dphi;
575  Vec4 pInSum;
576 
577 };
578 
579 //==========================================================================
580 
581 // A derived class for Les Houches events.
582 
583 class PhaseSpaceLHA : public PhaseSpace {
584 
585 public:
586 
587  // Constructor.
588  PhaseSpaceLHA() : strategy(), stratAbs(), nProc(), idProcSave(0),
589  xMaxAbsSum(), xSecSgnSum(), sigmaSgn() {}
590 
591  // Find maximal cross section for comparison with internal processes.
592  virtual bool setupSampling();
593 
594  // Construct the next process, by interface to Les Houches class.
595  virtual bool trialKin( bool , bool repeatSame = false);
596 
597  // Set scale, alpha_s and alpha_em if not done.
598  virtual bool finalKin() {sigmaProcessPtr->setScale(); return true;}
599 
600  // For Les Houches with negative event weight needs
601  virtual double sigmaSumSigned() const {return sigmaSgn;}
602 
603 private:
604 
605  // Constants.
606  static const double CONVERTPB2MB;
607 
608  // Local properties.
609  int strategy, stratAbs, nProc, idProcSave;
610  double xMaxAbsSum, xSecSgnSum, sigmaSgn;
611  vector<int> idProc;
612  vector<double> xMaxAbsProc;
613 
614 };
615 
616 //==========================================================================
617 
618 } // end namespace Pythia8
619 
620 #endif // Pythia8_PhaseSpace_H
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5