Line data Source code
1 : // PhaseSpace.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, 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/PythiaStdlib.h"
22 : #include "Pythia8/SigmaProcess.h"
23 : #include "Pythia8/SigmaTotal.h"
24 : #include "Pythia8/Settings.h"
25 : #include "Pythia8/StandardModel.h"
26 : #include "Pythia8/UserHooks.h"
27 :
28 : namespace Pythia8 {
29 :
30 : //==========================================================================
31 :
32 : // Forward reference to the UserHooks class.
33 : class UserHooks;
34 :
35 : //==========================================================================
36 :
37 : // PhaseSpace is a base class for phase space generators
38 : // used in the selection of hard-process kinematics.
39 :
40 : class PhaseSpace {
41 :
42 : public:
43 :
44 : // Destructor.
45 0 : virtual ~PhaseSpace() {}
46 :
47 : // Perform simple initialization and store pointers.
48 : void init(bool isFirst, SigmaProcess* sigmaProcessPtrIn,
49 : Info* infoPtrIn, Settings* settingsPtrIn, ParticleData* particleDataPtrIn,
50 : Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
51 : Couplings* couplingsPtrIn, SigmaTotal* sigmaTotPtrIn,
52 : UserHooks* userHooksPtrIn);
53 :
54 : // Update the CM energy of the event.
55 0 : void newECM(double eCMin) {eCM = eCMin; s = eCM * eCM;}
56 :
57 : // Store or replace Les Houches pointer.
58 0 : void setLHAPtr(LHAup* lhaUpPtrIn) {lhaUpPtr = lhaUpPtrIn;}
59 :
60 : // A pure virtual method, wherein an optimization procedure
61 : // is used to determine how phase space should be sampled.
62 : virtual bool setupSampling() = 0;
63 :
64 : // A pure virtual method, wherein a trial event kinematics
65 : // is to be selected in the derived class.
66 : virtual bool trialKin(bool inEvent = true, bool repeatSame = false) = 0;
67 :
68 : // A pure virtual method, wherein the accepted event kinematics
69 : // is to be constructed in the derived class.
70 : virtual bool finalKin() = 0;
71 :
72 : // Allow for nonisotropic decays when ME's available.
73 : void decayKinematics( Event& process);
74 :
75 : // Give back current or maximum cross section, or set latter.
76 0 : double sigmaNow() const {return sigmaNw;}
77 0 : double sigmaMax() const {return sigmaMx;}
78 0 : double biasSelectionWeight() const {return biasWt;}
79 0 : bool newSigmaMax() const {return newSigmaMx;}
80 0 : void setSigmaMax(double sigmaMaxIn) {sigmaMx = sigmaMaxIn;}
81 :
82 : // For Les Houches with negative event weight needs
83 0 : virtual double sigmaSumSigned() const {return sigmaMx;}
84 :
85 : // Give back constructed four-vectors and known masses.
86 0 : Vec4 p(int i) const {return pH[i];}
87 0 : double m(int i) const {return mH[i];}
88 :
89 : // Give back other event properties.
90 0 : double ecm() const {return eCM;}
91 0 : double x1() const {return x1H;}
92 0 : double x2() const {return x2H;}
93 0 : double sHat() const {return sH;}
94 0 : double tHat() const {return tH;}
95 0 : double uHat() const {return uH;}
96 0 : double pTHat() const {return pTH;}
97 0 : double thetaHat() const {return theta;}
98 0 : double phiHat() const {return phi;}
99 : double runBW3() const {return runBW3H;}
100 : double runBW4() const {return runBW4H;}
101 : double runBW5() const {return runBW5H;}
102 :
103 : // Inform whether beam particles are resolved in partons or scatter directly.
104 0 : virtual bool isResolved() const {return true;}
105 :
106 : protected:
107 :
108 : // Constructor.
109 0 : PhaseSpace() {}
110 :
111 : // Constants: could only be changed in the code itself.
112 : static const int NMAXTRY, NTRY3BODY;
113 : static const double SAFETYMARGIN, TINY, EVENFRAC, SAMESIGMA, WIDTHMARGIN,
114 : SAMEMASS, MASSMARGIN, EXTRABWWTMAX, THRESHOLDSIZE,
115 : THRESHOLDSTEP, YRANGEMARGIN, LEPTONXMIN, LEPTONXMAX,
116 : LEPTONXLOGMIN, LEPTONXLOGMAX, LEPTONTAUMIN,
117 : SHATMINZ, PT2RATMINZ, WTCORRECTION[11];
118 :
119 : // MBR constants: form factor appoximation with two exponents.
120 : static const double FFA1, FFA2,FFB1, FFB2;
121 :
122 : // Pointer to cross section.
123 : SigmaProcess* sigmaProcessPtr;
124 :
125 : // Pointer to various information on the generation.
126 : Info* infoPtr;
127 :
128 : // Pointer to the settings database.
129 : Settings* settingsPtr;
130 :
131 : // Pointer to the particle data table.
132 : ParticleData* particleDataPtr;
133 :
134 : // Pointer to the random number generator.
135 : Rndm* rndmPtr;
136 :
137 : // Pointers to incoming beams.
138 : BeamParticle* beamAPtr;
139 : BeamParticle* beamBPtr;
140 :
141 : // Pointer to Standard Model couplings.
142 : Couplings* couplingsPtr;
143 :
144 : // Pointer to the total/elastic/diffractive cross section object.
145 : SigmaTotal* sigmaTotPtr;
146 :
147 : // Pointer to userHooks object for user interaction with program.
148 : UserHooks* userHooksPtr;
149 :
150 : // Pointer to LHAup for generating external events.
151 : LHAup* lhaUpPtr;
152 :
153 : // Initialization data, normally only set once.
154 : bool useBreitWigners, doEnergySpread, showSearch, showViolation,
155 : increaseMaximum;
156 : int gmZmodeGlobal;
157 : double mHatGlobalMin, mHatGlobalMax, pTHatGlobalMin, pTHatGlobalMax,
158 : pTHatMinDiverge, minWidthBreitWigners;
159 :
160 : // Information on incoming beams.
161 : int idA, idB;
162 : double mA, mB, eCM, s;
163 : bool hasLeptonBeamA, hasLeptonBeamB, hasOneLeptonBeam,
164 : hasTwoLeptonBeams, hasOnePointLepton, hasTwoPointLeptons;
165 :
166 : // Cross section information.
167 : bool newSigmaMx, canModifySigma, canBiasSelection, canBias2Sel;
168 : int gmZmode;
169 : double bias2SelPow, bias2SelRef, wtBW, sigmaNw, sigmaMx, sigmaPos,
170 : sigmaNeg, biasWt;
171 :
172 : // Process-specific kinematics properties, almost always available.
173 : double mHatMin, mHatMax, sHatMin, sHatMax, pTHatMin, pTHatMax,
174 : pT2HatMin, pT2HatMax;
175 :
176 : // Event-specific kinematics properties, almost always available.
177 : double x1H, x2H, m3, m4, m5, s3, s4, s5, mHat, sH, tH, uH, pAbs, p2Abs,
178 : pTH, theta, phi, betaZ;
179 : Vec4 pH[12];
180 : double mH[12];
181 :
182 : // Reselect decay products momenta isotropically in phase space.
183 : void decayKinematicsStep( Event& process, int iRes);
184 :
185 : // Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
186 :
187 : // Determine how phase space should be sampled.
188 : void setup3Body();
189 : bool setupSampling123(bool is2, bool is3, ostream& os = cout);
190 :
191 : // Select a trial kinematics phase space point.
192 : bool trialKin123(bool is2, bool is3, bool inEvent = true,
193 : ostream& os = cout);
194 :
195 : // Presence and properties of any s-channel resonances.
196 : int idResA, idResB;
197 : double mResA, mResB, GammaResA, GammaResB, tauResA, tauResB, widResA,
198 : widResB;
199 : bool sameResMass;
200 :
201 : // Kinematics properties specific to 2 -> 1/2/3.
202 : bool useMirrorWeight;
203 : double tau, y, z, tauMin, tauMax, yMax, zMin, zMax, ratio34, unity34,
204 : zNeg, zPos, wtTau, wtY, wtZ, wt3Body, runBW3H, runBW4H, runBW5H,
205 : intTau0, intTau1, intTau2, intTau3, intTau4, intTau5, intTau6,
206 : intY0, intY12, intY34, intY56, mTchan1, sTchan1, mTchan2, sTchan2,
207 : frac3Flat, frac3Pow1, frac3Pow2;
208 : Vec4 p3cm, p4cm, p5cm;
209 :
210 : // Coefficients for optimized selection in 2 -> 1/2/3.
211 : int nTau, nY, nZ;
212 : double tauCoef[8], yCoef[8], zCoef[8], tauCoefSum[8], yCoefSum[8],
213 : zCoefSum[8];
214 :
215 : // Calculate kinematical limits for 2 -> 1/2/3.
216 : bool limitTau(bool is2, bool is3);
217 : bool limitY();
218 : bool limitZ();
219 :
220 : // Select kinematical variable between defined limits for 2 -> 1/2/3.
221 : void selectTau(int iTau, double tauVal, bool is2);
222 : void selectY(int iY, double yVal);
223 : void selectZ(int iZ, double zVal);
224 : bool select3Body();
225 :
226 : // Solve equation system for better phase space coefficients in 2 -> 1/2/3.
227 : void solveSys( int n, int bin[8], double vec[8], double mat[8][8],
228 : double coef[8], ostream& os = cout);
229 :
230 : // Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
231 : bool useBW[6];
232 : int idMass[6];
233 : double mPeak[6], sPeak[6], mWidth[6], mMin[6], mMax[6], mw[6], wmRat[6],
234 : mLower[6], mUpper[6], sLower[6], sUpper[6], fracFlat[6], fracInv[6],
235 : fracInv2[6], atanLower[6], atanUpper[6], intBW[6], intFlat[6],
236 : intInv[6], intInv2[6];
237 :
238 : // Setup mass selection for one resonance at a time. Split in two parts.
239 : void setupMass1(int iM);
240 : void setupMass2(int iM, double distToThresh);
241 :
242 : // Do mass selection and find the associated weight.
243 : void trialMass(int iM);
244 : double weightMass(int iM);
245 :
246 : // The error function erf(x) should normally be in your math library,
247 : // but if not uncomment this simple parametrization by Sergei Winitzki.
248 : //double erf(double x) { double x2 = x * x; double kx2 = 0.147 * x2;
249 : // double tmp = sqrt(1. - exp(-x2 * (4./M_PI + kx2) / (1. + kx2)));
250 : // return ((x >= 0.) ? tmp : -tmp); }
251 :
252 : };
253 :
254 : //==========================================================================
255 :
256 : // A derived class with 2 -> 1 kinematics set up in tau, y.
257 :
258 0 : class PhaseSpace2to1tauy : public PhaseSpace {
259 :
260 : public:
261 :
262 : // Constructor.
263 0 : PhaseSpace2to1tauy() {}
264 :
265 : // Optimize subsequent kinematics selection.
266 0 : virtual bool setupSampling() {if (!setupMass()) return false;
267 0 : return setupSampling123(false, false);}
268 :
269 : // Construct the trial kinematics.
270 0 : virtual bool trialKin(bool inEvent = true, bool = false) {wtBW = 1.;
271 0 : return trialKin123(false, false, inEvent);}
272 :
273 : // Construct the final event kinematics.
274 : virtual bool finalKin();
275 :
276 : private:
277 :
278 : // Set up allowed mass range.
279 : bool setupMass();
280 :
281 : };
282 :
283 : //==========================================================================
284 :
285 : // A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
286 :
287 0 : class PhaseSpace2to2tauyz : public PhaseSpace {
288 :
289 : public:
290 :
291 : // Constructor.
292 0 : PhaseSpace2to2tauyz() {}
293 :
294 : // Optimize subsequent kinematics selection.
295 0 : virtual bool setupSampling() {if (!setupMasses()) return false;
296 0 : return setupSampling123(true, false);}
297 :
298 : // Construct the trial kinematics.
299 : virtual bool trialKin(bool inEvent = true, bool = false) {
300 0 : if (!trialMasses()) return false;
301 0 : return trialKin123(true, false, inEvent);}
302 :
303 : // Construct the final event kinematics.
304 : virtual bool finalKin();
305 :
306 : private:
307 :
308 : // Set up for fixed or Breit-Wigner mass selection.
309 : bool setupMasses();
310 :
311 : // Select fixed or Breit-Wigner-distributed masses.
312 : bool trialMasses();
313 :
314 : // Pick off-shell initialization masses when on-shell not allowed.
315 : bool constrainedM3M4();
316 : bool constrainedM3();
317 : bool constrainedM4();
318 :
319 : };
320 :
321 : //==========================================================================
322 :
323 : // A derived class with 2 -> 2 kinematics set up for elastic scattering.
324 :
325 0 : class PhaseSpace2to2elastic : public PhaseSpace {
326 :
327 : public:
328 :
329 : // Constructor.
330 0 : PhaseSpace2to2elastic() {}
331 :
332 : // Construct the trial or final event kinematics.
333 : virtual bool setupSampling();
334 : virtual bool trialKin(bool inEvent = true, bool = false);
335 : virtual bool finalKin();
336 :
337 : // Are beam particles resolved in partons or scatter directly?
338 0 : virtual bool isResolved() const {return false;}
339 :
340 : private:
341 :
342 : // Constants: could only be changed in the code itself.
343 : static const double EXPMAX, CONVERTEL;
344 :
345 : // Kinematics properties specific to 2 -> 2 elastic.
346 : bool useCoulomb;
347 : double s1, s2, bSlope, lambda12S, tLow, tUpp, tAux, sigmaTot, rho,
348 : lambda, tAbsMin, phaseCst, alphaEM0, sigmaNuc, sigmaCou, signCou;
349 :
350 : // Calculation of alphaElectromagnetic.
351 : AlphaEM alphaEM;
352 :
353 : };
354 :
355 : //==========================================================================
356 :
357 : // A derived class with 2 -> 2 kinematics set up for diffractive scattering.
358 :
359 0 : class PhaseSpace2to2diffractive : public PhaseSpace {
360 :
361 : public:
362 :
363 : // Constructor.
364 0 : PhaseSpace2to2diffractive(bool isDiffAin = false, bool isDiffBin = false)
365 0 : : isDiffA(isDiffAin), isDiffB(isDiffBin) {}
366 :
367 : // Construct the trial or final event kinematics.
368 : virtual bool setupSampling();
369 : virtual bool trialKin(bool inEvent = true, bool = false);
370 : virtual bool finalKin();
371 :
372 : // Are beam particles resolved in partons or scatter directly?
373 0 : virtual bool isResolved() const {return false;}
374 :
375 : private:
376 :
377 : // Constants: could only be changed in the code itself.
378 : static const int NTRY;
379 : static const double EXPMAX, DIFFMASSMARGIN;
380 :
381 : // Initialization data, in constructor or read from Settings.
382 : bool isDiffA, isDiffB;
383 : int PomFlux;
384 : double epsilonPF, alphaPrimePF;
385 :
386 : // Initialization: kinematics properties specific to 2 -> 2 diffractive.
387 : double m3ElDiff, m4ElDiff, s1, s2, lambda12, lambda34, tLow, tUpp,
388 : cRes, sResXB, sResAX, sProton, bMin, bSlope, bSlope1, bSlope2,
389 : probSlope1, xIntPF, xtCorPF, mp24DL, coefDL, tAux, tAux1, tAux2;
390 :
391 : // Parameters for MBR model.
392 : double sdpmax, ddpmax, dymin0, dymax, amx, am1, am2, t;
393 : double eps, alph, alph2, m2min, dyminSD, dyminDD, dyminSigSD, dyminSigDD;
394 :
395 : };
396 :
397 : //==========================================================================
398 :
399 : // A derived class with 2 -> 3 kinematics set up for central diffractive
400 : // scattering.
401 :
402 0 : class PhaseSpace2to3diffractive : public PhaseSpace {
403 :
404 : public:
405 :
406 : // Constructor.
407 0 : PhaseSpace2to3diffractive() {}
408 :
409 : // Construct the trial or final event kinematics.
410 : virtual bool setupSampling();
411 : virtual bool trialKin(bool inEvent = true, bool = false);
412 : virtual bool finalKin();
413 :
414 : // Are beam particles resolved in partons or scatter directly?
415 0 : virtual bool isResolved() const {return false;}
416 :
417 : private:
418 :
419 : // Constants: could only be changed in the code itself.
420 : static const int NTRY, NINTEG2;
421 : static const double EXPMAX, DIFFMASSMIN, DIFFMASSMARGIN;
422 :
423 : // Local variables to calculate DPE kinematics.
424 : int PomFlux;
425 : double epsilonPF, alphaPrimePF, s1, s2, m5min, s5min, tLow[2], tUpp[2],
426 : bMin[2], tAux[2], bSlope1, bSlope2, probSlope1[2], tAux1[2],
427 : tAux2[2], bSlope, xIntPF, xIntInvPF, xtCorPF, mp24DL, coefDL,
428 : epsMBR, alphMBR, m2minMBR, dyminMBR, dyminSigMBR, dyminInvMBR,
429 : dpepmax, t1, t2;
430 : Vec4 p1, p2, p3, p4, p5;
431 :
432 : };
433 :
434 : //==========================================================================
435 :
436 : // A derived class for nondiffractive events. Hardly does anything, since
437 : // the real action is taken care of by the MultipartonInteractions class.
438 :
439 0 : class PhaseSpace2to2nondiffractive : public PhaseSpace {
440 :
441 : public:
442 :
443 : // Constructor.
444 0 : PhaseSpace2to2nondiffractive() {}
445 :
446 : // Construct the trial or final event kinematics.
447 0 : virtual bool setupSampling() {sigmaNw = sigmaProcessPtr->sigmaHat();
448 0 : sigmaMx = sigmaNw; return true;}
449 0 : virtual bool trialKin( bool , bool = false) {return true;}
450 0 : virtual bool finalKin() {return true;}
451 :
452 : private:
453 :
454 : };
455 :
456 : //==========================================================================
457 :
458 : // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
459 : // tau, y, pT2_4, pT2_5, phi_4, phi_5 and y_3 (partial cylindrical symmetry).
460 :
461 0 : class PhaseSpace2to3tauycyl : public PhaseSpace {
462 :
463 : public:
464 :
465 : // Constructor.
466 0 : PhaseSpace2to3tauycyl() {}
467 :
468 : // Optimize subsequent kinematics selection.
469 0 : virtual bool setupSampling() {if (!setupMasses()) return false;
470 0 : setup3Body(); return setupSampling123(false, true);}
471 :
472 : // Construct the trial kinematics.
473 : virtual bool trialKin(bool inEvent = true, bool = false) {
474 0 : if (!trialMasses()) return false;
475 0 : return trialKin123(false, true, inEvent);}
476 :
477 : // Construct the final event kinematics.
478 : virtual bool finalKin();
479 :
480 : private:
481 :
482 : // Constants: could only be changed in the code itself.
483 : static const int NITERNR;
484 :
485 : // Set up for fixed or Breit-Wigner mass selection.
486 : bool setupMasses();
487 :
488 : // Select fixed or Breit-Wigner-distributed masses.
489 : bool trialMasses();
490 :
491 : };
492 :
493 : //==========================================================================
494 :
495 : // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
496 : // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
497 : // Intended specifically for (essentially massless) 2 -> 3 QCD processes.
498 :
499 0 : class PhaseSpace2to3yyycyl : public PhaseSpace {
500 :
501 : public:
502 :
503 : // Constructor.
504 0 : PhaseSpace2to3yyycyl() {}
505 :
506 : // Optimize subsequent kinematics selection.
507 : virtual bool setupSampling();
508 :
509 : // Construct the trial kinematics.
510 : virtual bool trialKin(bool inEvent = true, bool = false);
511 :
512 : // Construct the final event kinematics.
513 : virtual bool finalKin();
514 :
515 : private:
516 :
517 : // Phase space cuts specifically for 2 -> 3 QCD processes.
518 : double pTHat3Min, pTHat3Max, pTHat5Min, pTHat5Max, RsepMin, R2sepMin;
519 : bool hasBaryonBeams;
520 :
521 : // Event kinematics choices.
522 : double pT3Min, pT3Max, pT5Min, pT5Max, y3Max, y4Max, y5Max,
523 : pT3, pT4, pT5, phi3, phi4, phi5, y3, y4, y5, dphi;
524 : Vec4 pInSum;
525 :
526 : };
527 :
528 : //==========================================================================
529 :
530 : // A derived class for Les Houches events.
531 :
532 0 : class PhaseSpaceLHA : public PhaseSpace {
533 :
534 : public:
535 :
536 : // Constructor.
537 0 : PhaseSpaceLHA() {idProcSave = 0;}
538 :
539 : // Find maximal cross section for comparison with internal processes.
540 : virtual bool setupSampling();
541 :
542 : // Construct the next process, by interface to Les Houches class.
543 : virtual bool trialKin( bool , bool repeatSame = false);
544 :
545 : // Set scale, alpha_s and alpha_em if not done.
546 0 : virtual bool finalKin() {sigmaProcessPtr->setScale(); return true;}
547 :
548 : // For Les Houches with negative event weight needs
549 0 : virtual double sigmaSumSigned() const {return sigmaSgn;}
550 :
551 : private:
552 :
553 : // Constants.
554 : static const double CONVERTPB2MB;
555 :
556 : // Local properties.
557 : int strategy, stratAbs, nProc, idProcSave;
558 : double xMaxAbsSum, xSecSgnSum, sigmaSgn;
559 : vector<int> idProc;
560 : vector<double> xMaxAbsProc;
561 :
562 : };
563 :
564 : //==========================================================================
565 :
566 : } // end namespace Pythia8
567 :
568 : #endif // Pythia8_PhaseSpace_H
|
