include/HepMC/HEPEVT_Wrapper.h

0001 //--------------------------------------------------------------------------
0002
0003 #ifndef HEPEVT_EntriesAllocation
0004 #define HEPEVT_EntriesAllocation 200000
0005 #endif  // HEPEVT_EntriesAllocation
0006
0007 //--------------------------------------------------------------------------
0008 #ifndef HEPMC_HEPEVT_COMMON_H
0009 #define HEPMC_HEPEVT_COMMON_H
0010 //////////////////////////////////////////////////////////////////////////
0011 //
0012 //      PARAMETER (NMXHEP=2000)
0013 //      COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0014 //     &        JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0015 /**********************************************************/
0016 /*           D E S C R I P T I O N :                      */
0017 /*--------------------------------------------------------*/
0018 /* NEVHEP          - event number (or some special meaning*/
0019 /*                    (see documentation for details)     */
0020 /* NHEP            - actual number of entries in current  */
0021 /*                    event.                              */
0022 /* ISTHEP[IHEP]    - status code for IHEP'th entry - see  */
0023 /*                    documentation for details           */
0024 /* IDHEP [IHEP]    - IHEP'th particle identifier according*/
0025 /*                    to PDG.                             */
0026 /* JMOHEP[IHEP][0] - pointer to position of 1st mother    */
0027 /* JMOHEP[IHEP][1] - pointer to position of 2nd mother    */
0028 /* JDAHEP[IHEP][0] - pointer to position of 1st daughter  */
0029 /* JDAHEP[IHEP][1] - pointer to position of 2nd daughter  */
0030 /* PHEP  [IHEP][0] - X momentum                           */
0031 /* PHEP  [IHEP][1] - Y momentum                           */
0032 /* PHEP  [IHEP][2] - Z momentum                           */
0033 /* PHEP  [IHEP][3] - Energy                               */
0034 /* PHEP  [IHEP][4] - Mass                                 */
0035 /* VHEP  [IHEP][0] - X vertex                             */
0036 /* VHEP  [IHEP][1] - Y vertex                             */
0037 /* VHEP  [IHEP][2] - Z vertex                             */
0038 /* VHEP  [IHEP][3] - production time                      */
0039 /*========================================================*/
0040 // Remember, array(1) is the first entry in a fortran array, array[0] is the
0041 //           first entry in a C array.
0042 //
0043 // This interface to HEPEVT common block treats the block as
0044 // an array of bytes --- the precision and number of entries
0045 // is determined "on the fly" by the wrapper and used to decode
0046 // each entry.
0047 //
0048 // HEPEVT_EntriesAllocation is the maximum size of the HEPEVT common block
0049 //   that can be interfaced.
0050 //   It is NOT the actual size of the HEPEVT common used in each
0051 //   individual application. The actual size can be changed on
0052 //   the fly using HEPEVT_Wrapper::set_max_number_entries().
0053 // Thus HEPEVT_EntriesAllocation should typically be set
0054 // to the maximum possible number of entries --- 10000 is a good choice
0055 // (and is the number used by ATLAS versions of Pythia).
0056 //
0057 // Note: a statement like    *( (int*)&hepevt.data[0] )
0058 //      takes the memory address of the first byte in HEPEVT,
0059 //      interprets it as an integer pointer,
0060 //      and dereferences the pointer.
0061 //      i.e. it returns an integer corresponding to nevhep
0062 //
0063
0064 #include <ctype.h>
0065
0066     const unsigned int hepevt_bytes_allocation =
0067                 sizeof(long int) * ( 2 + 6 * HEPEVT_EntriesAllocation )
0068                 + sizeof(double) * ( 9 * HEPEVT_EntriesAllocation );
0069
0070
0071 #ifdef _WIN32 // Platform: Windows MS Visual C++
0072 struct HEPEVT_DEF{
0073         char data[hepevt_bytes_allocation];
0074     };
0075 extern "C" HEPEVT_DEF HEPEVT;
0076 #define hepevt HEPEVT
0077
0078 #else
0079 extern "C" {
0080     extern struct {
0081     char data[hepevt_bytes_allocation];
0082     } hepevt_;
0083 }
0084 #define hepevt hepevt_
0085
0086 #endif // Platform
0087
0088 #endif  // HEPMC_HEPEVT_COMMON_H
0089
0090 //--------------------------------------------------------------------------
0091 #ifndef HEPMC_HEPEVT_WRAPPER_H
0092 #define HEPMC_HEPEVT_WRAPPER_H
0093
0094 //////////////////////////////////////////////////////////////////////////
0095 // Matt.Dobbs@Cern.CH, April 24, 2000, refer to:
0096 // M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
0097 // High Energy Physics", Computer Physics Communications (to be published).
0098 //
0099 // Generic Wrapper for the fortran HEPEVT common block
0100 // This class is intended for static use only - it makes no sense to
0101 // instantiate it.
0102 // Updated: June 30, 2000 (static initialization moved to separate .cxx file)
0103 //////////////////////////////////////////////////////////////////////////
0104 //
0105 // The index refers to the fortran style index:
0106 // i.e. index=1 refers to the first entry in the HEPEVT common block.
0107 // all indices must be >0
0108 // number_entries --> integer between 0 and max_number_entries() giving total
0109 //                    number of sequential particle indices
0110 // first_parent/child --> index of first mother/child if there is one,
0111 //                        zero otherwise
0112 // last_parent/child --> if number children is >1, address of last parent/child
0113 //                       if number of children is 1, same as first_parent/child
0114 //                       if there are no children, returns zero.
0115 // is_double_precision --> T or F depending if floating point variables
0116 //                         are 8 or 4 bytes
0117 //
0118
0119 #include <iostream>
0120 #include <cstdio>       // needed for formatted output using sprintf
0121
0122 namespace HepMC {
0123
0124     //! Generic Wrapper for the fortran HEPEVT common block
0125
0126     /// \class HEPEVT_Wrapper
0127     /// This class is intended for static use only - it makes no sense to
0128     /// instantiate it.
0129     ///
0130     class HEPEVT_Wrapper {
0131     public:
0132
0133         /// write information from HEPEVT common block
0134     static void print_hepevt( std::ostream& ostr = std::cout );
0135     /// write particle information to ostr
0136     static void print_hepevt_particle( int index,
0137                        std::ostream& ostr = std::cout );
0138         static bool is_double_precision();  //!< True if common block uses double
0139
0140         /// check for problems with HEPEVT common block
0141     static bool check_hepevt_consistency( std::ostream& ostr = std::cout );
0142
0143         /// set all entries in HEPEVT to zero
0144     static void zero_everything();
0145
0146     ////////////////////
0147     // Access Methods //
0148     ////////////////////
0149         static int    event_number();             //!< event number
0150         static int    number_entries();           //!< num entries in current evt
0151         static int    status( int index );        //!< status code
0152         static int    id( int index );            //!< PDG particle id
0153         static int    first_parent( int index );  //!< index of 1st mother
0154         static int    last_parent( int index );   //!< index of last mother
0155     static int    number_parents( int index ); //!< number of parents
0156         static int    first_child( int index );   //!< index of 1st daughter
0157         static int    last_child( int index );    //!< index of last daughter
0158     static int    number_children( int index ); //!< number of children
0159         static double px( int index );            //!< X momentum
0160         static double py( int index );            //!< Y momentum
0161         static double pz( int index );            //!< Z momentum
0162         static double e( int index );             //!< Energy
0163         static double m( int index );             //!< generated mass
0164         static double x( int index );             //!< X Production vertex
0165         static double y( int index );             //!< Y Production vertex
0166         static double z( int index );             //!< Z Production vertex
0167         static double t( int index );             //!< production time
0168
0169     ////////////////////
0170     // Set Methods    //
0171     ////////////////////
0172
0173     /// set event number
0174         static void set_event_number( int evtno );
0175     /// set number of entries in HEPEVT
0176         static void set_number_entries( int noentries );
0177     /// set particle status
0178         static void set_status( int index, int status );
0179     /// set particle ID
0180         static void set_id( int index, int id );
0181     /// define parents of a particle
0182         static void set_parents( int index, int firstparent, int lastparent );
0183     /// define children of a particle
0184         static void set_children( int index, int firstchild, int lastchild );
0185     /// set particle momentum
0186         static void set_momentum( int index, double px, double py,
0187                   double pz, double e );
0188     /// set particle mass
0189         static void set_mass( int index, double mass );
0190     /// set particle production vertex
0191         static void set_position( int index, double x, double y, double z,
0192                   double t );
0193     //////////////////////
0194     // HEPEVT Floorplan //
0195     //////////////////////
0196     static unsigned int sizeof_int();  //!< size of integer in bytes
0197     static unsigned int sizeof_real(); //!< size of real in bytes
0198         static int  max_number_entries();  //!< size of common block
0199     static void set_sizeof_int(unsigned int);  //!< define size of integer
0200     static void set_sizeof_real(unsigned int); //!< define size of real
0201     static void set_max_number_entries(unsigned int); //!< define size of common block
0202
0203     protected:
0204         /// navigate a byte array
0205     static double byte_num_to_double( unsigned int );
0206         /// navigate a byte array
0207     static int    byte_num_to_int( unsigned int );
0208         /// pretend common block is an array of bytes
0209     static void   write_byte_num( double, unsigned int );
0210         /// pretend common block is an array of bytes
0211     static void   write_byte_num( int, unsigned int );
0212     /// print output legend
0213     static void   print_legend( std::ostream& ostr = std::cout );
0214
0215     private:
0216     static unsigned int s_sizeof_int;
0217     static unsigned int s_sizeof_real;
0218     static unsigned int s_max_number_entries;
0219
0220     };
0221
0222     //////////////////////////////
0223     // HEPEVT Floorplan Inlines //
0224     //////////////////////////////
0225     inline unsigned int HEPEVT_Wrapper::sizeof_int(){ return s_sizeof_int; }
0226
0227     inline unsigned int HEPEVT_Wrapper::sizeof_real(){ return s_sizeof_real; }
0228
0229     inline int HEPEVT_Wrapper::max_number_entries()
0230     { return (int)s_max_number_entries; }
0231
0232     inline void HEPEVT_Wrapper::set_sizeof_int( unsigned int size )
0233     {
0234     if ( size != sizeof(short int) && size != sizeof(long int) && size != sizeof(int) ) {
0235         std::cerr << "HepMC is not able to handle integers "
0236               << " of size other than 2 or 4."
0237               << " You requested: " << size << std::endl;
0238     }
0239     s_sizeof_int = size;
0240     }
0241
0242     inline void HEPEVT_Wrapper::set_sizeof_real( unsigned int size ) {
0243     if ( size != sizeof(float) && size != sizeof(double) ) {
0244         std::cerr << "HepMC is not able to handle floating point numbers"
0245               << " of size other than 4 or 8."
0246               << " You requested: " << size << std::endl;
0247     }
0248     s_sizeof_real = size;
0249     }
0250
0251     inline void HEPEVT_Wrapper::set_max_number_entries( unsigned int size ) {
0252     s_max_number_entries = size;
0253     }
0254
0255     inline double HEPEVT_Wrapper::byte_num_to_double( unsigned int b ) {
0256     if ( b >= hepevt_bytes_allocation ) std::cerr
0257           << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
0258           << std::endl;
0259     if ( s_sizeof_real == sizeof(float) ) {
0260         float* myfloat = (float*)&hepevt.data[b];
0261         return (double)(*myfloat);
0262     } else if ( s_sizeof_real == sizeof(double) ) {
0263         double* mydouble = (double*)&hepevt.data[b];
0264         return (*mydouble);
0265     } else {
0266         std::cerr
0267         << "HEPEVT_Wrapper: illegal floating point number length."
0268         << s_sizeof_real << std::endl;
0269     }
0270     return 0;
0271     }
0272
0273     inline int HEPEVT_Wrapper::byte_num_to_int( unsigned int b ) {
0274     if ( b >= hepevt_bytes_allocation ) std::cerr
0275           << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
0276           << std::endl;
0277     if ( s_sizeof_int == sizeof(short int) ) {
0278         short int* myshortint = (short int*)&hepevt.data[b];
0279         return (int)(*myshortint);
0280     } else if ( s_sizeof_int == sizeof(long int) ) {
0281         long int* mylongint = (long int*)&hepevt.data[b];
0282         return (*mylongint);
0283        // on some 64 bit machines, int, short, and long are all different
0284     } else if ( s_sizeof_int == sizeof(int) ) {
0285         int* myint = (int*)&hepevt.data[b];
0286         return (*myint);
0287     } else {
0288         std::cerr
0289         << "HEPEVT_Wrapper: illegal integer number length."
0290         << s_sizeof_int << std::endl;
0291     }
0292     return 0;
0293     }
0294
0295     inline void HEPEVT_Wrapper::write_byte_num( double in, unsigned int b ) {
0296     if ( b >= hepevt_bytes_allocation ) std::cerr
0297           << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
0298           << std::endl;
0299     if ( s_sizeof_real == sizeof(float) ) {
0300         float* myfloat = (float*)&hepevt.data[b];
0301         (*myfloat) = (float)in;
0302     } else if ( s_sizeof_real == sizeof(double) ) {
0303         double* mydouble = (double*)&hepevt.data[b];
0304         (*mydouble) = (double)in;
0305     } else {
0306         std::cerr
0307         << "HEPEVT_Wrapper: illegal floating point number length."
0308         << s_sizeof_real << std::endl;
0309     }
0310     }
0311
0312     inline void HEPEVT_Wrapper::write_byte_num( int in, unsigned int b ) {
0313     if ( b >= hepevt_bytes_allocation ) std::cerr
0314           << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
0315           << std::endl;
0316     if ( s_sizeof_int == sizeof(short int) ) {
0317         short int* myshortint = (short int*)&hepevt.data[b];
0318         (*myshortint) = (short int)in;
0319     } else if ( s_sizeof_int == sizeof(long int) ) {
0320         long int* mylongint = (long int*)&hepevt.data[b];
0321         (*mylongint) = (int)in;
0322        // on some 64 bit machines, int, short, and long are all different
0323     } else if ( s_sizeof_int == sizeof(int) ) {
0324         int* myint = (int*)&hepevt.data[b];
0325         (*myint) = (int)in;
0326     } else {
0327         std::cerr
0328         << "HEPEVT_Wrapper: illegal integer number length."
0329         << s_sizeof_int << std::endl;
0330     }
0331     }
0332
0333     //////////////
0334     // INLINES  //
0335     //////////////
0336
0337     inline bool HEPEVT_Wrapper::is_double_precision()
0338     {
0339     // true if 8byte floating point numbers are used in the HepEVT common.
0340     return ( sizeof(double) == sizeof_real() );
0341     }
0342
0343     inline int HEPEVT_Wrapper::event_number()
0344     { return byte_num_to_int(0); }
0345
0346     inline int HEPEVT_Wrapper::number_entries()
0347     {
0348     int nhep = byte_num_to_int( 1*sizeof_int() );
0349     return ( nhep <= max_number_entries() ?
0350          nhep : max_number_entries() );
0351     }
0352
0353     inline int HEPEVT_Wrapper::status( int index )
0354     { return byte_num_to_int( (2+index-1) * sizeof_int() ); }
0355
0356     inline int HEPEVT_Wrapper::id( int index )
0357     {
0358     return byte_num_to_int( (2+max_number_entries()+index-1)
0359                 * sizeof_int() );
0360     }
0361
0362     inline int HEPEVT_Wrapper::first_parent( int index )
0363     {
0364     int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1))
0365                       * sizeof_int() );
0366     return ( parent > 0 && parent <= number_entries() ) ?
0367                      parent : 0;
0368     }
0369
0370     inline int HEPEVT_Wrapper::last_parent( int index )
0371     {
0372     // Returns the Index of the LAST parent in the HEPEVT record
0373     // for particle with Index index.
0374     // If there is only one parent, the last parent is forced to
0375     // be the same as the first parent.
0376     // If there are no parents for this particle, both the first_parent
0377     // and the last_parent with return 0.
0378     // Error checking is done to ensure the parent is always
0379     // within range ( 0 <= parent <= nhep )
0380     //
0381     int firstparent = first_parent(index);
0382     int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1)+1)
0383                       * sizeof_int() );
0384     return ( parent > firstparent && parent <= number_entries() )
0385                            ? parent : firstparent;
0386     }
0387
0388     inline int HEPEVT_Wrapper::number_parents( int index ) {
0389     int firstparent = first_parent(index);
0390     return ( firstparent>0 ) ?
0391         ( 1+last_parent(index)-firstparent ) : 0;
0392     }
0393
0394     inline int HEPEVT_Wrapper::first_child( int index )
0395     {
0396     int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1))
0397                      * sizeof_int() );
0398     return ( child > 0 && child <= number_entries() ) ?
0399                        child : 0;
0400     }
0401
0402     inline int HEPEVT_Wrapper::last_child( int index )
0403     {
0404     // Returns the Index of the LAST child in the HEPEVT record
0405     // for particle with Index index.
0406     // If there is only one child, the last child is forced to
0407     // be the same as the first child.
0408     // If there are no children for this particle, both the first_child
0409     // and the last_child with return 0.
0410     // Error checking is done to ensure the child is always
0411     // within range ( 0 <= parent <= nhep )
0412     //
0413     int firstchild = first_child(index);
0414     int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1)+1)
0415                      * sizeof_int() );
0416     return ( child > firstchild && child <= number_entries() )
0417                         ? child : firstchild;
0418     }
0419
0420     inline int HEPEVT_Wrapper::number_children( int index )
0421     {
0422     int firstchild = first_child(index);
0423     return ( firstchild>0 ) ?
0424         ( 1+last_child(index)-firstchild ) : 0;
0425     }
0426
0427     inline double HEPEVT_Wrapper::px( int index )
0428     {
0429     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0430                  + (5*(index-1)+0) *sizeof_real() );
0431     }
0432
0433     inline double HEPEVT_Wrapper::py( int index )
0434     {
0435     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0436                  + (5*(index-1)+1) *sizeof_real() );
0437     }
0438
0439
0440     inline double HEPEVT_Wrapper::pz( int index )
0441     {
0442     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0443                  + (5*(index-1)+2) *sizeof_real() );
0444     }
0445
0446     inline double HEPEVT_Wrapper::e( int index )
0447     {
0448     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0449                  + (5*(index-1)+3) *sizeof_real() );
0450     }
0451
0452     inline double HEPEVT_Wrapper::m( int index )
0453     {
0454     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0455                  + (5*(index-1)+4) *sizeof_real() );
0456     }
0457
0458     inline double HEPEVT_Wrapper::x( int index )
0459     {
0460     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0461                    + ( 5*max_number_entries()
0462                        + (4*(index-1)+0) ) *sizeof_real() );
0463     }
0464
0465     inline double HEPEVT_Wrapper::y( int index )
0466     {
0467     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0468                    + ( 5*max_number_entries()
0469                        + (4*(index-1)+1) ) *sizeof_real() );
0470     }
0471
0472     inline double HEPEVT_Wrapper::z( int index )
0473     {
0474     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0475                    + ( 5*max_number_entries()
0476                        + (4*(index-1)+2) ) *sizeof_real() );
0477     }
0478
0479     inline double HEPEVT_Wrapper::t( int index )
0480     {
0481     return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
0482                    + ( 5*max_number_entries()
0483                        + (4*(index-1)+3) ) *sizeof_real() );
0484     }
0485
0486     inline void HEPEVT_Wrapper::set_event_number( int evtno )
0487     { write_byte_num( evtno, 0 ); }
0488
0489     inline void HEPEVT_Wrapper::set_number_entries( int noentries )
0490     { write_byte_num( noentries, 1*sizeof_int() ); }
0491
0492     inline void HEPEVT_Wrapper::set_status( int index, int status )
0493     {
0494         if ( index <= 0 || index > max_number_entries() ) return;
0495     write_byte_num( status, (2+index-1) * sizeof_int() );
0496     }
0497
0498     inline void HEPEVT_Wrapper::set_id( int index, int id )
0499     {
0500         if ( index <= 0 || index > max_number_entries() ) return;
0501     write_byte_num( id, (2+max_number_entries()+index-1) *sizeof_int() );
0502     }
0503
0504     inline void HEPEVT_Wrapper::set_parents( int index, int firstparent,
0505                          int lastparent )
0506     {
0507         if ( index <= 0 || index > max_number_entries() ) return;
0508     write_byte_num( firstparent, (2+2*max_number_entries()+2*(index-1))
0509                          *sizeof_int() );
0510     write_byte_num( lastparent, (2+2*max_number_entries()+2*(index-1)+1)
0511                     * sizeof_int() );
0512     }
0513
0514     inline void HEPEVT_Wrapper::set_children( int index, int firstchild,
0515                           int lastchild )
0516     {
0517         if ( index <= 0 || index > max_number_entries() ) return;
0518     write_byte_num( firstchild, (2+4*max_number_entries()+2*(index-1))
0519                      *sizeof_int() );
0520     write_byte_num( lastchild, (2+4*max_number_entries()+2*(index-1)+1)
0521                     *sizeof_int() );
0522     }
0523
0524     inline void HEPEVT_Wrapper::set_momentum( int index, double px,
0525                           double py, double pz, double e )
0526     {
0527         if ( index <= 0 || index > max_number_entries() ) return;
0528     write_byte_num( px, (2+6*max_number_entries()) *sizeof_int()
0529                 + (5*(index-1)+0) *sizeof_real() );
0530     write_byte_num( py, (2+6*max_number_entries())*sizeof_int()
0531                 + (5*(index-1)+1) *sizeof_real() );
0532     write_byte_num( pz, (2+6*max_number_entries())*sizeof_int()
0533                 + (5*(index-1)+2) *sizeof_real() );
0534     write_byte_num( e,  (2+6*max_number_entries())*sizeof_int()
0535                 + (5*(index-1)+3) *sizeof_real() );
0536     }
0537
0538     inline void HEPEVT_Wrapper::set_mass( int index, double mass )
0539     {
0540         if ( index <= 0 || index > max_number_entries() ) return;
0541     write_byte_num( mass, (2+6*max_number_entries())*sizeof_int()
0542                   + (5*(index-1)+4) *sizeof_real() );
0543     }
0544
0545     inline void HEPEVT_Wrapper::set_position( int index, double x, double y,
0546                           double z, double t )
0547     {
0548         if ( index <= 0 || index > max_number_entries() ) return;
0549     write_byte_num( x, (2+6*max_number_entries())*sizeof_int()
0550                + ( 5*max_number_entries()
0551                    + (4*(index-1)+0) ) *sizeof_real() );
0552     write_byte_num( y, (2+6*max_number_entries())*sizeof_int()
0553                + ( 5*max_number_entries()
0554                    + (4*(index-1)+1) ) *sizeof_real() );
0555     write_byte_num( z, (2+6*max_number_entries())*sizeof_int()
0556                + ( 5*max_number_entries()
0557                    + (4*(index-1)+2) ) *sizeof_real() );
0558     write_byte_num( t, (2+6*max_number_entries())*sizeof_int()
0559                + ( 5*max_number_entries()
0560                    + (4*(index-1)+3) ) *sizeof_real() );
0561     }
0562
0563 } // HepMC
0564
0565 #endif  // HEPMC_HEPEVT_WRAPPER_H
0566 //--------------------------------------------------------------------------