using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...

[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / xmldoc / MultipleInteractions.xml

<chapter name="Multiple Interactions">

<h2>Multiple Interactions</h2>

The starting point for the multiple interactions physics scenario in
PYTHIA is provided by <ref>Sjo87</ref>. Recent developments have 
included a more careful study of flavour and colour correlations, 
junction topologies and the relationship to beam remnants 
<ref>Sjo04</ref>, and interleaving with initial-state radiation 
<ref>Sjo05</ref>, making use of transverse-momentum-ordered
initial- and final-state showers. 

<p/>
A big unsolved issue is how the colour of all these subsystems is 
correlated. For sure there is a correlation coming from the colour
singlet nature of the incoming beams, but in addition final-state 
colour rearrangements may change the picture. Indeed such extra
effects appear necessary to describe data, e.g. on 
<ei>&lt;pT&gt;(n_ch)</ei>. A simple implementation of colour 
rearrangement is found as part of the
<aloc href="BeamRemnants">beam remnants</aloc> description.

<h3>Main variables</h3>

The maximum <ei>pT</ei> to be allowed for multiple interactions is
related to the nature of the hard process itself. It involves a
delicate balance between not doublecounting and not leaving any
gaps in the coverage. The best procedure may depend on information 
only the user has: how the events were generated and mixed (e.g. with 
Les Houches Accord external input), and how they are intended to be 
used. Therefore a few options are available, with a sensible default 
behaviour.
<modepick name="MultipleInteractions:pTmaxMatch" default="0" min="0" max="2">
Way in which the maximum scale for multiple interactions is set
to match the scale of the hard process itself.
<option value="0"><b>(i)</b> if the final state of the hard process 
(not counting subsequent resonance decays) contains only quarks 
(<ei>u, d, s, c ,b</ei>), gluons and photons then <ei>pT_max</ei> 
is chosen to be the factorization scale for internal processes 
and the <code>scale</code> value for Les Houches input; 
<b>(ii)</b> if not, interactions are allowed to go all the way up 
to the kinematical limit. 
The reasoning is that the former kind of processes are generated by
the multiple-interactions machinery and so would doublecount hard
processes if allowed to overlap the same <ei>pT</ei> range, 
while no such danger exists in the latter case.
</option>
<option value="1">always use the factorization scale for an internal
process and the <code>scale</code> value for Les Houches input, 
i.e. the lower value. This should avoid doublecounting, but
may leave out some interactions that ought to have been simulated.
</option>
<option value="2">always allow multiple interactions up to the 
kinematical limit. This will simulate all possible event topologies, 
but may lead to doublecounting.
</option>
</modepick>

<p/>
The rate of interactions is determined by 
<parm name="MultipleInteractions:alphaSvalue" default="0.127" min="0.06" max="0.25">
The value of <ei>alpha_strong</ei> at <ei>m_Z</ei>. Default value is 
picked equal to the one used in CTEQ 5L. 
</parm> 

<p/>
The actual value is then regulated by the running to the scale 
<ei>pT^2</ei>, at which it is evaluated
<modepick name="MultipleInteractions:alphaSorder" default="1" min="0" max="2">
The order at which <ei>alpha_strong</ei> runs at scales away from 
<ei>m_Z</ei>.
<option value="0">zeroth order, i.e. <ei>alpha_strong</ei> is kept 
fixed.</option>
<option value="1">first order, which is the normal value.</option>
<option value="2">second order. Since other parts of the code do 
not go to second order there is no strong reason to use this option, 
but there is also nothing wrong with it.</option>
</modepick>

<p/>
QED interactions are regulated by the <ei>alpha_electromagnetic</ei>
value at the <ei>pT^2</ei> scale of an interaction.
 
<modepick name="MultipleInteractions:alphaEMorder" default="1" min="-1" max="1">
The running of <ei>alpha_em</ei> used in hard processes.
<option value="1">first-order running, constrained to agree with
<code>StandardModel:alphaEMmZ</code> at the <ei>Z^0</ei> mass.
</option>
<option value="0">zeroth order, i.e. <ei>alpha_em</ei> is kept 
fixed at its value at vanishing momentum transfer.</option>
<option value="-1">zeroth order, i.e. <ei>alpha_em</ei> is kept 
fixed, but at <code>StandardModel:alphaEMmZ</code>, i.e. its value
at the <ei>Z^0</ei> mass.
</option> 
</modepick>

<p/>
Note that the choices of <ei>alpha_strong</ei> and <ei>alpha_em</ei> 
made here override the ones implemented in the normal process machinery, 
but only for the interactions generated by the 
<code>MultipleInteractions</code> class.

<p/>
In addition there is the possibility of a global rescaling of 
cross sections (which could not easily be accommodated by a 
changed <ei>alpha_strong</ei>, since <ei>alpha_strong</ei> runs)
<parm name="MultipleInteractions:Kfactor" default="1.0" min="0.5" max="4.0">
Multiply all cross sections by this fix factor.
</parm>

<p/>
There are two complementary ways of regularizing the small-<ei>pT</ei> 
divergence, a sharp cutoff and a smooth dampening. These can be 
combined as desired, but it makes sense to coordinate with how the 
same issue is handled in <aloc href="SpacelikeShowers">spacelike 
showers</aloc>. Actually, by default, the parameters defined here are 
used also for the spacelike showers, but this can be overridden.

<p/>
Regularization of the divergence of the QCD cross section for 
<ei>pT -> 0</ei> is obtained by a factor <ei>pT^4 / (pT0^2 + pT^2)^2</ei>, 
and by using an <ei>alpha_s(pT0^2 + pT^2)</ei>. An energy dependence 
of the <ei>pT0</ei> choice is introduced by two further parameters, 
so that <ei>pT0Ref</ei> is the <ei>pT0</ei> value for the reference 
cm energy, <ei>pT0Ref = pT0(ecmRef)</ei>.   
<note>Warning:</note> if a large <ei>pT0</ei> is picked for multiple 
interactions, such that the integrated interaction cross section is 
below the nondiffractive inelastic one, this <ei>pT0</ei> will 
automatically be scaled down to cope.

<p/>
The actual pT0 parameter used at a given cm energy scale, <ei>ecmNow</ei>,
is obtained as
<eq>
     pT0 = pT0(ecmNow) = pT0Ref * (ecmNow / ecmRef)^ecmPow 
</eq>
where <ei>pT0Ref</ei>, <ei>ecmRef</ei> and <ei>ecmPow</ei> are the 
three parameters below.

<parm name="MultipleInteractions:pT0Ref" default="2.15" min="0.5" max="10.0">
The <ei>pT0Ref</ei> scale in the above formula.
<note>Note:</note> <ei>pT0Ref</ei> is one of the key parameters in a
complete PYTHIA tune. Its value is intimately tied to a number of other
choices, such as that of colour flow description, so unfortunately it is
difficult to give an independent meaning to <ei>pT0Ref</ei>.
</parm>

<parm name="MultipleInteractions:ecmRef" default="1800.0" min="1.">
The <ei>ecmRef</ei> reference energy scale introduced above.
</parm>

<parm name="MultipleInteractions:ecmPow" default="0.16" min="0.0" max="0.5">
The <ei>ecmPow</ei> energy rescaling pace introduced above.
</parm>

<p/>
Alternatively, or in combination, a sharp cut can be used.
<parm name="MultipleInteractions:pTmin" default="0.2" min="0.1" max="10.0">
Lower cutoff in <ei>pT</ei>, below which no further interactions 
are allowed. Normally <ei>pT0</ei> above would be used to provide 
the main regularization of the cross section for <ei>pT -> 0</ei>, 
in which case <ei>pTmin</ei> is used  mainly for technical reasons. 
It is possible, however, to set <ei>pT0Ref = 0</ei> and use 
<ei>pTmin</ei> to provide a step-function regularization, or to 
combine them in intermediate approaches. Currently <ei>pTmin</ei> 
is taken to be energy-independent.  
</parm>

<p/> 
The choice of impact-parameter dependence is regulated by several
parameters.

<modepick name="MultipleInteractions:bProfile" default="2" 
min="0" max="3">
Choice of impact parameter profile for the incoming hadron beams.
<option value="0">no impact parameter dependence at all.</option>
<option value="1">a simple Gaussian matter distribution; 
no free parameters.</option>
<option value="2">a double Gaussian matter distribution, 
with the two free parameters <ei>coreRadius</ei> and 
<ei>coreFraction</ei>.</option>
<option value="3">an overlap function, i.e. the convolution of 
the matter distributions of the two incoming hadrons, of the form
<ei>exp(- b^expPow)</ei>, where <ei>expPow</ei> is a free 
parameter.</option> 
</modepick>

<parm name="MultipleInteractions:coreRadius" default="0.4" min="0.1" max="1.">
When assuming a double Gaussian matter profile, <ei>bProfile = 2</ei>,
the inner core is assumed to have a radius that is a factor
<ei>coreRadius</ei> smaller than the rest.
</parm> 

<parm name="MultipleInteractions:coreFraction" default="0.5" min="0." max="1."> 
When assuming a double Gaussian matter profile, <ei>bProfile = 2</ei>,
the inner core is assumed to have a fraction <ei>coreFraction</ei> 
of the matter content of the hadron.
</parm> 

<parm name="MultipleInteractions:expPow" default="1." min="0.4" max="10.">
When <ei>bProfile = 3</ei> it gives the power of the assumed overlap 
shape <ei>exp(- b^expPow)</ei>. Default corresponds to a simple 
exponential drop, which is not too dissimilar from the overlap 
obtained with the standard double Gaussian parameters. For 
<ei>expPow = 2</ei> we reduce to the simple Gaussian, <ei>bProfile = 1</ei>, 
and for <ei>expPow -> infinity</ei> to no impact parameter dependence 
at all, <ei>bProfile = 0</ei>. For small <ei>expPow</ei> the program 
becomes slow and unstable, so the min limit must be respected.
</parm> 

<p/> 
It is possible to regulate the set of processes that are included in the
multiple-interactions framework.

<modepick name="MultipleInteractions:processLevel" default="3" 
min="0" max="3">
Set of processes included in the machinery.
<option value="0">only the simplest <ei>2 -> 2</ei> QCD processes
between quarks and gluons, giving no new flavours, i.e. dominated by
<ei>t</ei>-channel gluon exchange.</option>
<option value="1">also <ei>2 -> 2</ei> QCD processes giving new flavours
(including charm and bottom), i.e. proceeding through <ei>s</ei>-channel 
gluon exchange.</option>
<option value="2">also <ei>2 -> 2</ei> processes involving one or two
photons in the final state, <ei>s</ei>-channel <ei>gamma</ei>
boson exchange and <ei>t</ei>-channel <ei>gamma/Z^0/W^+-</ei>
boson exchange.</option>
<option value="3">also charmonium and bottomonium production, via
colour singlet and colour octet channels.</option> 
</modepick>

<h3>Further variables</h3>

These should normally not be touched. Their only function is for
cross-checks.

<modeopen name="MultipleInteractions:nQuarkIn" default="5" min="0" max="5">
Number of allowed incoming quark flavours in the beams; a change 
to 4 would thus exclude <ei>b</ei> and <ei>bbar</ei> as incoming 
partons, etc.
</modeopen>

<modeopen name="MultipleInteractions:nSample" default="1000" min="100"> 
The allowed <ei>pT</ei> range is split (unevenly) into 100 bins, 
and in each of these the interaction cross section is evaluated in 
<ei>nSample</ei> random phase space points. The full integral is used 
at initialization, and the differential one during the run as a
"Sudakov form factor" for the choice of the hardest interaction.
A larger number implies increased accuracy of the calculations.
</modeopen>

<h3>The process library</h3>

The processes used to generate multiple interactions form a subset
of the standard library of hard processes. The input is slightly
different from the standard hard-process machinery, however, 
since incoming flavours, the <ei>alpha_strong</ei> value and most
of the kinematics are aready fixed when the process is called.

<h3>Technical notes</h3>

Relative to the articles mentioned above, not much has happened.
The main news is a technical one, that the phase space of the 
<ei>2 -> 2</ei> (massless) QCD processes is now sampled in 
<ei>dy_3 dy_4 dpT^2</ei>, where <ei>y_3</ei> and <ei>y_4</ei> are 
the rapidities of the two produced partons. One can show that
<eq>
    (dx_1 / x_1) * (dx_2 / x_2) * d(tHat) = dy_3 * dy_4 * dpT^2
</eq>
Furthermore, since cross sections are dominated by the "Rutherford"
one of <ei>t</ei>-channel gluon exchange, which is enhanced by a 
factor of 9/4 for each incoming gluon, effective structure functions 
are defined as
<eq>
    F(x, pT2) = (9/4) * xg(x, pT2) + sum_i xq_i(x, pT2) 
</eq>
With this technical shift of factors 9/4 from cross sections to parton 
densities, a common upper estimate of 
<eq>
    d(sigmaHat)/d(pT2) &lt; pi * alpha_strong^2 / pT^4   
</eq>
is obtained. 

<p/>
In fact this estimate can be reduced by a factor of 1/2 for the 
following reason: for any configuration <ei>(y_3, y_4, pT2)</ei> also 
one with <ei>(y_4, y_3, pT2)</ei> lies in the phase space. Not both 
of those can enjoy being enhanced by the <ei>tHat -> 0</ei> 
singularity of 
<eq>
    d(sigmaHat) propto 1/tHat^2. 
</eq>
Or if they are, which is possible with identical partons like 
<ei>q q -> q q</ei> and <ei>g g -> g g</ei>, each singularity comes 
with half the strength. So, when integrating/averaging over the two 
configurations, the estimated <ei>d(sigmaHat)/d(pT2)</ei> drops. 
Actually, it drops even further, since the naive estimate above is 
based on
<eq>
    (4 /9) * (1 + (uHat/sHat)^2) &lt; 8/9 &lt; 1
</eq>
The 8/9 value would be approached for <ei>tHat -> 0</ei>, which 
implies <ei>sHat >> pT2</ei> and thus a heavy parton-distribution 
penalty, while parton distributions are largest for 
<ei>tHat = uHat = -sHat/2</ei>, where the above expression 
evaluates to 5/9. A fudge factor is therefore introduced to go the 
final step, so it can easily be modifed when further non-Rutherford 
processes are added, or should parton distributions change significantly.

<p/>
At initialization, it is assumed that  
<eq>
    d(sigma)/d(pT2) &lt; d(sigmaHat)/d(pT2) * F(x_T, pT2) * F(x_T, pT2)
       * (2 y_max(pT))^2
</eq>
where the first factor is the upper estimate as above, the second two
the parton density sum evaluated at <ei>y_3 = y_ 4 = 0</ei> so that 
<ei>x_1 = x_2 = x_T = 2 pT / E_cm</ei>, where the product is expected 
to be maximal, and the final is the phase space for
<ei>-y_max &lt; y_{3,4} &lt; y_max</ei>.
The right-hand side expression is scanned logarithmically in <ei>y</ei>, 
and a <ei>N</ei> is determined such that it always is below 
<ei>N/pT^4</ei>.

<p/>
To describe the dampening of the cross section at <ei>pT -> 0</ei> by
colour screening, the actual cross section is multiplied by a 
regularization factor <ei>(pT^2 / (pT^2 + pT0^2))^2</ei>, and the 
<ei>alpha_s</ei> is evaluated at a scale <ei>pT^2 + pT0^2</ei>, 
where <ei>pT0</ei> is a free parameter of the order of 2 - 4 GeV. 
Since <ei>pT0</ei> can be energy-dependent,  an ansatz
<eq>
    pT0(ecm) = pT0Ref * (ecm/ecmRef)^ecmPow
</eq>
is used, where <ei>ecm</ei> is the current cm frame energy, 
<ei>ecmRef</ei> is an arbitrary reference energy where <ei>pT0Ref</ei> 
is defined, and <ei>ecmPow</ei> gives the energy rescaling pace. For 
technical reasons, also an absolute lower <ei>pT</ei> scale <ei>pTmin</ei>, 
by default 0.2 GeV, is introduced. In principle, it is possible to 
recover older scenarios with a sharp <ei>pT</ei> cutoff by setting 
<ei>pT0 = 0</ei> and letting <ei>pTmin</ei> be a larger number. 

<p/>
The above scanning strategy is then slightly modified: instead of
an upper estimate <ei>N/pT^4</ei> one of the form 
<ei>N/(pT^2 + r * pT0^2)^2</ei> is used. At first glance, <ei>r = 1</ei> 
would seem to be fixed by the form of the regularization procedure, 
but this does not take into account the nontrivial dependence on 
<ei>alpha_s</ei>, parton distributions and phase space. A better 
Monte Carlo efficiency is obtained for <ei>r</ei> somewhat below unity, 
and currently <ei>r = 0.25</ei> is hardcoded.

In the generation a trial <ei>pT2</ei> is then selected according to
<eq>
    d(Prob)/d(pT2) = (1/sigma_ND) * N/(pT^2 + r * pT0^2)^2 * ("Sudakov")
</eq>
For the trial <ei>pT2</ei>, a <ei>y_3</ei> and a <ei>y_4</ei> are then 
selected, and incoming flavours according to the respective 
<ei>F(x_i, pT2)</ei>, and then the cross section is evaluated for this 
flavour combination. The ratio of trial/upper estimate gives the 
probability of survival.

<p/>
Actually, to profit from the factor 1/2 mentioned above, the cross
section for the combination with <ei>y_3</ei> and <ei>y_4</ei> 
interchanged is also tried, which corresponds to exchanging <ei>tHat</ei>
and <ei>uHat</ei>, and the average formed, while the final kinematics 
is given by the relative importance of the two.

<p/>
Furthermore, since large <ei>y</ei> values are disfavoured by dropping 
PDF's, a factor 
<eq>
   WT_y = (1 - (y_3/y_max)^2) * (1 - (y_4/y_max)^2) 
</eq>
is evaluated, and used as a survival probability before the more
time-consuming PDF+ME evaluation, with surviving events given a 
compensating weight <ei>1/WT_y</ei>. 

<p/>
An impact-parameter dependencs is also allowed. Based on the hard 
<ei>pT</ei> scale of the first interaction, and enhancement/depletion 
factor is picked, which multiplies the rate of subsequent interactions.

<p/>
Parton densities are rescaled and modified to take into account the 
energy-momentum and flavours kicked out by already-considered 
interactions.

</chapter>

<!-- Copyright (C) 2008 Torbjorn Sjostrand -->