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 You are free to write your own derived classes, using the rules and
22 methods to be described below. Normally, pointers to objects of such
23 derived classes should be handed in with the
24 <code><aloc href="ProgramFlow">Pythia::init( LHAup*)</aloc></code>
25 method. However, with the LHEF format a filename can replace the
26 pointer, see further below.
29 Let us now describe the methods at your disposal to do the job.
31 <method name="LHAup::LHAup( int strategy = 3)">
32 the base class constructor takes the choice of mixing/weighting
33 strategy as optional input argument, and calls <code>setStrategy</code>,
34 see below. It also reserves some space for processes and particles.
37 <method name="virtual LHAup::~LHAup()">
38 the destructor does not need to do anything.
41 <method name="void LHAup::setPtr(Info* infoPtr)">
42 this method only sets the pointer that allows some information
43 to be accessed, and is automatically called by
44 <code>Pythia::init(...)</code>.
47 <h3>Initialization</h3>
49 The <code>LHAup</code> class stores information equivalent to the
50 <code>/HEPRUP/</code> commonblock, as required to initialize the event
51 generation chain. The main difference is that the vector container
52 now allows a flexible number of subprocesses to be defined. For the
53 rest, names have been modified, since the 6-character-limit does not
54 apply, and variables have been regrouped for clarity, but nothing
55 fundamental is changed.
57 <method name="virtual bool LHAup::setInit()">
58 this pure virtual method has to be implemented in the derived class,
59 to set relevant information when called. It should return false if it
60 fails to set the info.
64 Inside <code>setInit()</code>, such information can be set by the following
66 <method name="void LHAup::setBeamA( int identity, double energy,
67 int pdfGroup, int pdfSet)">
69 <methodmore name="void LHAup::setBeamB( int identity, double energy,
70 int pdfGroup, int pdfSet)">
71 sets the properties of the first and second incoming beam, respectively
72 (cf. the Fortran <code>IDBMUP(1), EBMUP(i), PDFGUP(i), PDFSUP(i)</code>,
73 with <code>i</code> 1 or 2). The parton distribution information
74 defaults to zero. These numbers can be used to tell which PDF sets were
75 used when the hard process was generated, while the normal
76 <aloc href="PDFSelection">PDF Selection</aloc> is used for the further
77 event generation in PYTHIA.
80 <method name="void LHAup::setStrategy( int strategy)">
81 sets the event weighting and cross section strategy. The default,
82 provided in the class constructor, is 3, which is the natural value
84 <argument name="strategy">
85 chosen strategy (cf. <code>IDWTUP</code>; see <ref>Sjo06</ref>
86 section 9.9.1 for extensive comments).
87 <argoption value="1"> events come with non-negative weight, given in units
88 of pb, with an average that converges towards the cross section of the
89 process. PYTHIA is in charge of the event mixing, i.e. for each new
90 try decides which process should be generated, and then decides whether
91 is should be kept, based on a comparison with <code>xMax</code>.
92 Accepted events therefore have unit weight.</argoption>
93 <argoption value="-1"> as option 1, except that cross sections can now be
94 negative and events after unweighting have weight +-1. You can use
95 <code><aloc href="EventInformation">Info::weight()</aloc></code>
96 to find the weight of the current event. A correct event mixing requires
97 that a process that can take both signs should be split in two, one limited
98 to positive or zero and the other to negative or zero values, with
99 <code>xMax</code> chosen appropriately for the two.</argoption>
100 <argoption value="2"> events come with non-negative weight, in unspecified
101 units, but such that <code>xMax</code> can be used to unweight the events
102 to unit weight. Again PYTHIA is in charge of the event mixing.
103 The total cross section of a process is stored in
104 <code>xSec</code>.</argoption>
105 <argoption value="-2"> as option 2, except that cross sections can now be
106 negative and events after unweighting have weight +-1. As for option -1
107 processes with indeterminate sign should be split in two.</argoption>
108 <argoption value="3"> events come with unit weight, and are thus accepted
109 as is. The total cross section of the process is stored in
110 <code>xSec</code>.</argoption>
111 <argoption value="-3"> as option 3, except that events now come with weight
112 +-1. Unlike options -1 and -2 processes with indeterminate sign need not be
113 split in two, unless you intend to mix with internal PYTHIA processes
114 (see below).</argoption>
115 <argoption value="4"> events come with non-negative weight, given in units
116 of pb, with an average that converges towards the cross section of the
117 process, like for option 1. No attempt is made to unweight the events,
118 however, but all are generated in full, and retain their original weight.
119 For consistency with normal PYTHIA units, the weight stored in
120 <code>Info::weight()</code> has been converted to mb, however.
122 <argoption value="-4"> as option 4, except that events now can come
123 either with positive or negative weights.</argoption>
124 <note>Note 1</note>: if several processes have already been mixed and
125 stored in a common event file, either LHEF or some private format, it
126 would be problematical to read back events in a different order. Since it
127 is then not feasible to let PYTHIA pick the next process type, strategies
128 +-1 and +-2 would not work. Instead strategy 3 would be the recommended
129 choice, or -3 if negative-weight events are required.
130 <note>Note 2</note>: it is possible to switch on internally implemented
131 processes and have PYTHIA mix these with LHA ones according to their relative
132 cross sections for strategies +-1, +-2 and 3. It does not work for strategy
133 -3 unless the positive and negative sectors of the cross sections are in
134 separate subprocesses (as must always be the case for -1 and -2), since
135 otherwise the overall mixture of PYTHIA and LHA processes will be off.
136 Mixing is not possible for strategies +-4, since the weighting procedure
137 is not specified by the standard. (For instance, the intention may be to
138 have events biased towards larger <ei>pT</ei> values in some particular
143 <method name="void LHAup::addProcess( int idProcess, double xSec,
144 double xErr, double xMax)">
145 sets info on an allowed process (cf. <code>LPRUP, XSECUP, XERRUP,
147 Each new call will append one more entry to the list of processes.
148 The choice of strategy determines which quantities are mandatory:
149 <code>xSec</code> for strategies +-2 and +-3,
150 <code>xErr</code> never, and
151 <code>xMax</code> for strategies +-1 and +-2.
154 <note>Note</note>: PYTHIA does not make active use of the (optional)
155 <code>xErr</code> values, but calculates a statistical cross section
156 error based on the spread of event-to-event weights. This should work
157 fine for strategy options +-1, but not for the others. Specifically,
158 for options +-2 and +-3 the weight spread may well vanish, and anyway
159 is likely to be an underestimate of the true error. If the author of the
160 LHA input information does provide error information you may use that -
161 this information is displayed at initialization. If not, then a relative
162 error decreasing like <ei>1/sqrt(n_acc)</ei>, where <ei>n_acc</ei>
163 is the number of accepted events, should offer a reasonable estimate.
165 <method name="void LHAup::setXSec( int i, double xSec)">
166 update the <code>xSec</code> value of the <code>i</code>'th process
167 added with <code>addProcess</code> method (i.e. <code>i</code> runs
168 from 0 through <code>sizeProc() - 1</code>, see below).
171 <method name="void LHAup::setXErr( int i, double xErr)">
172 update the <code>xErr</code> value of the <code>i</code>'th process
173 added with <code>addProcess</code> method.
176 <method name="void LHAup::setXMax( int i, double xMax)">
177 update the <code>xMax</code> value of the <code>i</code>'th process
178 added with <code>addProcess</code> method.
182 Information is handed back by the following methods
183 (that normally you would not need to touch):
184 <method name="int LHAup::idBeamA()">
186 <methodmore name="int LHAup::idBeamB()">
188 <methodmore name="double LHAup::eBeamA()">
190 <methodmore name="double LHAup::eBeamB()">
192 <methodmore name="int LHAup::pdfGroupBeamA()">
194 <methodmore name="int LHAup::pdfGroupBeamB()">
196 <methodmore name="int LHAup::pdfSetBeamA()">
198 <methodmore name="int LHAup::pdfSetBeamB()">
199 for the beam properties.
201 <method name="int LHAup::strategy()">
202 for the strategy choice.
204 <method name="int LHAup::sizeProc()">
205 for the number of subprocesses.
207 <method name="int LHAup::idProcess(i)">
209 <methodmore name="double LHAup::xSec(i)">
211 <methodmore name="double LHAup::xErr(i)">
213 <methodmore name="double LHAup::xMax(i)">
214 for process <code>i</code> in the range <code>0 <= i <
218 <method name="void LHAup::listInit(ostream& os = cout)">
219 prints the above initialization information. This method is
220 automatically called from <code>Pythia::init(...)</code>,
221 so would normally not need to be called directly by the user.
226 The <code>LHAup</code> class also stores information equivalent to the
227 <code>/HEPEUP/</code> commonblock, as required to hand in the next
228 parton-level configuration for complete event generation. The main
229 difference is that the vector container now allows a flexible number
230 of partons to be defined. For the rest, names have been modified,
231 since the 6-character-limit does not apply, and variables have been
232 regrouped for clarity, but nothing fundamental is changed.
235 The LHA standard is based on Fortran arrays beginning with
236 index 1, and mother information is defined accordingly. In order to
237 be compatible with this convention, the zeroth line of the C++ particle
238 array is kept empty, so that index 1 also here corresponds to the first
239 particle. One small incompatibility is that the <code>sizePart()</code>
240 method returns the full size of the particle array, including the
241 empty zeroth line, and thus is one larger than the true number of
242 particles (<code>NUP</code>).
244 <method name="virtual bool LHAup::setEvent(int idProcess = 0)">
245 this pure virtual method has to be implemented in the derived class,
246 to set relevant information when called. For strategy options +-1
247 and +-2 the input <code>idProcess</code> value specifies which process
248 that should be generated, while <code>idProcess</code> is irrelevant
249 for strategies +-3 and +-4. The method should return false if it fails
250 to set the info, i.e. normally that the supply of events in a file is
251 exhausted. If so, no event is generated, and <code>Pythia::next()</code>
252 returns false. You can then interrogate
253 <code><aloc href="EventInformation">Info::atEndOfFile()</aloc></code>
254 to confirm that indeed the failure is caused in this method, and decide
255 to break out of the event generation loop.
258 Inside a normal <code>setEvent(...)</code> call, information can be set
259 by the following methods:
260 <method name="void LHAup::setProcess( int idProcess, double weight,
261 double scale, double alphaQED, double alphaQCD)">
262 tells which kind of process occured, with what weight, at what scale,
263 and which <ei>alpha_EM</ei> and <ei>alpha_strong</ei> were used
264 (cf. <code>IDPRUP, XWTGUP, SCALUP, AQEDUP, AQCDUP</code>). This method
265 also resets the size of the particle list, and adds the empty zeroth
266 line, so it has to be called before the <code>addParticle</code> method below.
268 <method name="void LHAup::addParticle( int id, int status, int mother1,
269 int mother2, int colourTag1, int colourTag2, double p_x, double p_y,
270 double p_z, double e, double m, double tau, double spin)">
271 gives the properties of the next particle handed in (cf. <code>IDUP, ISTUP,
272 MOTHUP(1,..), MOTHUP(2,..), ICOLUP(1,..), ICOLUP(2,..), PUP(J,..),
273 VTIMUP, SPINUP</code>) .
277 Information is handed back by the following methods:
278 <method name="int LHAup::idProcess()">
282 <method name="double LHAup::weight()">.
283 Note that the weight stored in <code>Info::weight()</code> as a rule
284 is not the same as the above <code>weight()</code>: the method here gives
285 the value before unweighting while the one in <code>info</code> gives
286 the one after unweighting and thus normally is 1 or -1. Only with strategy
287 options +-3 and +-4 would the value in <code>info</code> be the same as
288 here, except for a conversion from pb to mb for +-4.
291 <method name="double LHAup::scale()">
293 <methodmore name="double LHAup::alphaQED()">
295 <methodmore name="double LHAup::alphaQCD()">
296 scale and couplings at that scale.
299 <method name="int LHAup::sizePart()">
300 the size of the particle array, which is one larger than the number
301 of particles in the event, since the zeroth entry is kept empty
305 <method name="int LHAup::id(int i)">
307 <methodmore name="int LHAup::status(int i)">
309 <methodmore name="int LHAup::mother1(int i)">
311 <methodmore name="int LHAup::mother2(int i)">
313 <methodmore name="int LHAup::col1(int i)">
315 <methodmore name="int LHAup::col2(int i)">
317 <methodmore name="double LHAup::px(int i)">
319 <methodmore name="double LHAup::py(int i)">
321 <methodmore name="double LHAup::pz(int i)">
323 <methodmore name="double LHAup::e(int i)">
325 <methodmore name="double LHAup::m(int i)">
327 <methodmore name="double LHAup::tau(int i)">
329 <methodmore name="double LHAup::spin(int i)">
330 for particle <code>i</code> in the range
331 <code>0 <= i < sizePart()</code>. (But again note that
332 <code>i = 0</code> is an empty line, so the true range begins at 1.)
336 In the LHEF description <ref>Alw06</ref> an extension to
337 include information on the parton densities of the colliding partons
338 is suggested. This optional further information can be set by
339 <method name="void LHAup::setPdf( int id1, int id2, double x1, double x2,
340 double scalePDF, double xpdf1, double xpdf2)">
341 which gives the flavours , the <ei>x</ei> and the <ie>Q</ei> scale
342 (in GeV) at which the parton densities <ei>x*f_i(x, Q)</ei> have been
347 This information is returned by the methods
348 <method name="bool LHAup::pdfIsSet()">
350 <methodmore name="int LHAup::id1()">
352 <methodmore name="int LHAup::id2()">
354 <methodmore name="double LHAup::x1()">
356 <methodmore name="double LHAup::x2()">
358 <methodmore name="double LHAup::scalePDF()">
360 <methodmore name="double LHAup::xpdf1()">
362 <methodmore name="double LHAup::xpdf2()">
363 where the first one tells whether this optional information has been set
364 for the current event. (<code>setPdf(...)</code> must be called after the
365 <code>setProcess(...)</code> call of the event for this to work.)
369 <method name="void LHAup::listEvent(ostream& os = cout)">
370 prints the above information for the current event. In cases where the
371 <code>LHAup</code> object is not available to the user, the
372 <code>Pythia::LHAeventList(ostream& os = cout)</code> method can
373 be used, which is a wrapper for the above.
376 <method name="virtual bool LHAup::skipEvent(int nSkip)">
377 skip ahead <code>nSkip</code> events in the Les Houches generation
378 sequence, without doing anything further with them. Mainly
379 intended for debug purposes, e.g. when an event at a known
380 location in a Les Houches Event File is causing problems.
381 Will return false if operation fails, specifically if the
382 end of an LHEF has been reached. The implementation in the base class
383 simply executes <code>setEvent()</code> the requested number of times.
384 The derived <code>LHAupLHEF</code> class (see below) only uses the
385 <code>setNewEventLHEF(...)</code> part of its <code>setEvent()</code>
386 method, and other derived classes could choose other shortcuts.
390 The LHA expects the decay of resonances to be included as part of the
391 hard process, i.e. if unstable particles are produced in a process then
392 their decays are also described. This includes <ei>Z^0, W^+-, H^0</ei>
393 and other short-lived particles in models beyond the Standard Model.
394 Should this not be the case then PYTHIA will perform the decays of all
395 resonances it knows how to do, in the same way as for internal processes.
396 Note that you will be on slippery ground if you then restrict the decay of
397 these resonances to specific allowed channels since, if this is not what
398 was intended, you will obtain the wrong cross section and potentially the
399 wrong mix of different event types. (Since the original intention is
400 unknown, the cross section will not be corrected for the fraction of
401 open channels, i.e. the procedure used for internal processes is not
402 applied in this case.)
404 <h3>An interface to Les Houches Event Files</h3>
406 The LHEF standard <ref>Alw06</ref> specifies a format where a single file
407 packs initialization and event information. This has become the most
408 frequently used procedure to process external parton-level events in
409 Pythia. Therefore a special
410 <code><aloc href="ProgramFlow">Pythia::init(fileName)</aloc></code>
411 initialization option exists, where the LHEF name is provided as input.
412 Internally this name is then used to create an instance of the derived
413 class <code>LHAupLHEF</code>, which can do the job of reading an LHEF.
416 An example how to generate events from an LHEF is found in
417 <code>main12.cc</code>. Note the use of
418 <code>Info::atEndOfFile()</code> to find out when the whole
419 LHEF has been processed.
422 To allow the sequential use of several event files the
423 <code>Pythia::init(...)</code> method has an optional second argument:
424 <code>Pythia::init(fileName, bool skipInit = false)</code>.
425 If called with this argument <code>true</code> then there will be no
426 initialization, except that the existing <code>LHAupLHEF</code> class
427 instance will be deleted and replaced by ones pointing to the new file.
428 It is assumed (but never checked) that the initialization information is
429 identical, and that the new file simply contains further events of
430 exactly the same kind as the previous one. An example of this possibility,
431 and the option to mix with internal processes, is found in
432 <code>main13.cc</code>.
435 The workhorses of the <code>LHAupLHEF</code> class are three methods
436 found in the base class, so as to allow them to be reused in other
437 contexts. Specifically, it allows derived classes where one parton-level
438 configuration can be reused several times, e.g. in the context of
439 matrix-element-to-parton-shower matching (example in preparation).
440 To begin with also a small utility routine.
442 <method name="bool LHAup::fileFound()">
443 always returns true in the base class, but in <code>LHAupLHEF</code>
444 it returns false if the LHEF provided in the constructor is not
445 found and opened correctly.
448 <method name="bool LHAup::setInitLHEF(ifstream& is)">
449 read in and set all required initialization information from the
450 specified stream. Return false if it fails.
453 <method name="bool LHAup::setNewEventLHEF(ifstream& is)">
454 read in event information from the specified stream into a staging area
455 where it can be reused by <code>setOldEventLHEF</code>.
458 <method name="bool LHAup::setOldEventLHEF()">
459 store the event information from the staging area into the normal
460 location. Thus a single <code>setNewEventLHEF</code> call can be
461 followed by several <code>setOldEventLHEF</code> ones, so as to
462 process the same configuration several times. This method currently
463 only returns true, i.e. any errors should be caught by the preceding
464 <code>setNewEventLHEF</code> call.
468 <h3>A runtime Fortran interface</h3>
470 The runtime Fortran interface requires linking to an external Fortran
471 code. In order to avoid problems with unresolved external references
472 when this interface is not used, the code has been put in a separate
473 <code>LHAFortran.h</code> file, that is not included in any of the
474 other library files. Instead it should be included in the
475 user-supplied main program, together with the implementation of two
476 methods below that call the Fortran program to do its part of the job.
479 The <code>LHAupFortran</code> class derives from <code>LHAup</code>.
480 It reads initialization and event information from the LHA standard
481 Fortran commonblocks, assuming these commonblocks behave like two
482 <code>extern "C" struct</code> named <code>heprup_</code> and
483 <code>hepeup_</code>. (Note the final underscore, to match how the
484 gcc compiler internally names Fortran files.)
487 The instantiation does not require any arguments.
490 The user has to supply implementations of the <code>fillHepRup()</code>
491 and <code>fillHepEup()</code> methods, that is to do the actual calling
492 of the external Fortran routines that fill the <code>HEPRUP</code> and
493 <code>HEPEUP</code> commonblocks. The translation of this information to
494 the C++ structure is provided by the existing <code>setInit()</code> and
495 <code>setEvent()</code> code.
498 Up to and including version 8.125 the <code>LHAupFortran</code> class
499 was used to construct a runtime interface to PYTHIA 6.4. This was
500 convenient in the early days of PYTHIA 8 evolution, when this program
501 did not yet contain hard-process generation, and the LHEF standard
502 did not yet exist. Nowadays it is more of a bother, since a full
503 cross-platform support leads to many possible combinations. Therefore
504 the support has been reduced in the current version. Only the
505 <code>main51.cc</code> example remains as an illustration, where the
506 previously separate interface code
507 (<code>include/Pythia6Interface.h</code>) has been inserted in the
508 beginning. You also need to modify the <code>examples/Makefile</code>
509 to link <code>main51.cc</code> properly also to a PYTHIA 6.4 library
510 version, see commented-out section for ideas how to to this.
512 <h3>Methods for LHEF output</h3>
514 The main objective of the <code>LHAup</code> class is to feed information
515 from an external program into PYTHIA. It can be used to export information
516 as well, however. Specifically, there are four routines in the base class
517 that can be called to write a Les Houches Event File. These should be
518 called in sequence in order to build up the proper file structure.
520 <method name="bool LHAup::openLHEF(string filename)">
521 Opens a file with the filename indicated, and writes a header plus a brief
522 comment with date and time information.
525 <method name="bool LHAup::initLHEF()">
526 Writes initialization information to the file above. Such information should
527 already have been set with the methods described in the "Initialization"
531 <method name="bool LHAup::eventLHEF()">
532 Writes event information to the file above. Such information should
533 already have been set with the methods described in the "Event input"
534 section above. This call should be repeated once for each event to be
538 <method name="bool LHAup::closeLHEF(bool updateInit = false)">
539 Writes the closing tag and closes the file. Optionally, if
540 <code>updateInit = true</code>, this routine will reopen the file from
541 the beginning, rewrite the same header as <code>openLHEF()</code> did,
542 and then call <code>initLHEF()</code> again to overwrite the old
543 information. This is especially geared towards programs, such as PYTHIA
544 itself, where the cross section information is not available at the
545 beginning of the run, but only is obtained by Monte Carlo integration
546 in parallel with the event generation itself. Then the
547 <code>setXSec( i, xSec)</code>, <code>setXErr( i, xSec)</code> and
548 <code>setXMax( i, xSec)</code> can be used to update the relevant
549 information before <code>closeLHEF</code> is called.
550 <note>Warning:</note> overwriting the beginning of a file without
551 upsetting anything is a delicate operation. It only works when the new
552 lines require exactly as much space as the old ones did. Thus, if you add
553 another process in between, the file will be corrupted.
556 <h3>PYTHIA 8 output to an LHEF</h3>
558 The above methods could be used by any program to write an LHEF.
559 For PYTHIA 8 to do this, a derived class already exists,
560 <code>LHAupFromPYTHIA8</code>. In order for it to do its job,
561 it must gain access to the information produced by PYTHIA,
562 specifically the <code>process</code> event record and the
563 generic information stored in <code>info</code>. Therefore, if you
564 are working with an instance <code>pythia</code> of the
565 <code>Pythia</code> class, you have to instantiate
566 <code>LHAupFromPYTHIA8</code> with pointers to the
567 <code>process</code> and <code>info</code> objects of
569 <br/><code>LHAupFromPYTHIA8 myLHA(&pythia.process, &pythia.info);</code>
572 The method <code>setInit()</code> should be called to store the
573 <code>pythia</code> initialization information in the LHA object,
574 and <code>setEvent()</code> to store event information.
575 Furthermore, <code>updateSigma()</code> can be used at the end
576 of the run to update cross-section information, cf.
577 <code>closeLHEF(true)</code> above. An example how the
578 generation, translation and writing methods should be ordered is
579 found in <code>main20.cc</code>.
582 Currently there are some limitations, that could be overcome if
583 necessary. Firstly, you may mix many processes in the same run,
584 but the cross-section information stored in <code>info</code> only
585 refers to the sum of them all, and therefore they are all classified
586 as a common process 9999. Secondly, you should generate your events
587 in the CM frame of the collision, since this is the assumed frame of
588 stored Les Houches events, and no boosts have been implemented
589 for the case that <code>Pythia::process</code> is not in this frame.
593 <!-- Copyright (C) 2010 Torbjorn Sjostrand -->