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, from the respective Fortran
commonblocks, or from PYTHIA 8 itself.
You are free to write your own derived classes, using the rules and
methods to be described below. Normally, pointers to objects of such
derived classes should be handed in with the
Pythia::setLHAupPtr( LHAup*)
method. However, with the LHEF format a filename can replace the
pointer, see further below.
Let us now describe the methods at your disposal to do the job.
LHAup::LHAup( int strategy = 3) setStrategy,
see below. It also reserves some space for processes and particles.
virtual LHAup::~LHAup() Pythia::init().
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.
virtual bool LHAup::setInit() setInit(), such information can be set by the following
methods:
void LHAup::setBeamA( int identity, double energy, int pdfGroup = 0, int pdfSet = 0) IDBMUP(1), EBMUP(i), PDFGUP(i), PDFSUP(i),
with i 1 or 2). 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.
void LHAup::setStrategy( int strategy) 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
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
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.
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.
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.
void LHAup::setXSec( int i, double xSec) xSec value of the i'th process
added with addProcess method (i.e. i runs
from 0 through sizeProc() - 1, see below).
void LHAup::setXErr( int i, double xErr) xErr value of the i'th process
added with addProcess method.
void LHAup::setXMax( int i, double xMax) xMax value of the i'th process
added with addProcess method.
void LHAup::setInfoHeader(string &key, string &val) key to have value val.
This is a wrapper function to the
Info::setHeader function that
should be used in any classes derived from LHAup.
Information is handed back by the following methods
(that normally you would not need to touch):
int LHAup::idBeamA() i in the range 0 <= i <
sizeProc().
double LHAup::xSecSum() Pythia::init(),
so would normally not need to be called directly by the user.
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).
virtual bool LHAup::setEvent(int idProcess = 0) idProcess value specifies which process
that should be generated, while idProcess is irrelevant
for strategies +-3 and +-4. The method 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
Info::atEndOfFile()
to confirm that indeed the failure is caused in this method, and decide
to break out of the event generation loop.
Inside a normal setEvent(...) call, information can be set
by the following methods:
void LHAup::setProcess( int idProcess, double weight, double scale, double alphaQED, double alphaQCD) 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.
void LHAup::addParticle( int id, int status, int mother1, int mother2, int colourTag1, int colourTag2, double p_x, double p_y, double p_z, double e, double m, double tau, double spin, double scale) IDUP, ISTUP,
MOTHUP(1,..), MOTHUP(2,..), ICOLUP(1,..), ICOLUP(2,..), PUP(J,..),
VTIMUP, SPINUP; while scale is a new optional property,
see further below) .
Information is handed back by the following methods:
int LHAup::idProcess() 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.
double LHAup::scale() i in the range
0 <= i < sizePart(). (But again note that
i = 0 is an empty line, so the true range begins at 1.)
From the information in the event record it is possible to set
the flavour and x values of the initiators
void LHAup::setIdX(int id1, int id2, double x1, double x2) true.
This information is returned by the methods
bool LHAup::pdfIsSet() setPdf(...) must be called after the
setProcess(...) call of the event for this to work.)
Note that the flavour and x values usually but not always
agree with those obtained by the same methods without pdf
in their names, see explanation in the
Event Information description.
The maximum scale for parton-shower evolution of a Les Houches event is
regulated by the
TimeShower:pTmaxMatch
and
SpaceShower:pTmaxMatch
modes. If you want to guarantee that the input scale value
is respected, as is often the case in matching/merging procedures, you
should set both of these modes to 1. That only affects the hard process,
while resonance decays are still processed using the resonance mass to
set the upper limit. However, the optional
Beams:strictLHEFscale = on
setting restricts also resonance-decay emissions to be below the input
scale value.
As a further non-standard feature, it is also possible to read in the
separate scale values of all final particles. Such scale values could be used
e.g. to restrict the maximum scale for shower evolutions for each parton
separately. This reading will only be applied if the
Beams:setProductionScaleFromLHEF switch is true (see
Beam Parameters for details).
This information is returned by the method
double LHAup::scale(int i). When no such information
has been read from the LHEF, the scale defaults to -1.
void LHAup::listEvent() LHAup object is not available to the user, the
Pythia::LHAeventList() method can
be used, which is a wrapper for the above.
virtual bool LHAup::skipEvent(int nSkip) nSkip events in the Les Houches generation
sequence, without doing anything further with them. Mainly
intended for debug purposes, e.g. when an event at a known
location in a Les Houches Event File is causing problems.
Will return false if operation fails, specifically if the
end of an LHEF has been reached. The implementation in the base class
simply executes setEvent() the requested number of times.
The derived LHAupLHEF class (see below) only uses the
setNewEventLHEF(...) part of its setEvent()
method, and other derived classes could choose other shortcuts.
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.)
Even if PYTHIA can select resonance decay modes according to its
internal tables, there is normally no way for it to know which
decay angular correlations should exist in the simulated process.
Therefore almost all decays are isotropic. The exceptions are Higgs and
top decays, in the decay chains H → WW/ZZ → f fbar f' fbar'
and t → b W → b f fbar, where the process-independent
correlations implemented for internal processes are used. If part of
the decay chain has already been set, however (e.g. H → WW/ZZ
or t → b W), then decay is still isotropic.
process
event record, but leaves the LHA event record itself unchanged.
mode LesHouches:idRenameBeams
(default = 1000022; minimum = 0)mode LesHouches:setLifetime
(default = 1; minimum = 0; maximum = 2)VTIMUP
when the Les Houches event record is stored into the PYTHIA
process one. The reason is that some matrix-element
generators (like POWHEG) do not set decay times, so that it is up to
PYTHIA to make that selection. This is particularly important for
the tau lepton.
option 0 : all decay times are taken from the Les Houches input.
option 1 : the decay time of tau leptons is generated
like for internal PYTHIA taus, whereas all other decay times
are taken from the Les Houches input.
option 2 : all decay times are generated by PYTHIA, thus
completely disregarding the Les Houches values. This option could
go wrong in BSM scenarios with long-lived particles, if PYTHIA
has not been provided with the information to select those lifetimes
correctly.
mode LesHouches:setLeptonMass
(default = 1; minimum = 0; maximum = 2)option 0 : all lepton masses are taken from the Les Houches input.
option 1 : if the input lepton mass deviates by more than 10%
from the PYTHIA (data table) mass then its mass is reset according to the
PYTHIA value. This should catch weird masses, while allowing sensible
variations.
option 2 : each lepton mass is reset according to the PYTHIA value.
LesHouches:matchInOut = on, but not
always perfectly. One possibility then is to change the
tolerance to such errors.
mode LesHouches:setQuarkMass
(default = 1; minimum = 0; maximum = 2)option 0 : all quark masses are taken from the Les Houches input.
option 1 : if the input c or b mass is
more than 50% away from the PYTHIA (data table) mass then its mass is
reset according to the PYTHIA value.
option 2 : if the input mass, for all quarks except the top, is
more than 50% away from the PYTHIA (data table) mass then its mass is
reset according to the PYTHIA value.
LesHouches:matchInOut = true, but not always perfectly.
One possibility then is to change the
tolerance to such errors.
parm LesHouches:mRecalculate
(default = -1.)flag LesHouches:matchInOut
(default = on)Beams:frameType = 4 and Beams:LHEF to be the file
name, see Beam Parameters. Internally
this name is then used to create an instance of the derived class
LHAupLHEF, which can do the job of reading an LHEF.
As some information in a Les Houches Event File init block is only known
at the end of generation, some programs choose to output this as a
separate file. If so, the name of this file can be specified by
Beams:LHEFheader.
The two key compulsory parts of an LHEF is the initialization information
stored in an init block, enclosed by a matching <init>
- </init> pair of lines, and the event input, with each
event enclosed by a matching <event> -
</event> pair of lines. In the case of the no-beams
extension the init block may be empty, but the <init>
and </init> lines must be included for the file parsing
to work as expected. It is also possible to have a non-empty init block,
with the beams assigned code 0, and optionally a number of specified
"processes".
The latest update of the LHEF format [But14] introduced a
multitude of different optional features. This means that apart
from the <init> and <event>
tags, a plethora of new, optional information is available.
Furthermore, the inclusion of an arbitrary number of attributes into
the tags should be supported. The LHEF reader in Pythia adheres to
the updated LHEF format without any restriction. The new generation
information available through the updated LHEF format can be
retrieved by using Pythia's Info class. For a detailed
description, please consult the section "Les Houches Event File 3.0
information" in Event Information.
The LHEF reader can also read in and store header blocks. By default
this option is switched on, but may be controlled through the
Beams:readLHEFheaders
flag if necessary. The information can later be read out through the
Info class for further processing.
Due to the non-standard nature of the information in these blocks they
are stored whole, and PYTHIA itself makes no further attempt to process
their meaning.
Because Les Houches Event files tend not to adhere strictly to XML
conventions, to consistently read in header information, certain
choices must be made. The primary goal is to make as much information
available as possible. First, information sitting directly in the
<header> block is stored under the key "base". Second, the tags
starting and ending each sub block must be on their own line. Finally,
the contents of comment blocks, <!-- -->, are still stored. The
header keys are formed hierarchically from the names of the header
blocks. This behaviour is illustrated in the following example:
<header>
BaseA
<hblock1>
1A
<hblock11>
11A <hblock111>
</hblock111> 11B
</hblock11>
1B
</hblock1>
<hblock2>
2A
<!-- 2B -->
</hblock2>
BaseB
</header>
which would lead to the following information being stored in the
Info class:
| Key | Value |
|---|---|
| base | BaseA BaseB |
| hblock1 | 1A 1B |
| hblock1.hblock11 | 11A <hblock111> </hblock111> 11B |
| hblock2 | 2A <!-- 2B --> |
README file in the main directory for details
on how to do this.
An example how to generate events from an LHEF is found in
main11.cc. Note the use of
Info::atEndOfFile() to find out when the whole
LHEF has been processed.
To allow the sequential use of several event files, the
Beams:newLHEFsameInit
can be set true. Then there will be no
initialization, except that the existing LHAupLHEF class
instance will be deleted and replaced by one 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
main12.cc. A variant, based on input in a command file,
is given in main13.cc.
In C++, real numbers are printed with an 'E' to denote the exponent part,
e.g. 1.23E+04, and are read in accordingly. Other languages may use other
letters, e.g. Fortran allows either 'E' or 'D'. A file using
the latter convention would not be readable by the standard routines.
In case you have such an "incorrectly formatted" file, a conversion to
a new corrected file could be done e.g. using sed, as a
one-line command
sed -e 's/\([0-9]\.\{0,1\}\)[dD]\([+-]\{0,1\}[0-9]\)/\1E\2/g' old.lhe > new.lhe
This replaces a 'd' or 'D' with an 'E' only when it occurs in the combination
<header>...</header> block. For conversion only
inside the <init>...</init> and
<event>...</event> blocks, create a file
convert.sed containing
/<init>/,/<\/init>/bconv
/<event>/,/<\/event>/bconv
b
:conv
s/\([0-9]\.\{0,1\}\)[dD]\([+-]\{0,1\}[0-9]\)/\1E\2/g
and run it with
sed -f convert.sed old.lhe > new.lheThe workhorses of the
LHAupLHEF class are three methods
found in the base class, so as to allow them to be reused in other
contexts.
bool LHAup::setInitLHEF(ifstream& is, bool readHeaders = false) setOldEventLHEF.
bool LHAup::setOldEventLHEF() setNewEventLHEF call can be
followed by several setOldEventLHEF ones, so as to
process the same configuration several times. This method currently
only returns true, i.e. any errors should be caught by the preceding
setNewEventLHEF call.
These three main methods build on a number of container classes and a
generic LHEF reader class (called Reader) found in
LHEF3.h and LHEF3.cc. The Reader
handles all the parsing and storage necessary to adhere with
[But14]. (A matching Writer class is also
available; see documentation in LHEF3.h how it can be
used.) All parsing that is not strictly part of the LHEF format
(e.g. the reading of header information) is instead performed directly in
the LHAupLHEF methods.
Some other small utility routines are:
bool LHAup::fileFound() LHAupLHEF
it returns false if the LHEF provided in the constructor is not
found and opened correctly.
bool LHAup::useExternal() LHAupLHEF
it returns false if the LHAupLHEF instance is constructed to
work on an input LHE file, while it returns true if the LHAupLHEF
instance is constructed to use externally provided input streams instead.
For the latter, the LHAupLHEF instance should have been
constructed with the class constructor
LHAupLHEF(Info* infoPtrIn, istream* isIn, istream* isHeadIn,
bool readHeadersIn, bool setScalesFromLHEFIn).
void LHAup::setInfoHeader(const string &key, const string &val) Info class.
A few other methods, most of them derived from the base class,
streamlines file opening and closing, e.g. if several LHE files are
to be read consecutively, without the need for a complete
reinitialization. This presupposes that the events are of the same
kind, only split e.g. to limit file sizes.
bool LHAup::newEventFile(const char* fileIn) include/Pythia8Plugins/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, and used to create a derived class that
contains 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.
Up to and including version 8.125 the LHAupFortran class
was used to construct a runtime interface to PYTHIA 6.4. This was
convenient in the early days of PYTHIA 8 evolution, when this program
did not yet contain hard-process generation, and the LHEF standard
did not yet exist. Nowadays it is more of a bother, since a full
cross-platform support leads to many possible combinations. Therefore
this support has been removed, but can still be recuperated from
previous code versions, in a reduced form up to version 8.176.
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.
bool LHAup::openLHEF(string filename) verbose = false only
leaves a single blank between the information fields.
bool LHAup::closeLHEF(bool updateInit = false) 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.
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.
The LHEF standard is the agreed format to store the particles of a
hard process, as input to generators, whereas output of final states
is normally handled using the HepMC
standard. It is possible to use LHEF also here, however. It requires
that the above initialization is replaced by
LHAupFromPYTHIA8 myLHA(&pythia.event, &pythia.info);
process is replaced by event.
In addition, the PartonLevel:all = off command found in
main20.cc obviously must be removed if one wants to
obtain complete events.
Pythia::Writer of LHEF3.h and
LHEF3.cc) for this purpose. The functions of this
file writer are used in the LHEF3FromPYTHIA8.
This latter class allows users to output PYTHIA events
in LHEF3 format from a PYTHIA main program. An example of how to use
LHEF3FromPYTHIA8 is found in the
main20lhef3.cc example. Please note that, although
similar, the usage of LHEF3FromPYTHIA8 differs from
the usage of LHAupFromPYTHIA8, with LHEF3FromPYTHIA8
requiring fewer function calls.
To print a comprehensive LHE file, LHEF3FromPYTHIA8
is constructed with pointers to an Event object,
as well as pointers to instances of Settings,
Info and ParticleData, giving e.g.
a constructor call
LHEF3FromPYTHIA8 myLHEF3(&pythia.event, &pythia.settings,
&pythia.info, &pythia.particleData);
As a next step, you should open the output file by using the
LHAupFromPYTHIA8 member function
openLHEF(string name)
name is the output file name.
Then, the method setInit() should be called to store the
initialization information (read from settings and
info) and write the header and init blocks into the
output file. Note that at this stage, the cross section printed
in the init block is not sensible, as no integration has yet
taken place. The init block can be updated at the end of
the event generation (see below).
During event generation, you should use setEvent() to
write the event information (as read from info and
event) to the output file.
Finally, before leaving your main program, it is necessary to
close the output file by using the
LHAupFromPYTHIA8 member function
closeLHEF( bool doUpdate = false)
doUpdate is optional.
If doUpdate is used, and if
doUpdate = true, then the init block of the output
file will be updated with the latest cross section information.
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.