]>
Commit | Line | Data |
---|---|---|
63ba5337 | 1 | // MultipartonInteractions.h is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2012 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 multiparton interactions physics. | |
7 | // SigmaMultiparton stores allowed processes by in-flavour combination. | |
8 | // MultipartonInteractions: generates multiparton parton-parton interactions. | |
9 | ||
10 | #ifndef Pythia8_MultipartonInteractions_H | |
11 | #define Pythia8_MultipartonInteractions_H | |
12 | ||
13 | #include "Basics.h" | |
14 | #include "BeamParticle.h" | |
15 | #include "Event.h" | |
16 | #include "Info.h" | |
17 | #include "PartonSystems.h" | |
18 | #include "PythiaStdlib.h" | |
19 | #include "Settings.h" | |
20 | #include "SigmaTotal.h" | |
21 | #include "SigmaProcess.h" | |
22 | #include "StandardModel.h" | |
23 | #include "UserHooks.h" | |
24 | ||
25 | namespace Pythia8 { | |
26 | ||
27 | //========================================================================== | |
28 | ||
29 | // SigmaMultiparton is a helper class to MultipartonInteractions. | |
30 | // It packs pointers to the allowed processes for different | |
31 | // flavour combinations and levels of ambition. | |
32 | ||
33 | class SigmaMultiparton { | |
34 | ||
35 | public: | |
36 | ||
37 | // Constructor. | |
38 | SigmaMultiparton() {} | |
39 | ||
40 | // Destructor. | |
41 | ~SigmaMultiparton() { | |
42 | for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i]; | |
43 | for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];} | |
44 | ||
45 | // Initialize list of processes. | |
46 | bool init(int inState, int processLevel, Info* infoPtr, | |
47 | Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtrIn, | |
48 | BeamParticle* beamAPtr, BeamParticle* beamBPtr, Couplings* couplingsPtr); | |
49 | ||
50 | // Calculate cross section summed over possibilities. | |
51 | double sigma( int id1, int id2, double x1, double x2, double sHat, | |
52 | double tHat, double uHat, double alpS, double alpEM, | |
53 | bool restore = false, bool pickOtherIn = false); | |
54 | ||
55 | // Return whether the other, rare processes were selected. | |
56 | bool pickedOther() {return pickOther;} | |
57 | ||
58 | // Return one subprocess, picked according to relative cross sections. | |
59 | SigmaProcess* sigmaSel(); | |
60 | bool swapTU() {return pickedU;} | |
61 | ||
62 | // Return code or name of a specified process, for statistics table. | |
63 | int nProc() const {return nChan;} | |
64 | int codeProc(int iProc) const {return sigmaT[iProc]->code();} | |
65 | string nameProc(int iProc) const {return sigmaT[iProc]->name();} | |
66 | ||
67 | private: | |
68 | ||
69 | // Constants: could only be changed in the code itself. | |
70 | static const double MASSMARGIN, OTHERFRAC; | |
71 | ||
72 | // Number of processes. Some use massive matrix elements. | |
73 | int nChan; | |
74 | vector<bool> needMasses; | |
75 | vector<double> m3Fix, m4Fix, sHatMin; | |
76 | ||
77 | // Vector of process list, one for t-channel and one for u-channel. | |
78 | vector<SigmaProcess*> sigmaT, sigmaU; | |
79 | ||
80 | // Values of cross sections in process list above. | |
81 | vector<double> sigmaTval, sigmaUval; | |
82 | double sigmaTsum, sigmaUsum; | |
83 | bool pickOther, pickedU; | |
84 | ||
85 | // Pointer to the random number generator. | |
86 | Rndm* rndmPtr; | |
87 | ||
88 | }; | |
89 | ||
90 | //========================================================================== | |
91 | ||
92 | // The MultipartonInteractions class contains the main methods for the | |
93 | // generation of multiparton parton-parton interactions in hadronic collisions. | |
94 | ||
95 | class MultipartonInteractions { | |
96 | ||
97 | public: | |
98 | ||
99 | // Constructor. | |
100 | MultipartonInteractions() : bIsSet(false) {} | |
101 | ||
102 | // Initialize the generation process for given beams. | |
103 | bool init( bool doMPIinit, int iDiffSysIn, Info* infoPtrIn, | |
104 | Settings& settings, ParticleData* particleDataPtr, Rndm* rndmPtrIn, | |
105 | BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, | |
106 | Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn, | |
107 | SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn, | |
108 | ostream& os = cout); | |
109 | ||
110 | // Reset impact parameter choice and update the CM energy. | |
111 | void reset(); | |
112 | ||
113 | // Select first = hardest pT in minbias process. | |
114 | void pTfirst(); | |
115 | ||
116 | // Set up kinematics for first = hardest pT in minbias process. | |
117 | void setupFirstSys( Event& process); | |
118 | ||
119 | // Find whether to limit maximum scale of emissions. | |
120 | bool limitPTmax( Event& event); | |
121 | ||
122 | // Prepare system for evolution. | |
123 | void prepare(Event& event, double pTscale = 1000.) { | |
124 | if (!bSetInFirst) overlapNext(event, pTscale); } | |
125 | ||
126 | // Select next pT in downwards evolution. | |
127 | double pTnext( double pTbegAll, double pTendAll, Event& event); | |
128 | ||
129 | // Set up kinematics of acceptable interaction. | |
130 | bool scatter( Event& event); | |
131 | ||
132 | // Set "empty" values to avoid query of undefined quantities. | |
133 | void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac | |
134 | = xPDF1now = xPDF2now = 0.; bIsSet = false;} | |
135 | ||
136 | // Get some information on current interaction. | |
137 | double Q2Ren() const {return pT2Ren;} | |
138 | double alphaSH() const {return alpS;} | |
139 | double alphaEMH() const {return alpEM;} | |
140 | double x1H() const {return x1;} | |
141 | double x2H() const {return x2;} | |
142 | double Q2Fac() const {return pT2Fac;} | |
143 | double pdf1() const {return xPDF1now;} | |
144 | double pdf2() const {return xPDF2now;} | |
145 | double bMPI() const {return (bIsSet) ? bNow / bAvg : 0.;} | |
146 | double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;} | |
147 | ||
148 | // For x-dependent matter profile, return incoming valence/sea | |
149 | // decision from trial interactions. | |
150 | int getVSC1() const {return vsc1;} | |
151 | int getVSC2() const {return vsc2;} | |
152 | ||
153 | // Update and print statistics on number of processes. | |
154 | // Note: currently only valid for MB systems, not diffraction?? | |
155 | void accumulate() { int iBeg = (infoPtr->isMinBias()) ? 0 : 1; | |
156 | for (int i = iBeg; i < infoPtr->nMPI(); ++i) | |
157 | ++nGen[ infoPtr->codeMPI(i) ];} | |
158 | void statistics(bool resetStat = false, ostream& os = cout); | |
159 | void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin(); | |
160 | iter != nGen.end(); ++iter) iter->second = 0; } | |
161 | ||
162 | private: | |
163 | ||
164 | // Constants: could only be changed in the code itself. | |
165 | static const bool SHIFTFACSCALE, PREPICKRESCATTER; | |
166 | static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN, | |
167 | EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX, | |
168 | KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV, WTACCWARN; | |
169 | ||
170 | // Initialization data, read from Settings. | |
171 | bool allowRescatter, allowDoubleRes, canVetoMPI; | |
172 | int pTmaxMatch, alphaSorder, alphaEMorder, bProfile, processLevel, | |
173 | rescatterMode, nQuarkIn, nSample, enhanceScreening; | |
174 | double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius, | |
175 | coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP, | |
176 | mMaxPertDiff, mMinPertDiff; | |
177 | ||
178 | // x-dependent matter profile: | |
179 | // Constants. | |
180 | static const int XDEP_BBIN; | |
181 | static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC, | |
182 | XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM; | |
183 | ||
184 | // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast | |
185 | // integration. Only needed during initialisation. | |
186 | vector <double> sigmaIntWgt, sigmaSumWgt; | |
187 | ||
188 | // a1: value of a1 constant, taken from settings database. | |
189 | // a0now (a02now): tuned value of a0 (squared value). | |
190 | // bstepNow: current size of bins in b. | |
191 | // a2max: given an xMin, a maximal (squared) value of the | |
192 | // width, to be used in overestimate Omax(b). | |
193 | // enhanceBmax, retain enhanceB as enhancement factor for the hardest | |
194 | // enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2, | |
195 | // and enhanceBnow to store the value for the current | |
196 | // interaction. | |
197 | double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow; | |
198 | ||
199 | // Storage for trial interactions. | |
200 | int id1Save, id2Save; | |
201 | double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave, | |
202 | alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave, | |
203 | xPDF2nowSave; | |
204 | SigmaProcess *dSigmaDtSelSave; | |
205 | ||
206 | // vsc1, vsc2: for minimum-bias events with trial interaction, store | |
207 | // decision on whether hard interaction was valence or sea. | |
208 | int vsc1, vsc2; | |
209 | ||
210 | // Other initialization data. | |
211 | bool hasBaryonBeams, hasLowPow, globalRecoilFSR; | |
212 | int iDiffSys, nMaxGlobalRecoilFSR; | |
213 | double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR, | |
214 | pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax, | |
215 | pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101], | |
216 | zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv, | |
217 | probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh, | |
218 | fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax; | |
219 | ||
220 | // Properties specific to current system. | |
221 | bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel; | |
222 | int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel; | |
223 | double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2, | |
224 | tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now, | |
225 | dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel; | |
226 | ||
227 | // Stored values for mass interpolation for diffractive systems. | |
228 | int nStep, iStepFrom, iStepTo; | |
229 | double eCMsave, eStepSize, eStepSave, eStepFrom, eStepTo, pT0Save[5], | |
230 | pT4dSigmaMaxSave[5], pT4dProbMaxSave[5], sigmaIntSave[5], | |
231 | sudExpPTSave[5][101], zeroIntCorrSave[5], normOverlapSave[5], | |
232 | kNowSave[5], bAvgSave[5], bDivSave[5], probLowBSave[5], | |
233 | fracAhighSave[5], fracBhighSave[5], fracChighSave[5], | |
234 | fracABChighSave[5], cDivSave[5], cMaxSave[5]; | |
235 | ||
236 | // Pointer to various information on the generation. | |
237 | Info* infoPtr; | |
238 | ||
239 | // Pointer to the random number generator. | |
240 | Rndm* rndmPtr; | |
241 | ||
242 | // Pointers to the two incoming beams. | |
243 | BeamParticle* beamAPtr; | |
244 | BeamParticle* beamBPtr; | |
245 | ||
246 | // Pointers to Standard Model couplings. | |
247 | Couplings* couplingsPtr; | |
248 | ||
249 | // Pointer to information on subcollision parton locations. | |
250 | PartonSystems* partonSystemsPtr; | |
251 | ||
252 | // Pointer to total cross section parametrization. | |
253 | SigmaTotal* sigmaTotPtr; | |
254 | ||
255 | // Pointer to user hooks. | |
256 | UserHooks* userHooksPtr; | |
257 | ||
258 | // Collections of parton-level 2 -> 2 cross sections. Selected one. | |
259 | SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq; | |
260 | SigmaMultiparton* sigma2Sel; | |
261 | SigmaProcess* dSigmaDtSel; | |
262 | ||
263 | // Statistics on generated 2 -> 2 processes. | |
264 | map<int, int> nGen; | |
265 | ||
266 | // alphaStrong and alphaEM calculations. | |
267 | AlphaStrong alphaS; | |
268 | AlphaEM alphaEM; | |
269 | ||
270 | // Scattered partons. | |
271 | vector<int> scatteredA, scatteredB; | |
272 | ||
273 | // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2. | |
274 | void upperEnvelope(); | |
275 | ||
276 | // Integrate the parton-parton interaction cross section. | |
277 | void jetCrossSection(); | |
278 | ||
279 | // Evaluate "Sudakov form factor" for not having a harder interaction. | |
280 | double sudakov(double pT2sud, double enhance = 1.); | |
281 | ||
282 | // Do a quick evolution towards the next smaller pT. | |
283 | double fastPT2( double pT2beg); | |
284 | ||
285 | // Calculate the actual cross section, either for the first interaction | |
286 | // (including at initialization) or for any subsequent in the sequence. | |
287 | double sigmaPT2scatter(bool isFirst = false); | |
288 | ||
289 | // Find the partons that may rescatter. | |
290 | void findScatteredPartons( Event& event); | |
291 | ||
292 | // Calculate the actual cross section for a rescattering. | |
293 | double sigmaPT2rescatter( Event& event); | |
294 | ||
295 | // Calculate factor relating matter overlap and interaction rate. | |
296 | void overlapInit(); | |
297 | ||
298 | // Pick impact parameter and interaction rate enhancement, | |
299 | // either before the first interaction (for minbias) or after it. | |
300 | void overlapFirst(); | |
301 | void overlapNext(Event& event, double pTscale); | |
302 | ||
303 | }; | |
304 | ||
305 | //========================================================================== | |
306 | ||
307 | } // end namespace Pythia8 | |
308 | ||
309 | #endif // Pythia8_MultipartonInteractions_H |