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 $p_{\perp}$, and where it is interesting to trace this drop over several orders of magnitude. Experimentally this may be solved by prescaling events rates already at the trigger level, so that all high-$p_{\perp}$ events are saved but only a fraction of the lower-$p_{\perp}$ 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 $p_{\perp}$. 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 $p_{\perp}$ itself. (For production of a lepton pair by $s$-channel resonances, the invariant mass would be a better choice.) It is not possible to specify beforehand the jet $p_{\perp}$'s an event will contain, since this is a combination of the $\hat{p}_{\perp}$ 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 $\hat{p}_{\perp}$ scale and the jet $p_{\perp}$'s. Therefore the full $\hat{p}_{\perp}$ 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 $p_{\perp}$ is there a problem with this strategy: since the naive jet cross section is divergent for $\hat{p}_{\perp} \to 0$, 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-$p_{\perp}$ events; then a weight like WTXS $=(p_{\perp}/p_{\perp 0})^2$ (with $p_{\perp 0}$ 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)$\geq 1$, 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, $E_{\mathrm{cm}}$, $\hat{s}$, $\hat{t}$, $\hat{u}$, $\hat{p}_{\perp}$, $x_1$, $x_2$, $x_{\mathrm{F}}$, $\tau$, $y$, $\tau'$, $\cos\hat{\theta}$, 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-$p_{\perp}$, 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 $\hat{p}_{\perp}^6$, where four powers of $\hat{p}_{\perp}$ are motivated by the partonic cross section behaving like $1/\hat{p}_{\perp}^4$, 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
      END

Note that, in PYEVWT one cannot look for $\hat{p}_{\perp}$ in PARI(17), since this variable is only set at the end of the event generation. Instead the internal VINT(48) is used. The dummy copy of the PYEVWT routine found in the PYTHIA code shows what is available and how to access this.