next up previous contents
Next: How to Run with Up: The Process Generation Program Previous: General Event Information   Contents


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 event 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 naïve 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:


\fbox{\texttt{CALL PYEVWT(WTXS)}}

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.


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.


\fbox{\texttt{CALL UPVETO(IVETO)}}

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 $p_{\perp}$-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 $\mathrm{W}^{\pm}$, $\mathrm{Z}^0$, $\t $, 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 $(p_x, p_y, p_z, E, m)$ of the particle.
The rest of the 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 $\t\to \b\mathrm{W}^+ \to \b\u\overline{\mathrm{d}}$ decay sequence, the gluon radiation off the $\b $ is rather unambiguous, while $\u $ and $\overline{\mathrm{d}}$ form part of the same radiating system, although some sensible separation is provided by the program.


next up previous contents
Next: How to Run with Up: The Process Generation Program Previous: General Event Information   Contents
Stephen Mrenna 2007-10-30