StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MultipartonInteractions.cc
1 // MultipartonInteractions.cc 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 // Function definitions (not found in the header) for the
7 // SigmaMultiparton and MultipartonInteractions classes.
8 
9 #include "MultipartonInteractions.h"
10 
11 // Internal headers for special processes.
12 #include "SigmaQCD.h"
13 #include "SigmaEW.h"
14 #include "SigmaOnia.h"
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // The SigmaMultiparton class.
21 
22 //--------------------------------------------------------------------------
23 
24 // Constants: could be changed here if desired, but normally should not.
25 // These are of technical nature, as described for each.
26 
27 // The sum of outgoing masses must not be too close to the cm energy.
28 const double SigmaMultiparton::MASSMARGIN = 0.1;
29 
30 // Fraction of time not the dominant "gluon t-channel" process is picked.
31 const double SigmaMultiparton::OTHERFRAC = 0.2;
32 
33 //--------------------------------------------------------------------------
34 
35 // Initialize the generation process for given beams.
36 
37 bool SigmaMultiparton::init(int inState, int processLevel, Info* infoPtr,
38  Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtrIn,
39  BeamParticle* beamAPtr, BeamParticle* beamBPtr, Couplings* couplingsPtr) {
40 
41  // Store input pointer for future use.
42  rndmPtr = rndmPtrIn;
43 
44  // Reset vector sizes (necessary in case of re-initialization).
45  if (sigmaT.size() > 0) {
46  for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i];
47  sigmaT.resize(0);
48  }
49  if (sigmaU.size() > 0) {
50  for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];
51  sigmaU.resize(0);
52  }
53 
54  // Always store mimimal set of processes:QCD 2 -> 2 t-channel.
55 
56  // Gluon-gluon instate.
57  if (inState == 0) {
58  sigmaT.push_back( new Sigma2gg2gg() );
59  sigmaU.push_back( new Sigma2gg2gg() );
60 
61  // Quark-gluon instate.
62  } else if (inState == 1) {
63  sigmaT.push_back( new Sigma2qg2qg() );
64  sigmaU.push_back( new Sigma2qg2qg() );
65 
66  // Quark-(anti)quark instate.
67  } else {
68  sigmaT.push_back( new Sigma2qq2qq() );
69  sigmaU.push_back( new Sigma2qq2qq() );
70  }
71 
72  // Normally store QCD processes to new flavour.
73  if (processLevel > 0) {
74  if (inState == 0) {
75  sigmaT.push_back( new Sigma2gg2qqbar() );
76  sigmaU.push_back( new Sigma2gg2qqbar() );
77  sigmaT.push_back( new Sigma2gg2QQbar(4, 121) );
78  sigmaU.push_back( new Sigma2gg2QQbar(4, 121) );
79  sigmaT.push_back( new Sigma2gg2QQbar(5, 123) );
80  sigmaU.push_back( new Sigma2gg2QQbar(5, 123) );
81  } else if (inState == 2) {
82  sigmaT.push_back( new Sigma2qqbar2gg() );
83  sigmaU.push_back( new Sigma2qqbar2gg() );
84  sigmaT.push_back( new Sigma2qqbar2qqbarNew() );
85  sigmaU.push_back( new Sigma2qqbar2qqbarNew() );
86  sigmaT.push_back( new Sigma2qqbar2QQbar(4, 122) );
87  sigmaU.push_back( new Sigma2qqbar2QQbar(4, 122) );
88  sigmaT.push_back( new Sigma2qqbar2QQbar(5, 124) );
89  sigmaU.push_back( new Sigma2qqbar2QQbar(5, 124) );
90  }
91  }
92 
93  // Optionally store electroweak processes, mainly photon production.
94  if (processLevel > 1) {
95  if (inState == 0) {
96  sigmaT.push_back( new Sigma2gg2ggamma() );
97  sigmaU.push_back( new Sigma2gg2ggamma() );
98  sigmaT.push_back( new Sigma2gg2gammagamma() );
99  sigmaU.push_back( new Sigma2gg2gammagamma() );
100  } else if (inState == 1) {
101  sigmaT.push_back( new Sigma2qg2qgamma() );
102  sigmaU.push_back( new Sigma2qg2qgamma() );
103  } else if (inState == 2) {
104  sigmaT.push_back( new Sigma2qqbar2ggamma() );
105  sigmaU.push_back( new Sigma2qqbar2ggamma() );
106  sigmaT.push_back( new Sigma2ffbar2gammagamma() );
107  sigmaU.push_back( new Sigma2ffbar2gammagamma() );
108  sigmaT.push_back( new Sigma2ffbar2ffbarsgm() );
109  sigmaU.push_back( new Sigma2ffbar2ffbarsgm() );
110  }
111  if (inState >= 2) {
112  sigmaT.push_back( new Sigma2ff2fftgmZ() );
113  sigmaU.push_back( new Sigma2ff2fftgmZ() );
114  sigmaT.push_back( new Sigma2ff2fftW() );
115  sigmaU.push_back( new Sigma2ff2fftW() );
116  }
117  }
118 
119  // Optionally store charmonium and bottomonium production.
120  if (processLevel > 2) {
121  if (inState == 0) {
122  sigmaT.push_back( new Sigma2gg2QQbar3S11g(4, 401) );
123  sigmaU.push_back( new Sigma2gg2QQbar3S11g(4, 401) );
124  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 0, 402) );
125  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 0, 402) );
126  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 1, 403) );
127  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 1, 403) );
128  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 2, 404) );
129  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 2, 404) );
130  sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 0, 411) );
131  sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 0, 411) );
132  sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 1, 412) );
133  sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 1, 412) );
134  sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 2, 413) );
135  sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 2, 413) );
136  sigmaT.push_back( new Sigma2gg2QQbar3S11g(5, 501) );
137  sigmaU.push_back( new Sigma2gg2QQbar3S11g(5, 501) );
138  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 0, 502) );
139  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 0, 502) );
140  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 1, 503) );
141  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 1, 503) );
142  sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 2, 504) );
143  sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 2, 504) );
144  sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 0, 511) );
145  sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 0, 511) );
146  sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 1, 512) );
147  sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 1, 512) );
148  sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 2, 513) );
149  sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 2, 513) );
150  } else if (inState == 1) {
151  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 0, 405) );
152  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 0, 405) );
153  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 1, 406) );
154  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 1, 406) );
155  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 2, 407) );
156  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 2, 407) );
157  sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 0, 414) );
158  sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 0, 414) );
159  sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 1, 415) );
160  sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 1, 415) );
161  sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 2, 416) );
162  sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 2, 416) );
163  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 0, 505) );
164  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 0, 505) );
165  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 1, 506) );
166  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 1, 506) );
167  sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 2, 507) );
168  sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 2, 507) );
169  sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 0, 514) );
170  sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 0, 514) );
171  sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 1, 515) );
172  sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 1, 515) );
173  sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 2, 516) );
174  sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 2, 516) );
175  } else if (inState == 2) {
176  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 0, 408) );
177  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 0, 408) );
178  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 1, 409) );
179  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 1, 409) );
180  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 2, 410) );
181  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 2, 410) );
182  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 0, 417) );
183  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 0, 417) );
184  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 1, 418) );
185  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 1, 418) );
186  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 2, 419) );
187  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 2, 419) );
188  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 0, 508) );
189  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 0, 508) );
190  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 1, 509) );
191  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 1, 509) );
192  sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 2, 510) );
193  sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 2, 510) );
194  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 0, 517) );
195  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 0, 517) );
196  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 1, 518) );
197  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 1, 518) );
198  sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 2, 519) );
199  sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 2, 519) );
200  }
201  }
202 
203  // Resize arrays to match sizes above.
204  nChan = sigmaT.size();
205  needMasses.resize(nChan);
206  m3Fix.resize(nChan);
207  m4Fix.resize(nChan);
208  sHatMin.resize(nChan);
209  sigmaTval.resize(nChan);
210  sigmaUval.resize(nChan);
211 
212  // Initialize the processes.
213  for (int i = 0; i < nChan; ++i) {
214  sigmaT[i]->init( infoPtr, settingsPtr, particleDataPtr, rndmPtr,
215  beamAPtr, beamBPtr, couplingsPtr);
216  sigmaT[i]->initProc();
217  sigmaU[i]->init( infoPtr, settingsPtr, particleDataPtr, rndmPtr,
218  beamAPtr, beamBPtr, couplingsPtr);
219  sigmaU[i]->initProc();
220 
221  // Prepare for massive kinematics (but fixed masses!) where required.
222  needMasses[i] = false;
223  int id3Mass = sigmaT[i]->id3Mass();
224  int id4Mass = sigmaT[i]->id4Mass();
225  m3Fix[i] = 0.;
226  m4Fix[i] = 0.;
227  if (id3Mass > 0 || id4Mass > 0) {
228  needMasses[i] = true;
229  m3Fix[i] = particleDataPtr->m0(id3Mass);
230  m4Fix[i] = particleDataPtr->m0(id4Mass);
231  }
232  sHatMin[i] = pow2( m3Fix[i] + m4Fix[i] + MASSMARGIN);
233  }
234 
235  // Done.
236  return true;
237 
238 }
239 
240 //--------------------------------------------------------------------------
241 
242 // Calculate cross section summed over possibilities.
243 
244 double SigmaMultiparton::sigma( int id1, int id2, double x1, double x2,
245  double sHat, double tHat, double uHat, double alpS, double alpEM,
246  bool restore, bool pickOtherIn) {
247 
248  // Choose either the dominant process (in slot 0) or the rest of them.
249  if (restore) pickOther = pickOtherIn;
250  else pickOther = (rndmPtr->flat() < OTHERFRAC);
251 
252  // Iterate over all subprocesses.
253  sigmaTsum = 0.;
254  sigmaUsum = 0.;
255  for (int i = 0; i < nChan; ++i) {
256  sigmaTval[i] = 0.;
257  sigmaUval[i] = 0.;
258 
259  // Skip the not chosen processes.
260  if (i == 0 && pickOther) continue;
261  if (i > 0 && !pickOther) continue;
262 
263  // t-channel-sampling contribution.
264  if (sHat > sHatMin[i]) {
265  sigmaT[i]->set2KinMPI( x1, x2, sHat, tHat, uHat,
266  alpS, alpEM, needMasses[i], m3Fix[i], m4Fix[i]);
267  sigmaTval[i] = sigmaT[i]->sigmaHatWrap(id1, id2);
268  sigmaT[i]->pickInState(id1, id2);
269  // Correction factor for tHat rescaling in massive kinematics.
270  if (needMasses[i]) sigmaTval[i] *= sigmaT[i]->sHBetaMPI() / sHat;
271  sigmaTsum += sigmaTval[i];
272  }
273 
274  // u-channel-sampling contribution.
275  if (sHat > sHatMin[i]) {
276  sigmaU[i]->set2KinMPI( x1, x2, sHat, uHat, tHat,
277  alpS, alpEM, needMasses[i], m3Fix[i], m4Fix[i]);
278  sigmaUval[i] = sigmaU[i]->sigmaHatWrap( id1, id2);
279  sigmaU[i]->pickInState(id1, id2);
280  // Correction factor for tHat rescaling in massive kinematics.
281  if (needMasses[i]) sigmaUval[i] *= sigmaU[i]->sHBetaMPI() / sHat;
282  sigmaUsum += sigmaUval[i];
283  }
284 
285  // Average of t- and u-channel sampling; corrected for not selected channels.
286  }
287  double sigmaAvg = 0.5 * (sigmaTsum + sigmaUsum);
288  if (pickOther) sigmaAvg /= OTHERFRAC;
289  if (!pickOther) sigmaAvg /= (1. - OTHERFRAC);
290  return sigmaAvg;
291 
292 }
293 
294 //--------------------------------------------------------------------------
295 
296 // Return one subprocess, picked according to relative cross sections.
297 
298 SigmaProcess* SigmaMultiparton::sigmaSel() {
299 
300  // Decide between t- and u-channel-sampled kinematics.
301  pickedU = (rndmPtr->flat() * (sigmaTsum + sigmaUsum) < sigmaUsum);
302 
303  // Pick one of t-channel-sampled processes.
304  if (!pickedU) {
305  double sigmaRndm = sigmaTsum * rndmPtr->flat();
306  int iPick = -1;
307  do sigmaRndm -= sigmaTval[++iPick];
308  while (sigmaRndm > 0.);
309  return sigmaT[iPick];
310 
311  // Pick one of u-channel-sampled processes.
312  } else {
313  double sigmaRndm = sigmaUsum * rndmPtr->flat();
314  int iPick = -1;
315  do sigmaRndm -= sigmaUval[++iPick];
316  while (sigmaRndm > 0.);
317  return sigmaU[iPick];
318  }
319 
320 }
321 
322 //==========================================================================
323 
324 // The MultipartonInteractions class.
325 
326 //--------------------------------------------------------------------------
327 
328 // Constants: could be changed here if desired, but normally should not.
329 // These are of technical nature, as described for each.
330 
331 // Factorization scale pT2 by default, but could be shifted to pT2 + pT02.
332 // (A priori possible, but problems for flavour threshold interpretation.)
333 const bool MultipartonInteractions::SHIFTFACSCALE = false;
334 
335 // Pick one parton to represent rescattering cross section to speed up.
336 const bool MultipartonInteractions::PREPICKRESCATTER = true;
337 
338 // Naive upper estimate of cross section too pessimistic, so reduce by this.
339 const double MultipartonInteractions::SIGMAFUDGE = 0.7;
340 
341 // The r value above, picked to allow a flatter correct/trial cross section.
342 const double MultipartonInteractions::RPT20 = 0.25;
343 
344 // Reduce pT0 by factor pT0STEP if sigmaInt < SIGMASTEP * sigmaND.
345 const double MultipartonInteractions::PT0STEP = 0.9;
346 const double MultipartonInteractions::SIGMASTEP = 1.1;
347 
348 // Stop if pT0 or pTmin fall below this, or alpha_s blows up.
349 const double MultipartonInteractions::PT0MIN = 0.2;
350 
351 // Refuse too low expPow in impact parameter profile.
352 const double MultipartonInteractions::EXPPOWMIN = 0.4;
353 
354 // Define low-b region by interaction rate above given number.
355 const double MultipartonInteractions::PROBATLOWB = 0.6;
356 
357 // Basic step size for b integration; sometimes modified.
358 const double MultipartonInteractions::BSTEP = 0.01;
359 
360 // Stop b integration when integrand dropped enough.
361 const double MultipartonInteractions::BMAX = 1e-8;
362 
363 // Do not allow too large argument to exp function.
364 const double MultipartonInteractions::EXPMAX = 50.;
365 
366 // Convergence criterion for k iteration.
367 const double MultipartonInteractions::KCONVERGE = 1e-7;
368 
369 // Conversion of GeV^{-2} to mb for cross section.
370 const double MultipartonInteractions::CONVERT2MB = 0.389380;
371 
372 // Stay away from division by zero in Jacobian for tHat -> pT2.
373 const double MultipartonInteractions::ROOTMIN = 0.01;
374 
375 // No need to reinitialize parameters if energy close to previous.
376 const double MultipartonInteractions::ECMDEV = 0.01;
377 
378 // Settings for x-dependent matter profile:
379 // Number of bins in b (with these settings, no bStep increase and
380 // reintegration needed with a1 ~ 0.20 up to ECM ~ 40TeV).
381 const int MultipartonInteractions::XDEP_BBIN = 500;
382 // Initial value of a0.
383 const double MultipartonInteractions::XDEP_A0 = 1.0;
384 // Width of form ( XDEP_A1 + a1 * log(1 / x) ).
385 const double MultipartonInteractions::XDEP_A1 = 1.0;
386 // Initial step size in b and increment.
387 const double MultipartonInteractions::XDEP_BSTEP = 0.02;
388 const double MultipartonInteractions::XDEP_BSTEPINC = 0.01;
389 // Overlap-weighted cross section in last bin of b must be beneath
390 // XDEP_CUTOFF * sigmaInt.
391 const double MultipartonInteractions::XDEP_CUTOFF = 1e-4;
392 // a0 is calculated in units of sqrt(mb), so convert to fermi.
393 const double MultipartonInteractions::XDEP_SMB2FM = sqrt(0.1);
394 
395 //--------------------------------------------------------------------------
396 
397 // Initialize the generation process for given beams.
398 
399 bool MultipartonInteractions::init( bool doMPIinit, int diffractiveModeIn,
400  Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtr,
401  Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
402  Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
403  SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn, ostream& os) {
404 
405  // Store input pointers for future use. Done if no initialization.
406  diffractiveMode = diffractiveModeIn;
407  infoPtr = infoPtrIn;
408  rndmPtr = rndmPtrIn;
409  beamAPtr = beamAPtrIn;
410  beamBPtr = beamBPtrIn;
411  couplingsPtr = couplingsPtrIn;
412  partonSystemsPtr = partonSystemsPtrIn;
413  sigmaTotPtr = sigmaTotPtrIn;
414  userHooksPtr = userHooksPtrIn;
415  if (!doMPIinit) return false;
416 
417  // If both beams are baryons then softer PDF's than for mesons/Pomerons.
418  hasBaryonBeams = ( beamAPtr->isBaryon() && beamBPtr->isBaryon() );
419 
420  // Matching in pT of hard interaction to further interactions.
421  pTmaxMatch = settings.mode("MultipartonInteractions:pTmaxMatch");
422 
423  // Parameters of alphaStrong generation.
424  alphaSvalue = settings.parm("MultipartonInteractions:alphaSvalue");
425  alphaSorder = settings.mode("MultipartonInteractions:alphaSorder");
426 
427  // Parameters of alphaEM generation.
428  alphaEMorder = settings.mode("MultipartonInteractions:alphaEMorder");
429 
430  // Parameters of cross section generation.
431  Kfactor = settings.parm("MultipartonInteractions:Kfactor");
432 
433  // Regularization of QCD evolution for pT -> 0.
434  pT0Ref = settings.parm("MultipartonInteractions:pT0Ref");
435  ecmRef = settings.parm("MultipartonInteractions:ecmRef");
436  ecmPow = settings.parm("MultipartonInteractions:ecmPow");
437  pTmin = settings.parm("MultipartonInteractions:pTmin");
438 
439  // Impact parameter profile: nondiffractive topologies.
440  if (diffractiveMode == 0) {
441  bProfile = settings.mode("MultipartonInteractions:bProfile");
442  coreRadius = settings.parm("MultipartonInteractions:coreRadius");
443  coreFraction = settings.parm("MultipartonInteractions:coreFraction");
444  expPow = settings.parm("MultipartonInteractions:expPow");
445  expPow = max(EXPPOWMIN, expPow);
446  // x-dependent impact parameter profile.
447  a1 = settings.parm("MultipartonInteractions:a1");
448 
449  // Impact parameter profile: diffractive topologies.
450  } else {
451  bProfile = settings.mode("Diffraction:bProfile");
452  coreRadius = settings.parm("Diffraction:coreRadius");
453  coreFraction = settings.parm("Diffraction:coreFraction");
454  expPow = settings.parm("Diffraction:expPow");
455  expPow = max(EXPPOWMIN, expPow);
456  }
457 
458  // Process sets to include in machinery.
459  processLevel = settings.mode("MultipartonInteractions:processLevel");
460 
461  // Parameters of rescattering description.
462  allowRescatter = settings.flag("MultipartonInteractions:allowRescatter");
463  allowDoubleRes
464  = settings.flag("MultipartonInteractions:allowDoubleRescatter");
465  rescatterMode = settings.mode("MultipartonInteractions:rescatterMode");
466  ySepResc = settings.parm("MultipartonInteractions:ySepRescatter");
467  deltaYResc = settings.parm("MultipartonInteractions:deltaYRescatter");
468 
469  // Rescattering not yet implemented for x-dependent impact profile.
470  if (bProfile == 4) allowRescatter = false;
471 
472  // A global recoil FSR stategy restricts rescattering.
473  globalRecoilFSR = settings.flag("TimeShower:globalRecoil");
474  nMaxGlobalRecoilFSR = settings.mode("TimeShower:nMaxGlobalRecoil");
475 
476  // Various other parameters.
477  nQuarkIn = settings.mode("MultipartonInteractions:nQuarkIn");
478  nSample = settings.mode("MultipartonInteractions:nSample");
479 
480  // Optional dampening at small pT's when large multiplicities.
481  enhanceScreening = settings.mode("MultipartonInteractions:enhanceScreening");
482 
483  // Parameters for diffractive systems.
484  sigmaPomP = settings.parm("Diffraction:sigmaRefPomP");
485  mPomP = settings.parm("Diffraction:mRefPomP");
486  pPomP = settings.parm("Diffraction:mPowPomP");
487  mMinPertDiff = settings.parm("Diffraction:mMinPert");
488 
489  // Possibility to allow user veto of MPI
490  canVetoMPI = (userHooksPtr > 0) ? userHooksPtr->canVetoMPIEmission() : false;
491 
492  // Some common combinations for double Gaussian, as shorthand.
493  if (bProfile == 2) {
494  fracA = pow2(1. - coreFraction);
495  fracB = 2. * coreFraction * (1. - coreFraction);
496  fracC = pow2(coreFraction);
497  radius2B = 0.5 * (1. + pow2(coreRadius));
498  radius2C = pow2(coreRadius);
499 
500  // Some common combinations for exp(b^pow), as shorthand.
501  } else if (bProfile == 3) {
502  hasLowPow = (expPow < 2.);
503  expRev = 2. / expPow - 1.;
504  }
505 
506  // Initialize alpha_strong generation.
507  alphaS.init( alphaSvalue, alphaSorder);
508  double Lambda3 = alphaS.Lambda3();
509 
510  // Initialize alphaEM generation.
511  alphaEM.init( alphaEMorder, &settings);
512 
513  // Attach matrix-element calculation objects.
514  sigma2gg.init( 0, processLevel, infoPtr, &settings, particleDataPtr,
515  rndmPtr, beamAPtr, beamBPtr, couplingsPtr);
516  sigma2qg.init( 1, processLevel, infoPtr, &settings, particleDataPtr,
517  rndmPtr, beamAPtr, beamBPtr, couplingsPtr);
518  sigma2qqbarSame.init( 2, processLevel, infoPtr, &settings, particleDataPtr,
519  rndmPtr, beamAPtr, beamBPtr, couplingsPtr);
520  sigma2qq.init( 3, processLevel, infoPtr, &settings, particleDataPtr,
521  rndmPtr, beamAPtr, beamBPtr, couplingsPtr);
522 
523  // Calculate invariant mass of system.
524  eCM = infoPtr->eCM();
525  sCM = eCM * eCM;
526  mMaxPertDiff = eCM;
527  eCMsave = eCM;
528 
529  // Get the total inelastic and nondiffractive cross section.
530  if (!sigmaTotPtr->hasSigmaTot()) return false;
531  bool isNonDiff = (diffractiveMode == 0);
532  sigmaND = sigmaTotPtr->sigmaND();
533  double sigmaMaxViol = 0.;
534 
535  // Output initialization info - first part.
536  bool showMPI = settings.flag("Init:showMultipartonInteractions");
537  if (showMPI) {
538  os << "\n *------- PYTHIA Multiparton Interactions Initialization "
539  << "---------* \n"
540  << " | "
541  << " | \n";
542  if (isNonDiff)
543  os << " | minbias, sigmaNonDiffractive = " << fixed
544  << setprecision(2) << setw(7) << sigmaND << " mb | \n";
545  else if (diffractiveMode == 1)
546  os << " | diffraction XB "
547  << " | \n";
548  else if (diffractiveMode == 2)
549  os << " | diffraction AX "
550  << " | \n";
551  os << " | "
552  << " | \n";
553  }
554 
555  // For diffraction need to cover range of diffractive masses.
556  nStep = (diffractiveMode == 0) ? 1 : 5;
557  eStepSize = (nStep < 2) ? 1.
558  : log(mMaxPertDiff / mMinPertDiff) / (nStep - 1.);
559  for (int iStep = 0; iStep < nStep; ++iStep) {
560 
561  // Update and output current diffractive mass and
562  // fictitious Pomeron-proton cross section for normalization.
563  if (nStep > 1) {
564  eCM = mMinPertDiff * pow( mMaxPertDiff / mMinPertDiff,
565  iStep / (nStep - 1.) );
566  sCM = eCM * eCM;
567  sigmaND = sigmaPomP * pow( eCM / mPomP, pPomP);
568  if (showMPI) os << " | diffractive mass = " << scientific
569  << setprecision(3) << setw(9) << eCM << " GeV and sigmaNorm = "
570  << fixed << setw(6) << sigmaND << " mb | \n";
571  }
572 
573  // Set current pT0 scale.
574  pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
575 
576  // The pT0 value may need to be decreased, if sigmaInt < sigmaND.
577  double pT4dSigmaMaxBeg = 0.;
578  for ( ; ; ) {
579 
580  // Derived pT kinematics combinations.
581  pT20 = pT0*pT0;
582  pT2min = pTmin*pTmin;
583  pTmax = 0.5*eCM;
584  pT2max = pTmax*pTmax;
585  pT20R = RPT20 * pT20;
586  pT20minR = pT2min + pT20R;
587  pT20maxR = pT2max + pT20R;
588  pT20min0maxR = pT20minR * pT20maxR;
589  pT2maxmin = pT2max - pT2min;
590 
591  // Provide upper estimate of interaction rate d(Prob)/d(pT2).
592  upperEnvelope();
593 
594  // Setup binning in b for x-dependent matter profile.
595  if (bProfile == 4) {
596  sigmaIntWgt.resize(XDEP_BBIN);
597  sigmaSumWgt.resize(XDEP_BBIN);
598  bstepNow = XDEP_BSTEP;
599  }
600 
601  // Integrate the parton-parton interaction cross section.
602  pT4dSigmaMaxBeg = pT4dSigmaMax;
603  jetCrossSection();
604 
605  // If the overlap-weighted cross section has not fallen below
606  // cutoff, then increase bin size in b and reintegrate.
607  while (bProfile == 4
608  && sigmaIntWgt[XDEP_BBIN - 1] > XDEP_CUTOFF * sigmaInt) {
609  bstepNow += XDEP_BSTEPINC;
610  jetCrossSection();
611  }
612 
613  // Sufficiently big SigmaInt or reduce pT0; maybe also pTmin.
614  if (sigmaInt > SIGMASTEP * sigmaND) break;
615  if (showMPI) os << fixed << setprecision(2) << " | pT0 = "
616  << setw(5) << pT0 << " gives sigmaInteraction = " << setw(7)
617  << sigmaInt << " mb: rejected | \n";
618  if (pTmin > pT0) pTmin *= PT0STEP;
619  pT0 *= PT0STEP;
620 
621  // Give up if pT0 and pTmin fall too low.
622  if ( max(pT0, pTmin) < max(PT0MIN, Lambda3) ) {
623  infoPtr->errorMsg("Error in MultipartonInteractions::init:"
624  " failed to find acceptable pT0 and pTmin");
625  infoPtr->setTooLowPTmin(true);
626  return false;
627  }
628  }
629 
630  // Output for accepted pT0.
631  if (showMPI) os << fixed << setprecision(2) << " | pT0 = "
632  << setw(5) << pT0 << " gives sigmaInteraction = "<< setw(7)
633  << sigmaInt << " mb: accepted | \n";
634 
635  // Calculate factor relating matter overlap and interaction rate.
636  overlapInit();
637 
638  // Maximum violation relative to first estimate.
639  sigmaMaxViol = max( sigmaMaxViol, pT4dSigmaMax / pT4dSigmaMaxBeg);
640 
641  // Save values calculated.
642  if (nStep > 1) {
643  pT0Save[iStep] = pT0;
644  pT4dSigmaMaxSave[iStep] = pT4dSigmaMax;
645  pT4dProbMaxSave[iStep] = pT4dProbMax;
646  sigmaIntSave[iStep] = sigmaInt;
647  for (int j = 0; j <= 100; ++j) sudExpPTSave[iStep][j] = sudExpPT[j];
648  zeroIntCorrSave[iStep] = zeroIntCorr;
649  normOverlapSave[iStep] = normOverlap;
650  kNowSave[iStep] = kNow;
651  bAvgSave[iStep] = bAvg;
652  bDivSave[iStep] = bDiv;
653  probLowBSave[iStep] = probLowB;
654  fracAhighSave[iStep] = fracAhigh;
655  fracBhighSave[iStep] = fracBhigh;
656  fracChighSave[iStep] = fracBhigh;
657  fracABChighSave[iStep] = fracABChigh;
658  cDivSave[iStep] = cDiv;
659  cMaxSave[iStep] = cMax;
660  }
661 
662  // End of loop over diffractive masses.
663  }
664 
665  // Output details for x-dependent matter profile.
666  if (bProfile == 4 && showMPI)
667  os << " | "
668  << " | \n"
669  << fixed << setprecision(2)
670  << " | x-dependent matter profile: a1 = " << a1 << ", "
671  << "a0 = " << a0now * XDEP_SMB2FM << ", bStep = "
672  << bstepNow << " | \n";
673 
674  // End initialization printout.
675  if (showMPI) os << " | "
676  << " | \n"
677  << " *------- End PYTHIA Multiparton Interactions Initialization"
678  << " -----* " << endl;
679 
680  // Amount of violation from upperEnvelope to jetCrossSection.
681  if (sigmaMaxViol > 1.) {
682  ostringstream osWarn;
683  osWarn << "by factor " << fixed << setprecision(3) << sigmaMaxViol;
684  infoPtr->errorMsg("Warning in MultipartonInteractions::init:"
685  " maximum increased", osWarn.str());
686  }
687 
688  // Reset statistics.
689  SigmaMultiparton* dSigma;
690  for (int i = 0; i < 4; ++i) {
691  if (i == 0) dSigma = &sigma2gg;
692  else if (i == 1) dSigma = &sigma2qg;
693  else if (i == 2) dSigma = &sigma2qqbarSame;
694  else dSigma = &sigma2qq;
695  int nProc = dSigma->nProc();
696  for (int iProc = 0; iProc < nProc; ++iProc)
697  nGen[ dSigma->codeProc(iProc) ] = 0;
698  }
699 
700  // Additional setup for x-dependent matter profile.
701  if (bProfile == 4) {
702  sigmaIntWgt.clear();
703  sigmaSumWgt.clear();
704  }
705  // No preselection of sea/valence content and initialise a0.
706  vsc1 = 0;
707  vsc2 = 0;
708 
709  // Done.
710  return true;
711 }
712 
713 //--------------------------------------------------------------------------
714 
715 // Reset impact parameter choice and update the CM energy.
716 // For diffraction also interpolate parameters to current CM energy.
717 
718 void MultipartonInteractions::reset( ) {
719 
720  // Reset impact parameter choice.
721  bIsSet = false;
722  bSetInFirst = false;
723 
724  // Update CM energy. Done if not diffraction and not new energy.
725  eCM = infoPtr->eCM();
726  sCM = eCM * eCM;
727  if (nStep == 1 || abs( eCM / eCMsave - 1.) < ECMDEV) return;
728 
729  // Set fictitious Pomeron-proton cross section for diffractive system.
730  sigmaND = sigmaPomP * pow( eCM / mPomP, pPomP);
731 
732  // Current interpolation point.
733  eCMsave = eCM;
734  eStepSave = log(eCM / mMinPertDiff) / eStepSize;
735  iStepFrom = max( 0, min( nStep - 2, int( eStepSave) ) );
736  iStepTo = iStepFrom + 1;
737  eStepTo = max( 0., min( 1., eStepSave - iStepFrom) );
738  eStepFrom = 1. - eStepTo;
739 
740  // Update pT0 and combinations derived from it.
741  pT0 = eStepFrom * pT0Save[iStepFrom]
742  + eStepTo * pT0Save[iStepTo];
743  pT20 = pT0*pT0;
744  pT2min = pTmin*pTmin;
745  pTmax = 0.5*eCM;
746  pT2max = pTmax*pTmax;
747  pT20R = RPT20 * pT20;
748  pT20minR = pT2min + pT20R;
749  pT20maxR = pT2max + pT20R;
750  pT20min0maxR = pT20minR * pT20maxR;
751  pT2maxmin = pT2max - pT2min;
752 
753  // Update other parameters used in pT choice.
754  pT4dSigmaMax = eStepFrom * pT4dSigmaMaxSave[iStepFrom]
755  + eStepTo * pT4dSigmaMaxSave[iStepTo];
756  pT4dProbMax = eStepFrom * pT4dProbMaxSave[iStepFrom]
757  + eStepTo * pT4dProbMaxSave[iStepTo];
758  sigmaInt = eStepFrom * sigmaIntSave[iStepFrom]
759  + eStepTo * sigmaIntSave[iStepTo];
760  for (int j = 0; j <= 100; ++j)
761  sudExpPT[j] = eStepFrom * sudExpPTSave[iStepFrom][j]
762  + eStepTo * sudExpPTSave[iStepTo][j];
763 
764  // Update parameters related to the impact-parameter picture.
765  zeroIntCorr = eStepFrom * zeroIntCorrSave[iStepFrom]
766  + eStepTo * zeroIntCorrSave[iStepTo];
767  normOverlap = eStepFrom * normOverlapSave[iStepFrom]
768  + eStepTo * normOverlapSave[iStepTo];
769  kNow = eStepFrom * kNowSave[iStepFrom]
770  + eStepTo * kNowSave[iStepTo];
771  bAvg = eStepFrom * bAvgSave[iStepFrom]
772  + eStepTo * bAvgSave[iStepTo];
773  bDiv = eStepFrom * bDivSave[iStepFrom]
774  + eStepTo * bDivSave[iStepTo];
775  probLowB = eStepFrom * probLowBSave[iStepFrom]
776  + eStepTo * probLowBSave[iStepTo];
777  fracAhigh = eStepFrom * fracAhighSave[iStepFrom]
778  + eStepTo * fracAhighSave[iStepTo];
779  fracBhigh = eStepFrom * fracBhighSave[iStepFrom]
780  + eStepTo * fracBhighSave[iStepTo];
781  fracChigh = eStepFrom * fracChighSave[iStepFrom]
782  + eStepTo * fracChighSave[iStepTo];
783  fracABChigh = eStepFrom * fracABChighSave[iStepFrom]
784  + eStepTo * fracABChighSave[iStepTo];
785  cDiv = eStepFrom * cDivSave[iStepFrom]
786  + eStepTo * cDivSave[iStepTo];
787  cMax = eStepFrom * cMaxSave[iStepFrom]
788  + eStepTo * cMaxSave[iStepTo];
789 
790 }
791 
792 //--------------------------------------------------------------------------
793 
794 // Select first = hardest pT in minbias process.
795 // Requires separate treatment at low and high b values.
796 
797 void MultipartonInteractions::pTfirst() {
798  // Pick impact parameter and thereby interaction rate enhancement.
799  // This is not used for the x-dependent matter profile, which
800  // instead uses trial interactions.
801  if (bProfile == 4) isAtLowB = false;
802  else overlapFirst();
803  bSetInFirst = true;
804  double WTacc;
805 
806  // At low b values evolve downwards with Sudakov.
807  if (isAtLowB) {
808  pT2 = pT2max;
809  do {
810 
811  // Pick a pT using a quick-and-dirty cross section estimate.
812  pT2 = fastPT2(pT2);
813 
814  // If fallen below lower cutoff then need to restart at top.
815  if (pT2 < pT2min) {
816  pT2 = pT2max;
817  WTacc = 0.;
818 
819  // Else pick complete kinematics and evaluate cross-section correction.
820  } else {
821  WTacc = sigmaPT2scatter(true) / dSigmaApprox;
822  if (WTacc > 1.) infoPtr->errorMsg("Warning in "
823  "MultipartonInteractions::pTfirst: weight above unity");
824  }
825 
826  // Loop until acceptable pT and acceptable kinematics.
827  } while (WTacc < rndmPtr->flat() || !dSigmaDtSel->final2KinMPI());
828 
829  // At high b values make preliminary pT choice without Sudakov factor.
830  } else {
831 
832  while (true) {
833  do {
834  pT2 = pT20min0maxR / (pT20minR + rndmPtr->flat() * pT2maxmin) - pT20R;
835 
836  // Evaluate upper estimate of cross section for this pT2 choice.
837  dSigmaApprox = pT4dSigmaMax / pow2(pT2 + pT20R);
838 
839  // Pick complete kinematics and evaluate cross-section correction.
840  WTacc = sigmaPT2scatter(true) / dSigmaApprox;
841 
842  // Evaluate and include Sudakov factor.
843  if (bProfile != 4) WTacc *= sudakov( pT2, enhanceB);
844 
845  // Warn for weight above unity
846  if (WTacc > 1.) infoPtr->errorMsg("Warning in "
847  "MultipartonInteractions::pTfirst: weight above unity");
848 
849  // Loop until acceptable pT and acceptable kinematics.
850  } while (WTacc < rndmPtr->flat() || !dSigmaDtSel->final2KinMPI());
851 
852  // For x-dependent matter profile, use trial interactions to
853  // generate Sudakov, otherwise done.
854  if (bProfile != 4) break;
855  else {
856  // Save details of the original hard interaction.
857  pT2Save = pT2;
858  id1Save = id1;
859  id2Save = id2;
860  x1Save = x1;
861  x2Save = x2;
862  sHatSave = sHat;
863  tHatSave = tHat;
864  uHatSave = uHat;
865  alpSsave = alpS;
866  alpEMsave = alpEM;
867  pT2FacSave = pT2Fac;
868  pT2RenSave = pT2Ren;
869  xPDF1nowSave = xPDF1now;
870  xPDF2nowSave = xPDF2now;
871  // Save accepted kinematics and pointer to SigmaProcess.
872  dSigmaDtSel->saveKin();
873  dSigmaDtSelSave = dSigmaDtSel;
874 
875  // Put x1, x2 information into beam pointers to get correct
876  // PDF rescaling in trial interaction (for hard process, can
877  // be sea or valence, not companion).
878  beamAPtr->append( 0, id1, x1);
879  beamAPtr->xfISR( 0, id1, x1, pT2Fac * pT2Fac);
880  vsc1 = beamAPtr->pickValSeaComp();
881  beamBPtr->append( 0, id2, x2);
882  beamBPtr->xfISR( 0, id2, x2, pT2Fac * pT2Fac);
883  vsc2 = beamBPtr->pickValSeaComp();
884 
885  // Pick b according to O(b, x1, x2).
886  double w1 = XDEP_A1 + a1 * log(1. / x1);
887  double w2 = XDEP_A1 + a1 * log(1. / x2);
888  double fac = a02now * (w1 * w1 + w2 * w2);
889  double expb2 = rndmPtr->flat();
890  b2now = - fac * log(expb2);
891  bNow = sqrt(b2now);
892 
893  // Enhancement factor for the hard process and overestimate
894  // for fastPT2. Note that existing framework has a (1. / sigmaND)
895  // present.
896  enhanceB = sigmaND / M_PI / fac * expb2;
897  enhanceBmax = sigmaND / 2. / M_PI / a02now *
898  exp( -b2now / 2. / a2max );
899 
900  // Trial interaction with dummy event.
901  Event evDummy;
902  double pTtrial = pTnext(pTmax, pTmin, evDummy);
903 
904  // Restore beams.
905  beamAPtr->clear();
906  beamBPtr->clear();
907 
908  // Accept if fallen beneath factorisation scale.
909  if (pTtrial < sqrt(pT2FacSave)) {
910  // Restore previous values and original sigma.
911  swap(pT2, pT2Save);
912  swap(pT2Fac, pT2FacSave);
913  swap(pT2Ren, pT2RenSave);
914  swap(id1, id1Save);
915  swap(id2, id2Save);
916  swap(x1, x1Save);
917  swap(x2, x2Save);
918  swap(sHat, sHatSave);
919  swap(tHat, tHatSave);
920  swap(uHat, uHatSave);
921  swap(alpS, alpSsave);
922  swap(alpEM, alpEMsave);
923  swap(xPDF1now, xPDF1nowSave);
924  swap(xPDF2now, xPDF2nowSave);
925  if (dSigmaDtSel == dSigmaDtSelSave) dSigmaDtSel->swapKin();
926  else swap(dSigmaDtSel, dSigmaDtSelSave);
927 
928  // Accept.
929  bIsSet = true;
930  break;
931  }
932  } // if (bProfile == 4)
933  } // while (true)
934 
935  // End handling for high b.
936  }
937 
938 }
939 
940 //--------------------------------------------------------------------------
941 
942 // Set up kinematics for first = hardest pT in minbias process.
943 
944 void MultipartonInteractions::setupFirstSys( Event& process) {
945 
946  // Last beam-status particles. Offset relative to normal beam locations.
947  int sizeProc = process.size();
948  int nBeams = 3;
949  for (int i = 3; i < sizeProc; ++i)
950  if (process[i].statusAbs() < 20) nBeams = i + 1;
951  int nOffset = nBeams - 3;
952 
953  // Remove any partons of previous failed interactions.
954  if (sizeProc > nBeams) {
955  process.popBack( sizeProc - nBeams);
956  process.initColTag();
957  }
958 
959  // Entries 3 and 4, now to be added, come from 1 and 2.
960  process[1 + nOffset].daughter1(3 + nOffset);
961  process[2 + nOffset].daughter1(4 + nOffset);
962 
963  // Negate beam status, if not already done. (Case with offset beams.)
964  process[1 + nOffset].statusNeg();
965  process[2 + nOffset].statusNeg();
966 
967  // Loop over four partons and offset info relative to subprocess itself.
968  int colOffset = process.lastColTag();
969  for (int i = 1; i <= 4; ++i) {
970  Particle parton = dSigmaDtSel->getParton(i);
971  if (i <= 2 ) parton.mothers( i + nOffset, 0);
972  else parton.mothers( 3 + nOffset, 4 + nOffset);
973  if (i <= 2 ) parton.daughters( 5 + nOffset, 6 + nOffset);
974  else parton.daughters( 0, 0);
975  int col = parton.col();
976  if (col > 0) parton.col( col + colOffset);
977  int acol = parton.acol();
978  if (acol > 0) parton.acol( acol + colOffset);
979 
980  // Put the partons into the event record.
981  process.append(parton);
982  }
983 
984  // Set scale from which to begin evolution.
985  process.scale( sqrt(pT2Fac) );
986 
987  // Info on subprocess - specific to mimimum-bias events.
988  string nameSub = dSigmaDtSel->name();
989  int codeSub = dSigmaDtSel->code();
990  int nFinalSub = dSigmaDtSel->nFinal();
991  double pTMPI = dSigmaDtSel->pTMPIFin();
992  infoPtr->setSubType( nameSub, codeSub, nFinalSub);
993  infoPtr->setTypeMPI( codeSub, pTMPI, 0, 0, enhanceB / zeroIntCorr);
994 
995  // Further standard info on process.
996  infoPtr->setPDFalpha( id1, id2, xPDF1now, xPDF2now, pT2Fac, alpEM, alpS,
997  pT2Ren);
998  double m3 = dSigmaDtSel->m(3);
999  double m4 = dSigmaDtSel->m(4);
1000  double theta = dSigmaDtSel->thetaMPI();
1001  double phi = dSigmaDtSel->phiMPI();
1002  infoPtr->setKin( x1, x2, sHat, tHat, uHat, sqrt(pT2), m3, m4, theta, phi);
1003 
1004 }
1005 
1006 //--------------------------------------------------------------------------
1007 
1008 // Find whether to limit maximum scale of emissions.
1009 
1010 bool MultipartonInteractions::limitPTmax( Event& event) {
1011 
1012  // User-set cases.
1013  if (pTmaxMatch == 1) return true;
1014  if (pTmaxMatch == 2) return false;
1015 
1016  // Look if only quarks (u, d, s, c, b), gluons and photons in final state.
1017  bool onlyQGP = true;
1018  for (int i = 5; i < event.size(); ++i)
1019  if (event[i].status() != -21) {
1020  int idAbs = event[i].idAbs();
1021  if (idAbs > 5 && idAbs != 21 && idAbs != 22) onlyQGP = false;
1022  }
1023  return (onlyQGP);
1024 
1025 }
1026 
1027 //--------------------------------------------------------------------------
1028 
1029 // Select next pT in downwards evolution.
1030 
1031 double MultipartonInteractions::pTnext( double pTbegAll, double pTendAll,
1032  Event& event) {
1033 
1034  // Initial values.
1035  bool pickRescatter, acceptKin;
1036  double dSigmaScatter, dSigmaRescatter, WTacc;
1037  double pT2end = pow2( max(pTmin, pTendAll) );
1038 
1039  // With the x-dependent matter profile, it is possible to reuse
1040  // values that have been stored during trial interactions for a
1041  // slight speedup. bIsSet is false during trial interactions,
1042  // counter 21 in case partonLevel is retried and counter 22 for
1043  // the first pass through partonLevel.
1044  if (bProfile == 4 && bIsSet && infoPtr->getCounter(21) == 1
1045  && infoPtr->getCounter(22) == 1) {
1046 
1047  // Minimum bias.
1048  if (bSetInFirst) {
1049  if (pT2Save < pT2end) return 0.;
1050  pT2 = pT2Save;
1051  pT2Fac = pT2FacSave;
1052  pT2Ren = pT2RenSave;
1053  id1 = id1Save;
1054  id2 = id2Save;
1055  x1 = x1Save;
1056  x2 = x2Save;
1057  sHat = sHatSave;
1058  tHat = tHatSave;
1059  uHat = uHatSave;
1060  alpS = alpSsave;
1061  alpEM = alpEMsave;
1062  xPDF1now = xPDF1nowSave;
1063  xPDF2now = xPDF2nowSave;
1064  if (dSigmaDtSel == dSigmaDtSelSave) dSigmaDtSel->swapKin();
1065  else dSigmaDtSel = dSigmaDtSelSave;
1066  return sqrt(pT2);
1067 
1068  // Hard process.
1069  } else {
1070  return (pT2 < pT2end) ? 0. : sqrt(pT2);
1071  }
1072  }
1073 
1074  // Do not allow rescattering while sill FSR with global recoil.
1075  bool allowRescatterNow = allowRescatter;
1076  if (globalRecoilFSR && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoilFSR)
1077  allowRescatterNow = false;
1078 
1079  // Initial pT2 value.
1080  pT2 = pow2(pTbegAll);
1081 
1082  // Find the set of already scattered partons on the two sides.
1083  if (allowRescatterNow) findScatteredPartons( event);
1084 
1085  // Pick a pT2 using a quick-and-dirty cross section estimate.
1086  do {
1087  do {
1088  pT2 = fastPT2(pT2);
1089  if (pT2 < pT2end) return 0.;
1090 
1091  // Initial values: no rescattering.
1092  i1Sel = 0;
1093  i2Sel = 0;
1094  dSigmaSum = 0.;
1095  pickRescatter = false;
1096 
1097  // Pick complete kinematics and evaluate interaction cross-section.
1098  dSigmaScatter = sigmaPT2scatter(false);
1099 
1100  // Also cross section from rescattering if allowed.
1101  dSigmaRescatter = (allowRescatterNow) ? sigmaPT2rescatter( event) : 0.;
1102 
1103  // Normalize to dSigmaApprox, which was set in fastPT2 above.
1104  WTacc = (dSigmaScatter + dSigmaRescatter) / dSigmaApprox;
1105  if (WTacc > 1.) infoPtr->errorMsg("Warning in MultipartonInteractions::"
1106  "pTnext: weight above unity");
1107 
1108  // Idea suggested by Gosta Gustafson: increased screening in events
1109  // with large activity can be simulated by pT0_eff = sqrt(n) * pT0.
1110  if (enhanceScreening > 0) {
1111  int nSysNow = infoPtr->nMPI() + 1;
1112  if (enhanceScreening == 2) nSysNow += infoPtr->nISR();
1113  double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) );
1114  WTacc *= WTscreen;
1115  }
1116 
1117  // x-dependent matter profile overlap weighting.
1118  if (bProfile == 4) {
1119  double w1 = XDEP_A1 + a1 * log(1. / x1);
1120  double w2 = XDEP_A1 + a1 * log(1. / x2);
1121  double fac = a02now * (w1 * w1 + w2 * w2);
1122  // Correct enhancement factor and weight
1123  enhanceBnow = sigmaND / M_PI / fac * exp( - b2now / fac);
1124  double oWgt = enhanceBnow / enhanceBmax;
1125  if (oWgt > 1.) infoPtr->errorMsg("Warning in MultipartonInteractions::"
1126  "pTnext: overlap weight above unity");
1127  WTacc *= oWgt;
1128  }
1129 
1130  // Decide whether to keep the event based on weight.
1131  } while (WTacc < rndmPtr->flat());
1132 
1133  // When rescattering possible: new interaction or rescattering?
1134  if (allowRescatterNow) {
1135  pickRescatter = (i1Sel > 0 || i2Sel > 0);
1136 
1137  // Restore kinematics for selected scattering/rescattering.
1138  id1 = id1Sel;
1139  id2 = id2Sel;
1140  x1 = x1Sel;
1141  x2 = x2Sel;
1142  sHat = sHatSel;
1143  tHat = tHatSel;
1144  uHat = uHatSel;
1145  sigma2Sel->sigma( id1, id2, x1, x2, sHat, tHat, uHat, alpS, alpEM,
1146  true, pickOtherSel);
1147  }
1148 
1149  // Pick one of the possible channels summed above.
1150  dSigmaDtSel = sigma2Sel->sigmaSel();
1151  if (sigma2Sel->swapTU()) swap( tHat, uHat);
1152 
1153  // Decide to keep event based on kinematics (normally always OK).
1154  // Rescattering: need to provide incoming four-vectors and masses.
1155  if (pickRescatter) {
1156  Vec4 p1Res = (i1Sel == 0) ? 0.5 * eCM * x1Sel * Vec4( 0., 0., 1., 1.)
1157  : event[i1Sel].p();
1158  Vec4 p2Res = (i2Sel == 0) ? 0.5 * eCM * x2Sel * Vec4( 0., 0., -1., 1.)
1159  : event[i2Sel].p();
1160  double m1Res = (i1Sel == 0) ? 0. : event[i1Sel].m();
1161  double m2Res = (i2Sel == 0) ? 0. : event[i2Sel].m();
1162  acceptKin = dSigmaDtSel->final2KinMPI( i1Sel, i2Sel, p1Res, p2Res,
1163  m1Res, m2Res);
1164  // New interaction: already stored values enough.
1165  } else acceptKin = dSigmaDtSel->final2KinMPI();
1166  } while (!acceptKin);
1167 
1168  // Done.
1169  return sqrt(pT2);
1170 
1171 }
1172 
1173 //--------------------------------------------------------------------------
1174 
1175 // Set up the kinematics of the 2 -> 2 scattering process,
1176 // and store the scattering in the event record.
1177 
1178 bool MultipartonInteractions::scatter( Event& event) {
1179 
1180  // Last beam-status particles. Offset relative to normal beam locations.
1181  int sizeProc = event.size();
1182  int nBeams = 3;
1183  for (int i = 3; i < sizeProc; ++i)
1184  if (event[i].statusAbs() < 20) nBeams = i + 1;
1185  int nOffset = nBeams - 3;
1186 
1187  // Loop over four partons and offset info relative to subprocess itself.
1188  int motherOffset = event.size();
1189  int colOffset = event.lastColTag();
1190  for (int i = 1; i <= 4; ++i) {
1191  Particle parton = dSigmaDtSel->getParton(i);
1192  if (i <= 2 ) parton.mothers( i + nOffset, 0);
1193  else parton.mothers( motherOffset, motherOffset + 1);
1194  if (i <= 2 ) parton.daughters( motherOffset + 2, motherOffset + 3);
1195  else parton.daughters( 0, 0);
1196  int col = parton.col();
1197  if (col > 0) parton.col( col + colOffset);
1198  int acol = parton.acol();
1199  if (acol > 0) parton.acol( acol + colOffset);
1200 
1201  // Put the partons into the event record.
1202  event.append(parton);
1203  }
1204 
1205  // Allow veto of MPI. If so restore event record to before scatter.
1206  if (canVetoMPI && userHooksPtr->doVetoMPIEmission(sizeProc, event)) {
1207  event.popBack(event.size() - sizeProc);
1208  return false;
1209  }
1210 
1211  // Store participating partons as a new set in list of all systems.
1212  int iSys = partonSystemsPtr->addSys();
1213  partonSystemsPtr->setInA(iSys, motherOffset);
1214  partonSystemsPtr->setInB(iSys, motherOffset + 1);
1215  partonSystemsPtr->addOut(iSys, motherOffset + 2);
1216  partonSystemsPtr->addOut(iSys, motherOffset + 3);
1217  partonSystemsPtr->setSHat(iSys, sHat);
1218 
1219  // Tag double rescattering graphs that annihilate one initial colour.
1220  bool annihil1 = false;
1221  bool annihil2 = false;
1222  if (i1Sel > 0 && i2Sel > 0) {
1223  if (event[motherOffset].col() == event[motherOffset + 1].acol()
1224  && event[motherOffset].col() > 0) annihil1 = true;
1225  if (event[motherOffset].acol() == event[motherOffset + 1].col()
1226  && event[motherOffset].acol() > 0) annihil2 = true;
1227  }
1228 
1229  // Beam remnant A: add scattered partons to list.
1230  BeamParticle& beamA = *beamAPtr;
1231  int iA = beamA.append( motherOffset, id1, x1);
1232 
1233  // Find whether incoming partons are valence or sea, so prepared for ISR.
1234  if (i1Sel == 0) {
1235  beamA.xfISR( iA, id1, x1, pT2);
1236  beamA.pickValSeaComp();
1237 
1238  // Remove rescattered parton from final state and change history.
1239  // Propagate existing colour labels throught graph.
1240  } else {
1241  beamA[iA].companion(-10);
1242  event[i1Sel].statusNeg();
1243  event[i1Sel].daughters( motherOffset, motherOffset);
1244  event[motherOffset].mothers( i1Sel, i1Sel);
1245  int colOld = event[i1Sel].col();
1246  if (colOld > 0) {
1247  int colNew = event[motherOffset].col();
1248  for (int i = motherOffset; i < motherOffset + 4; ++i) {
1249  if (event[i].col() == colNew) event[i].col( colOld);
1250  if (event[i].acol() == colNew) event[i].acol( colOld);
1251  }
1252  }
1253  int acolOld = event[i1Sel].acol();
1254  if (acolOld > 0) {
1255  int acolNew = event[motherOffset].acol();
1256  for (int i = motherOffset; i < motherOffset + 4; ++i) {
1257  if (event[i].col() == acolNew) event[i].col( acolOld);
1258  if (event[i].acol() == acolNew) event[i].acol( acolOld);
1259  }
1260  }
1261  }
1262 
1263  // Beam remnant B: add scattered partons to list.
1264  BeamParticle& beamB = *beamBPtr;
1265  int iB = beamB.append( motherOffset + 1, id2, x2);
1266 
1267  // Find whether incoming partons are valence or sea, so prepared for ISR.
1268  if (i2Sel == 0) {
1269  beamB.xfISR( iB, id2, x2, pT2);
1270  beamB.pickValSeaComp();
1271 
1272  // Remove rescattered parton from final state and change history.
1273  // Propagate existing colour labels throught graph.
1274  } else {
1275  beamB[iB].companion(-10);
1276  event[i2Sel].statusNeg();
1277  event[i2Sel].daughters( motherOffset + 1, motherOffset + 1);
1278  event[motherOffset + 1].mothers( i2Sel, i2Sel);
1279  int colOld = event[i2Sel].col();
1280  if (colOld > 0) {
1281  int colNew = event[motherOffset + 1].col();
1282  for (int i = motherOffset; i < motherOffset + 4; ++i) {
1283  if (event[i].col() == colNew) event[i].col( colOld);
1284  if (event[i].acol() == colNew) event[i].acol( colOld);
1285  }
1286  }
1287  int acolOld = event[i2Sel].acol();
1288  if (acolOld > 0) {
1289  int acolNew = event[motherOffset + 1].acol();
1290  for (int i = motherOffset; i < motherOffset + 4; ++i) {
1291  if (event[i].col() == acolNew) event[i].col( acolOld);
1292  if (event[i].acol() == acolNew) event[i].acol( acolOld);
1293  }
1294  }
1295  }
1296 
1297  // Annihilating colour in double rescattering requires relabelling
1298  // of one colour into the other in the whole preceding event.
1299  if (annihil1 || annihil2) {
1300  int colLeft = (annihil1) ? event[motherOffset].col()
1301  : event[motherOffset].acol();
1302  int mother1 = event[motherOffset].mother1();
1303  int mother2 = event[motherOffset + 1].mother1();
1304  int colLost = (annihil1)
1305  ? event[mother1].col() + event[mother2].acol() - colLeft
1306  : event[mother1].acol() + event[mother2].col() - colLeft;
1307  for (int i = 0; i < motherOffset; ++i) {
1308  if (event[i].col() == colLost) event[i].col( colLeft );
1309  if (event[i].acol() == colLost) event[i].acol( colLeft );
1310  }
1311  }
1312 
1313  // Store info on subprocess code and rescattered partons.
1314  int codeMPI = dSigmaDtSel->code();
1315  double pTMPI = dSigmaDtSel->pTMPIFin();
1316  infoPtr->setTypeMPI( codeMPI, pTMPI, i1Sel, i2Sel,
1317  enhanceBnow / zeroIntCorr);
1318 
1319  // Done.
1320  return true;
1321 }
1322 
1323 //--------------------------------------------------------------------------
1324 
1325 // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
1326 
1327 void MultipartonInteractions::upperEnvelope() {
1328 
1329  // Initially determine constant in jet cross section upper estimate
1330  // d(sigma_approx)/d(pT2) < const / (pT2 + r * pT20)^2.
1331  pT4dSigmaMax = 0.;
1332 
1333  // Loop thorough allowed pT range logarithmically evenly.
1334  for (int iPT = 0; iPT < 100; ++iPT) {
1335  double pT = pTmin * pow( pTmax/pTmin, 0.01 * (iPT + 0.5) );
1336  pT2 = pT*pT;
1337  pT2shift = pT2 + pT20;
1338  pT2Ren = pT2shift;
1339  pT2Fac = (SHIFTFACSCALE) ? pT2shift : pT2;
1340  xT = 2. * pT / eCM;
1341 
1342  // Evaluate parton density sums at x1 = x2 = xT.
1343  double xPDF1sumMax = (9./4.) * beamAPtr->xf(21, xT, pT2Fac);
1344  for (int id = 1; id <= nQuarkIn; ++id)
1345  xPDF1sumMax += beamAPtr->xf( id, xT, pT2Fac)
1346  + beamAPtr->xf(-id, xT, pT2Fac);
1347  double xPDF2sumMax = (9./4.) * beamBPtr->xf(21, xT, pT2Fac);
1348  for (int id = 1; id <= nQuarkIn; ++id)
1349  xPDF2sumMax += beamBPtr->xf( id, xT, pT2Fac)
1350  + beamBPtr->xf(-id, xT, pT2Fac);
1351 
1352  // Evaluate alpha_strong and _EM, matrix element and phase space volume.
1353  alpS = alphaS.alphaS(pT2Ren);
1354  alpEM = alphaEM.alphaEM(pT2Ren);
1355  double dSigmaPartonApprox = CONVERT2MB * Kfactor * 0.5 * M_PI
1356  * pow2(alpS / pT2shift);
1357  double yMax = log(1./xT + sqrt(1./(xT*xT) - 1.));
1358  double volumePhSp = pow2(2. * yMax);
1359 
1360  // Final comparison to determine upper estimate.
1361  double dSigmaApproxNow = SIGMAFUDGE * xPDF1sumMax * xPDF2sumMax
1362  * dSigmaPartonApprox * volumePhSp;
1363  double pT4dSigmaNow = pow2(pT2 + pT20R) * dSigmaApproxNow;
1364  if ( pT4dSigmaNow > pT4dSigmaMax) pT4dSigmaMax = pT4dSigmaNow;
1365  }
1366 
1367  // Get wanted constant by dividing by the nondiffractive cross section.
1368  pT4dProbMax = pT4dSigmaMax / sigmaND;
1369 
1370 }
1371 
1372 //--------------------------------------------------------------------------
1373 
1374 // Integrate the parton-parton interaction cross section,
1375 // using stratified Monte Carlo sampling.
1376 // Store result in pT bins for use as Sudakov form factors.
1377 
1378 void MultipartonInteractions::jetCrossSection() {
1379 
1380  // Common factor from bin size in dpT2 / (pT2 + r * pT20)^2 and statistics.
1381  double sigmaFactor = (1. / pT20minR - 1. / pT20maxR) / (100. * nSample);
1382 
1383  // Reset overlap-weighted cross section for x-dependent matter profile.
1384  if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++)
1385  sigmaIntWgt[bBin] = 0.;
1386 
1387  // Loop through allowed pT range evenly in dpT2 / (pT2 + r * pT20)^2.
1388  sigmaInt = 0.;
1389  double dSigmaMax = 0.;
1390  sudExpPT[100] = 0.;
1391 
1392  for (int iPT = 99; iPT >= 0; --iPT) {
1393  double sigmaSum = 0.;
1394 
1395  // Reset pT-binned overlap-weigted integration.
1396  if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++)
1397  sigmaSumWgt[bBin] = 0.;
1398 
1399  // In each pT bin sample a number of random pT values.
1400  for (int iSample = 0; iSample < nSample; ++iSample) {
1401  double mappedPT2 = 1. - 0.01 * (iPT + rndmPtr->flat());
1402  pT2 = pT20min0maxR / (pT20minR + mappedPT2 * pT2maxmin) - pT20R;
1403 
1404  // Evaluate cross section dSigma/dpT2 in phase space point.
1405  double dSigma = sigmaPT2scatter(true);
1406 
1407  // Multiply by (pT2 + r * pT20)^2 to compensate for pT sampling. Sum.
1408  dSigma *= pow2(pT2 + pT20R);
1409  sigmaSum += dSigma;
1410  if (dSigma > dSigmaMax) dSigmaMax = dSigma;
1411 
1412  // Overlap-weighted cross section for x-dependent matter profile.
1413  // Note that dSigma can be 0. when points are rejected.
1414  if (bProfile == 4 && dSigma > 0.) {
1415  double w1 = XDEP_A1 + a1 * log(1. / x1);
1416  double w2 = XDEP_A1 + a1 * log(1. / x2);
1417  double fac = XDEP_A0 * XDEP_A0 * (w1 * w1 + w2 * w2);
1418  double b = 0.5 * bstepNow;
1419  for (int bBin = 0; bBin < XDEP_BBIN; bBin++) {
1420  double wgt = exp( - b * b / fac ) / fac / M_PI;
1421  sigmaSumWgt[bBin] += dSigma * wgt;
1422  b += bstepNow;
1423  }
1424  }
1425  }
1426 
1427  // Store total cross section and exponent of Sudakov.
1428  sigmaSum *= sigmaFactor;
1429  sigmaInt += sigmaSum;
1430  sudExpPT[iPT] = sudExpPT[iPT + 1] + sigmaSum / sigmaND;
1431 
1432  // Sum overlap-weighted cross section
1433  if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++) {
1434  sigmaSumWgt[bBin] *= sigmaFactor;
1435  sigmaIntWgt[bBin] += sigmaSumWgt[bBin];
1436  }
1437 
1438  // End of loop over pT values.
1439  }
1440 
1441  // Update upper estimate of differential cross section. Done.
1442  if (dSigmaMax > pT4dSigmaMax) {
1443  pT4dSigmaMax = dSigmaMax;
1444  pT4dProbMax = dSigmaMax / sigmaND;
1445  }
1446 
1447 }
1448 
1449 //--------------------------------------------------------------------------
1450 
1451 // Evaluate "Sudakov form factor" for not having a harder interaction
1452 // at the selected b value, given the pT scale of the event.
1453 
1454 double MultipartonInteractions::sudakov(double pT2sud, double enhance) {
1455 
1456  // Find bin the pT2 scale falls in.
1457  double xBin = (pT2sud - pT2min) * pT20maxR
1458  / (pT2maxmin * (pT2sud + pT20R));
1459  xBin = max(1e-6, min(100. - 1e-6, 100. * xBin) );
1460  int iBin = int(xBin);
1461 
1462  // Interpolate inside bin. Optionally include enhancement factor.
1463  double sudExp = sudExpPT[iBin] + (xBin - iBin)
1464  * (sudExpPT[iBin + 1] - sudExpPT[iBin]);
1465  return exp( -enhance * sudExp);
1466 
1467 }
1468 
1469 //--------------------------------------------------------------------------
1470 
1471 // Pick a trial next pT, based on a simple upper estimate of the
1472 // d(sigma)/d(pT2) spectrum.
1473 
1474 double MultipartonInteractions::fastPT2( double pT2beg) {
1475 
1476  // Use d(Prob)/d(pT2) < pT4dProbMax / (pT2 + r * pT20)^2.
1477  double pT20begR = pT2beg + pT20R;
1478  double pT4dProbMaxNow = pT4dProbMax * enhanceBmax;
1479  double pT2try = pT4dProbMaxNow * pT20begR
1480  / (pT4dProbMaxNow - pT20begR * log(rndmPtr->flat())) - pT20R;
1481 
1482  // Save cross section associated with ansatz above. Done.
1483  dSigmaApprox = pT4dSigmaMax / pow2(pT2try + pT20R);
1484  return pT2try;
1485 
1486 }
1487 
1488 //--------------------------------------------------------------------------
1489 
1490 // Calculate the actual cross section to decide whether fast choice is OK.
1491 // Select flavours and kinematics for interaction at given pT.
1492 // Slightly different treatment for first interaction and subsequent ones.
1493 
1494 double MultipartonInteractions::sigmaPT2scatter(bool isFirst) {
1495 
1496  // Derive recormalization and factorization scales, amd alpha_strong/em.
1497  pT2shift = pT2 + pT20;
1498  pT2Ren = pT2shift;
1499  pT2Fac = (SHIFTFACSCALE) ? pT2shift : pT2;
1500  alpS = alphaS.alphaS(pT2Ren);
1501  alpEM = alphaEM.alphaEM(pT2Ren);
1502 
1503  // Derive rapidity limits from chosen pT2.
1504  xT = 2. * sqrt(pT2) / eCM;
1505  if (xT >= 1.) return 0.;
1506  xT2 = xT*xT;
1507  double yMax = log(1./xT + sqrt(1./xT2 - 1.));
1508 
1509  // Select rapidities y3 and y4 of the two produced partons.
1510  double y3 = yMax * (2. * rndmPtr->flat() - 1.);
1511  double y4 = yMax * (2. * rndmPtr->flat() - 1.);
1512  y = 0.5 * (y3 + y4);
1513 
1514  // Reject some events at large rapidities to improve efficiency.
1515  // (Works for baryons, not pions or Pomerons if they have hard PDF's.)
1516  double WTy = (hasBaryonBeams)
1517  ? (1. - pow2(y3/yMax)) * (1. - pow2(y4/yMax)) : 1.;
1518  if (WTy < rndmPtr->flat()) return 0.;
1519 
1520  // Failure if x1 or x2 exceed what is left in respective beam.
1521  x1 = 0.5 * xT * (exp(y3) + exp(y4));
1522  x2 = 0.5 * xT * (exp(-y3) + exp(-y4));
1523  if (isFirst) {
1524  if (x1 > 1. || x2 > 1.) return 0.;
1525  } else {
1526  if (x1 > beamAPtr->xMax() || x2 > beamBPtr->xMax()) return 0.;
1527  }
1528  tau = x1 * x2;
1529 
1530  // Begin evaluate parton densities at actual x1 and x2.
1531  double xPDF1[21];
1532  double xPDF1sum = 0.;
1533  double xPDF2[21];
1534  double xPDF2sum = 0.;
1535 
1536  // For first interaction use normal densities.
1537  if (isFirst) {
1538  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1539  if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xf(21, x1, pT2Fac);
1540  else xPDF1[id+10] = beamAPtr->xf(id, x1, pT2Fac);
1541  xPDF1sum += xPDF1[id+10];
1542  }
1543  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1544  if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xf(21, x2, pT2Fac);
1545  else xPDF2[id+10] = beamBPtr->xf(id, x2, pT2Fac);
1546  xPDF2sum += xPDF2[id+10];
1547  }
1548 
1549  // For subsequent interactions use rescaled densities.
1550  } else {
1551  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1552  if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xfMPI(21, x1, pT2Fac);
1553  else xPDF1[id+10] = beamAPtr->xfMPI(id, x1, pT2Fac);
1554  xPDF1sum += xPDF1[id+10];
1555  }
1556  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1557  if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xfMPI(21, x2, pT2Fac);
1558  else xPDF2[id+10] = beamBPtr->xfMPI(id, x2, pT2Fac);
1559  xPDF2sum += xPDF2[id+10];
1560  }
1561  }
1562 
1563  // Select incoming flavours according to actual PDF's.
1564  id1 = -nQuarkIn - 1;
1565  double temp = xPDF1sum * rndmPtr->flat();
1566  do { xPDF1now = xPDF1[(++id1) + 10]; temp -= xPDF1now; }
1567  while (temp > 0. && id1 < nQuarkIn);
1568  if (id1 == 0) id1 = 21;
1569  id2 = -nQuarkIn-1;
1570  temp = xPDF2sum * rndmPtr->flat();
1571  do { xPDF2now = xPDF2[(++id2) + 10]; temp -= xPDF2now;}
1572  while (temp > 0. && id2 < nQuarkIn);
1573  if (id2 == 0) id2 = 21;
1574 
1575  // Assign pointers to processes relevant for incoming flavour choice:
1576  // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest).
1577  // Factor 4./9. per incoming gluon to compensate for preweighting.
1578  SigmaMultiparton* sigma2Tmp;
1579  double gluFac = 1.;
1580  if (id1 == 21 && id2 == 21) {
1581  sigma2Tmp = &sigma2gg;
1582  gluFac = 16. / 81.;
1583  } else if (id1 == 21 || id2 == 21) {
1584  sigma2Tmp = &sigma2qg;
1585  gluFac = 4. / 9.;
1586  } else if (id1 == -id2) sigma2Tmp = &sigma2qqbarSame;
1587  else sigma2Tmp = &sigma2qq;
1588 
1589  // Prepare to generate differential cross sections.
1590  sHat = tau * sCM;
1591  double root = sqrtpos(1. - xT2 / tau);
1592  tHat = -0.5 * sHat * (1. - root);
1593  uHat = -0.5 * sHat * (1. + root);
1594 
1595  // Evaluate cross sections, include possibility of K factor.
1596  // Use kinematics for two symmetrical configurations (tHat <-> uHat).
1597  // (Not necessary, but makes for better MC efficiency.)
1598  double dSigmaPartonCorr = Kfactor * gluFac
1599  * sigma2Tmp->sigma( id1, id2, x1, x2, sHat, tHat, uHat, alpS, alpEM);
1600 
1601  // Combine cross section, pdf's and phase space integral.
1602  double volumePhSp = pow2(2. * yMax) / WTy;
1603  double dSigmaScat = dSigmaPartonCorr * xPDF1sum * xPDF2sum * volumePhSp;
1604 
1605  // Dampen cross section at small pT values; part of formalism.
1606  // Use pT2 corrected for massive kinematics at this step.??
1607  // double pT2massive = dSigmaDtSel->pT2MPI();
1608  double pT2massive = pT2;
1609  dSigmaScat *= pow2( pT2massive / (pT2massive + pT20) );
1610 
1611  // Sum up total contribution for all scatterings and rescatterings.
1612  dSigmaSum += dSigmaScat;
1613 
1614  // Save values for comparison with rescattering processes.
1615  i1Sel = 0;
1616  i2Sel = 0;
1617  id1Sel = id1;
1618  id2Sel = id2;
1619  x1Sel = x1;
1620  x2Sel = x2;
1621  sHatSel = sHat;
1622  tHatSel = tHat;
1623  uHatSel = uHat;
1624  sigma2Sel = sigma2Tmp;
1625  pickOtherSel = sigma2Tmp->pickedOther();
1626 
1627  // For first interaction: pick one of the possible channels summed above.
1628  if (isFirst) {
1629  dSigmaDtSel = sigma2Tmp->sigmaSel();
1630  if (sigma2Tmp->swapTU()) swap( tHat, uHat);
1631  }
1632 
1633  // Done.
1634  return dSigmaScat;
1635 }
1636 
1637 //--------------------------------------------------------------------------
1638 
1639 // Find the partons that are allowed to rescatter.
1640 
1641 void MultipartonInteractions::findScatteredPartons( Event& event) {
1642 
1643  // Reset arrays.
1644  scatteredA.resize(0);
1645  scatteredB.resize(0);
1646  double yTmp, probA, probB;
1647 
1648  // Loop though the event record and catch "final" partons.
1649  for (int i = 0; i < event.size(); ++i)
1650  if (event[i].isFinal() && (event[i].idAbs() <= nQuarkIn
1651  || event[i].id() == 21)) {
1652  yTmp = event[i].y();
1653 
1654  // Different strategies to determine which partons may rescatter.
1655  switch(rescatterMode) {
1656 
1657  // Case 0: step function at origin
1658  case 0:
1659  if ( yTmp > 0.) scatteredA.push_back( i);
1660  if (-yTmp > 0.) scatteredB.push_back( i);
1661  break;
1662 
1663  // Case 1: step function as ySepResc.
1664  case 1:
1665  if ( yTmp > ySepResc) scatteredA.push_back( i);
1666  if (-yTmp > ySepResc) scatteredB.push_back( i);
1667  break;
1668 
1669  // Case 2: linear rise from ySep - deltaY to ySep + deltaY.
1670  case 2:
1671  probA = 0.5 * (1. + ( yTmp - ySepResc) / deltaYResc);
1672  if (probA > rndmPtr->flat()) scatteredA.push_back( i);
1673  probB = 0.5 * (1. + (-yTmp - ySepResc) / deltaYResc);
1674  if (probB > rndmPtr->flat()) scatteredB.push_back( i);
1675  break;
1676 
1677  // Case 3: rise like (1/2) * ( 1 + tanh((y - ySep) / deltaY) ).
1678  // Use that (1/2) (1 + tanh(x)) = 1 / (1 + exp(-2x)).
1679  case 3:
1680  probA = 1. / (1. + exp(-2. * ( yTmp - ySepResc) / deltaYResc));
1681  if (probA > rndmPtr->flat()) scatteredA.push_back( i);
1682  probB = 1. / (1. + exp(-2. * (-yTmp - ySepResc) / deltaYResc));
1683  if (probB > rndmPtr->flat()) scatteredB.push_back( i);
1684  break;
1685 
1686  // Case 4 and undefined values: all partons can rescatter.
1687  default:
1688  scatteredA.push_back( i);
1689  scatteredB.push_back( i);
1690  break;
1691 
1692  // End of loop over partons. Done.
1693  }
1694  }
1695 
1696 }
1697 
1698 //--------------------------------------------------------------------------
1699 
1700 // Rescattering contribution summed over all already scattered partons.
1701 // Calculate the actual cross section to decide whether fast choice is OK.
1702 // Select flavours and kinematics for interaction at given pT.
1703 
1704 double MultipartonInteractions::sigmaPT2rescatter( Event& event) {
1705 
1706  // Derive xT scale from chosen pT2.
1707  xT = 2. * sqrt(pT2) / eCM;
1708  if (xT >= 1.) return 0.;
1709  xT2 = xT*xT;
1710 
1711  // Pointer to cross section and total rescatter contribution.
1712  SigmaMultiparton* sigma2Tmp;
1713  double dSigmaResc = 0.;
1714 
1715  // Normally save time by picking one random scattered parton from side A.
1716  int nScatA = scatteredA.size();
1717  int iScatA = -1;
1718  if (PREPICKRESCATTER && nScatA > 0) iScatA = min( nScatA - 1,
1719  int( rndmPtr->flat() * double(nScatA) ) );
1720 
1721  // Loop over all already scattered partons from side A.
1722  for (int iScat = 0; iScat < nScatA; ++iScat) {
1723  if (PREPICKRESCATTER && iScat != iScatA) continue;
1724  int iA = scatteredA[iScat];
1725  int id1Tmp = event[iA].id();
1726 
1727  // Calculate maximum possible sHat and check whether big enough.
1728  double x1Tmp = event[iA].pPos() / eCM;
1729  double sHatMax = x1Tmp * beamBPtr->xMax() * sCM;
1730  if (sHatMax < 4. * pT2) continue;
1731 
1732  // Select rapidity separation between two produced partons.
1733  double dyMax = log( sqrt(0.25 * sHatMax / pT2)
1734  * (1. + sqrtpos(1. - 4. * pT2 / sHatMax)) );
1735  double dy = dyMax * (2. * rndmPtr->flat() - 1.);
1736 
1737  // Reconstruct kinematical variables, especially x2.
1738  // Incoming c/b masses handled better if tau != x1 * x2.
1739  double m2Tmp = event[iA].m2();
1740  double sHatTmp = m2Tmp + 4. * pT2 * pow2(cosh(dy));
1741  double x2Tmp = (sHatTmp - m2Tmp) /(x1Tmp * sCM);
1742  double tauTmp = sHatTmp / sCM;
1743  double root = sqrtpos(1. - xT2 / tauTmp);
1744  double tHatTmp = -0.5 * sHatTmp * (1. - root);
1745  double uHatTmp = -0.5 * sHatTmp * (1. + root);
1746 
1747  // Begin evaluate parton densities at actual x2.
1748  double xPDF2[21];
1749  double xPDF2sum = 0.;
1750 
1751  // Use rescaled densities, with preweighting 9/4 for gluons.
1752  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1753  if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xfMPI(21, x2Tmp, pT2Fac);
1754  else xPDF2[id+10] = beamBPtr->xfMPI(id, x2Tmp, pT2Fac);
1755  xPDF2sum += xPDF2[id+10];
1756  }
1757 
1758  // Select incoming flavour according to actual PDF's.
1759  int id2Tmp = -nQuarkIn - 1;
1760  double temp = xPDF2sum * rndmPtr->flat();
1761  do { xPDF2now = xPDF2[(++id2Tmp) + 10]; temp -= xPDF2now;}
1762  while (temp > 0. && id2Tmp < nQuarkIn);
1763  if (id2Tmp == 0) id2Tmp = 21;
1764 
1765  // Assign pointers to processes relevant for incoming flavour choice:
1766  // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest).
1767  // Factor 4./9. for incoming gluon 2 to compensate for preweighting.
1768  if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg;
1769  else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg;
1770  else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame;
1771  else sigma2Tmp = &sigma2qq;
1772  double gluFac = (id2Tmp == 21) ? 4. / 9. : 1.;
1773 
1774  // Evaluate cross sections, include possibility of K factor.
1775  // Use kinematics for two symmetrical configurations (tHat <-> uHat).
1776  // (Not necessary, but makes for better MC efficiency.)
1777  double dSigmaPartonCorr = Kfactor * gluFac
1778  * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp,
1779  uHatTmp, alpS, alpEM);
1780 
1781  // Combine cross section, pdf's and phase space integral.
1782  double volumePhSp = 4. *dyMax;
1783  double dSigmaCorr = dSigmaPartonCorr * xPDF2sum * volumePhSp;
1784 
1785  // Dampen cross section at small pT values; part of formalism.
1786  // Use pT2 corrected for massive kinematics at this step.
1787  //?? double pT2massive = dSigmaDtSel->pT2MPI();
1788  double pT2massive = pT2;
1789  dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) );
1790 
1791  // Compensate for prepicked rescattering if appropriate.
1792  if (PREPICKRESCATTER) dSigmaCorr *= nScatA;
1793 
1794  // Sum up total contribution for all scatterings or rescattering only.
1795  dSigmaSum += dSigmaCorr;
1796  dSigmaResc += dSigmaCorr;
1797 
1798  // Determine whether current rescattering should be selected.
1799  if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) {
1800  i1Sel = iA;
1801  i2Sel = 0;
1802  id1Sel = id1Tmp;
1803  id2Sel = id2Tmp;
1804  x1Sel = x1Tmp;
1805  x2Sel = x2Tmp;
1806  sHatSel = sHatTmp;
1807  tHatSel = tHatTmp;
1808  uHatSel = uHatTmp;
1809  sigma2Sel = sigma2Tmp;
1810  pickOtherSel = sigma2Tmp->pickedOther();
1811  }
1812  }
1813 
1814  // Normally save time by picking one random scattered parton from side B.
1815  int nScatB = scatteredB.size();
1816  int iScatB = -1;
1817  if (PREPICKRESCATTER && nScatB > 0) iScatB = min( nScatB - 1,
1818  int( rndmPtr->flat() * double(nScatB) ) );
1819 
1820  // Loop over all already scattered partons from side B.
1821  for (int iScat = 0; iScat < nScatB; ++iScat) {
1822  if (PREPICKRESCATTER && iScat != iScatB) continue;
1823  int iB = scatteredB[iScat];
1824  int id2Tmp = event[iB].id();
1825 
1826  // Calculate maximum possible sHat and check whether big enough.
1827  double x2Tmp = event[iB].pNeg() / eCM;
1828  double sHatMax = beamAPtr->xMax() * x2Tmp * sCM;
1829  if (sHatMax < 4. * pT2) continue;
1830 
1831  // Select rapidity separation between two produced partons.
1832  double dyMax = log( sqrt(0.25 * sHatMax / pT2)
1833  * (1. + sqrtpos(1. - 4. * pT2 / sHatMax)) );
1834  double dy = dyMax * (2. * rndmPtr->flat() - 1.);
1835 
1836  // Reconstruct kinematical variables, especially x1.
1837  // Incoming c/b masses handled better if tau != x1 * x2.
1838  double m2Tmp = event[iB].m2();
1839  double sHatTmp = m2Tmp + 4. * pT2 * pow2(cosh(dy));
1840  double x1Tmp = (sHatTmp - m2Tmp) /(x2Tmp * sCM);
1841  double tauTmp = sHatTmp / sCM;
1842  double root = sqrtpos(1. - xT2 / tauTmp);
1843  double tHatTmp = -0.5 * sHatTmp * (1. - root);
1844  double uHatTmp = -0.5 * sHatTmp * (1. + root);
1845 
1846  // Begin evaluate parton densities at actual x1.
1847  double xPDF1[21];
1848  double xPDF1sum = 0.;
1849 
1850  // Use rescaled densities, with preweighting 9/4 for gluons.
1851  for (int id = -nQuarkIn; id <= nQuarkIn; ++id) {
1852  if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xfMPI(21, x1Tmp, pT2Fac);
1853  else xPDF1[id+10] = beamAPtr->xfMPI(id, x1Tmp, pT2Fac);
1854  xPDF1sum += xPDF1[id+10];
1855  }
1856 
1857  // Select incoming flavour according to actual PDF's.
1858  int id1Tmp = -nQuarkIn - 1;
1859  double temp = xPDF1sum * rndmPtr->flat();
1860  do { xPDF1now = xPDF1[(++id1Tmp) + 10]; temp -= xPDF1now;}
1861  while (temp > 0. && id1Tmp < nQuarkIn);
1862  if (id1Tmp == 0) id1Tmp = 21;
1863 
1864  // Assign pointers to processes relevant for incoming flavour choice:
1865  // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest).
1866  // Factor 4./9. for incoming gluon 2 to compensate for preweighting.
1867  if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg;
1868  else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg;
1869  else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame;
1870  else sigma2Tmp = &sigma2qq;
1871  double gluFac = (id1Tmp == 21) ? 4. / 9. : 1.;
1872 
1873  // Evaluate cross sections, include possibility of K factor.
1874  // Use kinematics for two symmetrical configurations (tHat <-> uHat).
1875  // (Not necessary, but makes for better MC efficiency.)
1876  double dSigmaPartonCorr = Kfactor * gluFac
1877  * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp,
1878  uHatTmp, alpS, alpEM);
1879 
1880  // Combine cross section, pdf's and phase space integral.
1881  double volumePhSp = 4. *dyMax;
1882  double dSigmaCorr = dSigmaPartonCorr * xPDF1sum * volumePhSp;
1883 
1884  // Dampen cross section at small pT values; part of formalism.
1885  // Use pT2 corrected for massive kinematics at this step.
1886  //?? double pT2massive = dSigmaDtSel->pT2MPI();
1887  double pT2massive = pT2;
1888  dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) );
1889 
1890  // Compensate for prepicked rescattering if appropriate.
1891  if (PREPICKRESCATTER) dSigmaCorr *= nScatB;
1892 
1893  // Sum up total contribution for all scatterings or rescattering only.
1894  dSigmaSum += dSigmaCorr;
1895  dSigmaResc += dSigmaCorr;
1896 
1897  // Determine whether current rescattering should be selected.
1898  if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) {
1899  i1Sel = 0;
1900  i2Sel = iB;
1901  id1Sel = id1Tmp;
1902  id2Sel = id2Tmp;
1903  x1Sel = x1Tmp;
1904  x2Sel = x2Tmp;
1905  sHatSel = sHatTmp;
1906  tHatSel = tHatTmp;
1907  uHatSel = uHatTmp;
1908  sigma2Sel = sigma2Tmp;
1909  pickOtherSel = sigma2Tmp->pickedOther();
1910  }
1911  }
1912 
1913  // Double rescatter: loop over already scattered partons from both sides.
1914  // As before, allow to use only one parton per side to speed up.
1915  if (allowDoubleRes) {
1916  for (int iScat1 = 0; iScat1 < nScatA; ++iScat1) {
1917  if (PREPICKRESCATTER && iScat1 != iScatA) continue;
1918  int iA = scatteredA[iScat1];
1919  int id1Tmp = event[iA].id();
1920  for (int iScat2 = 0; iScat2 < nScatB; ++iScat2) {
1921  if (PREPICKRESCATTER && iScat2 != iScatB) continue;
1922  int iB = scatteredB[iScat2];
1923  int id2Tmp = event[iB].id();
1924 
1925  // Calculate current sHat and check whether big enough.
1926  double sHatTmp = (event[iA].p() + event[iB].p()).m2Calc();
1927  if (sHatTmp < 4. * pT2) continue;
1928 
1929  // Reconstruct other kinematical variables. (x values not needed.)
1930  double x1Tmp = event[iA].pPos() / eCM;
1931  double x2Tmp = event[iB].pNeg() / eCM;
1932  double tauTmp = sHatTmp / sCM;
1933  double root = sqrtpos(1. - xT2 / tauTmp);
1934  double tHatTmp = -0.5 * sHatTmp * (1. - root);
1935  double uHatTmp = -0.5 * sHatTmp * (1. + root);
1936 
1937  // Assign pointers to processes relevant for incoming flavour choice:
1938  // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest).
1939  if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg;
1940  else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg;
1941  else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame;
1942  else sigma2Tmp = &sigma2qq;
1943 
1944  // Evaluate cross sections, include possibility of K factor.
1945  // Use kinematics for two symmetrical configurations (tHat <-> uHat).
1946  // (Not necessary, but makes for better MC efficiency.)
1947  double dSigmaPartonCorr = Kfactor
1948  * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp,
1949  uHatTmp, alpS, alpEM);
1950 
1951  // Combine cross section and Jacobian tHat -> pT2
1952  // (with safety minvalue).
1953  double dSigmaCorr = dSigmaPartonCorr / max(ROOTMIN, root);
1954 
1955  // Dampen cross section at small pT values; part of formalism.
1956  // Use pT2 corrected for massive kinematics at this step.
1957  //?? double pT2massive = dSigmaDtSel->pT2MPI();
1958  double pT2massive = pT2;
1959  dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) );
1960 
1961  // Compensate for prepicked rescattering if appropriate.
1962  if (PREPICKRESCATTER) dSigmaCorr *= nScatA * nScatB;
1963 
1964  // Sum up total contribution for all scatterings or rescattering only.
1965  dSigmaSum += dSigmaCorr;
1966  dSigmaResc += dSigmaCorr;
1967 
1968  // Determine whether current rescattering should be selected.
1969  if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) {
1970  i1Sel = iA;
1971  i2Sel = iB;
1972  id1Sel = id1Tmp;
1973  id2Sel = id2Tmp;
1974  x1Sel = x1Tmp;
1975  x2Sel = x2Tmp;
1976  sHatSel = sHatTmp;
1977  tHatSel = tHatTmp;
1978  uHatSel = uHatTmp;
1979  sigma2Sel = sigma2Tmp;
1980  pickOtherSel = sigma2Tmp->pickedOther();
1981  }
1982  }
1983  }
1984  }
1985 
1986  // Done.
1987  return dSigmaResc;
1988 }
1989 
1990 
1991 //--------------------------------------------------------------------------
1992 
1993 // Calculate factor relating matter overlap and interaction rate,
1994 // i.e. k in <n_interaction(b)> = k * overlap(b) (neglecting
1995 // n_int = 0 corrections and energy-momentum conservation effects).
1996 
1997 void MultipartonInteractions::overlapInit() {
1998 
1999  // Initial values for iteration. Step size of b integration.
2000  nAvg = sigmaInt / sigmaND;
2001  kNow = 0.5;
2002  int stepDir = 1;
2003  double deltaB = BSTEP;
2004  if (bProfile == 2) deltaB *= min( 0.5, 2.5 * coreRadius);
2005  if (bProfile == 3) deltaB *= max(1., pow(2. / expPow, 1. / expPow));
2006 
2007  // Further variables, with dummy initial values.
2008  double nNow = 0.;
2009  double kLow = 0.;
2010  double nLow = 0.;
2011  double kHigh = 0.;
2012  double nHigh = 0.;
2013  double overlapNow = 0.;
2014  double probNow = 0.;
2015  double overlapInt = 0.5;
2016  double probInt = 0.;
2017  double probOverlapInt = 0.;
2018  double bProbInt = 0.;
2019  normPi = 1. / (2. * M_PI);
2020 
2021  // Subdivision into low-b and high-b region by interaction rate.
2022  bool pastBDiv = false;
2023  double overlapHighB = 0.;
2024 
2025  // For x-dependent matter profile, try to find a0 rather than k.
2026  // Existing framework is used, but with substitutions:
2027  // a0 tuned according to Int( Pint(b), d^2b ) = sigmaND,
2028  // nAvg = sigmaND, kNow = a0now, kLow = a0low, kHigh = a0high,
2029  // nNow = probInt, nLow = probIntLow, nHigh = probIntHigh.
2030  double rescale2 = 1.;
2031  if (bProfile == 4) {
2032  nAvg = sigmaND;
2033  kNow = XDEP_A0 / 2.0;
2034  }
2035 
2036  // First close k into an interval by binary steps,
2037  // then find k by successive interpolation.
2038  do {
2039  if (stepDir == 1) kNow *= 2.;
2040  else if (stepDir == -1) kNow *= 0.5;
2041  else kNow = kLow + (nAvg - nLow) * (kHigh - kLow) / (nHigh - nLow);
2042 
2043  // Overlap trivial if no impact parameter dependence.
2044  if (bProfile <= 0 || bProfile > 4) {
2045  probInt = 0.5 * M_PI * (1. - exp(-kNow));
2046  probOverlapInt = probInt / M_PI;
2047  bProbInt = probInt;
2048 
2049  // Ratio of b-integrated k * overlap / (1 - exp( - k * overlap)).
2050  nNow = M_PI * kNow * overlapInt / probInt;
2051 
2052  // Else integrate overlap over impact parameter.
2053  } else if (bProfile < 4) {
2054 
2055  // Reset integrals.
2056  overlapInt = (bProfile == 3) ? 0. : 0.5;
2057  probInt = 0.;
2058  probOverlapInt = 0.;
2059  bProbInt = 0.;
2060  pastBDiv = false;
2061  overlapHighB = 0.;
2062 
2063  // Step in b space.
2064  double b = -0.5 * deltaB;
2065  double bArea = 0.;
2066  do {
2067  b += deltaB;
2068  bArea = 2. * M_PI * b * deltaB;
2069 
2070  // Evaluate overlap at current b value.
2071  if (bProfile == 1) {
2072  overlapNow = normPi * exp( -b*b);
2073  } else if (bProfile == 2) {
2074  overlapNow = normPi * ( fracA * exp( -min(EXPMAX, b*b))
2075  + fracB * exp( -min(EXPMAX, b*b / radius2B)) / radius2B
2076  + fracC * exp( -min(EXPMAX, b*b / radius2C)) / radius2C );
2077  } else {
2078  overlapNow = normPi * exp( -pow( b, expPow));
2079  overlapInt += bArea * overlapNow;
2080  }
2081  if (pastBDiv) overlapHighB += bArea * overlapNow;
2082 
2083  // Calculate interaction probability and integrate.
2084  probNow = 1. - exp( -min(EXPMAX, M_PI * kNow * overlapNow));
2085  probInt += bArea * probNow;
2086  probOverlapInt += bArea * overlapNow * probNow;
2087  bProbInt += b * bArea * probNow;
2088 
2089  // Check when interaction probability has dropped sufficiently.
2090  if (!pastBDiv && probNow < PROBATLOWB) {
2091  bDiv = b + 0.5 * deltaB;
2092  pastBDiv = true;
2093  }
2094 
2095  // Continue out in b until overlap too small.
2096  } while (b < 1. || b * probNow > BMAX);
2097 
2098  // Ratio of b-integrated k * overlap / (1 - exp( - k * overlap)).
2099  nNow = M_PI * kNow * overlapInt / probInt;
2100 
2101  // x-dependent matter profile.
2102  } else if (bProfile == 4) {
2103  rescale2 = pow2(kNow / XDEP_A0);
2104  probInt = 0.;
2105  double b = 0.5 * bstepNow;
2106  for (int bBin = 0; bBin < XDEP_BBIN; bBin++) {
2107  double bArea = 2. * M_PI * b * bstepNow;
2108  double pIntNow = 1 - exp( -min(EXPMAX, sigmaIntWgt[bBin] / rescale2) );
2109  probInt += bArea * rescale2 * pIntNow;
2110  b += bstepNow;
2111  }
2112  nNow = probInt;
2113  }
2114 
2115  // Replace lower or upper limit of k.
2116  if (nNow < nAvg) {
2117  kLow = kNow;
2118  nLow = nNow;
2119  if (stepDir == -1) stepDir = 0;
2120  } else {
2121  kHigh = kNow;
2122  nHigh = nNow;
2123  if (stepDir == 1) stepDir = -1;
2124  }
2125 
2126  // Continue iteration until convergence.
2127  } while (abs(nNow - nAvg) > KCONVERGE * nAvg);
2128 
2129  // Save relevant final numbers for overlap values.
2130  if (bProfile >= 0 && bProfile < 4) {
2131  double avgOverlap = probOverlapInt / probInt;
2132  zeroIntCorr = probOverlapInt / overlapInt;
2133  normOverlap = normPi * zeroIntCorr / avgOverlap;
2134  bAvg = bProbInt / probInt;
2135 
2136  // Values for x-dependent matter profile.
2137  } else if (bProfile == 4) {
2138  // bAvg = Int(b * Pint(b), d2b) / sigmaND.
2139  // zeroIntCorr = Int(<n(b)> * Pint(b), d2b) / sigmaInt.
2140  bAvg = 0.;
2141  zeroIntCorr = 0.;
2142  double b = 0.5 * bstepNow;
2143  for (int bBin = 0; bBin < XDEP_BBIN; bBin++) {
2144  double bArea = 2. * M_PI * b * bstepNow;
2145  double pIntNow = 1. - exp( -min(EXPMAX, sigmaIntWgt[bBin] / rescale2) );
2146  bAvg += sqrt(rescale2) * b * bArea * rescale2 * pIntNow;
2147  zeroIntCorr += bArea * sigmaIntWgt[bBin] * pIntNow;
2148  b += bstepNow;
2149  }
2150  bAvg /= nNow;
2151  zeroIntCorr /= sigmaInt;
2152 
2153  // Other required values.
2154  a0now = kNow;
2155  infoPtr->a0MPI(a0now * XDEP_SMB2FM);
2156  a02now = a0now * a0now;
2157  double xMin = 2. * pTmin / infoPtr->eCM();
2158  a2max = a0now * (XDEP_A1 + a1 * log(1. / xMin));
2159  a2max *= a2max;
2160  }
2161 
2162  // Relative rates for preselection of low-b and high-b region.
2163  // Other useful combinations for subsequent selection.
2164  if (bProfile > 0 && bProfile <= 3) {
2165  probLowB = M_PI * bDiv*bDiv;
2166  double probHighB = M_PI * kNow * overlapHighB;
2167  if (bProfile == 1) probHighB = M_PI * kNow * 0.5 * exp( -bDiv*bDiv);
2168  else if (bProfile == 2) {
2169  fracAhigh = fracA * exp( -bDiv*bDiv);
2170  fracBhigh = fracB * exp( -bDiv*bDiv / radius2B);
2171  fracChigh = fracC * exp( -bDiv*bDiv / radius2C);
2172  fracABChigh = fracAhigh + fracBhigh + fracChigh;
2173  probHighB = M_PI * kNow * 0.5 * fracABChigh;
2174  } else {
2175  cDiv = pow( bDiv, expPow);
2176  cMax = max(2. * expRev, cDiv);
2177  }
2178  probLowB /= (probLowB + probHighB);
2179  }
2180 
2181 }
2182 
2183 //--------------------------------------------------------------------------
2184 
2185 // Pick impact parameter and interaction rate enhancement beforehand,
2186 // i.e. before even the hardest interaction for minimum-bias events.
2187 
2188 void MultipartonInteractions::overlapFirst() {
2189 
2190  // Trivial values if no impact parameter dependence.
2191  if (bProfile <= 0 || bProfile > 4) {
2192  bNow = bAvg;
2193  enhanceB = zeroIntCorr;
2194  bIsSet = true;
2195  isAtLowB = true;
2196  return;
2197  }
2198 
2199  // Preliminary choice between and inside low-b and high-b regions.
2200  double overlapNow = 0.;
2201  double probAccept = 0.;
2202  do {
2203 
2204  // Treatment in low-b region: pick b flat in area.
2205  if (rndmPtr->flat() < probLowB) {
2206  isAtLowB = true;
2207  bNow = bDiv * sqrt(rndmPtr->flat());
2208 
2209  // Evaluate overlap and from that acceptance probability.
2210  if (bProfile == 1) overlapNow = normPi * exp( -bNow*bNow);
2211  else if (bProfile == 2) overlapNow = normPi *
2212  ( fracA * exp( -bNow*bNow)
2213  + fracB * exp( -bNow*bNow / radius2B) / radius2B
2214  + fracC * exp( -bNow*bNow / radius2C) / radius2C );
2215  else overlapNow = normPi * exp( -pow( bNow, expPow));
2216  probAccept = 1. - exp( -min(EXPMAX, M_PI * kNow * overlapNow));
2217 
2218  // Treatment in high-b region: pick b according to overlap.
2219  } else {
2220  isAtLowB = false;
2221 
2222  // For simple and double Gaussian pick b according to exp(-b^2 / r^2).
2223  if (bProfile == 1) {
2224  bNow = sqrt(bDiv*bDiv - log(rndmPtr->flat()));
2225  overlapNow = normPi * exp( -min(EXPMAX, bNow*bNow));
2226  } else if (bProfile == 2) {
2227  double pickFrac = rndmPtr->flat() * fracABChigh;
2228  if (pickFrac < fracAhigh)
2229  bNow = sqrt(bDiv*bDiv - log(rndmPtr->flat()));
2230  else if (pickFrac < fracAhigh + fracBhigh)
2231  bNow = sqrt(bDiv*bDiv - radius2B * log(rndmPtr->flat()));
2232  else bNow = sqrt(bDiv*bDiv - radius2C * log(rndmPtr->flat()));
2233  overlapNow = normPi * ( fracA * exp( -min(EXPMAX, bNow*bNow))
2234  + fracB * exp( -min(EXPMAX, bNow*bNow / radius2B)) / radius2B
2235  + fracC * exp( -min(EXPMAX, bNow*bNow / radius2C)) / radius2C );
2236 
2237  // For exp( - b^expPow) transform to variable c = b^expPow so that
2238  // f(b) = b * exp( - b^expPow) -> f(c) = c^r * exp(-c) with r = expRev.
2239  // case hasLowPow: expPow < 2 <=> r > 0: preselect according to
2240  // f(c) < N exp(-c/2) and then accept with N' * c^r * exp(-c/2).
2241  } else if (hasLowPow) {
2242  double cNow, acceptC;
2243  do {
2244  cNow = cDiv - 2. * log(rndmPtr->flat());
2245  acceptC = pow(cNow / cMax, expRev) * exp( -0.5 * (cNow - cMax));
2246  } while (acceptC < rndmPtr->flat());
2247  bNow = pow( cNow, 1. / expPow);
2248  overlapNow = normPi * exp( -cNow);
2249 
2250  // case !hasLowPow: expPow >= 2 <=> - 1 < r < 0: preselect according to
2251  // f(c) < N exp(-c) and then accept with N' * c^r.
2252  } else {
2253  double cNow, acceptC;
2254  do {
2255  cNow = cDiv - log(rndmPtr->flat());
2256  acceptC = pow(cNow / cDiv, expRev);
2257  } while (acceptC < rndmPtr->flat());
2258  bNow = pow( cNow, 1. / expPow);
2259  overlapNow = normPi * exp( -cNow);
2260  }
2261  double temp = M_PI * kNow *overlapNow;
2262  probAccept = (1. - exp( -min(EXPMAX, temp))) / temp;
2263  }
2264 
2265  // Confirm choice of b value. Derive enhancement factor.
2266  } while (probAccept < rndmPtr->flat());
2267 
2268  // Same enhancement for hardest process and all subsequent MPI
2269  enhanceB = enhanceBmax = enhanceBnow = (normOverlap / normPi) * overlapNow;
2270 
2271  // Done.
2272  bIsSet = true;
2273 
2274 }
2275 
2276 //--------------------------------------------------------------------------
2277 
2278 // Pick impact parameter and interaction rate enhancement afterwards,
2279 // i.e. after a hard interaction is known but before rest of MPI treatment.
2280 
2281 void MultipartonInteractions::overlapNext(Event& event, double pTscale) {
2282 
2283  // Default, valid for bProfile = 0. Also initial Sudakov.
2284  enhanceB = zeroIntCorr;
2285  if (bProfile <= 0 || bProfile > 4) return;
2286  double pT2scale = pTscale*pTscale;
2287 
2288  // Use trial interaction for x-dependent matter profile.
2289  if (bProfile == 4) {
2290  double pTtrial = 0.;
2291  do {
2292  // Pick b according to wanted O(b, x1, x2).
2293  double expb2 = rndmPtr->flat();
2294  double w1 = XDEP_A1 + a1 * log(1. / infoPtr->x1());
2295  double w2 = XDEP_A1 + a1 * log(1. / infoPtr->x2());
2296  double fac = a02now * (w1 * w1 + w2 * w2);
2297  b2now = - fac * log(expb2);
2298  bNow = sqrt(b2now);
2299 
2300  // Enhancement factor for the hard process and overestimate
2301  // for fastPT2. Note that existing framework has a (1. / sigmaND)
2302  // present.
2303  enhanceB = sigmaND / M_PI / fac * expb2;
2304  enhanceBmax = sigmaND / 2. / M_PI / a02now
2305  * exp( -b2now / 2. / a2max );
2306 
2307  // Trial interaction. Keep going until pTtrial < pTscale.
2308  pTtrial = pTnext(pTmax, pTmin, event);
2309  } while (pTtrial > pTscale);
2310  bIsSet = true;
2311  return;
2312  }
2313 
2314  // Begin loop over pT-dependent rejection of b value.
2315  do {
2316 
2317  // Flat enhancement distribution for simple Gaussian.
2318  if (bProfile == 1) {
2319  double expb2 = rndmPtr->flat();
2320  // Same enhancement for hardest process and all subsequent MPI.
2321  enhanceB = enhanceBmax = enhanceBnow = normOverlap * expb2;
2322  bNow = sqrt( -log(expb2));
2323 
2324  // For double Gaussian go the way via b, according to each Gaussian.
2325  } else if (bProfile == 2) {
2326  double bType = rndmPtr->flat();
2327  double b2 = -log( rndmPtr->flat() );
2328  if (bType < fracA) ;
2329  else if (bType < fracA + fracB) b2 *= radius2B;
2330  else b2 *= radius2C;
2331  // Same enhancement for hardest process and all subsequent MPI.
2332  enhanceB = enhanceBmax = enhanceBnow = normOverlap *
2333  ( fracA * exp( -min(EXPMAX, b2))
2334  + fracB * exp( -min(EXPMAX, b2 / radius2B)) / radius2B
2335  + fracC * exp( -min(EXPMAX, b2 / radius2C)) / radius2C );
2336  bNow = sqrt(b2);
2337 
2338  // For exp( - b^expPow) transform to variable c = b^expPow so that
2339  // f(b) = b * exp( - b^expPow) -> f(c) = c^r * exp(-c) with r = expRev.
2340  // case hasLowPow: expPow < 2 <=> r > 0:
2341  // f(c) < r^r exp(-r) for c < 2r, < (2r)^r exp(-r-c/2) for c > 2r.
2342  } else if (bProfile == 3 && hasLowPow) {
2343  double cNow, acceptC;
2344  double probLowC = expRev / (expRev + pow(2., expRev) * exp( - expRev));
2345  do {
2346  if (rndmPtr->flat() < probLowC) {
2347  cNow = 2. * expRev * rndmPtr->flat();
2348  acceptC = pow( cNow / expRev, expRev) * exp(expRev - cNow);
2349  } else {
2350  cNow = 2. * (expRev - log( rndmPtr->flat() ));
2351  acceptC = pow(0.5 * cNow / expRev, expRev)
2352  * exp(expRev - 0.5 * cNow);
2353  }
2354  } while (acceptC < rndmPtr->flat());
2355  // Same enhancement for hardest process and all subsequent MPI.
2356  enhanceB = enhanceBmax = enhanceBnow = normOverlap *exp(-cNow);
2357  bNow = pow( cNow, 1. / expPow);
2358 
2359  // case !hasLowPow: expPow >= 2 <=> - 1 < r < 0:
2360  // f(c) < c^r for c < 1, f(c) < exp(-c) for c > 1.
2361  } else if (bProfile == 3 && !hasLowPow) {
2362  double cNow, acceptC;
2363  double probLowC = expPow / (2. * exp(-1.) + expPow);
2364  do {
2365  if (rndmPtr->flat() < probLowC) {
2366  cNow = pow( rndmPtr->flat(), 0.5 * expPow);
2367  acceptC = exp(-cNow);
2368  } else {
2369  cNow = 1. - log( rndmPtr->flat() );
2370  acceptC = pow( cNow, expRev);
2371  }
2372  } while (acceptC < rndmPtr->flat());
2373  // Same enhancement for hardest process and all subsequent MPI.
2374  enhanceB = enhanceBmax = enhanceBnow = normOverlap * exp(-cNow);
2375  bNow = pow( cNow, 1. / expPow);
2376  }
2377 
2378  // Evaluate "Sudakov form factor" for not having a harder interaction.
2379  } while (sudakov(pT2scale, enhanceB) < rndmPtr->flat());
2380 
2381  // Done.
2382  bIsSet = true;
2383 }
2384 
2385 //--------------------------------------------------------------------------
2386 
2387 // Printe statistics on number of multiparton-interactions processes.
2388 
2389 void MultipartonInteractions::statistics(bool resetStat, ostream& os) {
2390 
2391  // Header.
2392  os << "\n *------- PYTHIA Multiparton Interactions Statistics -----"
2393  << "---*\n"
2394  << " | "
2395  << " |\n"
2396  << " | Note: excludes hardest subprocess if already listed above "
2397  << " |\n"
2398  << " | "
2399  << " |\n"
2400  << " | Subprocess Code | Times"
2401  << " |\n"
2402  << " | | "
2403  << " |\n"
2404  << " |------------------------------------------------------------"
2405  << "-|\n"
2406  << " | | "
2407  << " |\n";
2408 
2409  // Loop over existing processes. Sum of all subprocesses.
2410  int numberSum = 0;
2411  for ( map<int, int>::iterator iter = nGen.begin(); iter != nGen.end();
2412  ++iter) {
2413  int code = iter -> first;
2414  int number = iter->second;
2415  numberSum += number;
2416 
2417  // Find process name that matches code.
2418  string name = " ";
2419  bool foundName = false;
2420  SigmaMultiparton* dSigma;
2421  for (int i = 0; i < 4; ++i) {
2422  if (i == 0) dSigma = &sigma2gg;
2423  else if (i == 1) dSigma = &sigma2qg;
2424  else if (i == 2) dSigma = &sigma2qqbarSame;
2425  else dSigma = &sigma2qq;
2426  int nProc = dSigma->nProc();
2427  for (int iProc = 0; iProc < nProc; ++iProc)
2428  if (dSigma->codeProc(iProc) == code) {
2429  name = dSigma->nameProc(iProc);
2430  foundName = true;
2431  }
2432  if (foundName) break;
2433  }
2434 
2435  // Print individual process info.
2436  os << " | " << left << setw(40) << name << right << setw(5) << code
2437  << " | " << setw(11) << number << " |\n";
2438  }
2439 
2440  // Print summed process info.
2441  os << " | "
2442  << " |\n"
2443  << " | " << left << setw(45) << "sum" << right << " | " << setw(11)
2444  << numberSum << " |\n";
2445 
2446  // Listing finished.
2447  os << " | | "
2448  << " |\n"
2449  << " *------- End PYTHIA Multiparton Interactions Statistics ----"
2450  << "-*" << endl;
2451 
2452  // Optionally reset statistics contents.
2453  if (resetStat) resetStatistics();
2454 
2455 }
2456 
2457 //==========================================================================
2458 
2459 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5