Matrix Element Merging
CKKW-L merging [Lon01] allows for a consistent merging of parton
showers with
matrix elements to include multiple well-separated partons. The
algorithm implemented in PYTHIA is described in [Lon11]. To
perform matrix element merging, the user has to supply LHE
files [Alw07] for the hard process and the corresponding
process with up to N additional jets.
The usage of the merging procedure is illustrated in a few
example main programs
(main81.cc
, main82.cc
,
main83.cc
and main84.cc
, together with the
input files main81.cmnd
, main82.cmnd
and
main84.cmnd
). These examples should of course only serve
as an illustration, and as such will not make use of the merging in
all possible ways. For full generality, the example programs link to
LHAPDF, FastJet and HepMC. Of course the user is welcome to remove
these dependencies. To remove the FastJet dependence, the functions
calculating example observables have to be deleted. Removing the
LHAPDF dependence requires changing the cmnd input files to choose an
inbuilt PDF, as outlined in the
PDF documentation. The
HepMC dependence can be removed by erasing the code allowing for HepMC
output.
Three very short LHE files (w+_production_lhc_0.lhe
,
w+_production_lhc_1.lhe
, w+_production_lhc_2.lhe
)
are included in the distribution. These files are not intended for
physics studies, but only serve as input for the example main
programs. For realistic studies, the user has to supply LHE files.
In the generation of LHE files, the value of the factorisation
scale used in the PDFs is not important, since the cross section will
be multiplied by ratios of PDFs to adjust to the PYTHIA starting
scales. The same is true for the renormalisation scale (and starting
value αs(MZ)) used to evaluate
αs. Coupling and scale choices by the user
will be transferred to the merging routines.
LHE files should be regularised with a jet measure, and
additional cuts on the partons in the LHEF generation avoided as much
as possible. This means that the merging scale is always a more
stringent cut than all other cuts on the partons. Of course, if the
hard process itself is divergent, a cut needs to be made. However,
this should be chosen in such a way as to not exclude regions that
will be available to the matrix elements with additional jets. An
example is QCD di-jet production with additional jets: Say the
2 -> 2 process is regularised with a pTmin cut
of pTminCut = 100
GeV, and the
2 - >3 sample is regularised with a kTmin-cut of
kTminCut = 50
GeV. This would mean that when reclustering
the emission in the 2 -> 3 sample, we could end up with a
pT = pTminNow of the 2 -> 2 configuration with
pTminCut > pTminNow, which is excluded in the
2 -> 2 sample. Thus, the 2 -> 3 sample will include a
Sudakov factor not included in the 2 -> 2 sample, resulting
in merging scale dependencies. Such dependencies can be avoided if
the cuts on the hard process are minimal.
Of course, additional cuts on electroweak particles are
allowed. These should be the same for all samples with
0 <= n <= N additional partons.
If it is not possible to generate LHE files with minimal cuts,
the user can choose to use the MergingHooks
structures in
order to decide how much influence to attribute to parton shower
histories in which the reclustered lowest multiplicity process does
not pass the matrix element cuts. This is described below.
When generating LHE files, please refrain from explicitly putting
unstable particles (e.g. massive gauge bosons) in the final state.
Rather, specify a resonance by its decay products, e.g. if Les Houches
events for the pp → Z + jets → e+e- + jets process
are desired, generate the matrix element events with the Z decay
included. From a physical point of view, on-shell final massive gauge
bosons should not be considered part of a hard process, since only
the boson decay products will be detectable. Furthermore,
non-narrow-width approximation contributions are not present if the
ME generator only produces on-shell bosons. Interference effects
between different production channels for the decay products would
also be neglected. These points seem an unnecessary restriction on
the accuracy of the ME calculation. In addition, there is a
technical reason for this strategy. Since some matrix element
generators choose to put additional information on intermediate
bosons into Les Houches events, depending on if they pass a certain
criterion (e.g. being close to the mass shell), without exact
knowledge of this criterion, the only feasible way of bookkeeping the
hard process is by identifying outgoing decay products. The code
implemented in PYTHIA8 uses this strategy and may become unreliable
otherwise.
For all merging purposes, processes with different charge of
outgoing leptons are considered different processes. That means
e.g. that e+νe+ jets and
e-ν̄e + jets
are considered independent processes. If the user wishes to generate
distributions including effects of more than one process, merged
samples for all independent processes should be generated separately
and added afterwards.
When the matrix element merging is used to produce HepMC
[Dob01] files to be analysed with RIVET [Buc10],
special care needs to taken in how the cross section is read by RIVET
(see below).
To specify the merging conditions, additionally information on
the merging scale value and the functional definition of the merging
scale is needed. A few standard definitions of merging scales are
available. This makes the user interface less cumbersome.
Different choices intrinsic to the CKKW-L merging procedure might
be relevant for the user as well. The base
class MergingHooks
gives the user the opportunity to
define the functional form of the merging scale. In the following,
the usage of the merging machinery to consistently include LHE files
with additional jets into PYTHIA will be discussed.
Merging with merging scale defined in kT
The quickest way to include processes with additional jets is to
produce LHE files with one of the standard ways to define the merging
scale. The following switches are provided:
flag
Merging:doKTMerging
(default = off
)
If the additional jets in the LHE files have been regulated by
a kT cut, the user can supply the merging scale definition by
setting this flag to on. kT here and below means cutting on
Durham kT for e+e- collisions, and cutting on
longitudinally invariant kT for hadronic collisions.
Since currently, a few slightly different definitions of
longitudinally invariant kT are considered, a specific form
can be chosen by setting the switch
mode
Merging:ktType
(default = 1
; minimum = 1
; maximum = 3
)
Precise functional definition of longitudinally
invariant kT. For e+e- collisions, Durham kT is
always defined by the square root of min{ 2*min[
Ei2, Ej2] * [ 1 -
cosθij] }, so that this switch will have no effect.
option
1 : Longitudinally invariant kT is defined by
the square root of the minimum of minimal jet kinematic pT
(pTkin,min = min{ pT,i } ) and
pTlon,min = min{ min[ pT,i2,
pT,j2] * [ (Δyij)2 +
(Δφij)2 ] / D2 } , i.e.
kT = min{ √pTkin,min, √pTlon,min }
for hadronic collisions. Note that the true rapidity of partons is used.
option
2 : Longitudinally invariant kT is defined by
the square root of the minimum of minimal jet kinematic pT
(pTkin,min = min{ pT,i } ) and
pTlon,min = min{ min[ pT,i2,
pT,j2] * [
(Δηij)2 +
(Δφij)2 ] / D2 }, i.e.
kT = min{ √pTkin,min, √pTlon,min }
for hadronic collisions. Note that the pseudorapidity of partons is used.
option
3 : Longitudinally invariant kT is defined by
the square root of the minimum of minimal jet kinematic pT
(pTkin,min = min{ pT,i } ) and
pTlon,min = min{ min[ pT,i2,
pT,j2] * [ cosh(Δηij) -
cos(Δφij) ] / D2 } , i.e.
kT = min{ √pTkin,min, √pTlon,min }
for hadronic collisions.
mode
Merging:nJetMax
(default = 0
; minimum = 0
)
Maximal number of additional jets in the matrix element
parm
Merging:TMS
(default = 0.0
)
The value of the merging scale. The name is inspired by the scale in
evolution equations, which is often called 't', and the suffix 'MS' stands
for merging scale. In the particular case of kT-merging, this
would be the value of the kT-cut in GeV. For any merging scale
definition, this input is considered the actual value of the merging
scale.
word
Merging:Process
(default = void
)
The string specifying the hard core process, in MG/ME notation. If
e.g. W + jets merging should be performed, set this to
pp>e+ve
(without white spaces or quotation marks).
This string may contain resonances in the MG/ME notation, e.g. for merging
pp→Z W+→q q̄ e+νe + jets,
the string pp>(z>jj)(w+>e+ve)
would be applicable.
flag
Merging:doMGMerging
(default = off
)
Even easier, but highly non-general, is to perform the merging with
MadGraph/MadEvent-produced LHE files, with a merging scale defined by
a kT cut. For this, set this switch to on. The merging scale
value will be read from the +1 jet LHE file by searching for the
string ktdurham
, and extracting the value from
value = ktdurham
. Also, the hard process will be read from
the +0 jet LHE file, from the line containing the string @1
(the tag specifying the first process in the MadGraph process card).
For this to work, PYTHIA should be initialised on LHE files called
NameOfYourLesHouchesFile_0.lhe
(+0 jet sample) and
NameOfYourLesHouchesFile_1.lhe
(+1 jet sample) and the
same naming convention for LHE files with two or more additional jets.
Since for this option, the merging scale value is read from the
LHEF, no merging scale value needs to be supplied by setting
Merging:TMS . Also, the hard process is read from LHEF, the
input Merging::Process does not have to be defined.
However, the maximal number of merged jets still has to be supplied by
setting Merging:nJetMax.
Histogramming the events
After the event has been processed,
histograms for observables of interest need to be filled. In order to
achieve good statistical accuracy for all jet multiplicities and all
subprocesses contributing to one jet multiplicity, generally a fixed
number of unit-weighted events is read from each Les Houches Event
file. To then arrive at the correct prediction, for each of these
events, histogram bins should be filled with the corresponding cross
section, or weighted with unit weight and normalised at the end to
the generated cross section for each jet multiplicity separately.
Still another, even more important, event weight that has to
applied on an event-by-event basis is the CKKW-L-weight. This
corrective weight is the main outcome of the merging procedure and
includes the correct no-emission probabilities, PDF weights and
αs factors. This means that the merging
implementation will generate weighted events. The CKKW-L-weight can be
accessed by the following function:
double Info::mergingWeight()
Returns the CKKW-L weight for the current event.
Note that to avoid confusion, this function does not include the
the weight of a phase space point (given
by Info::weight()). This weight will differ from
unity when reading in weighted Les Houches events. In this case, the
full weight with which to fill histogram bins is
Info::mergingWeight() * Info::weight().
Finally, to arrive at a correct relative normalisation of the
contributions from different number of additional jets in the matrix
element, each histogram should be rescaled with the accepted cross
section given by
Info::sigmaGen(). The accepted cross section includes
the effect of vetoes generating Sudakov form factors for the matrix
elements, and is in general only known after the run.
This final step can of course be skipped if the accepted cross
section had been estimated before the histgramming run, and histogram
bins had instead been filled with the weight
Info::mergingWeight() * σest(number of
additional jets in current ME sample). This is the way HepMC
events should be weighted to produce correct relative weights of
events (see below, and particularly examine the example program
main84.cc
).
Examples how to use these options are given in main81.cc
(kT merging) and main84.cc
(automatic MG/ME merging
for RIVET usage).
Merging with user-defined merging scale function
For all other merging scale definitions, the procedure is
slightly more complicated, since the user has to write a small piece
of code defining the merging scale. To allow for a user defined
procedure, set the input
flag
Merging:doUserMerging
(default = off
)
General user defined merging on/off.
Then, set
the Merging:nJetMax, Merging:TMS
and Merging:Process input as before.
Since during execution, PYTHIA needs to evaluate the merging
scale with the definition of the user, the user interface is designed
in a way similar to the
UserHooks
strategy. The class controlling the merging
scale definition is called MergingHooks
.
Initialisation
To initialise the merging with user-defined merging scale, we
should construct a class derived from MergingHooks
, with
a constructor and destructor
MergingHooks::MergingHooks()
virtual MergingHooks::~MergingHooks()
The constructor and destructor do not need to do anything.
For the class to be called during execution, a pointer to an
object of the class should be handed in with the
Pythia::setMergingHooksPtr( MergingHooks*)
method.
An examples of this procedure are given in main82.cc
.
Defining a merging scale
Then, in the spirit of the UserHooks
class, the user
needs to supply the process to be merged by defining a methods to
evaluate the merging scale variable.
virtual double MergingHooks::tmsDefinition(const Event& event)
This method will have to calculate the value of the merging scale
defined in some variable from the input event record. An example of
such a function is given in main82.cc
.
The base class MergingHooks
contains many functions
giving information on the hard process, to make the definition of the
merging scale as easy as possible:
int MergingHooks::nMaxJets()
Return the maximum number of additional jets to be merged.
int MergingHooks::nHardOutPartons()
Returns the number of outgoing partons in the hard core process.
int MergingHooks::nHardOutLeptons()
Returns the number of outgoing leptons in the hard core process.
int MergingHooks::nHardInPartons()
Returns the number of incoming partons in the hard core process.
int MergingHooks::nHardInLeptons()
Returns the number of incoming leptons in the hard core process.
int MergingHooks::nResInCurrent()
The number of resonances in the hard process reconstructed from the
current event. If e.g. the ME configuration was
pp -> (w+->e+ve)(z -> mu+mu-)jj, and the ME generator put
both intermediate bosons into the LHE file, this will return 2.
double MergingHooks::tms()
Returns the value used as the merging scale.
Filling output histograms for the event then proceeds along the
lines described above in "Histogramming the events".
The full procedure is outlined in main82.cc
. Special
care needs to be taken when the output is stored in the form of HepMC
files for RIVET usage.
Defining a cut on lowest jet multiplicity events
It can sometimes happen that when generating LHE files, a fairly
restrictive cut has been used when generating the lowest multiplicity
matrix element configurations. Then, it can happen that states that
are (in the generation of a parton shower history) constructed by
reclustering from higher multiplicity configurations, do not pass
this matrix element cut.
Consider as an example pure QCD dijet merging, when up to one
additional jet should be merged. Three-jet matrix element
configurations for which the reclustered two-jet state does not pass
the cuts applied to the two-jet matrix element would never have been
produced by showering the two-jet matrix element. This means that the
three-jet matrix element includes regions of phase space that would
never have been populated by the parton shower. Thus, since the
matrix element phase space is larger than the shower phase space,
merging scale dependencies are expected. A priori, this is not
troublesome, since the aim of matrix element merging is to include
regions of phase space outside the range of the parton shower
approximation into the shower. An example is the inclusion of
configurations with only unordered histories.
Clearly, if the parton shower phase space is very constrained by
applying stringent cuts to the two-jet matrix element, merging scale
dependencies can become sizable, as was e.g. seen in [Lon11]
when forcing shower emissions to be ordered both in the evolution
variable and in rapidity. To influence the effect of large phase
space differences for shower emissions and matrix element
configurations due to LHEF generation cuts, the user has to write a
small piece of code overwriting method
virtual double MergingHooks::dampenIfFailCuts(const Event&event)
multiplicity reclustered state as an input Event. From this input
event, the user can then check if matrix element cuts are
fulfilled. The return value will be internally multiplied to the
CKKW-L weight of the current event. Thus, if the user wishes to
suppress contributions not passing particular cuts, a number smaller
than unity can be returned.
Note that this method gives the user access to the lowest
multiplicity state, which ( e.g. in the case of incomplete histories)
does not have to be a 2 → 2 configuration. Also, changing the
weight of the current event by hand is of course a major intervention
in the algorithm, and should be considered very carefully. Generally,
if this facility would have to be used extensively, it is certainly
preferable to be less restrictive when applying additional,
non-merging-scale-related cuts to the matrix element.
An example how to force a cut on lowest multiplicity reclustered
states for pure QCD matrix element configurations is given by
main83.cc
(to be used with e.g. main82.cmnd
).
Influencing the construction of all possible histories
Even more powerful - and dangerous - is influencing the construction
of histories directly. This should only be attempted by expert users. If you
believe manipulations completely unavoidable, we advise you to take great care
when redefining the following functions.
virtual bool MergingHooks::canCutOnRecState()
In the base class this method returns false. If you redefine it
to return true then the method doCutOnRecState(...)
will be called for each reclustered state encountered in the generation of
all possible histories of the matrix element state.
virtual bool MergingHooks::doCutOnRecState(const Event&event)
This routine will be supplied internally with every possible reclustered
event that can be reached by reclustering any number of partons in
the matrix element input state. The new, reclustered, states can then be
analysed. If the method returns false, no further clusterings of the
reclustered state will be attempted, thus disallowing all history branches
which contain the disallowed state.
Clearly, these methods are highly intrusive. It could e.g. happen that no
history is allowed, which would make merging impossible. One example where
this method could be useful is if cuts on the core 2 -> 2 processes
have to be checked, and the method
MergingHooks::dampenIfFailCuts(const Event& event)
is not
sufficiently effective.
Matrix element merging and HepMC output for RIVET
An example how to produce matrix element merged events to be analysed
with RIVET is given by main84.cc
.
The main issue is that the output of separate RIVET runs can not
in general be combined. To perform a matrix element merging, we
however need to runs over different LHE files. The solution to this
problem (so far) is to only perform one RIVET run for all matrix
elements, i.e. print the events for all ME parton multiplicities,
with the correct weights, to a single HepMC file. Since the correct
weight includes the cross section of the different samples after
Sudakov vetoes --- which is not a priori known --- the cross sections
have to be estimated in a test run, before the actual production run
is performed. Finally, the cross section of the last event in the
HepMC file has to be taken as the full merged cross section
sigma_merge = Sum_{i=0}^N Sum_{j=0}*^{nEvents}
sigma_est(i)*wckkwl(j).
This procedure is outlined in main84.cc
.
Further variables
For more advanced manipulations of the merging machinery, all
parameter changes that were investigated in [Lon11] are
supplied. Please check [Lon11] for a detailed discussion of
the switches.
These switches allow enthusiastic users to perform a systematic
assessment of the merging prescription. Apart from this, we advise the
non-expert user to keep the default values.
flag
Merging:includeMassive
(default = on
)
If on, use the correct massive evolution variable and massive
splitting kernels in the reconstruction and picking of parton shower
histories of the matrix element. If off, reconstruct evolution
scales, kinematics and splitting kernels as if all partons were
massless.
flag
Merging:enforceStrongOrdering
(default = off
)
If on, preferably pick parton shower histories of the matrix element
which have strongly ordered consecutive splittings, i.e. paths in
which consecutive reclustered evolution scales are separated by a
user-defined factor.
parm
Merging:scaleSeparationFactor
(default = 1.0
; minimum = 1.0
; maximum = 10.0
)
The factor by which scales should differ to be classified as strongly
ordered.
flag
Merging:orderInRapidity
(default = off
)
If on, preferably pick parton shower histories of the matrix element
with consecutive splittings ordered in rapidity and pT.
flag
Merging:pickByFullP
(default = on
)
If on, pick parton shower histories of the matrix element by the full
shower splitting kernels, including potential ME corrections and
Jacobians from joined evolution measures.
flag
Merging:pickByPoPT2
(default = off
)
If on, pick parton shower histories of the matrix element by the
shower splitting kernels divided by the evolution pT.
flag
Merging:pickBySumPT
(default = off
)
If on, exclusively pick parton shower histories of the matrix element
for which have the smallest sum of scalar evolution pT for consecutive
splittings has been calculated.
flag
Merging:includeRedundant
(default = off
)
If on, then also include PDF ratios and αs
factors in the splitting probabilities used for picking a parton shower
history of the matrix element, when picking histories by the full shower
splitting probability. As argued in [Lon11], this should not
be done since a reweighting with PDF ratios and αs
factors will be performed. However, it can give useful insight in how
sensitive the results are to the prescription on how to choose PS
histories.
parm
Merging:nonJoinedNorm
(default = 1.0
; minimum = 0.0
; maximum = 10.0
)
Normalisation factor with which to multiply splitting probability for
splittings without joined evolution equation.
parm
Merging:fsrInRecNorm
(default = 1.0
; minimum = 0.0
; maximum = 10.0
)
Normalisation factor with which to multiply splitting probability for
final state splittings with an initial state recoiler.
parm
Merging:aCollFSR
(default = 1.0
; minimum = 0.0
; maximum = 10.0
)
Factor with which to multiply the scalar pT of a final state
splitting, when choosing the history by the smallest sum of scalar
pT. Default value taken from Herwig++ [Tul09].
parm
Merging:aCollISR
(default = 0.9
; minimum = 0.0
; maximum = 10.0
)
Factor with which to multiply the scalar pT of an initial state
splitting, when choosing the history by the smallest sum of scalar
pT. Default value taken from Herwig++ [Tul09].
mode
Merging:unorderedScalePrescrip
(default = 0
; minimum = 0
; maximum = 1
)
When the parton shower history of the matrix element contains a
sequence of splittings which are not ordered in evolution pT
(called an unordered history), this sequence is interpreted as a combined
emission. Then, a decision on which starting scale for trial emissions
off reconstructed states in this sequence of unordered splittings has
to be made. Two options are available:
option
0 : Use larger of the two reconstructed (unordered)
scales as starting scale.
option
1 : Use smaller of the two reconstructed (unordered)
scales as starting scale.
mode
Merging:unorderedASscalePrescrip
(default = 1
; minimum = 0
; maximum = 1
)
Prescription which scale to use to evaluate αs
weight for splittings in a sequence of splittings which are not ordered
in evolution pT.
option
0 : Use the combined splitting scale as argument in
αs, for both splittings.
option
1 : Use the true reconstructed scale as as argument in
αs, for each splitting separately.
mode
Merging:incompleteScalePrescrip
(default = 0
; minimum = 0
; maximum = 2
)
When no complete parton shower history (i.e. starting from a
2 → 2 process) for a matrix element with additional jets
can be found, such a configuration is said to have an incomplete history.
Since in incomplete histories, not all shower starting scales are
determined by clusterings, a prescription for setting the starting scale
of trial showers in incomplete histories is needed. Three options are
provided.
option
0 : Use factorisation scale as shower starting scale
for incomplete histories.
option
1 : Use sHat as shower starting scale for
incomplete histories.
option
2 : Use s as shower starting scale for
incomplete histories.
flag
Merging:allowColourShuffling
(default = off
)
If on, this will allow the algorithm to swap one colour index in the state,
when trying to find all possible clusterings, if no clustering has been
found, but more clusterings had been requested. In this way, some incomplete
histories can be avoided. Generally, we advise the non-expert user to not
touch this switch, because a slight change in the colour structure can change
the radiation pattern. To however study the sensitivity of the predictions on
these effects, allowing for colour reshuffling can be useful.