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