using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / include / PhaseSpace.h
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