2 <chapter name="Matrix Element Merging">
4 <h2>Matrix Element Merging</h2>
6 CKKW-L merging <ref>Lon01</ref> allows for a consistent merging of parton
8 matrix elements to include multiple well-separated partons. The
9 algorithm implemented in PYTHIA is described in <ref>Lon11</ref>. To
10 perform matrix element merging, the user has to supply LHE
11 files <ref>Alw07</ref> for the hard process and the corresponding
12 process with up to N additional jets.
14 <p/> The usage of the merging procedure is illustrated in a few
16 (<code>main81.cc</code>, <code>main82.cc</code>,
17 <code>main83.cc</code> and <code>main84.cc</code>, together with the
18 input files <code>main81.cmnd</code>, <code>main82.cmnd</code> and
19 <code>main84.cmnd</code>). These examples should of course only serve
20 as an illustration, and as such will not make use of the merging in
21 all possible ways. For full generality, the example programs link to
22 LHAPDF, FastJet and HepMC. Of course the user is welcome to remove
23 these dependencies. To remove the FastJet dependence, the functions
24 calculating example observables have to be deleted. Removing the
25 LHAPDF dependence requires changing the cmnd input files to choose an
26 inbuilt PDF, as outlined in the
27 <a href="PDFSelection.html" target="page">PDF documentation</a>. The
28 HepMC dependence can be removed by erasing the code allowing for HepMC
31 <p/> Please note that a detailed tutorial on merging in Pythia is available
33 <strong>http://home.thep.lu.se/~torbjorn/pythia8/mergingworksheet8160.pdf
36 <p/> Three very short LHE files (<code>w+_production_lhc_0.lhe</code>,
37 <code>w+_production_lhc_1.lhe</code>, <code>w+_production_lhc_2.lhe</code>)
38 are included in the distribution. These files are not intended for
39 physics studies, but only serve as input for the example main
40 programs. For realistic studies, the user has to supply LHE files.
42 <p/> In the generation of LHE files, the value of the factorisation
43 scale used in the PDFs is not important, since the cross section will
44 be multiplied by ratios of PDFs to adjust to the PYTHIA starting
45 scales. The same is true for the renormalisation scale (and starting
46 value <ei>α<sub>s</sub>(M<sub>Z</sub>)</ei>) used to evaluate
47 <ei>α<sub>s</sub></ei>. Coupling and scale choices by the user
48 will be transferred to the merging routines.
50 <p/> Multi-jet events can suffer from infrared divergences in the
51 calculation. Sensible matrix element generator (MEG) outputs should not
52 contain phase space points in which the calculation is ill-defined, meaning
53 infrared regions need to be removed by cuts. This is most conveniently done
54 by disallowing the MEG to produce partons below a
55 minimal parton-parton separation in a certain jet algorithm. Using
56 dedicated cuts to regularise MEG output is of course possible as well. Any
57 regularisation criterion defines the matrix element region: The parts of
58 phase space in which the fixed order calculation is considered valid and
59 preferable to the parton shower. Matrix element merging is combining
60 MEG events in the matrix element region with parton shower events in regions
61 outside the regularisation cut (often called parton shower region). Because
62 the regularisation cut defines a boundary between the matrix element
63 and parton shower regions, i.e. the regions to be merged into one inclusive
64 sample, it is usually called <ei> merging scale </ei>. Since many different
65 cut choices may regularise the MEG calculation, many different merging scale
66 definitions are possible. A few standard choices are listed below, as well as
67 documentation on how to use a user-defined cut criterion. In combining matrix
68 element and parton shower regions, the CKKW-L prescription tries to minimise
69 the dependence on the merging scale. This can only be achieved if the
70 combination of MEG events and parton shower populates the whole phase space.
71 Additional cuts on the partons in the LHEF generation should hence be
72 avoided as much as possible, meaning that the merging scale cut should always
73 pose a more stringent cut than all other cuts on the partons. Of course, if
74 the hard process itself is divergent, cuts need to be made. However, this
75 should be chosen in such a way as to not exclude regions that will be
76 available to the matrix elements with additional jets. An example is QCD
77 di-jet production with additional jets: Say the <ei>2 -> 2</ei> process is
78 regularised with a <ei>pTmin</ei> cut of pTminCut = 100 GeV, and
79 the <ei>2 - >3</ei> sample is regularised with a <ei>kTmin</ei>-cut of
80 kTminCut = 50 GeV. This would mean that when reclustering
81 the emission in the 2 -> 3 sample, we could end up with a
82 pT value <ei>pTminNow</ei> of the 2 -> 2 configuration with
83 <ei>pTminCut > pTminNow</ei>, which is excluded in the
84 2 -> 2 sample. Thus, the 2 -> 3 sample will include a
85 Sudakov factor not included in the 2 -> 2 sample, resulting
86 in merging scale dependencies. Such dependencies can be avoided if
87 the additional cuts on the hard process are minimal.
89 <p/> Of course, additional cuts on electroweak particles are
90 allowed. These should be the same for all samples with
91 <ei>0 <= n <= N</ei> additional partons.
93 <p/> If it is not possible to generate LHE files with minimal cuts,
94 the user can choose to use the <code>MergingHooks</code> structures in
95 order to decide how much influence to attribute to parton shower
96 histories in which the reclustered lowest multiplicity process does
97 not pass the matrix element cuts. This is described below.
99 <p/> When generating LHE files, we advise against putting
100 unstable particles (e.g. massive gauge bosons) in the final state.
101 Rather, specify a resonance by its decay products, e.g. if Les Houches
102 events for the <ei>pp → Z + jets → e+e- + jets</ei> process
103 are desired, generate the matrix element events with the <ei>Z</ei> decay
104 included. From a physical point of view, on-shell final massive gauge
105 bosons should not be considered part of a hard process, since only
106 the boson decay products will be detectable. Furthermore,
107 non-narrow-width approximation contributions are not present if the
108 ME generator only produces on-shell bosons. Interference effects
109 between different production channels for the decay products would
110 also be neglected. These points seem an unnecessary restriction on
111 the accuracy of the ME calculation. In addition, there is a
112 technical reason for this strategy. Since some matrix element
113 generators choose to put additional information on intermediate
114 bosons into Les Houches events, depending on if they pass a certain
115 criterion (e.g. being close to the mass shell), without exact
116 knowledge of this criterion, the only feasible way of bookkeeping the
117 hard process is by identifying outgoing decay products.
119 <p/> Despite these considerations, (massive) gauge bosons in the final state
120 are allowed in the hard process definition. This is useful particularly for
121 Higgs physics, when different decays of the Higgs boson need to be simulated
122 after the LHEF generation.
124 <p/> For all merging purposes, processes with different charge of
125 outgoing leptons are considered different processes. That means
126 e.g. that <ei>e+ν<sub>e</sub>+ jets</ei> and
127 <ei>e-ν̄<sub>e</sub> + jets</ei>
128 are considered independent processes. If the user wishes to generate
129 distributions including effects of more than one process, merged
130 samples for all independent processes should be generated separately
131 and added afterwards. Alternatively, to combine simple processes,
132 combined LHE files can be used in conjunction with flexible containers (see
135 <p/> When the matrix element merging is used to produce HepMC
136 <ref>Dob01</ref> files to be analysed with RIVET <ref>Buc10</ref>,
137 special care needs to taken in how the cross section is read by RIVET
140 <p/> To specify the merging conditions, additionally information on
141 the merging scale value and the functional definition of the merging
142 scale is needed. A few standard definitions of merging scales are
143 available. We hope this makes the user interface less cumbersome.
145 <p/> Different choices intrinsic to the CKKW-L merging procedure might
146 be relevant for the user as well. The base
147 class <code>MergingHooks</code> gives the user the opportunity to
148 define the functional form of the merging scale. In the following,
149 the usage of the merging machinery to consistently include LHE files
150 with additional jets into PYTHIA will be discussed.
153 <h3>Merging scale definitions</h3>
155 <p/> The quickest way to include processes with additional jets is to
156 produce LHE files with one of the standard ways to define the merging
157 scale. Three standard ways to define a merging scale (minimal <ei>kT</ei>,
158 minimal evolution <ei>pT</ei> and by three cuts) are implemented. All of these
159 prescriptions are equivalent - different definitions have only been introduced
160 for the convenience of users, who might be limited by which cuts can be used
161 in the generation of LHE files. Below, we describe how to switch on and use
162 these different merging scale definitions.
164 <h4>Merging with merging scale defined in kT:</h4>
166 <flag name="Merging:doKTMerging" default="off">
167 If the additional jets in the LHE files have been regulated by
168 a <ei>kT</ei> cut, the user can supply the merging scale definition by
169 setting this flag to on. <ei>kT</ei> here and below means cutting on
170 Durham <ei>kT</ei> for <ei>e+e-</ei> collisions, and cutting on
171 longitudinally invariant <ei>kT</ei> for hadronic collisions. Please note
172 that this particular merging scale definition will check <ei>kT</ei> between
173 all pairs of <ei>u,d,c,s,b,g</ei> partons.
177 Currently, the name longitudinally invariant <ei>kT</ei> is used
178 for a few jet recombination algorithms with slightly different
179 jet measures. A specific form can be chosen by setting the switch
181 <modepick name="Merging:ktType" default="1" min="1" max="3">
182 Precise functional definition of longitudinally
183 invariant <ei>kT</ei>. For e+e- collisions, <ei>Durham kT</ei> is
184 always defined by the square root of <ei>min{ 2*min[ </sub>
185 E<sub>i</sub><sup>2</sup>, E<sub>j</sub><sup>2</sup>] * [ 1 -
186 cosθ<sub>ij</sub>] }</ei>, so that this switch will have no effect.
187 <option value="1"> Longitudinally invariant <ei>kT</ei> is defined by
188 the square root of the minimum of minimal jet kinematic <ei>pT</ei>
189 (<ei>p<sub>Tkin,min</sub><sup>2</sup> = min{ p<sub>T,i</sub><sup>2</sup> }
190 </ei>) and <ei>p<sub>Tlon,min</sub><sup>2</sup> =
191 min{ min[ p<sub>T,i</sub><sup>2</sup>, p<sub>T,j</sub><sup>2</sup>] *
192 [ (Δy<sub>ij</sub>)<sup>2</sup> +
193 (Δφ<sub>ij</sub>)<sup>2</sup> ] / D<sup>2</sup> }</ei> , i.e.
194 <ei>kT = min{ √p<sub>Tkin,min</sub><sup>2</sup>,
195 √p<sub>Tlon,min</sub><sup>2</sup> }</ei> for hadronic collisions. Note
196 that the true rapidity of partons is used.
198 <option value="2">Longitudinally invariant <ei>kT</ei> is defined by
199 the square root of the minimum of minimal jet kinematic <ei>pT</ei>
200 (<ei>p<sub>Tkin,min</sub><sup>2</sup> = min{ p<sub>T,i</sub><sup>2</sup>
201 } </ei>) and <ei>p<sub>Tlon,min</sub><sup>2</sup> =
202 min{ min[ p<sub>T,i</sub><sup>2</sup>,
203 p<sub>T,j</sub><sup>2</sup>] * [
204 (Δη<sub>ij</sub>)<sup>2</sup> +
205 (Δφ<sub>ij</sub>)<sup>2</sup> ] / D<sup>2</sup> }</ei>, i.e.
206 <ei>kT = min{ √p<sub>Tkin,min</sub><sup>2</sup>,
207 √p<sub>Tlon,min</sub><sup>2</sup> }</ei>
208 for hadronic collisions. Note that the pseudorapidity of partons is used.
210 <option value="3"> Longitudinally invariant <ei>kT</ei> is defined by
211 the square root of the minimum of minimal jet kinematic <ei>pT</ei>
212 (<ei>p<sub>Tkin,min</sub><sup>2</sup> = min{ p<sub>T,i</sub><sup>2</sup>
214 <ei>p<sub>Tlon,min</sub><sup>2</sup> =
215 min{ min[ p<sub>T,i</sub><sup>2</sup>,
216 p<sub>T,j</sub><sup>2</sup>] * [ cosh(Δη<sub>ij</sub>) -
217 cos(Δφ<sub>ij</sub>) ] / D<sup>2</sup> } </ei>, i.e.
218 <ei>kT = min{ √p<sub>Tkin,min</sub><sup>2</sup>
219 , √p<sub>Tlon,min</sub><sup>2</sup> }</ei>
220 for hadronic collisions.
224 <modeopen name="Merging:nJetMax" default="0" min="0">
225 Maximal number of additional jets in the matrix element
228 <parm name="Merging:TMS" default="0.0">
229 The value of the merging scale. The name is inspired by the scale in
230 evolution equations, which is often called 't', and the suffix 'MS' stands
231 for merging scale. In the particular case of <ei>kT</ei>-merging, this
232 would be the value of the <ei>kT</ei>-cut in GeV. For any merging scale
233 definition, this input is considered the actual value of the merging
237 <word name="Merging:Process" default="void">
238 The string specifying the hard core process, in MG4/ME notation.
242 If e.g. <ei>W + jets</ei> merging should be performed, set this to
243 <code>pp>e+ve</code> (<i>without white spaces or quotation marks</i>).
244 This string may contain resonances in the MG/ME notation, e.g. for merging
245 <ei>pp→Z W<sup>+</sup>→q q̄ e+ν<sub>e</sub> + jets</ei>,
246 the string <code>pp>(z>jj)(w+>e+ve)</code> would be applicable.
249 A lot more flexible hard process definitions are possible. To not dwell too
250 much on these details here, we will come back to the process string at the end
254 <flag name="Merging:doMGMerging" default="off">
255 Even easier, but highly non-general, is to perform the merging with
256 MadGraph/MadEvent-produced LHE files, with a merging scale defined by
257 a <ei>kT</ei> cut. For this, set this switch to on. The merging scale
258 value will be read from the +1 jet LHE file by searching for the
259 string <code>ktdurham</code>, and extracting the value from <code>
260 value = ktdurham</code>. Also, the hard process will be read from
261 the +0 jet LHE file, from the line containing the string <code>@1</code>
262 (the tag specifying the first process in the MadGraph process card).
263 For this to work, PYTHIA should be initialised on LHE files called
264 <code>NameOfYourLesHouchesFile_0.lhe</code> (+0 jet sample) and
265 <code>NameOfYourLesHouchesFile_1.lhe</code> (+1 jet sample) and the
266 same naming convention for LHE files with two or more additional jets.
267 Since for this option, the merging scale value is read from the
268 LHEF, no merging scale value needs to be supplied by setting <code>
269 Merging:TMS </code>. Also, the hard process is read from LHEF, the
270 input <code>Merging::Process</code> does not have to be defined.
271 However, the maximal number of merged jets still has to be supplied by
272 setting <code>Merging:nJetMax</code>.
276 <h4>Merging with merging scale defined in Pythia evolution pT:</h4>
278 If the LHE file has been regularised by cutting on the minimal Pythia evolution
279 pT in the state, this can also be used as a merging scale right away. For this,
282 <flag name="Merging:doPTLundMerging" default="off">
283 The merging scale is then defined by finding the minimal Pythia evolution pT
284 between sets of radiator, emitted an recoiler partons. For this
285 particular merging scale definition, <ei>u,d,c,s,b,g</ei> are considered
286 partons. The Pythia evolution pT of a single three-parton set is defined by
288 <ei>pT<sub>evol</sub> = z<sub>ijk</sub>(1-z<sub>ijk</sub>)
289 Q<sub>ij</sub><sup>2</sup></ei> for FSR, where <ei>i</ei> is the radiating
290 parton, <ei>j</ei> is the emitted parton and <ei>k</ei> is the recoiler, and
291 <ei> Q<sub>ij</sub><sup>2</sup> =
292 (p<sub>i</sub> + p<sub>j</sub>)<sup>2</sup> </ei>, and
293 <ei>z<sub>ijk</sub> =
294 x<sub>i,jk</sub> / (x<sub>i,jk</sub> + x<sub>j,ik</sub>)</ei> with
295 <ei>x<sub>i,jk</sub> =
296 2 p<sub>i</sub> (p<sub>i</sub> + p<sub>j</sub> + p<sub>k</sub>)
297 / (p<sub>i</sub> + p<sub>j</sub> + p<sub>k</sub>)<sup>2</sup> </ei>
299 <ei>pT<sub>evol</sub> = (1-z<sub>ijk</sub>)
300 Q<sub>ij</sub><sup>2</sup></ei> for ISR, where <ei>i</ei> is the radiating
301 parton, <ei>j</ei> is the emitted parton and <ei>k</ei> is the second initial
303 <ei> Q<sub>ij</sub><sup>2</sup> =
304 -(p<sub>i</sub> - p<sub>j</sub>)<sup>2</sup> </ei>, and
305 <ei>z<sub>ijk</sub> =
306 (p<sub>i</sub> - p<sub>j</sub> + p<sub>k</sub>)<sup>2</sup>
307 / (p<sub>i</sub> + p<sub>k</sub>)<sup>2</sup> </ei>.
310 When using this option, the merging scale is defined by the minimum
311 <ei>pT<sub>evol</sub></ei> for all combinations of three partons in the event,
312 irrespectively of flavour or colour-connections. The merging scale value will
313 be read from the <code>Merging:TMS</code> parameter, so that this needs to be
314 set just as in the case of the <ei>kT</ei>-merging prescription. Of course you
315 will also need to set <code>Merging:Process</code> and the maximal number of
316 additional matrix element jets <code>Merging:nJetMax</code>.
319 <h4>Merging with merging scale defined by a combination of cuts:</h4>
321 It is possible to regularise QCD divergences in a LHE file by applying cuts
322 to the kinematical pT of jets (<ei>pT<sub>i</sub></ei>), combined with a cut on
323 <ei> ΔR<sub>ij</sub></ei> between jets and a cut on invariant mass
324 <ei> Q<sub>ij</sub></ei> of jet pairs. The combination of these standard cuts
325 can also serve as a merging scale. For this, use this setting
327 <flag name="Merging:doCutBasedMerging" default="off">
328 This switch will use cuts on (<ei>pT<sub>i</sub></ei>),
329 <ei> ΔR<sub>ij</sub></ei> and
330 <ei> Q<sub>ij</sub></ei> to define when parton shower emissions are allowed.
331 Please note for this particular merging scale definition, only light jets
332 (<ei>u,d,c,s,g</ei>) are checked.
336 The values of the cuts will then be read from
338 <parm name="Merging:QijMS" default="0.0">
339 The value of the invariant mass cut <ei> Q<sub>ij</sub></ei> of pairs of final
340 state partons used in the matrix element generation.
343 <parm name="Merging:pTiMS" default="0.0">
344 The value of the minimal transverse momentum cut <ei> pT<sub>i</sub></ei> on
345 final state partons, as used in the matrix element generation.
348 <parm name="Merging:dRijMS" default="0.0">
349 The value of the minimal <ei> ΔR<sub>ij</sub></ei> separation between
350 pairs of final state partons used in the matrix element generation, where
351 <ei> ΔR<sub>ij</sub><sup>2</sup> = (Δy<sub>ij</sub>)<sup>2</sup> +
352 (Δφ<sub>ij</sub>)<sup>2</sup></ei>.
356 With knowledge of these values, and <code>Merging:doCutBasedMerging</code>,
357 Pythia will use these cuts as a separation between matrix element phase space
358 and parton shower region. If e.g. the Les Houches Events have been generated
359 with the cuts <ei> ΔR<sub>ij</sub> = 0.1 </ei>,
360 <ei>pT<sub>i</sub>= 20 GeV</ei> and <ei> Q<sub>ij</sub> = 40 GeV</ei>, set
361 <code>Merging:QijMS=40.</code>,
362 <code>Merging:pTjMS=20.</code>,
363 <code>Merging:dRijMS=0.1</code> to perform a cut-based merging. Of course
364 you will also need to set <code>Merging:Process</code> and the
365 maximal number of additional matrix element jets
366 <code>Merging:nJetMax</code>.
368 <h4>Les Houches events outside the matrix element region</h4>
371 Before continuing, we would like to point out that in general, the user
372 should make sure that the events in the Les Houches file are actually
373 calculated using the regularisation cut definition and value(s) supplied to
374 Pythia as merging scale definition and value(s). However, if LH files with
375 a large number of events and loose merging scale cuts are available, it
376 might be useful to choose a higher merging scale value, e.g. for merging
377 scale variations as part of uncertainty assessments. If CKKW-L merging is
378 enabled, Pythia will by default check if events read from Les Houches file
379 are in the matrix element region as defined by the merging scale definition
380 and merging scale value. Events outside the matrix element region will be
381 discarded. To change this default behaviour, use the flag
383 <flag name="Merging:enforceCutOnLHE" default="on">
384 This will check if the events read from LHE file are in the matrix element
385 region as defined by the merging scale definition and value(s). If on, LHE
386 input outside the matrix element region will be rejected. If off, every
387 event is assumed to pass the merging scale cut.
390 <h4>Defining the hard process</h4>
392 To perform CKKW-L matrix element merging, the user has to decide on a hard
393 process, and communicate this choice to Pythia. This is done by setting the
394 input <strong>Merging:Process</strong>
397 For single processes in the Standard Model or the MSSM, MG4/ME notation is
398 applicable. However, for some purposes, using a single simple process
399 string is not satisfactory. Mixed W<sup>+</sup> and W<sup>-</sup> events
400 in a single LHE file is a common example. For this case, it would of course
401 be perfectly allowed to perform twice, once for W<sup>+</sup> production and
402 once for W<sup>-</sup> production, and then add the results. Nevertheless, it
403 seems reasonable to alleviate difficulties by allowing for less restrictive
404 hard process definitions. Two generalisations of the process tag are
405 available: Containers and user-defined particle tags. The syntax of these
406 settings is described below.
409 In case you want multiple processes in a LHE file to be treated on equal
410 footing (e.g. <ei>W<sup>+</sup> + jets</ei> and <ei>W<sup>-</sup> + jets</ei>),
411 you should use flexible containers do specify the hard process. So far, we
412 allow the use of the containers <code>LEPTONS</code>, <code>NEUTRINOS</code>,
413 <code>BQUARKS</code>. If you use these containers, the hard process definition
414 is relatively flexible, meaning that Pythia will attempt a merging of QCD jets
415 for each event in the LHE file, and assume that all particles matching one of
416 the containers are products of the hard process. This is best explained by
417 examples. If you want to have both <ei>pp → e+ ν<sub>e</sub> + jets</ei>
418 and <ei>pp → e- ν̄<sub>e</sub> + jets</ei> events in a single
419 file, you can set <code>Merging:Process=pp>LEPTONS,NEUTRINOS</code> as hard
420 process (note that for readability, you are allowed to use commata to separate
421 container names). Combining e.g. the processes
422 <ei>pp → e+ ν<sub>e</sub></ei> and
423 <ei>pp → μ+ ν<sub>μ</sub></ei> is possible with the hard process
424 definition <code>pp>LEPTONS,NEUTRINOS</code>.
427 For maximal flexibility, the definition of the hard process by these containers
428 does not mean that each Les Houches event needs to contain particles to match
429 each container. It is enough if one container is matched. This means that
430 with the string <code>pp>LEPTONS,NEUTRINOS</code>, you can immediately process
431 <ei>pp → e+ e- </ei> events mixed with
432 <ei>pp → e+ ν<sub>e</sub></ei> events, since particles matching at
433 least one container can be found in both cases. Another example for the usage
434 of containers is mixing <ei>pp → e+ ν<sub>e</sub></ei> and
435 <ei>pp → tt̄ → e+ ν<sub>e</sub> e- ν̄<sub>e</sub>
436 bb̄</ei>. This can be accommodated by the hard process string
437 <code>Merging:Process=pp>LEPTONS,NEUTRINOS,BQUARKS</code>.
440 There is however a conceptual limitation to containers: The hard process
441 definition is necessary to ensure that when constructing lower multiplicity
442 states (that will be used to calculate the correct merging weight), the
443 structure of the hard process will be preserved. If e.g. we want the hard process
444 to be <ei>pp → Z → bb̄ </ei>, we should ensure that the lowest
445 multiplicity state contains a colour-singlet bb̄-pair. When
446 reconstructing intermediate lower multiplicity states from multi-jet matrix
447 elements, we should thus always be able to find at least one bb̄-pair. By
448 mixing different processes in a LHE file, this requirement might already be
449 violated at the level of Les Houches events. Flexible containers cannot give
450 strong conditions which flavours should be preserved in the construction of
451 the hard process. In order to avoid non-sensible results, it is hence <ei>
452 assumed that all</ei> particles matching any of the containers will be part
453 of the lowest multiplicity process. This implies that if you decide to use the
454 <code>BQUARKS</code> container, all b-quarks in the LHE file will be
455 interpreted as hard process particles, and never as additional radiation.
458 Another way to specify the hard process particles is to explicitly define the
459 particle names and identifiers. This is necessary if the matrix element merging
460 in Pythia does not contain the particles of interest. To make sure that the
461 hard process is still treated correctly, it is possible to define particles
462 in the process string. If you e.g. want the hard process to contain a particle
463 "zeta~" with PDG identifier "12345", produced in proton collisions, you
464 have to include a user-defined particle tag by setting the process string to
465 <code>pp>{zeta~,12345}</code>. The user-defined particle is enclosed in
466 curly brackets, with syntax
467 <code>{particle_name,particle_identifier}</code>, where "particle_name"
468 and "particle_identifier" are the particle name and particle identifier used
469 for this particle in the input LHE file. User-defined particles are only
470 allowed in the final state. You are free to fix user-defined particles with
471 more common ones, as long as user-defined particles are put before more common
472 particles in the process string. This means that if you e.g. wanted the hard
473 process to contain a graviton in association with a positron and an
474 electron-neutrino, you have to define the hard process as
475 <code>pp>{G,39}e+ve</code>.
478 Below you can find a list of particles predefined in the merging. If you wish
479 to include a hard process with different final state particles, you may use
480 the "curly bracket notation" outlined above.
483 The set of incoming particles us limited to:
484 <code>e-</code> (electron), <code>e+</code> (positron), <code>mu-</code> (muon),
485 <code>mu+</code> (antimuon), <code>p</code> (proton, container to hold all
486 initial state coloured particles), <code>p~</code> (identical to
487 <code>p</code> container).
490 The following intermediate particles are allowed:
491 <code>a</code> (photon), <code>z</code> (Z boson),
492 <code>w-</code> (W<sup>-</sup> boson), <code>w+</code> (W<sup>+</sup> boson),
493 <code>h</code> (scalar Higgs boson), <code>W</code> (container to hold both
494 W<sup>-</sup> and W<sup>+</sup> boson), <code>t</code> (top quark),
495 <code>t~</code> (anti-top),
496 <code>dl</code>, <code>dl~</code>, <code>ul</code>, <code>ul~</code>,
497 <code>sl</code>, <code>sl~</code>, <code>cl</code>, <code>cl~</code>,
498 <code>b1</code>, <code>b1~</code>, <code>t1</code>, <code>t1~</code>,
499 <code>dr</code>, <code>dr~</code>, <code>ur</code>, <code>ur~</code>,
500 <code>sr</code>, <code>sr~</code>, <code>cr</code>, <code>cr~</code>,
501 <code>b2</code>, <code>b2~</code>, <code>t2</code>, <code>t2~</code> (all MSSM
505 We have pre-defined the outgoing particles:
506 <code>e+</code>, <code>e-</code>, <code>ve~</code>,
507 <code>ve</code>, <code>mu+</code>, <code>mu-</code>,
508 <code>vm~</code>, <code>vm</code>, <code>ta+</code>, <code>ta-</code>,
509 <code>vt~</code>, <code>vt</code> (all SM leptons and neutrinos),
510 <code>j~</code> (container to hold all final state coloured particles),
511 <code>j</code> (container to hold all final state coloured particles),
512 <code>NEUTRINOS</code> (container to hold all final state neutrinos and
513 anti-neutrinos), <code>LEPTONS</code> (container to hold all final state
514 leptons and anti-leptons), <code>BQUARKS</code> (container to hold final
515 state b-quarks), <code>d~</code>, <code>d</code>, <code>u~</code>,
516 <code>u</code>, <code>s~</code>, <code>s</code>, <code>c~</code>,
517 <code>c</code>, <code>b~</code>, <code>b</code>, <code>t~</code>,
518 <code>t</code> (all SM quarks), <code>a</code>, <code>z</code>,
519 <code>w-</code>, <code>w+</code> (all SM electro-weak bosons),
520 <code>h</code> (scalar Higgs boson), <code>W</code> (container to hold both
521 W<sup>-</sup> and W<sup>+</sup> boson), <code>n1</code> (MSSM neutralino),
522 <code>dl~</code>, <code>dl</code>, <code>ul~</code>, <code>ul</code>,
523 <code>sl~</code>, <code>sl</code>, <code>cl~</code>, <code>cl</code>,
524 <code>b1~</code>, <code>b1</code>, <code>t1~</code>, <code>t1</code>,
525 <code>dr~</code>, <code>dr</code>, <code>ur~</code>, <code>ur</code>,
526 <code>sr~</code>, <code>sr</code>, <code>cr~</code>, <code>cr</code>,
527 <code>b2~</code>, <code>b2</code>, <code>t2~</code>, <code>t2</code>
528 (all MSSM squarks). Other outgoing particles are possible if you use the
529 "curly bracket notation" described earlier.
532 <h3>Histogramming the events</h3>
533 After the event has been processed, histograms for observables of interest
534 need to be filled. In order to achieve good statistical accuracy for all jet
535 multiplicities and all subprocesses contributing to one jet multiplicity,
536 generally a fixed number of unit-weighted events is read from each Les
537 Houches Event file. To then arrive at the correct prediction, for each of these
538 events, histogram bins should be filled with the corresponding cross
539 section, or weighted with unit weight and normalised at the end to
540 the generated cross section for each jet multiplicity separately.
542 <p/> Still another, even more important, event weight that has to
543 applied on an event-by-event basis is the CKKW-L-weight. This
544 corrective weight is the main outcome of the merging procedure and
545 includes the correct no-emission probabilities, PDF weights and
546 α<sub>s</sub> factors. This means that the merging
547 implementation will generate weighted events. The CKKW-L-weight can be
548 accessed by the following function:
550 <p/><strong> double Info::mergingWeight() </strong> <br/>
551 Returns the CKKW-L weight for the current event.
553 <p/> Note that to avoid confusion, this function does not include the
554 the weight of a phase space point (given
555 by <strong>Info::weight()</strong>). This weight will differ from
556 unity when reading in weighted Les Houches events. In this case, the
557 full weight with which to fill histogram bins is
558 <strong>Info::mergingWeight() * Info::weight()</strong>.
560 <p/> Finally, to arrive at a correct relative normalisation of the
561 contributions from different number of additional jets in the matrix
562 element, each histogram should be rescaled with the accepted cross
564 <strong>Info::sigmaGen()</strong>. The accepted cross section includes
565 the effect of vetoes generating Sudakov form factors for the matrix
566 elements, and is in general only known after the run.
568 <p/> This final step can of course be skipped if the accepted cross
569 section had been estimated before the histogramming run, and histogram
570 bins had instead been filled with the weight
571 <strong>Info::mergingWeight() * σ<sub>est</sub>(number of
572 additional jets in current ME sample)</strong>. This is the way HepMC
573 events should be weighted to produce correct relative weights of
574 events (see below, and particularly examine the example program
575 <code>main84.cc</code>).
577 <p/> Examples how to use these options are given in <code>main81.cc</code>
578 (<ei>kT</ei> merging) and <code>main84.cc</code> (automatic MG/ME merging
582 <h3>Merging with user-defined merging scale function</h3>
584 <p/> For all other merging scale definitions, the procedure is
585 slightly more complicated, since the user has to write a small piece
586 of code defining the merging scale. To allow for a user defined
587 procedure, set the input
589 <flag name="Merging:doUserMerging" default="off">
590 General user defined merging on/off.
595 the <strong>Merging:nJetMax</strong>, <strong>Merging:TMS</strong>
596 and <strong>Merging:Process</strong> input as before.
598 <p/> Since during execution, PYTHIA needs to evaluate the merging
599 scale with the definition of the user, the user interface is designed
600 in a way similar to the
601 <code>UserHooks</code> strategy. The class controlling the merging
602 scale definition is called <code>MergingHooks</code>.
604 <h4>Initialisation</h4>
606 <p/> To initialise the merging with user-defined merging scale, we
607 should construct a class derived from <code>MergingHooks</code>, with
608 a constructor and destructor
611 <method name="MergingHooks::MergingHooks()">
613 <method name="virtual MergingHooks::~MergingHooks()">
615 The constructor and destructor do not need to do anything.
617 <p/> For the class to be called during execution, a pointer to an
618 object of the class should be handed in with the
619 <br/><code><a href="ProgramFlow.html" target="page">
620 Pythia::setMergingHooksPtr( MergingHooks*)</a></code> method.
622 An examples of this procedure are given in <code>main82.cc</code>.
624 <h4>Defining a merging scale</h4>
626 <p/> Then, in the spirit of the <code>UserHooks</code> class, the user
627 needs to supply the process to be merged by defining a methods to
628 evaluate the merging scale variable.
630 <method name="virtual double MergingHooks::tmsDefinition(const Event& event)">
631 This method will have to calculate the value of the merging scale
632 defined in some variable from the input event record. An example of
633 such a function is given in <code>main82.cc</code>.
636 <p/> The base class <code>MergingHooks</code> contains many functions
637 giving information on the hard process, to make the definition of the
638 merging scale as easy as possible:
640 <method name="int MergingHooks::nMaxJets()">
641 Return the maximum number of additional jets to be merged.
644 <method name="int MergingHooks::nHardOutPartons()">
645 Returns the number of outgoing partons in the hard core process.
648 <method name="int MergingHooks::nHardOutLeptons()">
649 Returns the number of outgoing leptons in the hard core process.
652 <method name="int MergingHooks::nHardInPartons()">
653 Returns the number of incoming partons in the hard core process.
656 <method name="int MergingHooks::nHardInLeptons()">
657 Returns the number of incoming leptons in the hard core process.
660 <method name="int MergingHooks::nResInCurrent()">
661 The number of resonances in the hard process reconstructed from the
662 current event. If e.g. the ME configuration was
663 <ei>pp -> (w+->e+ve)(z -> mu+mu-)jj</ei>, and the ME generator put
664 both intermediate bosons into the LHE file, this will return 2.
667 <method name="double MergingHooks::tms()">
668 Returns the value used as the merging scale.
671 <p/> Filling output histograms for the event then proceeds along the
672 lines described above in "Histogramming the events".
674 <p/> The full procedure is outlined in <code>main82.cc</code>. Special
675 care needs to be taken when the output is stored in the form of HepMC
676 files for RIVET usage.
678 <h4>Defining a cut on lowest jet multiplicity events</h4>
680 <p/> It can sometimes happen that when generating LHE files, a fairly
681 restrictive cut has been used when generating the lowest multiplicity
682 matrix element configurations. Then, it can happen that states that
683 are (in the generation of a parton shower history) constructed by
684 reclustering from higher multiplicity configurations, do not pass
685 this matrix element cut.
687 <p/> Consider as an example pure QCD dijet merging, when up to one
688 additional jet should be merged. Three-jet matrix element
689 configurations for which the reclustered two-jet state does not pass
690 the cuts applied to the two-jet matrix element would never have been
691 produced by showering the two-jet matrix element. This means that the
692 three-jet matrix element includes regions of phase space that would
693 never have been populated by the parton shower. Thus, since the
694 matrix element phase space is larger than the shower phase space,
695 merging scale dependencies are expected. A priori, this is not
696 troublesome, since the aim of matrix element merging is to include
697 regions of phase space outside the range of the parton shower
698 approximation into the shower. An example is the inclusion of
699 configurations with only unordered histories.
701 <p/> Clearly, if the parton shower phase space is very constrained by
702 applying stringent cuts to the two-jet matrix element, merging scale
703 dependencies can become sizable, as was e.g. seen in <ref>Lon11</ref>
704 when forcing shower emissions to be ordered both in the evolution
705 variable and in rapidity. To influence the effect of large phase
706 space differences for shower emissions and matrix element
707 configurations due to LHEF generation cuts, the user has to write a
708 small piece of code overwriting method
711 name="virtual double MergingHooks::dampenIfFailCuts(const Event&
712 event)"> This routine will be supplied internally with the lowest
713 multiplicity reclustered state as an input Event. From this input
714 event, the user can then check if matrix element cuts are
715 fulfilled. The return value will be internally multiplied to the
716 CKKW-L weight of the current event. Thus, if the user wishes to
717 suppress contributions not passing particular cuts, a number smaller
718 than unity can be returned.
721 <p/> Note that this method gives the user access to the lowest
722 multiplicity state, which ( e.g. in the case of incomplete histories)
723 does not have to be a <ei>2 → 2</ei> configuration. Also, changing the
724 weight of the current event by hand is of course a major intervention
725 in the algorithm, and should be considered very carefully. Generally,
726 if this facility would have to be used extensively, it is certainly
727 preferable to be less restrictive when applying additional,
728 non-merging-scale-related cuts to the matrix element.
730 <p/> An example how to force a cut on lowest multiplicity reclustered
731 states for pure QCD matrix element configurations is given by
732 <code>main83.cc</code> (to be used with e.g. <code>main82.cmnd</code>).
734 <h4>Influencing the construction of all possible histories</h4>
736 <p/> Even more powerful - and dangerous - is influencing the construction
737 of histories directly. This should only be attempted by expert users. If you
738 believe manipulations completely unavoidable, we advise you to take great care
739 when redefining the following functions.
741 <method name="virtual bool MergingHooks::canCutOnRecState()">
742 In the base class this method returns false. If you redefine it
743 to return true then the method <code>doCutOnRecState(...)</code>
744 will be called for each reclustered state encountered in the generation of
745 all possible histories of the matrix element state.
748 <method name="virtual bool MergingHooks::doCutOnRecState(const Event&
750 This routine will be supplied internally with every possible reclustered
751 event that can be reached by reclustering any number of partons in
752 the matrix element input state. The new, reclustered, states can then be
753 analysed. If the method returns false, the history to which the state belongs
754 will be treated as if it were unordered, i.e. this path will only be chosen
755 if no other histories are available. In this way, the number of histories
756 not fulfilling the user criterion will be minimised.
760 Clearly, these methods are highly intrusive. It could e.g. happen that no
761 history is allowed, which would make merging impossible. One example where
762 this method could be useful is if cuts on the core <ei>2 -> 2</ei> processes
763 have to be checked, and the method
764 <code>MergingHooks::dampenIfFailCuts(const Event& event)</code> is not
765 sufficiently effective.
769 <h3>Matrix element merging and HepMC output for RIVET</h3>
771 An example how to produce matrix element merged events to be analysed
772 with RIVET is given by <code>main84.cc</code>.
774 <p/> The main issue is that the output of separate RIVET runs can not
775 in general be combined. To perform a matrix element merging, we
776 however need to runs over different LHE files. The solution to this
777 problem (so far) is to only perform one RIVET run for all matrix
778 elements, i.e. print the events for all ME parton multiplicities,
779 with the correct weights, to a single HepMC file. Since the correct
780 weight includes the cross section of the different samples after
781 Sudakov vetoes --- which is not a priori known --- the cross sections
782 have to be estimated in a test run, before the actual production run
783 is performed. Finally, the cross section of the last event in the
784 HepMC file has to be taken as the full merged cross section
785 <ei>sigma_merge = Sum_{i=0}^N Sum_{j=0}*^{nEvents}
786 sigma_est(i)*wckkwl(j)</ei>.
788 <p/> This procedure is outlined in <code>main84.cc</code>.
791 <h3>Further variables</h3>
793 For more advanced manipulations of the merging machinery, all
794 parameter changes that were investigated in <ref>Lon11</ref> are
795 supplied. Please check <ref>Lon11</ref> for a detailed discussion of
798 <p/> These switches allow enthusiastic users to perform a systematic
799 assessment of the merging prescription. Apart from this, we advise the
800 non-expert user to keep the default values.
802 <flag name="Merging:includeMassive" default="on">
803 If on, use the correct massive evolution variable and massive
804 splitting kernels in the reconstruction and picking of parton shower
805 histories of the matrix element. If off, reconstruct evolution
806 scales, kinematics and splitting kernels as if all partons were
810 <flag name="Merging:enforceStrongOrdering" default="off">
811 If on, preferably pick parton shower histories of the matrix element
812 which have strongly ordered consecutive splittings, i.e. paths in
813 which consecutive reclustered evolution scales are separated by a
817 <parm name="Merging:scaleSeparationFactor" default="1.0" min="1.0" max="10.0">
818 The factor by which scales should differ to be classified as strongly
822 <flag name="Merging:orderInRapidity" default="off">
823 If on, preferably pick parton shower histories of the matrix element
824 with consecutive splittings ordered in rapidity and <ei>pT</ei>.
827 <flag name="Merging:pickByFullP" default="on">
828 If on, pick parton shower histories of the matrix element by the full
829 shower splitting kernels, including potential ME corrections and
830 Jacobians from joined evolution measures.
833 <flag name="Merging:pickByPoPT2" default="off">
834 If on, pick parton shower histories of the matrix element by the
835 shower splitting kernels divided by the evolution <ei>pT</ei>.
838 <flag name="Merging:pickBySumPT" default="off">
839 If on, exclusively pick parton shower histories of the matrix element
840 for which have the smallest sum of scalar evolution <ei>pT</ei> for consecutive
841 splittings has been calculated.
844 <flag name="Merging:includeRedundant" default="off">
845 If on, then also include PDF ratios and <ei>α<sub>s</sub></ei>
846 factors in the splitting probabilities used for picking a parton shower
847 history of the matrix element, when picking histories by the full shower
848 splitting probability. As argued in <ref>Lon11</ref>, this should not
849 be done since a reweighting with PDF ratios and <ei>α<sub>s</sub></ei>
850 factors will be performed. However, it can give useful insight in how
851 sensitive the results are to the prescription on how to choose PS
855 <parm name="Merging:nonJoinedNorm" default="1.0" min="0.0" max="10.0">
856 Normalisation factor with which to multiply splitting probability for
857 splittings without joined evolution equation.
860 <parm name="Merging:fsrInRecNorm" default="1.0" min="0.0" max="10.0">
861 Normalisation factor with which to multiply splitting probability for
862 final state splittings with an initial state recoiler.
865 <parm name="Merging:aCollFSR" default="1.0" min="0.0" max="10.0">
866 Factor with which to multiply the scalar <ei>pT</ei> of a final state
867 splitting, when choosing the history by the smallest sum of scalar
868 <ei>pT</ei>. Default value taken from Herwig++ <ref>Tul09</ref>.
871 <parm name="Merging:aCollISR" default="0.9" min="0.0" max="10.0">
872 Factor with which to multiply the scalar <ei>pT</ei> of an initial state
873 splitting, when choosing the history by the smallest sum of scalar
874 <ei>pT</ei>. Default value taken from Herwig++ <ref>Tul09</ref>.
877 <modepick name="Merging:unorderedScalePrescrip" default="0" min="0" max="1">
878 When the parton shower history of the matrix element contains a
879 sequence of splittings which are not ordered in evolution <ei>pT</ei>
880 (called an unordered history), this sequence is interpreted as a combined
881 emission. Then, a decision on which starting scale for trial emissions
882 off reconstructed states in this sequence of unordered splittings has
883 to be made. Two options are available:
884 <option value="0"> Use larger of the two reconstructed (unordered)
885 scales as starting scale.
887 <option value="1"> Use smaller of the two reconstructed (unordered)
888 scales as starting scale.
892 <modepick name="Merging:unorderedASscalePrescrip" default="1" min="0" max="1">
893 Prescription which scale to use to evaluate <ei>α<sub>s</sub></ei>
894 weight for splittings in a sequence of splittings which are not ordered
895 in evolution <ei>pT</ei>.
896 <option value="0"> Use the combined splitting scale as argument in
897 <ei>α<sub>s</sub></ei>, for both splittings.
899 <option value="1"> Use the true reconstructed scale as as argument in
900 <ei>α<sub>s</sub></ei>, for each splitting separately.
904 <modepick name="Merging:unorderedPDFscalePrescrip" default="0" min="0" max="1">
905 Prescription which scale to use to evaluate ratios of parton distibutions
906 for splittings in a sequence of splittings which are not ordered
907 in evolution <ei>pT</ei>.
908 <option value="0"> Use the combined splitting scale as argument in PDF ratios,
911 <option value="1"> Use the true reconstructed scale as argument in PDF
912 ratios, for each splitting separately.
916 <modepick name="Merging:incompleteScalePrescrip" default="0" min="0" max="2">
917 When no complete parton shower history (i.e. starting from a
918 <ei>2 → 2</ei> process) for a matrix element with additional jets
919 can be found, such a configuration is said to have an incomplete history.
920 Since in incomplete histories, not all shower starting scales are
921 determined by clusterings, a prescription for setting the starting scale
922 of trial showers in incomplete histories is needed. Three options are
924 <option value="0"> Use factorisation scale as shower starting scale
925 for incomplete histories.
927 <option value="1"> Use <ei>sHat</ei> as shower starting scale for
928 incomplete histories.
930 <option value="2"> Use <ei>s</ei> as shower starting scale for
931 incomplete histories.
935 <flag name="Merging:allowColourShuffling" default="off">
936 If on, this will allow the algorithm to swap one colour index in the state,
937 when trying to find all possible clusterings, if no clustering has been
938 found, but more clusterings had been requested. In this way, some incomplete
939 histories can be avoided. Generally, we advise the non-expert user to not
940 touch this switch, because a slight change in the colour structure can change
941 the radiation pattern. To however study the sensitivity of the predictions on
942 these effects, allowing for colour reshuffling can be useful.
945 <flag name="Merging:usePythiaQRenHard" default="on">
946 If on, this will allow the algorithm to use a dynamical renormalisation scale
947 to evaluate the strong couplings of the core hard process in dijet and
948 prompt photon events.
949 This means that the value of <ei>α<sub>s</sub></ei> used as coupling
950 of the hard process in the matrix element generation will be replaced with
951 a running coupling evaluated at the geometric mean of the squared transverse
952 masses of the two outgoing particles, as is the default prescription in
956 <flag name="Merging:usePythiaQFacHard" default="on">
957 If on, this will allow the algorithm to use a dynamical factorisation scale
958 to evaluate parton distributions associated with the hadronic cross section
959 of the core hard process in dijet and prompt photon events.
960 In the calculation of PDF ratios as part of the CKKW-L weight of an event,
961 parton distributions that should be evaluated at the scale of the core
962 2 - >2 process will be evaluated using the dynamical factorisation scale
963 Pythia would attribute to this process. This means that the hard process
964 factorisation scale is set to the smaller of the squared transverse masses
965 of the two outgoing particles.
970 <!-- Copyright (C) 2012 Torbjorn Sjostrand -->