pythia8130 distributed with AliRoot

[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / xmldoc / ProgramFlow.xml

<chapter name="Program Flow">

<h2>Program Flow</h2>

Recall that, to first order, the event generation process can be 
subdivided into three stages:
<ol>
<li>Initialization.</li>
<li>The event loop.</li>
<li>Finishing.</li>
</ol>
This is reflected in how the top-level <code>Pythia</code> class should 
be used in the user-supplied main program, further outlined in the 
following. Since the nature of the run is defined at the initialization
stage, this is where most of the PYTHIA user code has to be written. 
So as not to confuse the reader unduly, the description of initialization 
options has been subdivided into what would normally be used and what is 
intended for more special applications.

<h3>Initialization - normal usage</h3>

<ol>

<li>
Already at the top of the main program file, you need to include the proper 
header file
<pre>
    #include "Pythia.h"
</pre>
To simplify typing, it also makes sense to declare
<pre>
    using namespace Pythia8; 
</pre>
</li>

<p/>
<li>
The first step is to create a generator object, 
e.g. with
<pre>
     Pythia pythia;
</pre>
It is this object that we will use from now on. Normally a run
will only contain one <code>Pythia</code> object. (Hypothetically 
you  could use several <code>Pythia</code> objects sequentially, 
but if done in parallel the <code>static</code> character of some
program elements is likely not to give the desired behaviour.)<br/>
By default all output from <code>Pythia</code> will be on the 
<code>cout</code> stream, but the <code>list</code> methods below do
allow output to alternative streams or files (by an optional
last argument, a reference to an <code>ostream</code>, usually not
explicitly written out here).
</li>  

<p/>
<li> 
You next want to set up the character of the run. 
The pages under the "Setup Run Tasks" heading in the index
describe all the options available (with some very few exceptions,
found on the other pages).  
The default values and your modifications are stored in two databases, 
one for <aloc href="SettingsScheme">generic settings</aloc>
and one for <aloc href="ParticleDataScheme">particle data</aloc>. 
Both of these are static classes, and are 
initialized with their default values by the <code>Pythia</code> 
constructor. The default values can then be changed, primarily
by one of the two ways below, or by a combination of them.<br/>  

<p/>
a) You can use the dedicated methods of each class to change the 
database default values to fit the needs of your current run. However, 
the
<pre>
    pythia.readString(string);
</pre>
method provides a covenient uniform interface to all of them.
The information in the string is case-insensitive, but upper- and
lowercase can be combined for clarity. The rules are that<br/>
(i) if the first nonblank character of the string is a letter
it is assumed to contain a setting, and is sent on to 
<code>pythia.settings.readString(string)</code>;<br/> 
(ii) if instead the string begins with a digit it is assumed to 
contain particle data updates, and so sent on to 
<code>pythia.particleData.readString(string)</code>;<br/>
(iii) if none of the above, the string is assumed to be a comment,
i.e. nothing will be done.<br/>
In the former two cases, a warning is issued whenever a string
cannot be recognized (maybe because of a spelling mistake),
unless an optional second argument <code>false</code> is used to 
switch off warnings.<br/>
Some examples would be
<pre>
    pythia.readString("TimeShower:pTmin = 1.0");
    pythia.readString("111:mayDecay = false");
</pre>
The <code>readString(string)</code> method is intended primarily for 
a few changes. It can also be useful if you want to construct a
parser of input files that contain commands to several different 
libraries.<br/>

<p/>
b) You can read in a file containing a list of those variables 
you want to see changed, with a 
<pre>
    pythia.readFile(fileName);
</pre>
Each line in this file with be processes by the 
<code>readString(string)</code> method introduced above. You can thus 
freely mix comment lines and lines handed on to <code>Settings</code> 
or to <code>ParticleDataTable</code>.   
Again, an optional second argument <code>false</code> allows 
you to switch off warning messages for unknown variables.<br/>
This approach is better suited for more extensive changes than a direct
usage of <code>readString(string)</code>, and can also avoid having to
recompile and relink your main program between runs.
</li>

<p/>
<li>
Next comes the initialization stage, where all 
remaining details of the generation are to be specified. The 
<code>init(...)</code> method allows a few different input formats, 
so you can pick the one convenient for you:

<p/>
a) <code>pythia.init( idA, idB, eCM);</code><br/>
lets you specify the identities and the CM energy of the two incoming
beam particles, with A (B) assumed moving in the <ei>+z (-z)</ei> 
direction.

<p/>
b) <code>pythia.init( idA, idB, eA, eB);</code><br/>
is similar, but the two beam energies can be different, so the
collisions do not occur in the CM frame. If one of the beam energies 
is below the particle mass you obtain a fixed-target topology.

<p/>
c) <code>pythia.init( idA, idB, pxA, pyA, pzA, pxB, pyB, pzB);</code><br/>
is similar, but here you provide the three-momenta 
<ei>(p_x, p_y, p_z)</ei> of the two incoming particles,
to allow for arbitrary beam directions.

<p/>
d) <code>pythia.init(fileName);</code> <br/> 
assumes a file in the <aloc href="LesHouchesAccord">Les Houches 
Event File</aloc> format <ref>Alw06</ref> is provided.

<p/>
e) <code>pythia.init();</code><br/>
with no arguments will read the beam parameters from the 
<aloc href="MainProgramSettings"><code>Main</code></aloc> 
group of variables, which provides you with the same possibilities as
the above options a, b, c and d. If you don't change any of those you will 
default to proton-proton collisions at 14 TeV, i.e. the nominal LHC 
values.

<p/>
f) <code>pythia.init( LHAup*);</code> <br/>
assumes <aloc href="LesHouchesAccord">Les Houches Accord</aloc> 
<ref>Boo01</ref> initialization and event information is available 
in an <code>LHAup</code> class object, and that a pointer to this object 
is handed in.

<p/>
<li>
If you want to have a list of the generator and particle data used, 
either only what has been changed or everything, you can use 
<pre>
    pythia.settings.listChanged();
    pythia.settings.listAll();
    pythia.particleData.listChanged(); 
    pythia.particleData.listAll(); 
</pre>
</li>

</ol>

<h3>The event loop</h3>

<ol>

<li>
Inside the event generation loop you generate the 
next event using the <code>next()</code> method,
<pre>
    pythia.next();
</pre>
This method takes no arguments; everything has already been specified. 
It does return a bool value, however, <code>false</code> when the
generation failed. This can be a "programmed death" when the
supply of input parton-level configurations on file is exhausted,
but also caused by a failure of <code>Pythia</code> to generate an event,
or that an event was generated but something strange was detected
in it. It makes sense to allow a few <code>false</code> 
values before a run is aborted, so long as the related faulty
events are skipped.
</li>  
 
<p/>
<li>
The generated event is now stored in the <code>event</code> 
object, of type <aloc href="EventRecord"><code>Event</code></aloc>, 
which is a public member of <code>pythia</code>. You therefore have 
access to all the tools described on the pages under the "Study Output" 
header in the index. For instance, an event can be listed with 
<code>pythia.event.list()</code>, the identity of the <ei>i</ei>'th 
<aloc href="ParticleProperties">particle</aloc> is given by 
<code>pythia.event[i].id()</code>, and so on.<br/> 
The hard process - roughly the information normally stored in the 
Les Houches Accord event record - is available as a second object, 
<code>process</code>, also of type <code>Event</code>.<br/> 
A third public object is 
<aloc href="EventInformation"><code>info</code></aloc>, which offers 
a set of one-of-a kind pieces of information about the most recent
event.
</li> 

</ol>

<h3>Finishing</h3>

<ol>

<li>At the end of the generation process, you can call
<pre>
    pythia.statistics(); 
</pre>
to get some run statistics, on cross sections and the number of errors 
and warnings encountered. With optional argument <code>true</code>
also further statistics is printed. Currently this means the number of
different subprocesses generated in the multiple-interactions
framework. 
</li> 

</ol>

<h3>Advanced usage, mainly for initialization</h3>

A) Necessary data are automatically loaded when you use the 
default PYTHIA installation directory structure and run the main 
programs in the <code>examples</code> subdirectory. However, in the 
general case, you must provide the path to the <code>.xml</code> files, 
originally stored in the <code>xmldoc</code> directory, where default 
settings and particle data are found. This can be done in two ways.

<ol>

<li>
You can set the environment variable <code>PYTHIA8DATA</code> to
contain the location of the <code>xmldoc</code> directory. In the
<code>csh</code> and <code>tcsh</code> shells this could e.g. be 
<pre>
     setenv PYTHIA8DATA /home/myname/pythia81xx/xmldoc
</pre>
while in other shells it could be
<pre>
     export PYTHIA8DATA=/home/myname/pythia81xx/xmldoc
</pre>
where xx is the subversion number.<br/>
Recall that environment variables set locally are only defined in the 
current instance of the shell. The above lines should go into your 
<code>.cshrc</code> and <code>.bashrc</code> files, respectively, 
if you want a more permanant assignment.
</li>

<p/>
<li>
You can provide the path as argument to the <code>Pythia</code>
constructor, e.g.
<pre>
     Pythia pythia("/home/myname/pythia81xx/xmldoc");
</pre>
</li>
</ol>
where again xx is the subversion number.<br/>
When <code>PYTHIA8DATA</code> is set it takes precedence, else 
the path in the constructor is used, else one defaults to the 
<code>../xmldoc</code> directory.

<p/>
B) You can override the default behaviour of PYTHIA not only by the 
settings and particle data, but also by replacing some of the 
PYTHIA standard routines by ones of your own. Of course, this is only
possible if your routines fit into the general PYTHIA framework.
Therefore they must be coded according to the the rules relevant
in each case, as a derived class of a PYTHIA base class, and a pointer 
to such an object must be handed in by one of the methods below.
These calls must be made before the <code>pythia.init(...)</code> call.

<ol>

<li>
If you are not satisfied with the list of parton density functions that 
are implemented internally or available via the LHAPDF interface
(see the <aloc href="PDFSelection">PDF Selection</aloc> page), you 
can suppy your own by a call to the <code>setPDFPtr(...)</code> method
<pre>
      pythia.setPDFptr( pdfAPtr, pdfBPtr); 
</pre>
where <code>pdfAPtr</code> and <code>pdfBPtr</code> are pointers to 
two <code>Pythia</code> <aloc href="PartonDistributions">PDF 
objects</aloc>. Note that <code>pdfAPtr</code> and <code>pdfBPtr</code> 
cannot point to the same object; even if the PDF set is the same, 
two copies are needed to keep track of two separate sets of <ei>x</ei>
and density values.<br/>
If you further wish to use separate PDF's for the hard process of an
event than the ones being used for everything else, the extended form
<pre>
      pythia.setPDFptr( pdfAPtr, pdfBPtr, pdfHardAPtr, pdfHardBPtr); 
</pre>
allows you to specify those separately, and then the first two sets 
would only be used for the showers and for multiple interactions.
</li>

<p/>
<li>
If you want to perform some particle decays with an
external generator, you can call the <code>setDecayPtr(...)</code> 
method
<pre>
      pythia.setDecayPtr( decayHandlePtr, particles);
</pre>
where the <code>decayHandlePtr</code> derives from the 
<aloc href="ExternalDecays"><code>DecayHandler</code></aloc> base 
class and <code>particles</code> is a vector of particle codes to be 
handled. 
</li>

<p/>
<li>
If you want to use an external random number generator, 
you can call the <code>setRndmEnginePtr(...)</code> method
<pre>
      pythia.setRndmEnginePtr( rndmEnginePtr); 
</pre>
where <code>rndmEnginePtr</code> derives from the 
<aloc href="RandomNumbers"><code>RndmEngine</code></aloc> base class. 
The <code>Pythia</code> default random number generator is perfectly 
good, so this is only intended for consistency in bigger frameworks.
</li>

<p/>
<li>
If you want to interrupt the evolution at various stages, 
to interrogate the event and possibly veto it, or you want to
reweight the cross section, you can use   
<pre>
      pythia.setUserHooksPtr( userHooksPtr); 
</pre>
where <code>userHooksPtr</code> derives from the 
<aloc href="UserHooks"><code>UserHooks</code></aloc> base class.
</li>

<p/>
<li>
If you want to use your own parametrization of beam momentum spread and
interaction vertex, rather than the provided simple Gaussian 
parametrization (off by default), you can call
<pre>
      pythia.setBeamShapePtr( beamShapePtr); 
</pre>
where <code>beamShapePtr</code> derives from the 
<aloc href="BeamShape"><code>BeamShape</code></aloc> base class.
</li>

<p/>
<li>
If you want to implement a cross section of your own, but still make use
of the built-in phase space selection machinery, you can use
<pre>
      pythia.setSigmaPtr( sigmaPtr);
</pre>
where <code>sigmaPtr</code> of type <code>SigmaProcess*</code> is an
instance of a class derived from one of the <code>Sigma1Process</code>,
<code>Sigma2Process</code> and  <code>Sigma3Process</code> base classes
in their turn derived from 
<aloc href="SemiInternalProcesses"><code>SigmaProcess</code></aloc>. 
This call can be used repeatedly to hand in several different processes.
</li>

<p/>
<li>
If your cross section contains the production of a new resonance
with known analytical expression for all the relevant partial widths,
you can make this resonance available to the program with 
<pre>
      pythia.setResonancePtr( resonancePtr);
</pre>
where <code>resonancePtr</code> of type <code>ResonanceWidths*</code> 
is an instance of a class derived from the 
<aloc href="SemiInternalResonances"><code>ResonanceWidths</code></aloc> 
base class. In addition you need to add the particle to the normal 
<aloc href="ParticleDataScheme">particle and decay database</aloc>.
This procedure can be used repeatedly to hand in several different 
resonances.
</li>

<p/>
<li>
If you are a real expert and want to replace the PYTHIA initial- and 
final-state showers, you can use
<pre>
      pythia.setShowerPtr( timesDecPtr, timesPtr, spacePtr);
</pre>
where <code>timesDecPtr</code> and <code>timesPtr</code>
derive from the <code>TimeShower</code> base class, and 
<code>spacePtr</code> from <code>SpaceShower</code> 
(<aloc href="ImplementNewShowers">further instructions</aloc>). 
</li>

</ol>

<p/>
C) Some comments on collecting several tasks in the same run.
<ol>

<li>
PYTHIA has not been written for threadsafe execution. As a rule, you 
should not have more than one <code>Pythia</code> object in an executable
at any time. For multicore processors, if you want to use all cores, 
the most efficient way presumably is to start correspondingly many jobs, 
with different random number seeds, and add the statistics at the end.

<p/>
In some cases it is possible to use more than one <code>Pythia</code> 
object. The key example would be the simultaneous generation of signal 
and pileup events, see <code>main19.cc</code>. Here all signal processes 
must be switched on before the first initialization, and then switched 
off and replaced by the background ones before the second initialization. 
Also most other settings can be changed consistently in between the two
initializations, but in a few cases the last value will be used. Particle 
data is always based on the latest information. As a rule, however, it is 
safer to use two separate runs to store events on disk, in two separate 
files, and mix afterwards.
</li>

<p/>
<li>
When time is not an issue, it may be that you want to perform several 
separate subruns sequentially inside a run, e.g. to combine results for
several kinematical regions or to compare results for some different 
tunes of the underlying event. One way to go is to create and destroy a 
<code>pythia</code> object once for each subrun, in which case they are 
completely separate. You can also use the same <code>pythia</code> object, 
only doing a new <code>init(...)</code> call for each subrun. In that 
case, the settings and particle databases remain as they were in the  
previous subrun, only affected by the specific changes you introduced in 
the meantime. You can put those changes in the main program, with
<code>pythia.readString(string)</code>, using your own logic to decide
which ones to execute in which subrun. A corresponding possibility 
exists with <code>pythia.readFile(fileName, subrun)</code>, which 
as second argument can take a non-negative subrun number. (Or, 
alternatively use the longer form 
<code>pythia.readFile(fileName, warn, subrun)</code>.) Then only those 
sections of the file before any <code>Main:subrun = ...</code> line
or with matching <code>subrun</code> number will be read. That is, the
file could have a structure like
<pre>
    ( lines always read, i.e. "default values" always (re)set )
    Main:subrun = 1
    ( lines only read with readFile(fileName, 1) )
    Main:subrun = 2
    ( lines only read with readFile(fileName, 2) )
</pre>
Both of these possibilities are illustrated in <code>main08.cc</code>.
</li>

<p/>
<li>
When working with Les Houches Event Files, it may well be that your 
intended input event sample is spread over several files, that you all 
want to turn into complete events in one and the same run. There is no
problem with looping over several subruns, where each new subrun 
is initialized with a new file. However, in that case you will do
a complete re-initialization each time around. If you want to avoid
this, note that there is an optional second argument for LHEF
initialization: <code>pythia.init(fileName, skipInit)</code>.
Alternatively, the tag <code>Main:LHEFskipInit</code> can be put 
in a file of commands to obtain the same effect.
Here <code>skipInit</code> defaults to <code>false</code>,
but if set <code>true</code> then the new file will be simulated 
with the same initialization data as already set in a previous
<code>pythia.init(...)</code> call. The burden rests on you to ensure 
that this is indeed correct, e.g. that the two event samples have not 
been generated for different beam energies.
</li>

</ol>
   
</chapter>

<!-- Copyright (C) 2008 Torbjorn Sjostrand -->