Adding a new Generator

Under:
Event generators are responsible for creating the final state particles which are fed out to GEANT for simulation.  They can be as simple as particle guns, shooting individual particles along well defined trajectories, or complex hydrodynamical models of heavy-ion collisions.  Regardless of the complexities of the underlying physical model, the job of an event generator in the STAR framework is to add particles to an event record.  In this document we will describe the steps needed to add a new event generator to the STAR framework.  The document will be divided into three sections: (1)  An overview, presenting the general steps which are required; (2) A FORtran-specific HOWTO, providing guidance specific to the problem of interfacing FORtran with a C++ application; and (3) A document describing the STAR event record.

Contents:
1.0 Integrating Event Generators
2.0 Integrating FORtran Event Generators
3.0 The STAR Event Record


1.0 Integrating Event Generators

The STAR Event Generator Framework implements several C++ classes which facilitate the integration of FORtran and C++ event generators with the STAR simulation code.  The code is available in the CVS repository and can be checked out as
$ cvs co StRoot/StarGenerator
After checking out the generator area you will note that the code is organized into several directories, containing both CORE packages and concrete event generators.  Specifically:

StarGenerator/BASE  -- contains the classes implementing the STAR interface to event generators
StarGenerator/EVENT -- contains the classes implementing the STAR event record
StarGenerator/UTIL  -- contains random number generator base class and particle data
StarGenerator/TEST  -- contains test makers used for validating the event generators


The concrete event generators (at the time this document was assembled) include

StarGenerator/Hijing1_383
StarGenerator/Pepsi
StarGenerator/Pythia6_4_23
StarGenerator/Pythia8_1_62


1.1 Compiling your Generator

Your first task in integrating a new event generator is to create a directory for it under StarGenerator, and get your code to compile.   You should select a name for your directory which includes the name and version of your event generator.  It should also be CamelCased...  MyGenerator1_2_3, for example.  (Do not select a name which is ALL CAPS, as this has a special meaning for compilation).  Once you have your directory, you can begin moving your source files into the build area.  In general, we would like to minimize the number of edits to the source code to make it compile.  But you may find that you need to reorganize the directory structure of your code to get it to compile under cons.  (For certain, if your FORtran source ends in ".f" you will need to rename the file to ".F", and if your C++ files end in ".cpp" or ".cc", you may need to rename to ".cxx".)

1.2 Creating your Interface

Ok.  So the code compiles.  Now we need to interface the event generation machinery with the STAR framework.  This entails several things.  First, we need to expose the configuration of the event generator so that the end user can generate the desired event sample.  We must then initialize the concrete event generator at the start of the run, and then exercise the event generation machinery on each and every event.  Finally, we need to loop over all of the particles which were created by the event generator and push them onto the event record so that they are persistent (i.e. the full event can be analyzed at a later date) and so that the particles are made available to the Monte Carlo application for simulation.

The base class for all event generator interfaces is  StarGenerator

Taking a quick look at the code, we see that there are several "standard" methods defined for configuring an event generator:
These methods have been defined in order to establish a common interface amongst all event generators in STAR.  These methods set variables defined within the class, from which you will initialize your concrete event generator.
You may need to implement additional methods in order to expose the configuration of your event generator.  You should, of course, do this.

The two methods which StarGenerator requires you to implement are Init() and Generate().  These methods will respectively be called at the start of each run, and during each event.

Init() is responsible for initializing the event generator.  In this method, you should pass any of the configuration information on to your concrete event generator.  This may be through calls to subroutines in your event generator, or by setting values in common blocks.  However this is done, this is the place to do it.

Generate() will be called on every single event during the run.   This is where you should exercise the generation machinery of your event generator.  Every event generator handles this differently, so you will need to consult your manual to figure out the details.

Once Generate() has been called, you are ready to fill the event record.  The event record consists of two parts: (1) the particle record, and (2) the event-wise information describing the physical interaction being simulated.  At a minimum, you will need to fill the particle-wise information.  For more details, see The STAR Event Record  below.


2.0 Integrating FORtran Event Generators

Interfacing a FORtran event generator with ROOT involves (three) steps:

1. Interface the event generator's common blocks (at least the ones which we need to use) to C++
2. Map those common blocks onto C++ structures
3. Expose the C++ structures representing the common blocks to ROOT so that users may modify / access their contents

Let's look at the pythia event generator for a concrete example. 

If you examine the code in StRoot/StarGenerator/Pythia6_4_23/ there is a FORtran file named address.F.  Open that up in your favorite editor and have a look... You'll see several functions defined.  The first one is address_of_pyjets.  In it we declare the PYJETS common block, essentially just cutting and pasting the delcaration from the pythia source code in the same directory.
      Function address_of_pyjets()
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
Integer address_of_pyjets
Integer LOC
COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
address_of_pyjets = loc(N)
return
End
We use the intrinsic function LOC to return the address (i.e. pointer) to the first variable in the common block.  We have just created a FORtran function which returns a pointer to the data stored in this common block.  The remaining functions in address.F simply expose the remaining common blocks in pythia which we want access to.  By calling this function on the C++ side, we will obtain a pointer to the memory address where the common block resides.

Next we need to describe the memory layout to C++.  This is done in the file Pythia6.h.  Each common block setup in address.F has a corresponding structure defined in this header file.  So, let's take a look at the setup for the PyJets common block:

//                                                                              
// Interface to the PYJETS common block
//

#define address_of_pyjets F77_NAME( address_of_pyjets, ADDRESS_OF_PYJETS )
struct PyJets_t {
/* Layout of the memory. */
Int_t n;
Int_t npad;
Int_t _k[5][4000];
Double_t _p[5][4000];
Double_t _v[5][4000];
/* Add access methods which mimic fortran arrays */
Int_t &k( Int_t i, Int_t j ){ return _k[j-1][i-1]; }
Double_t &p( Int_t i, Int_t j ){ return _p[j-1][i-1]; }
Double_t &v( Int_t i, Int_t j ){ return _v[j-1][i-1]; }
};
extern "C" PyJets_t *address_of_pyjets();



First, notice the first line where we call the c-preprocessor macro "F77_NAME".  This line handles, in a portable way, the different conventions between FORtran and C++ compilers, when linking together object codes. 

Next, let's discuss "memory layout".  In this section of the code we map the FORtran memory onto a C-structure.   Every variable in the common block should be declared in the same order in the C struct as it was declared in FORtran, and with the corresponding C data type.  These are:
INTEGER          --> Int_t
REAL             --> Float_t
REAL *4          --> Float_t 
REAL *8          --> Double_t 
DOUBLE PRECISION --> Double_t
You probably noticed that there are two differences with the way we have declared the arrays.  First, the arrays all were declared with an "_" in front of their name.  This was a choice on my part, which I will explain in a moment.  The important thing to notice right now is that the indicies on the arrays are reversed, compared to their declarion in FORtran.  "INTEGER K(4000,5)" in FORtran becomes "Int_t _k[5][4000]" in C++.  The reason for this is that C++ and FORtran represent arrays differently in memory.  It is important to keep these differences in mind when mapping the memory of a FORtran common block --

1) The indices in the arrays will always be reversed between FORtran and C --   A(10,20,30) in FORtran becomes A[30][20][10] in C.
2) FORtran indices (by default) start from 1, C++ (always) from 0 --  i.e. B(1) in FORtran would be B[0] in C.
3) FORtran indices may start from any value.  An array declared as D(-10:10) would be declared in C as D[21], and D(-10) in FORtran is D[0] in C.

What about the underscore?

We need to make some design choices at this point.  Specifically, how do we expose the common blocks to the end user?  Do we want the end user to deal with the differences in C++ and FORtran, or do we want to provide a mechanism by which the FORtran behavior (i.e. count from 1, preserve the order of indices) can be emulated.

My preference is to do the latter -- provide the end user functions which emulate the behavior of the FORtran arrays, because these arrays are what is documented in the event generator's manual.   This will minimize the likelyhood that the end user will make mistakes in configuring t he event generator.


So we have created a c-struct which describes how the common block's memory is laid out, and we have defined a function in the FORtran library which returns the location of the memory address of the common block.  Now we need to expose that function to C++.  To do that, we need to declare a prototype of the function.  There are two things we need to do.

First, we need to define the name of the subroutine in a portable way.  This is done using a macro defined in #include "StarCallf77.h" --

#define address_of_pyjets F77_NAME( address_of_pyjets, ADDRESS_OF_PYJETS )

Next we need to declare to C++ that address_of_pyjets can be found in an external library, and will return a pointer to the PyJets_t structure

extern "C" PyJets_t *address_of_pyjets();

Now we are almost done.  We need to add a function in our generator class which returns a pointer (or reference) to the common blocks, and we need to add PyJets_t to the ROOT dictionary... In MyGeneratorLinkDef.h, add the line

#pragma link C++ struct PyJets_t+;

Finally, you need to expose the FORtran subroutines to C++.  Again, take a look at the code in Pythia6.h and Pythia6.cxx.  In the header we declare wrapper functions around the FORtran subroutines, and in the implementation file we expose the FORtran subroutines to C++.

Our first step is declaring the prototypes of the subroutines and implementing the C++ infterface.  Consider the SUBROUTINE PYINIT in pythia, which  initializes the event generator.  In FORtran it is declared as

      SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
...
      CHARACTER*(*) FRAME,BEAM,TARGET


So the variables FRAME, BEAM and TARGET are declared as character variables, and WIN is implicitly a double precision variable.

There are several webpages which show how to interface fortran and c++, e.g. http://www.yolinux.com/TUTORIALS/LinuxTutorialMixingFortranAndC.html

It is really a system-dependent thing.  We're keeping it simple and only supporting Linux. 

#define pyinit F77_NAME(pyinit,PYINIT) /* pythia initialization */
extern "C" void   type_of_call  pyinit( const char *frame,
                                        const char *beam,
                                        const char *targ,
                                        double *ener,
                                        int nframe,
                                        int nbeam,
                                        int ntarg );

So there's three character varaibles declared: frame, beam and targ.  These correspond to the character variables on the FORtran side.  Also there's a double precision variable ener... this is WIN.   But then there are three integer variables nframe, nbeam  and ntarg.  FORtran expects to get the size of the character variables when the subroutine is called, and it gets them after the last arguements in the list.

Again, we would like to hide this from the end user... so I like to define wrapper functions such as

#include <string>

void PyInit( string frame, string blue, string yellow, double energy )
{
   pyinit( frame.c_str(), blue.c_str(), yellow.c_str(), &energy,
           frame.size(),
           blue.size(),
           yellow.size() );
}


3.0 The STAR Event Record