How to Run with Varying Energies

It is possible to use PYTHIA in a mode where the energy can be varied from one event to the next, without the need to re-initialize with a new PYINIT call. This allows a significant speed-up of execution, although it is not as fast as running at a fixed energy. It can not be used for everything -- we will come to the fine print at the end -- but it should be applicable for most tasks.

The master switch to access this possibility is in MSTP(171). By default it is off, so you must set MSTP(171)=1 before initialization. There are two submodes of running, with MSTP(172) being 1 or 2. In the former mode, PYTHIA will generate an event at the requested energy. This means that you have to know which energy you want beforehand. In the latter mode, PYTHIA will often return without having generated an event -- with flag MSTI(61)=1 to signal that -- and you are then requested to give a new energy. The energy spectrum of accepted events will then, in the end, be your naive input spectrum weighted with the cross-section of the processes you study. We will come back to this.

The energy can be varied, whichever frame is given in the PYINIT call. (Except for 'USER', where such information is fed in via the HEPEUP common block and thus beyond the control of PYTHIA.) When the frame is 'CMS', PARP(171) should be filled with the fractional energy of each event, i.e. $E_{\mathrm{cm}} =$PARP(171)$\times$WIN, where WIN is the nominal c.m. energy of the PYINIT call. Here PARP(171) should normally be smaller than unity, i.e. initialization should be done at the maximum energy to be encountered. For the 'FIXT' frame, PARP(171) should be filled by the fractional beam energy of that one, i.e. $E_{\mathrm{beam}} =$PARP(171)$\times$WIN. For the '3MOM', '4MOM' and '5MOM' options, the two four-momenta are given in for each event in the same format as used for the PYINIT call. Note that there is a minimum c.m. energy allowed, PARP(2). If you give in values below this, the program will stop for MSTP(172)=1, and will return with MSTI(61)=1 for MSTP(172)=1.

To illustrate the use of the MSTP(172)=2 facility, consider the case of beamstrahlung in $\mathrm{e}^+\mathrm{e}^-$ linear colliders. This is just for convenience; what is said here can be translated easily into other situations. Assume that the beam spectrum is given by $D(z)$, where $z$ is the fraction retained by the original $\mathrm{e}$ after beamstrahlung. Therefore $0 \leq z \leq 1$ and the integral of $D(z)$ is unity. This is not perfectly general; one could imagine branchings $\mathrm{e}^- \to \mathrm{e}^- \gamma \to \mathrm{e}^-\mathrm{e}^+\mathrm{e}^-$, which gives a multiplication in the number of beam particles. This could either be expressed in terms of a $D(z)$ with integral larger than unity or in terms of an increased luminosity. We will assume the latter, and use $D(z)$ properly normalized. Given a nominal $s = 4 E_{\mathrm{beam}}^2$, the actual $s'$ after beamstrahlung is given by $s' = z_1 z_2 s$. For a process with a cross section $\sigma(s)$ the total cross section is then

$$\sigma_{\mathrm{tot}} = \int_0^1 \int_0^1 D(z_1) \, D(z_2) \sigma(z_1 z_2 s) \, \d z_1 \, \d z_2~.$$ (156)

The cross section $\sigma$ may in itself be an integral over a number of additional phase space variables. If the maximum of the differential cross section is known, a correct procedure to generate events is

81.

pick $z_1$ and $z_2$ according to $D(z_1) \, \d z_1$ and $D(z_2) \, \d z_2$, respectively;

82.

pick a set of phase space variables of the process, for the given $s'$ of the event;

83.

evaluate $\sigma(s')$ and compare with $\sigma_{\mathrm{max}}$;

84.

if event is rejected, then return to step 1 to generate new variables;

85.

else continue the generation to give a complete event.

You as a user are assumed to take care of step 1, and present the resulting kinematics with incoming $\mathrm{e}^+$ and $\mathrm{e}^-$ of varying energy. Thereafter PYTHIA will do steps 2-5, and either return an event or put MSTI(61)=1 to signal failure in step 4.

The maximization procedure does search in phase space to find $\sigma_{\mathrm{max}}$, but it does not vary the $s'$ energy in this process. Therefore the maximum search in the PYINIT call should be performed where the cross section is largest. For processes with increasing cross section as a function of energy this means at the largest energy that will ever be encountered, i.e. $s' = s$ in the case above. This is the `standard' case, but often one encounters other behaviours, where more complicated procedures are needed. One such case would be the process $\mathrm{e}^+\mathrm{e}^-\to \mathrm{Z}^{*0} \to \mathrm{Z}^0 \mathrm{h}^0$, which is known to have a cross section that increases near the threshold but is decreasing asymptotically. If one already knows that the maximum, for a given Higgs mass, appears at 300 GeV, say, then the PYINIT call should be made with that energy, even if subsequently one will be generating events for a 500 GeV collider.

In general, it may be necessary to modify the selection of $z_1$ and $z_2$ and assign a compensating event weight. For instance, consider a process with a cross section behaving roughly like $1/s$. Then the $\sigma_{\mathrm{tot}}$ expression above may be rewritten as

$$\sigma_{\mathrm{tot}} = \int_0^1 \int_0^1 \frac{D(z_1)}{z_1}... ..._2)}{z_2} \, z_1 z_2 \sigma(z_1 z_2 s) \, \d z_1 \, \d z_2 ~.$$ (157)

The expression $z_1 z_2 \sigma(s')$ is now essentially flat in $s'$, i.e. not only can $\sigma_{\mathrm{max}}$ be found at a convenient energy such as the maximum one, but additionally the PYTHIA generation efficiency (the likelihood of surviving step 4) is greatly enhanced. The price to be paid is that $z$ has to be selected according to $D(z)/z$ rather than according to $D(z)$. Note that $D(z)/z$ is not normalized to unity. One therefore needs to define

$${\cal I}_D = \int_0^1 \frac{D(z)}{z} \, \d z ~,$$ (158)

and a properly normalized

$$D'(z) = \frac{1}{{\cal I}_D} \, \frac{D(z)}{z} ~.$$ (159)

Then

$$\sigma_{\mathrm{tot}} = \int_0^1 \int_0^1 D'(z_1) \, D'(z_2)... ...cal I}_D^2 \, z_1 z_2 \sigma(z_1 z_2 s) \, \d z_1 \, \d z_2 ~.$$ (160)

Therefore the proper event weight is ${\cal I}_D^2 \, z_1 z_2$. This weight should be stored by you, for each event, in PARP(173). The maximum weight that will be encountered should be stored in PARP(174) before the PYINIT call, and not changed afterwards. It is not necessary to know the precise maximum; any value larger than the true maximum will do, but the inefficiency will be larger the cruder the approximation. Additionally you must put MSTP(173)=1 for the program to make use of weights at all. Often $D(z)$ is not known analytically; therefore ${\cal I}_D$ is also not known beforehand, but may have to be evaluated (by you) during the course of the run. Then you should just use the weight $z_1 z_2$ in PARP(173) and do the overall normalization yourself in the end. Since PARP(174)=1 by default, in this case you need not set this variable specially. Only the cross sections are affected by the procedure selected for overall normalization, the events themselves still are properly distributed in $s'$ and internal phase space.

Above it has been assumed tacitly that $D(z) \to 0$ for $z \to 0$. If not, $D(z)/z$ is divergent, and it is not possible to define a properly normalized $D'(z) = D(z)/z$. If the cross section is truly diverging like $1/s$, then a $D(z)$ which is nonvanishing for $z \to 0$ does imply an infinite total cross section, whichever way things are considered. In cases like that, it is necessary to impose a lower cut on $z$, based on some physics or detector consideration. Some such cut is anyway needed to keep away from the minimum c.m. energy required for PYTHIA events, see above.

The most difficult cases are those with a very narrow and high peak, such as the $\mathrm{Z}^0$. One could initialize at the energy of maximum cross section and use $D(z)$ as is, but efficiency might turn out to be very low. One might then be tempted to do more complicated transforms of the kind illustrated above. As a rule it is then convenient to work in the variables $\tau_z = z_1 z_2$ and $y_z = (1/2) \ln (z_1/z_2)$, cf. section .

Clearly, the better the behaviour of the cross section can be modelled in the choice of $z_1$ and $z_2$, the better the overall event generation efficiency. Even under the best of circumstances, the efficiency will still be lower than for runs with fix energy. There is also a non-negligible time overhead for using variable energies in the first place, from kinematics reconstruction and (in part) from the phase space selection. One should therefore not use variable energies when not needed, and not use a large range of energies $\sqrt{s'}$ if in the end only a smaller range is of experimental interest.

This facility may be combined with most other aspects of the program. For instance, it is possible to simulate beamstrahlung as above and still include bremsstrahlung with MSTP(11)=1. Further, one may multiply the overall event weight of PARP(173) with a kinematics-dependent weight given by PYEVWT, although it is not recommended (since the chances of making a mistake are also multiplied). However, a few things do not work.

$\bullet$

It is not possible to use pile-up events, i.e. you must have MSTP(131)=0.

$\bullet$

The possibility of giving in your own cross-section optimization coefficients, option MSTP(121)=2, would require more input than with fixed energies, and this option should therefore not be used. You can still use MSTP(121)=1, however.

$\bullet$

The multiple interactions scenario with MSTP(82)$\geq 2$ only works approximately for energies different from the initialization one. If the c.m. energy spread is smaller than a factor 2, say, the approximation should be reasonable, but if the spread is larger one may have to subdivide into subruns of different energy bins. The initialization should be made at the largest energy to be encountered -- whenever multiple interactions are possible (i.e. for incoming hadrons and resolved photons) this is where the cross sections are largest anyway, and so this is no further constraint. There is no simple possibility to change PARP(82) during the course of the run, i.e. an energy-independent $p_{\perp 0}$ must be assumed. The default option MSTP(82)=1 works fine, i.e. does not suffer from the constraints above. If so desired, $p_{\perp\mathrm{min}}=$PARP(81) can be set differently for each event, as a function of c.m. energy. Initialization should then be done with PARP(81) as low as it is ever supposed to become.