1 <chapter name="Les Houches Accord">
3 <h2>Les Houches Accord</h2>
5 The Les Houches Accord (LHA) for user processes <ref>Boo01</ref> is the
6 standard way to input parton-level information from a
7 matrix-elements-based generator into PYTHIA. The conventions for
8 which information should be stored has been defined in a Fortran context,
9 as two commonblocks. Here a C++ equivalent is defined, as a single class.
12 The <code>LHAup</code> class is a base class, containing reading and
13 printout functions, plus two pure virtual functions, one to set
14 initialization information and one to set information on each new event.
15 Derived classes have to provide these two virtual functions to do
16 the actual work. The existing derived classes are for reading information
17 from a Les Houches Event File (LHEF), from the respective Fortran
18 commonblocks, or from PYTHIA 8 itself.
21 Normally, pointers to objects of the derived classes should be handed in
22 with the <aloc href="ProgramFlow"><code>pythia.init( LHAup*)</code></aloc>
23 method. However, with the LHEF format a filename can replace the pointer,
26 <h3>Initialization</h3>
28 The <code>LHAup</code> class stores information equivalent to the
29 <code>/HEPRUP/</code> commonblock, as required to initialize the event
30 generation chain. The main difference is that the vector container
31 now allows a flexible number of subprocesses to be defined. For the
32 rest, names have been modified, since the 6-character-limit does not
33 apply, and variables have been regrouped for clarity, but nothing
34 fundamental is changed.
37 The pure virtual function <code>setInit()</code> has to be implemented in
38 the derived class, to set relevant information when called. It should
39 return <code>false</code> if it fails to set the info.
42 Inside <code>setInit()</code>, such information can be set by the following
44 <method name="setBeamA( identity, energy, pdfGroup, pdfSet)">
45 sets the properties of the first incoming beam (cf. the Fortran
46 <code>IDBMUP(1), EBMUP(1), PDFGUP(1), PDFSUP(1)</code>), and similarly
47 a <code>setBeamB</code> method exists. The parton distribution information
48 defaults to zero. These numbers can be used to tell which PDF sets were
49 used when the hard process was generated, while the normal
50 <aloc href="PDFSelection">PDF Selection</aloc> is used for the further
51 event generation in PYTHIA.
54 <method name="setStrategy( strategy)">
55 sets the event weighting and cross section strategy. The default,
56 provided in the class constructor, is 3, which is the natural value
58 <argument name="strategy">
59 chosen strategy (cf. <code>IDWTUP</code>; see <ref>Sjo06</ref>
60 section 9.9.1 for extensive comments).
61 <argoption value="1"> events come with non-negative weight, given in units
62 of pb, with an average that converges towards the cross section of the
63 process. PYTHIA is in charge of the event mixing, i.e. for each new
64 try decides which process should be generated, and then decides whether
65 is should be kept, based on a comparison with <code>xMax</code>.
66 Accepted events therefore have unit weight.</argoption>
67 <argoption value="-1"> as option 1, except that cross sections can now be
68 negative and events after unweighting have weight +-1. You can use
69 <aloc href="EventInformation"><code>pythia.info.weight()</code></aloc>
70 to find the weight of the current event. A correct event mixing requires
71 that a process that can take both signs should be split in two, one limited
72 to positive or zero and the other to negative or zero values, with
73 <code>xMax</code> chosen appropriately for the two.</argoption>
74 <argoption value="2"> events come with non-negative weight, in unspecified
75 units, but such that <code>xMax</code> can be used to unweight the events
76 to unit weight. Again PYTHIA is in charge of the event mixing.
77 The total cross section of a process is stored in
78 <code>xSec</code>.</argoption>
79 <argoption value="-2"> as option 2, except that cross sections can now be
80 negative and events after unweighting have weight +-1. As for option -1
81 processes with indeterminate sign should be split in two.</argoption>
82 <argoption value="3"> events come with unit weight, and are thus accepted
83 as is. The total cross section of the process is stored in
84 <code>xSec</code>.</argoption>
85 <argoption value="-3"> as option 3, except that events now come with weight
86 +-1. Unlike options -1 and -2 processes with indeterminate sign need not be
87 split in two, unless you intend to mix with internal PYTHIA processes
88 (see below).</argoption>
89 <argoption value="4"> events come with non-negative weight, given in units
90 of pb, with an average that converges towards the cross section of the
91 process, like for option 1. No attempt is made to unweight the events,
92 however, but all are generated in full, and retain their original weight.
93 For consistency with normal PYTHIA units, the weight stored in
94 <code>pythia.info.weight()</code> has been converted to mb, however.
96 <argoption value="-4"> as option 4, except that events now can come
97 either with positive or negative weights.</argoption>
98 <note>Note 1</note>: if several processes have already been mixed and
99 stored in a common event file, either LHEF or some private format, it
100 would be problematical to read back events in a different order. Since it
101 is then not feasible to let PYTHIA pick the next process type, strategies
102 +-1 and +-2 would not work. Instead strategy 3 would be the recommended
103 choice, or -3 if negative-weight events are required.
104 <note>Note 2</note>: it is possible to switch on internally implemented
105 processes and have PYTHIA mix these with LHA ones according to their relative
106 cross sections for strategies +-1, +-2 and 3. It does not work for strategy
107 -3 unless the positive and negative sectors of the cross sections are in
108 separate subprocesses (as must always be the case for -1 and -2), since
109 otherwise the overall mixture of PYTHIA and LHA processes will be off.
110 Mixing is not possible for strategies +-4, since the weighting procedure
111 is not specified by the standard. (For instance, the intention may be to
112 have events biased towards larger <ei>pT</ei> values in some particular
117 <method name="addProcess( idProcess, xSec, xErr, xMax)">
118 sets info on an allowed process (cf. <code>LPRUP, XSECUP, XERRUP,
120 Each new call will append one more entry to the list of processes.
121 The choice of strategy determines which quantities are mandatory:
122 <code>xSec</code> for strategies +-2 and +-3,
123 <code>xErr</code> never, and
124 <code>xMax</code> for strategies +-1 and +-2.
127 <note>Note</note>: PYTHIA does not make active use of the (optional)
128 <code>xErr</code> values, but calculates a statistical cross section
129 error based on the spread of event-to-event weights. This should work
130 fine for strategy options +-1, but not for the others. Specifically,
131 for options +-2 and +-3 the weight spread may well vanish, and anyway
132 is likely to be an underestimate of the true error. If the author of the
133 LHA input information does provide error information you may use that -
134 this information is displayed at initialization. If not, then a relative
135 error decreasing like <ei>1/sqrt(n_acc)</ei>, where <ei>n_acc</ei>
136 is the number of accepted events, should offer a reasonable estimate.
138 <method name="setXSec( i, xSec)">
139 update the <code>xSec</code> value of the <code>i</code>'th process
140 added with <code>addProcess</code> method (i.e. <code>i</code> runs
141 from 0 through <code>sizeProc() - 1</code>, see below).
144 <method name="setXErr( i, xErr)">
145 update the <code>xErr</code> value of the <code>i</code>'th process
146 added with <code>addProcess</code> method.
149 <method name="setXMax( i, xMax)">
150 update the <code>xMax</code> value of the <code>i</code>'th process
151 added with <code>addProcess</code> method.
155 Information is handed back by the following methods:
156 <method name="idBeamA(), eBeamA(), pdfGroupBeamA(), pdfSetBeamA()">
157 and similarly with <ei>A -> B</ei>, for the two beam particles.
159 <method name="strategy()">
160 for the strategy choice.
162 <method name="sizeProc()">
163 for the number of subprocesses.
165 <method name="idProcess(i), xSec(i), xErr(i), xMax(i)">
166 for process <code>i</code> in the range <code>0 <= i <
171 The information can also be printed using the <code>listInit()</code>
172 method, e.g. <code>LHAupObject->listInit()</code>.
173 This is automatically done by the <code>pythia.init</code> call.
177 The <code>LHAup</code> class also stores information equivalent to the
178 <code>/HEPEUP/</code> commonblock, as required to hand in the next
179 parton-level configuration for complete event generation. The main
180 difference is that the vector container now allows a flexible number
181 of partons to be defined. For the rest, names have been modified,
182 since the 6-character-limit does not apply, and variables have been
183 regrouped for clarity, but nothing fundamental is changed.
186 The LHA standard is based on Fortran arrays beginning with
187 index 1, and mother information is defined accordingly. In order to
188 be compatible with this convention, the zeroth line of the C++ particle
189 array is kept empty, so that index 1 also here corresponds to the first
190 particle. One small incompatibility is that the <code>sizePart()</code>
191 method returns the full size of the particle array, including the
192 empty zeroth line, and thus is one larger than the true number of
193 particles (<code>NUP</code>).
196 The pure virtual function <code>setEvent(idProcess)</code> has to be
197 implemented in the derived class, to set relevant information when
198 called. For strategy options +-1 and +-2 the input
199 <code>idProcess</code> value specifies which process that should be
200 generated by <code>setEvent(idProcess)</code>, while
201 <code>idProcess</code> is irrelevant for strategies +-3 and +-4.
202 The <code>setEvent(idProcess)</code> function should return
203 <code>false</code> if it fails to set the info, i.e. normally that the
204 supply of events in a file is exhausted. If so, no event is generated,
205 and <code>pythia.next()</code> returns <code>false</code>. You can then
207 <aloc href="EventInformation"><code>pythia.info.atEndOfFile()</code></aloc>
208 to confirm that indeed the failure is caused in the
209 <code>setEvent(idProcess)</code> function, and decide to break out of
210 the event generation loop.
213 Inside a normal <code>setEvent(...)</code> call, information can be set
214 by the following methods:
215 <method name="setProcess( idProcess, weight, scale, alphaQED, alphaQCD)">
216 tells which kind of process occured, with what weight, at what scale,
217 and which <ei>alpha_EM</ei> and <ei>alpha_strong</ei> were used
218 (cf. <code>IDPRUP, XWTGUP, SCALUP, AQEDUP, AQCDUP</code>). This method
219 also resets the size of the particle list, and adds the empty zeroth
220 line, so it has to be called before the <code>addParticle</code> method below.
222 <method name="addParticle( id, status, mother1, mother2, colourTag1,
223 colourTag2, p_x, p_y, p_z, e, m, tau, spin)">
224 gives the properties of the next particle handed in (cf. <code>IDUP, ISTUP,
225 MOTHUP(1,..), MOTHUP(2,..), ICOLUP(1,..), ICOLUP(2,..), PUP(J,..),
226 VTIMUP, SPINUP</code>) .
230 Information is handed back by the following methods:
231 <method name="idProcess(), weight(), scale(), alphaQED(), alphaQCD()">.
232 Note that the weight stored in <code>pythia.info.weight()</code> as a rule
233 is not the same as the above <code>weight()</code>: the method here gives
234 the value before unweighting while the one in <code>info</code> gives
235 the one after unweighting and thus normally is 1 or -1. Only with strategy
236 options +-3 and +-4 would the value in <code>info</code> be the same as
237 here, except for a conversion from pb to mb for +-4.
239 <method name="sizePart()">
240 for the size of the particle array, which is one larger than the number
241 of particles in the event, since the zeroth entry is kept empty
242 (see above). Thus a typical loop would be
243 <br/><code>for (int i = 1; i < sizePart(); ++i) {...}</code>
245 <method name="id(i), status(i), mother1(i), mother2(i), col1(i), col2(i),
246 px(i), py(i), pz(i), e(i), m(i), tau(i), spin(i)">
247 for particle <code>i</code> in the range
248 <code>0 <= i < size()</code>. (But again note that
249 <code>i = 0</code> is an empty line, so the true range begins at 1.)
253 In the LHEF description <ref>Alw06</ref> an extension to
254 include information on the parton densities of the colliding partons
255 is suggested. This optional further information can be set by
256 <method name="setPdf( id1, id2, x1, x2, scalePDF, xpdf1, xpdf2)">
257 which gives the flavours , the <ei>x</ei> and the <ie>Q</ei> scale
258 (in GeV) at which the parton densities <ei>x*f_i(x, Q)</ei> have been
263 This information is returned by the methods
264 <method name="pdfIsSet(), id1(), id2(), x1(), x2(), scalePDF(), xpdf1(), xpdf2()">
265 where the first one tells whether this optional information has been set
266 for the current event. (<code>setPdf(...)</code> must be called after the
267 <code>setProcess(...)</code> call of the event for this to work.)
271 The information can also be printed using the <code>listEvent()</code>
272 method, e.g. <code>LHAupObject->listEvent()</code>.
273 In cases where the <code>LHAupObject</code> is not available to the
274 user, the <code>pythia.LHAeventList()</code> method can be used, which
275 is a wrapper for the above.
278 The LHA expects the decay of resonances to be included as part of the
279 hard process, i.e. if unstable particles are produced in a process then
280 their decays are also described. This includes <ei>Z^0, W^+-, H^0</ei>
281 and other short-lived particles in models beyond the Standard Model.
282 Should this not be the case then PYTHIA will perform the decays of all
283 resonances it knows how to do, in the same way as for internal processes.
284 Note that you will be on slippery ground if you then restrict the decay of
285 these resonances to specific allowed channels since, if this is not what
286 was intended, you will obtain the wrong cross section and potentially the
287 wrong mix of different event types. (Since the original intention is
288 unknown, the cross section will not be corrected for the fraction of
289 open channels, i.e. the procedure used for internal processes is not
290 applied in this case.)
292 <h3>An interface to Les Houches Event Files</h3>
294 The LHEF standard <ref>Alw06</ref> specifies a format where a single file
295 packs initialization and event information. This has become the most
296 frequently used procedure to process external parton-level events in Pythia.
298 <aloc href="ProgramFlow"><code>pythia.init(fileName)</code></aloc>
299 initialization option exists, where the LHEF name is provided as input.
300 Internally this name is then used to create an instance of the derived
301 class <code>LHAupLHEF</code>, which can do the job of reading an LHEF.
304 An example how to generate events from an LHEF is found in
305 <code>main12.cc</code>. Note the use of
306 <code>pythia.info.atEndOfFile()</code> to find out when the whole
307 LHEF has been processed.
310 To allow the sequential use of several event files the
311 <code>init(...)</code> method has an optional second argument:
312 <code>pythia.init(fileName, bool skipInit = false)</code>.
313 If called with this argument <code>true</code> then there will be no
314 initialization, except that the existing <code>LHAupLHEF</code> class
315 instance will be deleted and replaced by ones pointing to the new file.
316 It is assumed (but never checked) that the initialization information is
317 identical, and that the new file simply contains further events of
318 exactly the same kind as the previous one. An example of this possibility,
319 and the option to mix with internal processes, is found in
320 <code>main13.cc</code>.
322 <h3>A runtime Fortran interface</h3>
324 The runtime Fortran interface requires linking to an external Fortran
325 code. In order to avoid problems with unresolved external references
326 when this interface is not used, the code has been put in a separate
327 <code>LHAFortran.h</code> file, that is not included in any of the
328 other library files. Instead it should be included in the
329 user-supplied main program, together with the implementation of two
330 methods below that call the Fortran program to do its part of the job.
333 The <code>LHAupFortran</code> class derives from <code>LHAup</code>.
334 It reads initialization and event information from the LHA standard
335 Fortran commonblocks, assuming these commonblocks behave like two
336 <code>extern "C" struct</code> named <code>heprup_</code> and
337 <code>hepeup_</code>. (Note the final underscore, to match how the
338 gcc compiler internally names Fortran files.)
341 The instantiation does not require any arguments.
344 The user has to supply implementations of the <code>fillHepRup()</code>
345 and <code>fillHepEup()</code> methods, that is to do the actual calling
346 of the external Fortran routines that fill the <code>HEPRUP</code> and
347 <code>HEPEUP</code> commonblocks. The translation of this information to
348 the C++ structure is provided by the existing <code>setInit()</code> and
349 <code>setEvent()</code> code.
353 <aloc href="AccessPYTHIA6Processes">here</aloc> for information how
354 PYTHIA 6.4 can be linked to make use of this facility.
356 <h3>Methods for LHEF output</h3>
358 The main objective of the <code>LHAup</code> class is to feed information
359 from an external program into PYTHIA. It can be used to export information
360 as well, however. Specifically, there are four routines in the base class
361 that can be called to write a Les Houches Event File. These should be
362 called in sequence in order to build up the proper file structure.
364 <method name="openLHEF(string filename)">
365 Opens a file with the filename indicated, and writes a header plus a brief
366 comment with date and time information.
369 <method name="initLHEF()">
370 Writes initialization information to the file above. Such information should
371 already have been set with the methods described in the "Initialization"
375 <method name="eventLHEF()">
376 Writes event information to the file above. Such information should
377 already have been set with the methods described in the "Event input"
378 section above. This call should be repeated once for each event to be
382 <method name="closeLHEF(bool updateInit = false)">
383 Writes the closing tag and closes the file. Optionally, if
384 <code>updateInit = true</code>, this routine will reopen the file from
385 the beginning, rewrite the same header as <code>openLHEF()</code> did,
386 and then call <code>initLHEF()</code> again to overwrite the old
387 information. This is especially geared towards programs, such as PYTHIA
388 itself, where the cross section information is not available at the
389 beginning of the run, but only is obtained by Monte Carlo integration
390 in parallel with the event generation itself. Then the
391 <code>setXSec( i, xSec)</code>, <code>setXErr( i, xSec)</code> and
392 <code>setXMax( i, xSec)</code> can be used to update the relevant
393 information before <code>closeLHEF</code> is called.
394 <note>Warning:</note> overwriting the beginning of a file without
395 upsetting anything is a delicate operation. It only works when the new
396 lines require exactly as much space as the old ones did. Thus, if you add
397 another process in between, the file will be corrupted.
400 <h3>PYTHIA 8 output to an LHEF</h3>
402 The above methods could be used by any program to write an LHEF.
403 For PYTHIA 8 to do this, a derived class already exists,
404 <code>LHAupFromPYTHIA8</code>. In order for it to do its job,
405 it must gain access to the information produced by PYTHIA,
406 specifically the <code>process</code> event record and the
407 generic information stored in <code>info</code>. Therefore, if you
408 are working with an instance <code>pythia</code> of the
409 <code>Pythia</code> class, you have to instantiate
410 <code>LHAupFromPYTHIA8<code> with pointers to the
411 <code>process</code> and <code>info</code> objects of
413 <br/>LHAupFromPYTHIA8 myLHA(&pythia.process, &pythia.info);
416 The method <code>setInit()</code> should be called to store the
417 <code>pythia</code> initialization information in the LHA object,
418 and <code>setEvent()</code> to store event information.
419 Furthermore, <code>updateSigma()</code> can be used at the end
420 of the run to update cross-section information, cf.
421 <code>closeLHEF(true)</code> above. An example how the
422 generation, translation and writing methods should be ordered is
423 found in <code>main20.cc</code>.
426 Currently there are some limitations, that could be overcome if
427 necessary. Firstly, you may mix many processes in the same run,
428 but the cross-section information stored in <code>info</code> only
429 refers to the sum of them all, and therefore they are all classified
430 as a common process 9999. Secondly, you should generate your events
431 in the CM frame of the collision, since this is the assumed frame of
432 stored Les Houches events, and no boosts have been implemented
433 for the case that <code>pythia.process</code> is not in this frame.
437 <!-- Copyright (C) 2008 Torbjorn Sjostrand -->