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
.