fd746d2ffc2a1cd41af3d6f0cbc01e43344bd7eb
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / include / MultipleInteractions.h
1 // MultipleInteractions.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 // This file contains the main classes for multiple interactions physics.
7 // SigmaMultiple stores allowed processes by in-flavour combination.
8 // MultipleInteractions: generates multiple parton-parton interactions.
9
10 #ifndef Pythia8_MultipleInteractions_H
11 #define Pythia8_MultipleInteractions_H
12
13 #include "Basics.h"
14 #include "BeamParticle.h"
15 #include "Event.h"
16 #include "Info.h"
17 #include "PythiaStdlib.h"
18 #include "Settings.h"
19 #include "SigmaTotal.h"
20 #include "SigmaProcess.h"
21 #include "StandardModel.h"
22
23 namespace Pythia8 {
24  
25 //**************************************************************************
26
27 // SigmaMultiple is a helper class to MultipleInteractions.
28 // It packs pointers to the allowed processes for different 
29 // flavour combinations and levels of ambition.
30
31 class SigmaMultiple {
32
33 public:
34
35   // Constructor.
36   SigmaMultiple() {}
37   
38   // Destructor.
39   ~SigmaMultiple() {
40     for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i];
41     for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];}   
42
43   // Initialize list of processes.
44   bool init(int inState, int processLevel);
45
46   // Calculate cross section summed over possibilities.
47   double sigma( int id1, int id2, double x1, double x2, double sHat, 
48     double tHat, double uHat, double alpS, double alpEM);
49
50   // Return one subprocess, picked according to relative cross sections.
51   SigmaProcess* sigmaSel();
52   bool swapTU() {return pickedU;}
53
54   // Return code or name of a specified process, for statistics table.
55   int    nProc() const {return nChan;}
56   int    codeProc(int iProc) const {return sigmaT[iProc]->code();}
57   string nameProc(int iProc) const {return sigmaT[iProc]->name();}
58
59 private:
60
61   // Constants: could only be changed in the code itself.
62   static const double MASSMARGIN, OTHERFRAC;
63
64   // Number of processes. Some use massive matrix elements.
65   int            nChan;
66   vector<bool>   needMasses;
67   vector<double> m3Fix, m4Fix, sHatMin;
68
69   // Vector of process list, one for t-channel and one for u-channel.
70   vector<SigmaProcess*> sigmaT, sigmaU;
71
72   // Values of cross sections in process list above.
73   vector<double> sigmaTval, sigmaUval;
74   double         sigmaTsum, sigmaUsum;
75   bool           pickedU;
76   
77 };
78  
79 //**************************************************************************
80
81 // The MultipleInteractions class contains the main methods for the 
82 // generation of multiple parton-parton interactions in hadronic collisions.
83
84 class MultipleInteractions {
85
86 public:
87
88   // Constructor.
89   MultipleInteractions() {sudExpPT.resize(NBINS+1);}
90
91   // Initialize the generation process for given beams.
92   bool init( bool doMIinit, Info* infoPtrIn, BeamParticle* beamAPtrIn, 
93     BeamParticle* beamBPtrIn, SigmaTotal* sigmaTotPtrIn, 
94     ostream& os = cout);
95
96   // Reset impact parameter choice and update the CM energy.
97   void clear() {bIsSet = false; bSetInFirst = false;
98     eCM = infoPtr->eCM(); sCM = eCM * eCM;}
99
100   // Select first = hardest pT in minbias process.
101   void pTfirst(); 
102
103   // Set up kinematics for first = hardest pT in minbias process.
104   void setupFirstSys( Event& process);
105
106   // Find whether to limit maximum scale of emissions.
107   bool limitPTmax( Event& event);
108
109   // Prepare system for evolution.
110   void prepare(double pTscale = 1000.) {
111     if (!bSetInFirst) overlapNext(pTscale);}
112
113   // Select next pT in downwards evolution.
114   double pTnext( double pTbegAll, double pTendAll);
115
116   // Set up kinematics of acceptable interaction.
117   void scatter( Event& event); 
118
119   // Get some information on current interaction.
120   double Q2Ren()     const {return pT2Ren;}
121   double alphaSH()   const {return alpS;}
122   double alphaEMH()  const {return alpEM;}
123   double x1H()       const {return x1;} 
124   double x2H()       const {return x2;} 
125   double Q2Fac()     const {return pT2Fac;}
126   double pdf1()      const {return xPDF1now;}
127   double pdf2()      const {return xPDF2now;}
128   double bMI()       const {return (bIsSet) ? bNow / bAvg : 0.;}
129   double enhanceMI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
130
131   // Update and print statistics on number of processes.
132   void accumulate() { int iBeg = (infoPtr->isMinBias()) ? 0 : 1; 
133     for (int i = iBeg; i < infoPtr->nMI(); ++i) ++nGen[ infoPtr->codeMI(i) ];}
134   void statistics(bool reset = false, ostream& os = cout);
135   
136 private: 
137
138   // Constants: could only be changed in the code itself.
139   static const bool   SHIFTFACSCALE;
140   static const int    NBINS;
141   static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, EXPPOWMIN,
142                       PROBATLOWB, BSTEP, BMAX, EXPMAX, KCONVERGE, 
143                       CONVERT2MB;
144
145   // Initialization data, read from Settings.
146   int    pTmaxMatch, alphaSorder, alphaEMorder, bProfile, 
147          processLevel, nQuarkIn, nSample;
148   double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, 
149          pTmin, coreRadius, coreFraction, expPow;
150
151   // Other initialization data.
152   bool   hasLowPow;
153   double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR, 
154          pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax, 
155          pT4dProbMax, dSigmaApprox, sigmaInt, zeroIntCorr, normOverlap, 
156          nAvg, kNow, normPi, bAvg, bDiv, probLowB, radius2B, radius2C, 
157          fracA, fracB, fracC, fracAhigh, fracBhigh, fracChigh, fracABChigh, 
158          expRev, cDiv, cMax;
159   vector<double> sudExpPT;
160
161   // Properties specific to current system.
162   int    id1, id2;
163   double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2, 
164          tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now;
165   bool   bIsSet, bSetInFirst, isAtLowB;
166
167   // Pointer to various information on the generation.
168   Info*  infoPtr;
169
170   // Pointers to the two incoming beams.
171   BeamParticle* beamAPtr;
172   BeamParticle* beamBPtr;
173
174   // Pointer to total cross section parametrization.
175   SigmaTotal* sigmaTotPtr;
176
177   // Collections of parton-level 2 -> 2 cross sections. Selected one.
178   SigmaMultiple sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
179   SigmaProcess* dSigmaDtSel;
180
181   // Statistics on generated 2 -> 2 processes.
182   map<int, int> nGen;
183
184   // alphaStrong and alphaEM calculations.
185   AlphaStrong alphaS;
186   AlphaEM alphaEM;
187
188   // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.  
189   void upperEnvelope();
190
191   // Integrate the parton-parton interaction cross section.
192   void jetCrossSection();
193
194   // Evaluate "Sudakov form factor" for not having a harder interaction.
195   double sudakov(double pT2sud, double enhance = 1.);
196
197   // Do a quick evolution towards the next smaller pT.
198   double fastPT2( double pT2beg);
199
200   // Calculate the actual cross section, either for the first interaction
201   // (including at initialization) or for any subsequent in the sequence. 
202   double sigmaPT2(bool isFirst = false);
203
204   // Calculate factor relating matter overlap and interaction rate.
205   void overlapInit();
206
207   // Pick impact parameter and interaction rate enhancement,
208   // either before the first interaction (for minbias) or after it.
209   void overlapFirst();
210   void overlapNext(double pTscale);
211
212 };
213  
214 //**************************************************************************
215
216 } // end namespace Pythia8
217
218 #endif // Pythia8_MultipleInteractions_H