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