How to Generate Weighted Events

By default PYTHIA generates unweighted events, i.e. all events in a run are on an equal footing. This means that corners of phase space with low cross sections are poorly populated, as it should be. However, sometimes one is interested in also exploring such corners, in order to gain a better understanding of physics. A typical example would be the jet cross section in hadron collisions, which is dropping rapidly with increasing jet , and where it is interesting to trace this drop over several orders of magnitude. Experimentally this may be solved by prescaling event rates already at the trigger level, so that all high- events are saved but only a fraction of the lower- ones. In this section we outline procedures to generate events in a similar manner.

Basically two approaches can be used. One is to piece together results from different subruns, where each subrun is restricted to some specific region of phase space. Within each subrun all events then have the same weight, but subruns have to be combined according to their relative cross sections. The other approach is to let each event come with an associated weight, that can vary smoothly as a function of . These two alternatives correspond to stepwise or smoothly varying prescaling factors in the experimental analogue. We describe them one after the other.

The phase space can be sliced in many different ways. However, for the
jet rate and many other processes, the most natural variable would be
itself. (For production of a lepton pair by -channel
resonances, the invariant mass would be a better choice.) It is not
possible to specify beforehand the jet 's an event will contain,
since this is a combination of the
of the hard
scattering process with additional showering activity, with
hadronization, with underlying event and with the jet clustering
approach actually used. However, one would expect a strong correlation
between the
scale and the jet 's. Therefore
the full
range can be subdivided into a set of
ranges by using the `CKIN(3)` and `CKIN(4)` variables as
lower and upper limits. This could be done e.g. for adjacent
non-overlapping bins 10-20,20-40,40-70, etc.

Only if one would like to cover also very small is there a
problem with this strategy: since the naïve jet cross section is
divergent for
, a unitarization procedure is
implied by setting `CKIN(3) = 0` (or some other low value). This
unitarization then disregards the actual `CKIN(3)` and `CKIN(4)`
values and generates events over the full phase space. In order not to
double-count, then events above the intended upper limit of the first
bin have to be removed by brute force.

A simple but complete example of a code performing this task (with some primitive histogramming) is the following:

C...All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Three Pythia functions return integers, so need declaring. INTEGER PYK,PYCHGE,PYCOMP C...EXTERNAL statement links PYDATA on most platforms. EXTERNAL PYDATA C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Selection of hard-scattering subprocesses. COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) C...Parameters. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Bins of pT. DIMENSION PTBIN(10) DATA PTBIN/0D0,10D0,20D0,40D0,70D0,110D0,170D0,250D0,350D0,1000D0/ C...Main parameters of run: c.m. energy and number of events per bin. ECM=2000D0 NEV=1000 C...Histograms. CALL PYBOOK(1,'dn_ev/dpThat',100,0D0,500D0) CALL PYBOOK(2,'dsigma/dpThat',100,0D0,500D0) CALL PYBOOK(3,'log10(dsigma/dpThat)',100,0D0,500D0) CALL PYBOOK(4,'dsigma/dpTjet',100,0D0,500D0) CALL PYBOOK(5,'log10(dsigma/dpTjet)',100,0D0,500D0) CALL PYBOOK(11,'dn_ev/dpThat, dummy',100,0D0,500D0) CALL PYBOOK(12,'dn/dpTjet, dummy',100,0D0,500D0) C...Loop over pT bins and initialize. DO 300 IBIN=1,9 CKIN(3)=PTBIN(IBIN) CKIN(4)=PTBIN(IBIN+1) CALL PYINIT('CMS','p','pbar',ECM) C...Loop over events. Remove unwanted ones in first pT bin. DO 200 IEV=1,NEV CALL PYEVNT PTHAT=PARI(17) IF(IBIN.EQ.1.AND.PTHAT.GT.PTBIN(IBIN+1)) GOTO 200 C...Store pThat. Cluster jets and store variable number of pTjet. CALL PYFILL(1,PTHAT,1D0) CALL PYFILL(11,PTHAT,1D0) CALL PYCELL(NJET) DO 100 IJET=1,NJET CALL PYFILL(12,P(N+IJET,5),1D0) 100 CONTINUE C...End of event loop. 200 CONTINUE C...Normalize cross section to pb/GeV and add up. FAC=1D9*PARI(1)/(DBLE(NEV)*5D0) CALL PYOPER(2,'+',11,2,1D0,FAC) CALL PYOPER(4,'+',12,4,1D0,FAC) C...End of loop over pT bins. 300 CONTINUE C...Take logarithm and plot. CALL PYOPER(2,'L',2,3,1D0,0D0) CALL PYOPER(4,'L',4,5,1D0,0D0) CALL PYNULL(11) CALL PYNULL(12) CALL PYHIST END

The alternative to slicing the phase space is to used weighted events.
This is possible by making use of the `PYEVWT` routine:

**Purpose:**- to allow you to reweight event cross sections,
by process type and kinematics of the hard scattering. There exists
two separate modes of usage, described in the following.

For`MSTP(142) = 1`, it is assumed that the cross section of the process is correctly given by default in PYTHIA, but that one wishes to generate events biased to a specific region of phase space. While the`WTXS`factor therefore multiplies the naïve cross section in the choice of subprocess type and kinematics, the produced event comes with a compensating weight`PARI(10) = 1./WTXS`, which should be used when filling histograms etc. In the`PYSTAT(1)`table, the cross sections are unchanged (up to statistical errors) compared with the standard cross sections, but the relative composition of events may be changed and need no longer be in proportion to relative cross sections. A typical example of this usage is if one wishes to enhance the production of high- events; then a weight like`WTXS`(with some fixed number) might be appropriate. See`PARI(2)`for a discussion of overall normalization issues.

For`MSTP(142) = 2`, on the other hand, it is assumed that the true cross section is really to be modified by the multiplicative factor`WTXS`. The generated events therefore come with unit weight, just as usual. This option is really equivalent to replacing the basic cross sections coded in PYTHIA, but allows more flexibility: no need to recompile the whole of PYTHIA.

The routine will not be called unless`MSTP(142)`, and never if `minimum-bias'-type events (including elastic and diffractive scattering) are to be generated as well. Further, cross sections for additional multiple interactions or pile-up events are never affected. A dummy routine`PYEVWT`is included in the program file, so as to avoid unresolved external references when the routine is not used. `WTXS:`- multiplication factor to ordinary event cross section;
to be set (by you) in
`PYEVWT`call. **Remark :**- at the time of selection, several variables in the
`MINT`and`VINT`arrays in the`PYINT1`common block contain information that can be used to make the decision. The routine provided in the program file explicitly reads the variables that have been defined at the time`PYEVWT`is called, and also calculates some derived quantities. The given list of information includes subprocess type`ISUB`, , , , , , , , , , , , , and a few more. Some of these may not be relevant for the process under study, and are then set to zero. **Warning:**- the weights only apply to the hard-scattering subprocesses. There is no way to reweight the shape of initial- and final-state showers, fragmentation, or other aspects of the event.

There are some limitations to the facility. `PYEVWT` is called
at an early stage of the generation process, when the hard kinematics
is selected, well before the full event is constructed. It then cannot
be used for low-, elastic or diffractive events, for which no
hard kinematics has been defined. If such processes are included,
the event weighting is switched off. Therefore it is no longer an
option to run with `CKIN(3) = 0`.

Which weight expression to use may take some trial and error. In the above case, a reasonable ansatz seems to be a weight behaving like , where four powers of are motivated by the partonic cross section behaving like , and the remaining two by the fall-off of parton densities. An example for the same task as above one would then be:

C...All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Three Pythia functions return integers, so need declaring. INTEGER PYK,PYCHGE,PYCOMP C...EXTERNAL statement links PYDATA on most platforms. EXTERNAL PYDATA C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Selection of hard-scattering subprocesses. COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) C...Parameters. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Main parameters of run: c.m. energy, pTmin and number of events. ECM=2000D0 CKIN(3)=5D0 NEV=10000 C...Histograms. CALL PYBOOK(1,'dn_ev/dpThat',100,0D0,500D0) CALL PYBOOK(2,'dsigma/dpThat',100,0D0,500D0) CALL PYBOOK(3,'log10(dsigma/dpThat)',100,0D0,500D0) CALL PYBOOK(4,'dsigma/dpTjet',100,0D0,500D0) CALL PYBOOK(5,'log10(dsigma/dpTjet)',100,0D0,500D0) C...Initialize with weighted events. MSTP(142)=1 CALL PYINIT('CMS','p','pbar',ECM) C...Loop over events; read out pThat and event weight. DO 200 IEV=1,NEV CALL PYEVNT PTHAT=PARI(17) WT=PARI(10) C...Store pThat. Cluster jets and store variable number of pTjet. CALL PYFILL(1,PTHAT,1D0) CALL PYFILL(2,PTHAT,WT) CALL PYCELL(NJET) DO 100 IJET=1,NJET CALL PYFILL(4,P(N+IJET,5),WT) 100 CONTINUE C...End of event loop. 200 CONTINUE C...Normalize cross section to pb/GeV, take logarithm and plot. FAC=1D9*PARI(2)/5D0 CALL PYFACT(2,FAC) CALL PYFACT(4,FAC) CALL PYOPER(2,'L',2,3,1D0,0D0) CALL PYOPER(4,'L',4,5,1D0,0D0) CALL PYHIST END C************************************************************* SUBROUTINE PYEVWT(WTXS) C...Double precision and integer declarations. IMPLICIT DOUBLE PRECISION(A-H, O-Z) IMPLICIT INTEGER(I-N) INTEGER PYK,PYCHGE,PYCOMP C...Common block. COMMON/PYINT1/MINT(400),VINT(400) C...Read out pThat^2 and set weight. PT2=VINT(48) WTXS=PT2**3 RETURN ENDNote that, in

The above solutions do not work if the weighting should also depend
on the shower evolution, i.e. not only on the hard process. Such a
rejection/acceptance `weight' is particularly convenient when
trying to combine matrix-element-generated events of different jet
multiplicities (in addition to whatever is the basic process) without
doublecounting, e.g. using the L-CKKW or MLM matching prescriptions
[Cat01]. To handle such situations, a routine `UPVETO` has
been introduced. The name is intended to emphasize the logical place
along with the `UPINIT` and `UPEVNT` routines for the handling
of external user processes (see section ), but it
can also be used for internal PYTHIA processes.

**Purpose:**- to allow the user the possibility to abort the
generation of an event immediately after parton showers but before
underlying event and handronization is considered. The routine is called
only if
`MSTP(143) = 1`, and only for the old, virtuality-ordered showers in`PYEVNT`, i.e. has not been implemented for the -ordered showers in`PYEVNW`. `IVETO :`- specifies choice made by user.
`= 0 :`- retain current event and generate it in full.
`= 1 :`- abort generation of current event and move to next.

**Note 1:**- if resonances like
, , , Higgs and
SUSY particles are handed undecayed from
`UPEVNT`, or are generated internally in PYTHIA, they will also be undecayed at this stage; if decayed their decay products will have been allowed to shower. **Note 2:**- all partons at the end of the shower phase are stored in
the
`HEPEVT`commonblock, see section . The interesting information is:`NHEP :`- the number of such partons, in entries 1 through
`NHEP`, `ISTHEP(I) :`- where all intermediate products have code 2 and
all the final ones (at the time
`UPVETO`is called) code 1, `IDHEP(I) :`- the particle identity code according to PDG conventions,
`JMOHEP(1,I) :`- the mother position,
`PHEP(J,I) :`- the of the particle.

`HEPEVT`variables are zeroed. **Note 3:**- the cross section of processes, as shown with
`CALL PYSTAT(1)`, is reduced when events are aborted. Such events are also counted in the "fraction of events that fail fragmentation cuts" in the last line of this table.

In the `HEPEVT` event record, all intermediate and `final'
particles of the hard process itself, i.e. of the matrix-element
calculation, are listed as documentation lines, with `ISTHEP(I) = 2`.
The `final'-state particles that actually are defined when `UPVETO`
is called, after shower evolution but before multiple interactions have
been added, have `ISTHEP(I) = 1`. These point back to the one of the
`ISTHEP(I) = 2` partons they originate from, as first mother. If a
system does not radiate, the same set of partons will be repeated twice,
once with `ISTHEP(I) = 2` and once with `ISTHEP(I) = 1`. A more
typical example would be that a set of partons with `ISTHEP(I) = 1`
point back to the same `mother' with `ISTHEP(I) = 2`. Note that
all the intermediate stages of the shower evolution are not shown,
only the original mother and all the final daughters. If the
mother index is zero for an `ISTHEP(I) = 1` parton, it comes
from initial-state radiation.

Of course, the separation of which radiation comes from where is often gauge-dependent, i.e. model-dependent, and so should not be over-stressed. For instance, in a decay sequence, the gluon radiation off the is rather unambiguous, while and form part of the same radiating system, although some sensible separation is provided by the program.