1 <chapter name="Multiple Interactions">
3 <h2>Multiple Interactions</h2>
5 The starting point for the multiple interactions physics scenario in
6 PYTHIA is provided by <ref>Sjo87</ref>. Recent developments have
7 included a more careful study of flavour and colour correlations,
8 junction topologies and the relationship to beam remnants
9 <ref>Sjo04</ref>, and interleaving with initial-state radiation
10 <ref>Sjo05</ref>, making use of transverse-momentum-ordered
11 initial- and final-state showers.
14 A big unsolved issue is how the colour of all these subsystems is
15 correlated. For sure there is a correlation coming from the colour
16 singlet nature of the incoming beams, but in addition final-state
17 colour rearrangements may change the picture. Indeed such extra
18 effects appear necessary to describe data, e.g. on
19 <ei><pT>(n_ch)</ei>. A simple implementation of colour
20 rearrangement is found as part of the
21 <aloc href="BeamRemnants">beam remnants</aloc> description.
23 <h3>Main variables</h3>
25 The maximum <ei>pT</ei> to be allowed for multiple interactions is
26 related to the nature of the hard process itself. It involves a
27 delicate balance between not doublecounting and not leaving any
28 gaps in the coverage. The best procedure may depend on information
29 only the user has: how the events were generated and mixed (e.g. with
30 Les Houches Accord external input), and how they are intended to be
31 used. Therefore a few options are available, with a sensible default
33 <modepick name="MultipleInteractions:pTmaxMatch" default="0" min="0" max="2">
34 Way in which the maximum scale for multiple interactions is set
35 to match the scale of the hard process itself.
36 <option value="0"><b>(i)</b> if the final state of the hard process
37 (not counting subsequent resonance decays) contains only quarks
38 (<ei>u, d, s, c ,b</ei>), gluons and photons then <ei>pT_max</ei>
39 is chosen to be the factorization scale for internal processes
40 and the <code>scale</code> value for Les Houches input;
41 <b>(ii)</b> if not, interactions are allowed to go all the way up
42 to the kinematical limit.
43 The reasoning is that the former kind of processes are generated by
44 the multiple-interactions machinery and so would doublecount hard
45 processes if allowed to overlap the same <ei>pT</ei> range,
46 while no such danger exists in the latter case.
48 <option value="1">always use the factorization scale for an internal
49 process and the <code>scale</code> value for Les Houches input,
50 i.e. the lower value. This should avoid doublecounting, but
51 may leave out some interactions that ought to have been simulated.
53 <option value="2">always allow multiple interactions up to the
54 kinematical limit. This will simulate all possible event topologies,
55 but may lead to doublecounting.
60 The rate of interactions is determined by
61 <parm name="MultipleInteractions:alphaSvalue" default="0.127" min="0.06" max="0.25">
62 The value of <ei>alpha_strong</ei> at <ei>m_Z</ei>. Default value is
63 picked equal to the one used in CTEQ 5L.
67 The actual value is then regulated by the running to the scale
68 <ei>pT^2</ei>, at which it is evaluated
69 <modepick name="MultipleInteractions:alphaSorder" default="1" min="0" max="2">
70 The order at which <ei>alpha_strong</ei> runs at scales away from
72 <option value="0">zeroth order, i.e. <ei>alpha_strong</ei> is kept
74 <option value="1">first order, which is the normal value.</option>
75 <option value="2">second order. Since other parts of the code do
76 not go to second order there is no strong reason to use this option,
77 but there is also nothing wrong with it.</option>
81 QED interactions are regulated by the <ei>alpha_electromagnetic</ei>
82 value at the <ei>pT^2</ei> scale of an interaction.
84 <modepick name="MultipleInteractions:alphaEMorder" default="1" min="-1" max="1">
85 The running of <ei>alpha_em</ei> used in hard processes.
86 <option value="1">first-order running, constrained to agree with
87 <code>StandardModel:alphaEMmZ</code> at the <ei>Z^0</ei> mass.
89 <option value="0">zeroth order, i.e. <ei>alpha_em</ei> is kept
90 fixed at its value at vanishing momentum transfer.</option>
91 <option value="-1">zeroth order, i.e. <ei>alpha_em</ei> is kept
92 fixed, but at <code>StandardModel:alphaEMmZ</code>, i.e. its value
93 at the <ei>Z^0</ei> mass.
98 Note that the choices of <ei>alpha_strong</ei> and <ei>alpha_em</ei>
99 made here override the ones implemented in the normal process machinery,
100 but only for the interactions generated by the
101 <code>MultipleInteractions</code> class.
104 In addition there is the possibility of a global rescaling of
105 cross sections (which could not easily be accommodated by a
106 changed <ei>alpha_strong</ei>, since <ei>alpha_strong</ei> runs)
107 <parm name="MultipleInteractions:Kfactor" default="1.0" min="0.5" max="4.0">
108 Multiply all cross sections by this fix factor.
112 There are two complementary ways of regularizing the small-<ei>pT</ei>
113 divergence, a sharp cutoff and a smooth dampening. These can be
114 combined as desired, but it makes sense to coordinate with how the
115 same issue is handled in <aloc href="SpacelikeShowers">spacelike
116 showers</aloc>. Actually, by default, the parameters defined here are
117 used also for the spacelike showers, but this can be overridden.
120 Regularization of the divergence of the QCD cross section for
121 <ei>pT -> 0</ei> is obtained by a factor <ei>pT^4 / (pT0^2 + pT^2)^2</ei>,
122 and by using an <ei>alpha_s(pT0^2 + pT^2)</ei>. An energy dependence
123 of the <ei>pT0</ei> choice is introduced by two further parameters,
124 so that <ei>pT0Ref</ei> is the <ei>pT0</ei> value for the reference
125 cm energy, <ei>pT0Ref = pT0(ecmRef)</ei>.
126 <note>Warning:</note> if a large <ei>pT0</ei> is picked for multiple
127 interactions, such that the integrated interaction cross section is
128 below the nondiffractive inelastic one, this <ei>pT0</ei> will
129 automatically be scaled down to cope.
132 The actual pT0 parameter used at a given cm energy scale, <ei>ecmNow</ei>,
135 pT0 = pT0(ecmNow) = pT0Ref * (ecmNow / ecmRef)^ecmPow
137 where <ei>pT0Ref</ei>, <ei>ecmRef</ei> and <ei>ecmPow</ei> are the
138 three parameters below.
140 <parm name="MultipleInteractions:pT0Ref" default="2.15" min="0.5" max="10.0">
141 The <ei>pT0Ref</ei> scale in the above formula.
142 <note>Note:</note> <ei>pT0Ref</ei> is one of the key parameters in a
143 complete PYTHIA tune. Its value is intimately tied to a number of other
144 choices, such as that of colour flow description, so unfortunately it is
145 difficult to give an independent meaning to <ei>pT0Ref</ei>.
148 <parm name="MultipleInteractions:ecmRef" default="1800.0" min="1.">
149 The <ei>ecmRef</ei> reference energy scale introduced above.
152 <parm name="MultipleInteractions:ecmPow" default="0.16" min="0.0" max="0.5">
153 The <ei>ecmPow</ei> energy rescaling pace introduced above.
157 Alternatively, or in combination, a sharp cut can be used.
158 <parm name="MultipleInteractions:pTmin" default="0.2" min="0.1" max="10.0">
159 Lower cutoff in <ei>pT</ei>, below which no further interactions
160 are allowed. Normally <ei>pT0</ei> above would be used to provide
161 the main regularization of the cross section for <ei>pT -> 0</ei>,
162 in which case <ei>pTmin</ei> is used mainly for technical reasons.
163 It is possible, however, to set <ei>pT0Ref = 0</ei> and use
164 <ei>pTmin</ei> to provide a step-function regularization, or to
165 combine them in intermediate approaches. Currently <ei>pTmin</ei>
166 is taken to be energy-independent.
170 The choice of impact-parameter dependence is regulated by several
173 <modepick name="MultipleInteractions:bProfile" default="2"
175 Choice of impact parameter profile for the incoming hadron beams.
176 <option value="0">no impact parameter dependence at all.</option>
177 <option value="1">a simple Gaussian matter distribution;
178 no free parameters.</option>
179 <option value="2">a double Gaussian matter distribution,
180 with the two free parameters <ei>coreRadius</ei> and
181 <ei>coreFraction</ei>.</option>
182 <option value="3">an overlap function, i.e. the convolution of
183 the matter distributions of the two incoming hadrons, of the form
184 <ei>exp(- b^expPow)</ei>, where <ei>expPow</ei> is a free
188 <parm name="MultipleInteractions:coreRadius" default="0.4" min="0.1" max="1.">
189 When assuming a double Gaussian matter profile, <ei>bProfile = 2</ei>,
190 the inner core is assumed to have a radius that is a factor
191 <ei>coreRadius</ei> smaller than the rest.
194 <parm name="MultipleInteractions:coreFraction" default="0.5" min="0." max="1.">
195 When assuming a double Gaussian matter profile, <ei>bProfile = 2</ei>,
196 the inner core is assumed to have a fraction <ei>coreFraction</ei>
197 of the matter content of the hadron.
200 <parm name="MultipleInteractions:expPow" default="1." min="0.4" max="10.">
201 When <ei>bProfile = 3</ei> it gives the power of the assumed overlap
202 shape <ei>exp(- b^expPow)</ei>. Default corresponds to a simple
203 exponential drop, which is not too dissimilar from the overlap
204 obtained with the standard double Gaussian parameters. For
205 <ei>expPow = 2</ei> we reduce to the simple Gaussian, <ei>bProfile = 1</ei>,
206 and for <ei>expPow -> infinity</ei> to no impact parameter dependence
207 at all, <ei>bProfile = 0</ei>. For small <ei>expPow</ei> the program
208 becomes slow and unstable, so the min limit must be respected.
212 It is possible to regulate the set of processes that are included in the
213 multiple-interactions framework.
215 <modepick name="MultipleInteractions:processLevel" default="3"
217 Set of processes included in the machinery.
218 <option value="0">only the simplest <ei>2 -> 2</ei> QCD processes
219 between quarks and gluons, giving no new flavours, i.e. dominated by
220 <ei>t</ei>-channel gluon exchange.</option>
221 <option value="1">also <ei>2 -> 2</ei> QCD processes giving new flavours
222 (including charm and bottom), i.e. proceeding through <ei>s</ei>-channel
223 gluon exchange.</option>
224 <option value="2">also <ei>2 -> 2</ei> processes involving one or two
225 photons in the final state, <ei>s</ei>-channel <ei>gamma</ei>
226 boson exchange and <ei>t</ei>-channel <ei>gamma/Z^0/W^+-</ei>
227 boson exchange.</option>
228 <option value="3">also charmonium and bottomonium production, via
229 colour singlet and colour octet channels.</option>
232 <h3>Further variables</h3>
234 These should normally not be touched. Their only function is for
237 <modeopen name="MultipleInteractions:nQuarkIn" default="5" min="0" max="5">
238 Number of allowed incoming quark flavours in the beams; a change
239 to 4 would thus exclude <ei>b</ei> and <ei>bbar</ei> as incoming
243 <modeopen name="MultipleInteractions:nSample" default="1000" min="100">
244 The allowed <ei>pT</ei> range is split (unevenly) into 100 bins,
245 and in each of these the interaction cross section is evaluated in
246 <ei>nSample</ei> random phase space points. The full integral is used
247 at initialization, and the differential one during the run as a
248 "Sudakov form factor" for the choice of the hardest interaction.
249 A larger number implies increased accuracy of the calculations.
252 <h3>The process library</h3>
254 The processes used to generate multiple interactions form a subset
255 of the standard library of hard processes. The input is slightly
256 different from the standard hard-process machinery, however,
257 since incoming flavours, the <ei>alpha_strong</ei> value and most
258 of the kinematics are aready fixed when the process is called.
260 <h3>Technical notes</h3>
262 Relative to the articles mentioned above, not much has happened.
263 The main news is a technical one, that the phase space of the
264 <ei>2 -> 2</ei> (massless) QCD processes is now sampled in
265 <ei>dy_3 dy_4 dpT^2</ei>, where <ei>y_3</ei> and <ei>y_4</ei> are
266 the rapidities of the two produced partons. One can show that
268 (dx_1 / x_1) * (dx_2 / x_2) * d(tHat) = dy_3 * dy_4 * dpT^2
270 Furthermore, since cross sections are dominated by the "Rutherford"
271 one of <ei>t</ei>-channel gluon exchange, which is enhanced by a
272 factor of 9/4 for each incoming gluon, effective structure functions
275 F(x, pT2) = (9/4) * xg(x, pT2) + sum_i xq_i(x, pT2)
277 With this technical shift of factors 9/4 from cross sections to parton
278 densities, a common upper estimate of
280 d(sigmaHat)/d(pT2) < pi * alpha_strong^2 / pT^4
285 In fact this estimate can be reduced by a factor of 1/2 for the
286 following reason: for any configuration <ei>(y_3, y_4, pT2)</ei> also
287 one with <ei>(y_4, y_3, pT2)</ei> lies in the phase space. Not both
288 of those can enjoy being enhanced by the <ei>tHat -> 0</ei>
291 d(sigmaHat) propto 1/tHat^2.
293 Or if they are, which is possible with identical partons like
294 <ei>q q -> q q</ei> and <ei>g g -> g g</ei>, each singularity comes
295 with half the strength. So, when integrating/averaging over the two
296 configurations, the estimated <ei>d(sigmaHat)/d(pT2)</ei> drops.
297 Actually, it drops even further, since the naive estimate above is
300 (4 /9) * (1 + (uHat/sHat)^2) < 8/9 < 1
302 The 8/9 value would be approached for <ei>tHat -> 0</ei>, which
303 implies <ei>sHat >> pT2</ei> and thus a heavy parton-distribution
304 penalty, while parton distributions are largest for
305 <ei>tHat = uHat = -sHat/2</ei>, where the above expression
306 evaluates to 5/9. A fudge factor is therefore introduced to go the
307 final step, so it can easily be modifed when further non-Rutherford
308 processes are added, or should parton distributions change significantly.
311 At initialization, it is assumed that
313 d(sigma)/d(pT2) < d(sigmaHat)/d(pT2) * F(x_T, pT2) * F(x_T, pT2)
316 where the first factor is the upper estimate as above, the second two
317 the parton density sum evaluated at <ei>y_3 = y_ 4 = 0</ei> so that
318 <ei>x_1 = x_2 = x_T = 2 pT / E_cm</ei>, where the product is expected
319 to be maximal, and the final is the phase space for
320 <ei>-y_max < y_{3,4} < y_max</ei>.
321 The right-hand side expression is scanned logarithmically in <ei>y</ei>,
322 and a <ei>N</ei> is determined such that it always is below
326 To describe the dampening of the cross section at <ei>pT -> 0</ei> by
327 colour screening, the actual cross section is multiplied by a
328 regularization factor <ei>(pT^2 / (pT^2 + pT0^2))^2</ei>, and the
329 <ei>alpha_s</ei> is evaluated at a scale <ei>pT^2 + pT0^2</ei>,
330 where <ei>pT0</ei> is a free parameter of the order of 2 - 4 GeV.
331 Since <ei>pT0</ei> can be energy-dependent, an ansatz
333 pT0(ecm) = pT0Ref * (ecm/ecmRef)^ecmPow
335 is used, where <ei>ecm</ei> is the current cm frame energy,
336 <ei>ecmRef</ei> is an arbitrary reference energy where <ei>pT0Ref</ei>
337 is defined, and <ei>ecmPow</ei> gives the energy rescaling pace. For
338 technical reasons, also an absolute lower <ei>pT</ei> scale <ei>pTmin</ei>,
339 by default 0.2 GeV, is introduced. In principle, it is possible to
340 recover older scenarios with a sharp <ei>pT</ei> cutoff by setting
341 <ei>pT0 = 0</ei> and letting <ei>pTmin</ei> be a larger number.
344 The above scanning strategy is then slightly modified: instead of
345 an upper estimate <ei>N/pT^4</ei> one of the form
346 <ei>N/(pT^2 + r * pT0^2)^2</ei> is used. At first glance, <ei>r = 1</ei>
347 would seem to be fixed by the form of the regularization procedure,
348 but this does not take into account the nontrivial dependence on
349 <ei>alpha_s</ei>, parton distributions and phase space. A better
350 Monte Carlo efficiency is obtained for <ei>r</ei> somewhat below unity,
351 and currently <ei>r = 0.25</ei> is hardcoded.
353 In the generation a trial <ei>pT2</ei> is then selected according to
355 d(Prob)/d(pT2) = (1/sigma_ND) * N/(pT^2 + r * pT0^2)^2 * ("Sudakov")
357 For the trial <ei>pT2</ei>, a <ei>y_3</ei> and a <ei>y_4</ei> are then
358 selected, and incoming flavours according to the respective
359 <ei>F(x_i, pT2)</ei>, and then the cross section is evaluated for this
360 flavour combination. The ratio of trial/upper estimate gives the
361 probability of survival.
364 Actually, to profit from the factor 1/2 mentioned above, the cross
365 section for the combination with <ei>y_3</ei> and <ei>y_4</ei>
366 interchanged is also tried, which corresponds to exchanging <ei>tHat</ei>
367 and <ei>uHat</ei>, and the average formed, while the final kinematics
368 is given by the relative importance of the two.
371 Furthermore, since large <ei>y</ei> values are disfavoured by dropping
374 WT_y = (1 - (y_3/y_max)^2) * (1 - (y_4/y_max)^2)
376 is evaluated, and used as a survival probability before the more
377 time-consuming PDF+ME evaluation, with surviving events given a
378 compensating weight <ei>1/WT_y</ei>.
381 An impact-parameter dependencs is also allowed. Based on the hard
382 <ei>pT</ei> scale of the first interaction, and enhancement/depletion
383 factor is picked, which multiplies the rate of subsequent interactions.
386 Parton densities are rescaled and modified to take into account the
387 energy-momentum and flavours kicked out by already-considered
392 <!-- Copyright (C) 2008 Torbjorn Sjostrand -->