using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / xmldoc / LesHouchesAccord.xml
1 <chapter name="Les Houches Accord">
2
3 <h2>Les Houches Accord</h2>
4
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.
10
11 <p/>
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.
19
20 <p/>
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, 
24 see below. 
25
26 <h3>Initialization</h3>
27
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.
35
36 <p/>
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.
40
41 <p/>
42 Inside <code>setInit()</code>, such information can be set by the following 
43 methods:
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.
52 </method>
53
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 
57 e.g. for an LHEF.
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.
95 </argoption>
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 
113 functional form.)
114 </argument>
115 </method>
116
117 <method name="addProcess( idProcess, xSec, xErr, xMax)"> 
118 sets info on an allowed process (cf. <code>LPRUP, XSECUP, XERRUP, 
119 XMAXUP</code>). 
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. 
125 </method>
126
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.
137
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).
142 </method>
143
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.
147 </method>
148
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.
152 </method>
153
154 <p/>
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.
158 </method>
159 <method name="strategy()">
160 for the strategy choice.
161 </method>
162 <method name="sizeProc()"> 
163 for the number of subprocesses.
164 </method>
165 <method name="idProcess(i), xSec(i), xErr(i), xMax(i)"> 
166 for process <code>i</code> in the range <code>0 &lt;= i &lt; 
167 sizeProc()</code>.   
168 </method>
169
170 <p/>
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.
174
175 <h3>Event input</h3>
176
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.
184
185 <p/>
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>). 
194
195 <p/>
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 
206 interrogate 
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.
211
212 <p/>
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.
221 </method>
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>) .
227 </method>
228
229 <p/>
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. 
238 </method>
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>
244 </method>
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 &lt;= i &lt; size()</code>. (But again note that 
249 <code>i = 0</code> is an empty line, so the true range begins at 1.)   
250 </method>
251
252 <p/>
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
259 evaluated.
260 </method>
261
262 <p/>
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.)
268 </method>
269
270 <p/>
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. 
276
277 <p/>
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.) 
291
292 <h3>An interface to Les Houches Event Files</h3>
293
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. 
297 Therefore a special 
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.
302
303 <p/>
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. 
308
309 <p/>
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>.
321
322 <h3>A runtime Fortran interface</h3>
323
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. 
331
332 <p/>
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.) 
339
340 <p/>
341 The instantiation does not require any arguments. 
342
343 <p/>
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. 
350
351 <p/>
352 See further 
353 <aloc href="AccessPYTHIA6Processes">here</aloc> for information how 
354 PYTHIA 6.4 can be linked to make use of this facility. 
355
356 <h3>Methods for LHEF output</h3>
357
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. 
363
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.
367 </method>
368
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"
372 section above.
373 </method>
374
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 
379 stored. 
380 </method>
381
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. 
398 </method>
399
400 <h3>PYTHIA 8 output to an LHEF</h3>
401
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 
412 <code>pythia</code>:
413 <br/>LHAupFromPYTHIA8 myLHA(&pythia.process, &pythia.info);
414
415 <p/>
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>.
424
425 <p/>
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. 
434  
435 </chapter>
436
437 <!-- Copyright (C) 2008 Torbjorn Sjostrand -->