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 naïve 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.
`PARP(171)``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.
`PARP(171)``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
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 , where
is the fraction retained by the original after beamstrahlung.
Therefore
and the integral of is unity.
This is not perfectly general; one could imagine branchings
, which gives a multiplication
in the number of beam particles. This could either be expressed in
terms of a with integral larger than unity or in terms of an
increased luminosity. We will assume the latter, and use
properly normalized. Given a nominal
,
the actual after beamstrahlung is given by
.
For a process with a cross section the total
cross section is then

(155) |

- 87.
- pick and according to and , respectively;
- 88.
- pick a set of phase space variables of the process, for the given of the event;
- 89.
- evaluate and compare with ;
- 90.
- if event is rejected, then return to step 1 to generate new variables;
- 91.
- else continue the generation to give a complete event.

The maximization procedure does search in phase space to find
, but it does not vary the 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. 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
,
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 and
and assign a compensating event weight. For instance, consider
a process with a cross section behaving roughly like . Then the
expression above may be rewritten as

(156) |

(157) |

(158) |

(159) |

Above it has been assumed tacitly that for . If not, is divergent, and it is not possible to define a properly normalized . If the cross section is truly diverging like , then a which is nonvanishing for 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 , 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 . One could initialize at the energy of maximum cross section and use 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 and , cf. section .

Clearly, the better the behaviour of the cross section can be modelled in the choice of and , 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 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.

- It is not possible to use pile-up events, i.e. you must have
`MSTP(131) = 0`. - 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. - The multiple interactions scenario with
`MSTP(82)`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 must be assumed. By contrast, for`MSTP(82) = 1``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.