User Hooks

Sometimes it may be convenient to step in during the generation process: to modify the built-in cross sections, to veto undesirable events or simply to collect statistics at various stages of the evolution. There is a base class UserHooks that gives you this access at a few selected places. This class in itself does nothing; the idea is that you should write your own derived class for your task. A few very simple derived classes come with the program, mainly as illustration.

For a derived class to be called during the execution, a pointer to an object of this class should be handed in with the
pythia.setUserHooksPtr( UserHooks*) method.

There are four distinct sets of routines. Ordered by increasing complexity, rather than by their appearance in the event-generation sequence, they are:
(i) Ones that gives you access to the event record in between the process-level and parton-level steps, or in between the parton-level and hadron-level ones. You can study the event record and decide whether to veto this event.
(ii) Ones that allow you to set a scale at with the combined parton-level MI+ISR+FSR downwards evolution in pT is temporarily interrupted, so the event can be studied and either vetoed or allowed to continue the evolution.
(iii) Ones that allow you to to study the event after the first few ISR/FSR emissions, so the event can be vetoed or allowed to continue the evolution.
(iv) Ones that gives you access to the properties of the trial hard process, so that you can modify the internal Pythia cross section by your own correction factors.
They are described further in the following.

Interrupt between the main generation levels

If your derived class redefines the
bool canVetoProcessLevel()
method to return true, then the method
bool doVetoProcessLevel(const Event& process)
will be called immediately after a hard process has been selected and stored in the process event record. You can study, but not modify, this record. Based on that you can decide whether to veto the event or let it continue to evolve. If you veto, then this event is not counted among the accepted ones, and do not contribute to the estimated cross section. The pytha.next() method will begin a completely new event, so the vetoed event will not appear in the output of pythia.next().

Note that this is different from setting the flag PartonLevel:all = off. Also in this case the event generation will stop after the process level, but an event generated up to this point is considered perfectly acceptable, and cross sections are not affected. That is, this option is intended for simple studies of hard processes, where one can save time by not generating the rest of the story. By contrast, the doVetoProcessLevel() allows you to throw away uninteresting events at an early stage to save time that way, but those events that do survive the veto are allowed to develop into complete final states (unless flags have been set otherwise).

The
bool canVetoPartonLevel()
and
bool doVetoPartonLevel(const Event& event)
are exact analogues to the above two methods, except that these ones are called after the parton level, i.e. when showers, multiple interactions and beam remnants have been set up, but hadronization and decays have not yet been performed. Information is now made available in the event event record. The difference relative to the HadronLevel:all = off flag setting follows the same pattern as above.

The effect of the vetoes can be studied in the output of the pythia.statistics() method. The "Selected" column represents the number of events that were found acceptable by the internal Pythia machinery, whereas the "Accepted" one are the events that also survived the user cuts.

Interrupt during the parton-level evolution, at a pT scale

During the parton-level evolution, multiple interactions (MI), initial-state radiation (ISR) and final-state radiation (FSR) are normally evolved downwards in one interleaved evolution sequence of decreasing pT values. For some applications, e.g matrix-element-matching approaches, it may be convenient to stop the evolution temporarily when the "hard" emissions have been considered, but before continuing with the more time-consuming soft activity. Based on these hard partons one can make a decision whether the event at all falls in the intended event class, e.g. has the "right" number of parton-level jets. If yes then, as for the methods above, the evolution will continue all the way up to a complete event. Also as above, if no, then the event will not be considered in the final cross section.

In this subsection we outline the possibility to interrupt at a given pT scale, in the next to interrupt after a given number of emissions.

To use this possibility you need to redefine
bool canVetoPT()
to return true and
double scaleVetoPT()
to return the pT scale at which you want to study the event.

The key routine, where you decide whether the event should be vetoed (return true) or not (false), is
bool doVetoPT(int iPos, const Event& event)
Here
iPos is the position/status when the routine is called:
= 0 when no MI, ISR or FSR occured above the veto scale;
= 1 when inside the interleaved MI + ISR + FSR evolution, after an MI process;
= 2 when inside the interleaved MI + ISR + FSR evolution, after an ISR emission;
= 3 when inside the interleaved MI + ISR + FSR evolution, after an FSR emission;
= 4 for the optional case where FSR is deferred from the interleaved evolution and only considered separately afterward (then alternative 3 would never occur);
= 5 is for subsequent resonance decays, and is called once for each decay in a chain such as t -> b W, W -> u dbar.

The event record contains a list of all partons generated so far, also including intermediate ones not part of the "current final state", and also those from further multiple interactions. This may not be desirable for comparisons with matrix-element calculations. The method
void subEvent(const Event& event, bool isHardest = true)
offers a simple recipe to extract a list of only the current partons from the hardest interaction, as relevant for iPos codes 0 - 4. With isHardest = false instead the latest "subprocess" is extracted, as relevant when iPos is 5, where it corresponds to the partons in the currently considered decay.

The result is stored in the class member Event workEvent. The daughter1() and daughter2() both return the position in the original event record (process or event), so you can trace the full history, if of interest. The workEvent can e.g. be sent on to a jet clustering algorithm. You are free to edit workEvent as you desire, e.g. boost to its rest frame before analysis, or remove particles that should not be analyzed.

Interrupt during the parton-level evolution, after a step

This option is closely related to the one above, so we do not repeat the introduction, nor the possibilities to study the event record, also by using subEvent(...). What is different is that this method gives access to the event as it looks like after each of the first few steps in the downwards evolution, irrespectively of the pT scales of these branchings. Furthermore, it is here assumed that the focus is on the hardest subprocess, so that ISR/FSR emissions associated with additional MI's are not considered.

To use the possibility to study the event after the first steps you need to redefine
bool canVetoStep()
to return true and
int numberVetoStep()
to return up to how many steps each of ISR and FSR (for the hardest interaction) that you want to be able to study. The number of steps defaults to the first one only.

The key routine, where you decide whether the event should be vetoed (return true) or not (false), is
bool doVetoStep( int iPos, int nISR, int nFSR, const Event& event)
Here
iPos is the position from where the routine has been called, options 2 - 5 of the doVetoPT(...) routine above, while options 0 and 1 are not relevant here;
nISR is the number of ISR emissions in the hardest process so far; and
nFSR is the number of FSR emissions in the hardest process so far.
For resonance decays, iPos = 5, the nISR is set 0 and nFSR refers to the number of emissions in the currently studied system.

Modify cross-sections

If you want to modify a cross section you need to redefine
bool canModifySigma()
to return true and
double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr, const PhaseSpace* phaseSpacePtr, bool inEvent)
to provide the factor by which you want to see the cross section modified. If you return unity then the normal cross section is obtained. Note that, unlike the methods above, these modifications do not lead to a difference between the number of "selected" events and the number of "accepted" ones, since the modifications occur already before the "selected" level. The integrated cross section of a process is modified, of course.

What makes the multiplySigmaBy(...) routine somewhat tricky to write is that the hard-process event has not yet been constructed, so one is restricted to use the information available in the phase-space and cross-section objects currently being accessed. Which of their methods are applicable depends on the process, in particular the number of final-state particles. The UserHooks code contains explicit instructions about which methods provide meaningful information.

The inEvent flag is true when this method is called from within the event-generation machinery and false when it is called at the initialization stage of the run, when the cross section is explored to find a maximum for later Monte Carlo usage. Cross-section modifications should be independent of this flag, for consistency, but if multiplySigmaBy(...) is used to collect statistics on the original kinematics distributions before cuts, then it is important to be able to exclude the initialization stage from comparisons.

Note that the cross section is only modifiable for normal hard processes. It does not affect the cross section in further multiple interactions, nor in elastic/diffractive/minimum-bias events.

One derived class is supplied as an example how this facility can be used to reweight cross sections in the same spirit as is done with QCD cross sections for the minimum-bias/underlying-event description:

class  SuppressSmallPT( pT0timesMI, numberAlphaS, useSameAlphaSasMI)  
suppress small-pT production for 2 -> 2 processes only, while leaving other processes unaffected. The basic suppression factor is pT^4 / ((k*pT0)^2 + pT^2)^2, where pT refers to the current hard subprocess and pT0 is the same energy-dependent dampening scale as used for multiple interactions. The optional arguments provide further variability.
argument pT0timesMI : corresponds to the additional factor k in the above formula. It is by default equal to 1 but can be used to explore deviations from the expected value.
argument numberAlphaS : if this number n is bigger than the default 0, the corresponding number of alpha_strong factors is also reweighted from the normal renormalization scale to a modified one, i.e. a further suppression factor ( alpha_s((k*pT0)^2 + Q^2_ren) / alpha_s(Q^2_ren) )^n is introduced.
argument useSameAlphaSasMI : regulates which kind of new alpha_strong value is evaluated for the numerator in the above expression. It is by default the same as set for multiple interactions (i.e. same starting value at M_Z and same order of running), but if false instead the one for hard subprocesses. The denominator alpha_s(Q^2_ren) is always the value used for the "original", unweighted cross section.

Final comments

All the possibilities above can be combined freely and also be combined with the standard flags. An event would then survive only if it survived each of the possible veto methods. There are no hidden interdependencies in this game, but of course some combinations may not be particularly meaningful. For instance, if you set PartonLevel:all = off then the doVetoPT(...) and doVetoPartonLevel(...) locations in the code are not even reached, so they would never be called.

An example how the above methods can be used for toy studies is found in main10.cc.