Les Houches Accord
The Les Houches Accord (LHA) for user processes [Boo01] is the
standard way to input parton-level information from a
matrix-elements-based generator into PYTHIA. The conventions for
which information should be stored has been defined in a Fortran context,
as two commonblocks. Here a C++ equivalent is defined, as a single class.
The LHAup
class is a base class, containing reading and
printout functions, plus two pure virtual functions, one to set
initialization information and one to set information on each new event.
Derived classes have to provide these two virtual functions to do
the actual work. The existing derived classes are for reading information
from a Les Houches Event File (LHEF), from the respective Fortran
commonblocks, or from PYTHIA 8 itself.
Normally, pointers to objects of the derived classes should be handed in
with the pythia.init( LHAup*)
method. However, with the LHEF format a filename can replace the pointer,
see below.
Initialization
The LHAup
class stores information equivalent to the
/HEPRUP/
commonblock, as required to initialize the event
generation chain. The main difference is that the vector container
now allows a flexible number of subprocesses to be defined. For the
rest, names have been modified, since the 6-character-limit does not
apply, and variables have been regrouped for clarity, but nothing
fundamental is changed.
The pure virtual function setInit()
has to be implemented in
the derived class, to set relevant information when called. It should
return false
if it fails to set the info.
Inside setInit()
, such information can be set by the following
methods:
method
setBeamA( identity, energy, pdfGroup, pdfSet)
sets the properties of the first incoming beam (cf. the Fortran
IDBMUP(1), EBMUP(1), PDFGUP(1), PDFSUP(1)
), and similarly
a setBeamB
method exists. The parton distribution information
defaults to zero. These numbers can be used to tell which PDF sets were
used when the hard process was generated, while the normal
PDF Selection is used for the further
event generation in PYTHIA.
method
setStrategy( strategy)
sets the event weighting and cross section strategy. The default,
provided in the class constructor, is 3, which is the natural value
e.g. for an LHEF.
argument
strategy :
chosen strategy (cf. IDWTUP
; see [Sjo06]
section 9.9.1 for extensive comments).
argumentoption
1 : events come with non-negative weight, given in units
of pb, with an average that converges towards the cross section of the
process. PYTHIA is in charge of the event mixing, i.e. for each new
try decides which process should be generated, and then decides whether
is should be kept, based on a comparison with xMax
.
Accepted events therefore have unit weight.
argumentoption
-1 : as option 1, except that cross sections can now be
negative and events after unweighting have weight +-1. You can use
pythia.info.weight()
to find the weight of the current event. A correct event mixing requires
that a process that can take both signs should be split in two, one limited
to positive or zero and the other to negative or zero values, with
xMax
chosen appropriately for the two.
argumentoption
2 : events come with non-negative weight, in unspecified
units, but such that xMax
can be used to unweight the events
to unit weight. Again PYTHIA is in charge of the event mixing.
The total cross section of a process is stored in
xSec
.
argumentoption
-2 : as option 2, except that cross sections can now be
negative and events after unweighting have weight +-1. As for option -1
processes with indeterminate sign should be split in two.
argumentoption
3 : events come with unit weight, and are thus accepted
as is. The total cross section of the process is stored in
xSec
.
argumentoption
-3 : as option 3, except that events now come with weight
+-1. Unlike options -1 and -2 processes with indeterminate sign need not be
split in two, unless you intend to mix with internal PYTHIA processes
(see below).
argumentoption
4 : events come with non-negative weight, given in units
of pb, with an average that converges towards the cross section of the
process, like for option 1. No attempt is made to unweight the events,
however, but all are generated in full, and retain their original weight.
For consistency with normal PYTHIA units, the weight stored in
pythia.info.weight()
has been converted to mb, however.
argumentoption
-4 : as option 4, except that events now can come
either with positive or negative weights.
Note 1: if several processes have already been mixed and
stored in a common event file, either LHEF or some private format, it
would be problematical to read back events in a different order. Since it
is then not feasible to let PYTHIA pick the next process type, strategies
+-1 and +-2 would not work. Instead strategy 3 would be the recommended
choice, or -3 if negative-weight events are required.
Note 2: it is possible to switch on internally implemented
processes and have PYTHIA mix these with LHA ones according to their relative
cross sections for strategies +-1, +-2 and 3. It does not work for strategy
-3 unless the positive and negative sectors of the cross sections are in
separate subprocesses (as must always be the case for -1 and -2), since
otherwise the overall mixture of PYTHIA and LHA processes will be off.
Mixing is not possible for strategies +-4, since the weighting procedure
is not specified by the standard. (For instance, the intention may be to
have events biased towards larger pT values in some particular
functional form.)
method
addProcess( idProcess, xSec, xErr, xMax)
sets info on an allowed process (cf. LPRUP, XSECUP, XERRUP,
XMAXUP
).
Each new call will append one more entry to the list of processes.
The choice of strategy determines which quantities are mandatory:
xSec
for strategies +-2 and +-3,
xErr
never, and
xMax
for strategies +-1 and +-2.
Note: PYTHIA does not make active use of the (optional)
xErr
values, but calculates a statistical cross section
error based on the spread of event-to-event weights. This should work
fine for strategy options +-1, but not for the others. Specifically,
for options +-2 and +-3 the weight spread may well vanish, and anyway
is likely to be an underestimate of the true error. If the author of the
LHA input information does provide error information you may use that -
this information is displayed at initialization. If not, then a relative
error decreasing like 1/sqrt(n_acc), where n_acc
is the number of accepted events, should offer a reasonable estimate.
method
setXSec( i, xSec)
update the xSec
value of the i
'th process
added with addProcess
method (i.e. i
runs
from 0 through sizeProc() - 1
, see below).
method
setXErr( i, xErr)
update the xErr
value of the i
'th process
added with addProcess
method.
method
setXMax( i, xMax)
update the xMax
value of the i
'th process
added with addProcess
method.
Information is handed back by the following methods:
method
idBeamA(), eBeamA(), pdfGroupBeamA(), pdfSetBeamA()
and similarly with A -> B, for the two beam particles.
method
strategy()
for the strategy choice.
method
sizeProc()
for the number of subprocesses.
method
idProcess(i), xSec(i), xErr(i), xMax(i)
for process i
in the range 0 <= i <
sizeProc()
.
The information can also be printed using the listInit()
method, e.g. LHAupObject->listInit()
.
This is automatically done by the pythia.init
call.
Event input
The LHAup
class also stores information equivalent to the
/HEPEUP/
commonblock, as required to hand in the next
parton-level configuration for complete event generation. The main
difference is that the vector container now allows a flexible number
of partons to be defined. For the rest, names have been modified,
since the 6-character-limit does not apply, and variables have been
regrouped for clarity, but nothing fundamental is changed.
The LHA standard is based on Fortran arrays beginning with
index 1, and mother information is defined accordingly. In order to
be compatible with this convention, the zeroth line of the C++ particle
array is kept empty, so that index 1 also here corresponds to the first
particle. One small incompatibility is that the sizePart()
method returns the full size of the particle array, including the
empty zeroth line, and thus is one larger than the true number of
particles (NUP
).
The pure virtual function setEvent(idProcess)
has to be
implemented in the derived class, to set relevant information when
called. For strategy options +-1 and +-2 the input
idProcess
value specifies which process that should be
generated by setEvent(idProcess)
, while
idProcess
is irrelevant for strategies +-3 and +-4.
The setEvent(idProcess)
function should return
false
if it fails to set the info, i.e. normally that the
supply of events in a file is exhausted. If so, no event is generated,
and pythia.next()
returns false
. You can then
interrogate
pythia.info.atEndOfFile()
to confirm that indeed the failure is caused in the
setEvent(idProcess)
function, and decide to break out of
the event generation loop.
Inside a normal setEvent(...)
call, information can be set
by the following methods:
method
setProcess( idProcess, weight, scale, alphaQED, alphaQCD)
tells which kind of process occured, with what weight, at what scale,
and which alpha_EM and alpha_strong were used
(cf. IDPRUP, XWTGUP, SCALUP, AQEDUP, AQCDUP
). This method
also resets the size of the particle list, and adds the empty zeroth
line, so it has to be called before the addParticle
method below.
method
addParticle( id, status, mother1, mother2, colourTag1, colourTag2, p_x, p_y, p_z, e, m, tau, spin)
gives the properties of the next particle handed in (cf. IDUP, ISTUP,
MOTHUP(1,..), MOTHUP(2,..), ICOLUP(1,..), ICOLUP(2,..), PUP(J,..),
VTIMUP, SPINUP
) .
Information is handed back by the following methods:
method
idProcess(), weight(), scale(), alphaQED(), alphaQCD()
Note that the weight stored in pythia.info.weight()
as a rule
is not the same as the above weight()
: the method here gives
the value before unweighting while the one in info
gives
the one after unweighting and thus normally is 1 or -1. Only with strategy
options +-3 and +-4 would the value in info
be the same as
here, except for a conversion from pb to mb for +-4.
method
sizePart()
for the size of the particle array, which is one larger than the number
of particles in the event, since the zeroth entry is kept empty
(see above). Thus a typical loop would be
for (int i = 1; i < sizePart(); ++i) {...}
method
id(i), status(i), mother1(i), mother2(i), col1(i), col2(i),px(i), py(i), pz(i), e(i), m(i), tau(i), spin(i)
for particle i
in the range
0 <= i < size()
. (But again note that
i = 0
is an empty line, so the true range begins at 1.)
In the LHEF description [Alw06] an extension to
include information on the parton densities of the colliding partons
is suggested. This optional further information can be set by
method
setPdf( id1, id2, x1, x2, scalePDF, xpdf1, xpdf2)
which gives the flavours , the x and the Q scale
(in GeV) at which the parton densities x*f_i(x, Q) have been
evaluated.
This information is returned by the methods
method
pdfIsSet(), id1(), id2(), x1(), x2(), scalePDF(), xpdf1(), xpdf2()
where the first one tells whether this optional information has been set
for the current event. (setPdf(...)
must be called after the
setProcess(...)
call of the event for this to work.)
The information can also be printed using the listEvent()
method, e.g. LHAupObject->listEvent()
.
In cases where the LHAupObject
is not available to the
user, the pythia.LHAeventList()
method can be used, which
is a wrapper for the above.
The LHA expects the decay of resonances to be included as part of the
hard process, i.e. if unstable particles are produced in a process then
their decays are also described. This includes Z^0, W^+-, H^0
and other short-lived particles in models beyond the Standard Model.
Should this not be the case then PYTHIA will perform the decays of all
resonances it knows how to do, in the same way as for internal processes.
Note that you will be on slippery ground if you then restrict the decay of
these resonances to specific allowed channels since, if this is not what
was intended, you will obtain the wrong cross section and potentially the
wrong mix of different event types. (Since the original intention is
unknown, the cross section will not be corrected for the fraction of
open channels, i.e. the procedure used for internal processes is not
applied in this case.)
An interface to Les Houches Event Files
The LHEF standard [Alw06] specifies a format where a single file
packs initialization and event information. This has become the most
frequently used procedure to process external parton-level events in Pythia.
Therefore a special
pythia.init(fileName)
initialization option exists, where the LHEF name is provided as input.
Internally this name is then used to create an instance of the derived
class LHAupLHEF
, which can do the job of reading an LHEF.
An example how to generate events from an LHEF is found in
main12.cc
. Note the use of
pythia.info.atEndOfFile()
to find out when the whole
LHEF has been processed.
To allow the sequential use of several event files the
init(...)
method has an optional second argument:
pythia.init(fileName, bool skipInit = false)
.
If called with this argument true
then there will be no
initialization, except that the existing LHAupLHEF
class
instance will be deleted and replaced by ones pointing to the new file.
It is assumed (but never checked) that the initialization information is
identical, and that the new file simply contains further events of
exactly the same kind as the previous one. An example of this possibility,
and the option to mix with internal processes, is found in
main13.cc
.
A runtime Fortran interface
The runtime Fortran interface requires linking to an external Fortran
code. In order to avoid problems with unresolved external references
when this interface is not used, the code has been put in a separate
LHAFortran.h
file, that is not included in any of the
other library files. Instead it should be included in the
user-supplied main program, together with the implementation of two
methods below that call the Fortran program to do its part of the job.
The LHAupFortran
class derives from LHAup
.
It reads initialization and event information from the LHA standard
Fortran commonblocks, assuming these commonblocks behave like two
extern "C" struct
named heprup_
and
hepeup_
. (Note the final underscore, to match how the
gcc compiler internally names Fortran files.)
The instantiation does not require any arguments.
The user has to supply implementations of the fillHepRup()
and fillHepEup()
methods, that is to do the actual calling
of the external Fortran routines that fill the HEPRUP
and
HEPEUP
commonblocks. The translation of this information to
the C++ structure is provided by the existing setInit()
and
setEvent()
code.
See further
here for information how
PYTHIA 6.4 can be linked to make use of this facility.
Methods for LHEF output
The main objective of the LHAup
class is to feed information
from an external program into PYTHIA. It can be used to export information
as well, however. Specifically, there are four routines in the base class
that can be called to write a Les Houches Event File. These should be
called in sequence in order to build up the proper file structure.
method
openLHEF(string filename)
Opens a file with the filename indicated, and writes a header plus a brief
comment with date and time information.
method
initLHEF()
Writes initialization information to the file above. Such information should
already have been set with the methods described in the "Initialization"
section above.
method
eventLHEF()
Writes event information to the file above. Such information should
already have been set with the methods described in the "Event input"
section above. This call should be repeated once for each event to be
stored.
method
closeLHEF(bool updateInit = false)
Writes the closing tag and closes the file. Optionally, if
updateInit = true
, this routine will reopen the file from
the beginning, rewrite the same header as openLHEF()
did,
and then call initLHEF()
again to overwrite the old
information. This is especially geared towards programs, such as PYTHIA
itself, where the cross section information is not available at the
beginning of the run, but only is obtained by Monte Carlo integration
in parallel with the event generation itself. Then the
setXSec( i, xSec)
, setXErr( i, xSec)
and
setXMax( i, xSec)
can be used to update the relevant
information before closeLHEF
is called.
Warning: overwriting the beginning of a file without
upsetting anything is a delicate operation. It only works when the new
lines require exactly as much space as the old ones did. Thus, if you add
another process in between, the file will be corrupted.
PYTHIA 8 output to an LHEF
The above methods could be used by any program to write an LHEF.
For PYTHIA 8 to do this, a derived class already exists,
LHAupFromPYTHIA8
. In order for it to do its job,
it must gain access to the information produced by PYTHIA,
specifically the process
event record and the
generic information stored in info
. Therefore, if you
are working with an instance pythia
of the
Pythia
class, you have to instantiate
LHAupFromPYTHIA8 with pointers to the
process
and info
objects of
pythia
:
LHAupFromPYTHIA8 myLHA(&pythia.process, &pythia.info);
The method setInit()
should be called to store the
pythia
initialization information in the LHA object,
and setEvent()
to store event information.
Furthermore, updateSigma()
can be used at the end
of the run to update cross-section information, cf.
closeLHEF(true)
above. An example how the
generation, translation and writing methods should be ordered is
found in main20.cc
.
Currently there are some limitations, that could be overcome if
necessary. Firstly, you may mix many processes in the same run,
but the cross-section information stored in info
only
refers to the sum of them all, and therefore they are all classified
as a common process 9999. Secondly, you should generate your events
in the CM frame of the collision, since this is the assumed frame of
stored Les Houches events, and no boosts have been implemented
for the case that pythia.process
is not in this frame.