[u/mrichter/AliRoot.git] / TEvtGen / HepMC / HEPEVT_Wrapper.h

//--------------------------------------------------------------------------

#ifndef HEPEVT_EntriesAllocation
#define HEPEVT_EntriesAllocation 20000
#endif  // HEPEVT_EntriesAllocation

//--------------------------------------------------------------------------
#ifndef HEPMC_HEPEVT_COMMON_H
#define HEPMC_HEPEVT_COMMON_H
//////////////////////////////////////////////////////////////////////////
//
//      PARAMETER (NMXHEP=2000) 
//      COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), 
//     &        JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
/**********************************************************/
/*           D E S C R I P T I O N :                      */
/*--------------------------------------------------------*/
/* NEVHEP          - event number (or some special meaning*/
/*                    (see documentation for details)     */
/* NHEP            - actual number of entries in current  */
/*                    event.                              */
/* ISTHEP[IHEP]    - status code for IHEP'th entry - see  */
/*                    documentation for details           */
/* IDHEP [IHEP]    - IHEP'th particle identifier according*/
/*                    to PDG.                             */
/* JMOHEP[IHEP][0] - pointer to position of 1st mother    */
/* JMOHEP[IHEP][1] - pointer to position of 2nd mother    */
/* JDAHEP[IHEP][0] - pointer to position of 1st daughter  */
/* JDAHEP[IHEP][1] - pointer to position of 2nd daughter  */
/* PHEP  [IHEP][0] - X momentum                           */
/* PHEP  [IHEP][1] - Y momentum                           */
/* PHEP  [IHEP][2] - Z momentum                           */
/* PHEP  [IHEP][3] - Energy                               */
/* PHEP  [IHEP][4] - Mass                                 */
/* VHEP  [IHEP][0] - X vertex                             */
/* VHEP  [IHEP][1] - Y vertex                             */
/* VHEP  [IHEP][2] - Z vertex                             */
/* VHEP  [IHEP][3] - production time                      */
/*========================================================*/
// Remember, array(1) is the first entry in a fortran array, array[0] is the
//           first entry in a C array.
//
// This interface to HEPEVT common block treats the block as
// an array of bytes --- the precision and number of entries 
// is determined "on the fly" by the wrapper and used to decode
// each entry.
//
// HEPEVT_EntriesAllocation is the maximum size of the HEPEVT common block 
//   that can be interfaced.
//   It is NOT the actual size of the HEPEVT common used in each
//   individual application. The actual size can be changed on
//   the fly using HEPEVT_Wrapper::set_max_number_entries().
// Thus HEPEVT_EntriesAllocation should typically be set
// to the maximum possible number of entries --- 10000 is a good choice
// (and is the number used by ATLAS versions of Pythia).
//
// Note: a statement like    *( (int*)&hepevt.data[0] )
//      takes the memory address of the first byte in HEPEVT,
//      interprets it as an integer pointer, 
//      and dereferences the pointer.
//      i.e. it returns an integer corresponding to nevhep
//

#include <ctype.h>

    const unsigned int hepevt_bytes_allocation = 
                sizeof(long int) * ( 2 + 6 * HEPEVT_EntriesAllocation )
                + sizeof(double) * ( 9 * HEPEVT_EntriesAllocation );


#ifdef _WIN32 // Platform: Windows MS Visual C++
struct HEPEVT_DEF{
        char data[hepevt_bytes_allocation];
    };
extern "C" HEPEVT_DEF HEPEVT;
#define hepevt HEPEVT

#else
extern "C" {
    extern struct {
	char data[hepevt_bytes_allocation];
    } hepevt_;
}
#define hepevt hepevt_

#endif // Platform

#endif  // HEPMC_HEPEVT_COMMON_H

//--------------------------------------------------------------------------
#ifndef HEPMC_HEPEVT_WRAPPER_H
#define HEPMC_HEPEVT_WRAPPER_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, April 24, 2000, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Generic Wrapper for the fortran HEPEVT common block
// This class is intended for static use only - it makes no sense to 
// instantiate it.
// Updated: June 30, 2000 (static initialization moved to separate .cxx file)
//////////////////////////////////////////////////////////////////////////
//
// The index refers to the fortran style index: 
// i.e. index=1 refers to the first entry in the HEPEVT common block.
// all indices must be >0
// number_entries --> integer between 0 and max_number_entries() giving total
//                    number of sequential particle indices
// first_parent/child --> index of first mother/child if there is one, 
//                        zero otherwise
// last_parent/child --> if number children is >1, address of last parent/child
//                       if number of children is 1, same as first_parent/child
//                       if there are no children, returns zero.
// is_double_precision --> T or F depending if floating point variables 
//                         are 8 or 4 bytes
//

#include <iostream>
#include <cstdio>       // needed for formatted output using sprintf 

namespace HepMC {

    //! Generic Wrapper for the fortran HEPEVT common block
    
    /// \class HEPEVT_Wrapper
    /// This class is intended for static use only - it makes no sense to 
    /// instantiate it.
    ///
    class HEPEVT_Wrapper {
    public:

        /// write information from HEPEVT common block
	static void print_hepevt( std::ostream& ostr = std::cout );
	/// write particle information to ostr
	static void print_hepevt_particle( int index, 
					   std::ostream& ostr = std::cout );
        static bool is_double_precision();  //!< True if common block uses double

        /// check for problems with HEPEVT common block
	static bool check_hepevt_consistency( std::ostream& ostr = std::cout );

        /// set all entries in HEPEVT to zero
	static void zero_everything();

	////////////////////
	// Access Methods //
	////////////////////
        static int    event_number();             //!< event number
        static int    number_entries();           //!< num entries in current evt
        static int    status( int index );        //!< status code
        static int    id( int index );            //!< PDG particle id
        static int    first_parent( int index );  //!< index of 1st mother
        static int    last_parent( int index );   //!< index of last mother
	static int    number_parents( int index ); //!< number of parents 
        static int    first_child( int index );   //!< index of 1st daughter
        static int    last_child( int index );    //!< index of last daughter
	static int    number_children( int index ); //!< number of children
        static double px( int index );            //!< X momentum       
        static double py( int index );            //!< Y momentum   
        static double pz( int index );            //!< Z momentum   
        static double e( int index );             //!< Energy
        static double m( int index );             //!< generated mass
        static double x( int index );             //!< X Production vertex
        static double y( int index );             //!< Y Production vertex
        static double z( int index );             //!< Z Production vertex
        static double t( int index );             //!< production time

	////////////////////
	// Set Methods    //
	////////////////////

	/// set event number
        static void set_event_number( int evtno );
	/// set number of entries in HEPEVT
        static void set_number_entries( int noentries );
	/// set particle status
        static void set_status( int index, int status );
	/// set particle ID
        static void set_id( int index, int id );
	/// define parents of a particle
        static void set_parents( int index, int firstparent, int lastparent );
	/// define children of a particle
        static void set_children( int index, int firstchild, int lastchild );
	/// set particle momentum
        static void set_momentum( int index, double px, double py,
				  double pz, double e );
	/// set particle mass
        static void set_mass( int index, double mass );
	/// set particle production vertex
        static void set_position( int index, double x, double y, double z, 
				  double t );
	//////////////////////
	// HEPEVT Floorplan //
	//////////////////////
	static unsigned int sizeof_int();  //!< size of integer in bytes
	static unsigned int sizeof_real(); //!< size of real in bytes
        static int  max_number_entries();  //!< size of common block
	static void set_sizeof_int(unsigned int);  //!< define size of integer
	static void set_sizeof_real(unsigned int); //!< define size of real
	static void set_max_number_entries(unsigned int); //!< define size of common block

    protected:
        /// navigate a byte array
	static double byte_num_to_double( unsigned int );
        /// navigate a byte array
	static int    byte_num_to_int( unsigned int );
        /// pretend common block is an array of bytes
	static void   write_byte_num( double, unsigned int );
        /// pretend common block is an array of bytes
	static void   write_byte_num( int, unsigned int );
	/// print output legend
	static void   print_legend( std::ostream& ostr = std::cout );

    private:
	static unsigned int s_sizeof_int;
	static unsigned int s_sizeof_real;
	static unsigned int s_max_number_entries;

    }; 

    //////////////////////////////
    // HEPEVT Floorplan Inlines //
    //////////////////////////////
    inline unsigned int HEPEVT_Wrapper::sizeof_int(){ return s_sizeof_int; }

    inline unsigned int HEPEVT_Wrapper::sizeof_real(){ return s_sizeof_real; }

    inline int HEPEVT_Wrapper::max_number_entries() 
    { return (int)s_max_number_entries; }

    inline void HEPEVT_Wrapper::set_sizeof_int( unsigned int size ) 
    {
	if ( size != sizeof(short int) && size != sizeof(long int) && size != sizeof(int) ) {
	    std::cerr << "HepMC is not able to handle integers "
		      << " of size other than 2 or 4."
		      << " You requested: " << size << std::endl;
	}
	s_sizeof_int = size;
    }

    inline void HEPEVT_Wrapper::set_sizeof_real( unsigned int size ) {
	if ( size != sizeof(float) && size != sizeof(double) ) {
	    std::cerr << "HepMC is not able to handle floating point numbers"
		      << " of size other than 4 or 8."
		      << " You requested: " << size << std::endl;
	}
	s_sizeof_real = size;
    }

    inline void HEPEVT_Wrapper::set_max_number_entries( unsigned int size ) {
	s_max_number_entries = size;
    }

    inline double HEPEVT_Wrapper::byte_num_to_double( unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_real == sizeof(float) ) {
	    float* myfloat = (float*)&hepevt.data[b];
	    return (double)(*myfloat);
	} else if ( s_sizeof_real == sizeof(double) ) {
	    double* mydouble = (double*)&hepevt.data[b];
	    return (*mydouble);
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal floating point number length." 
		<< s_sizeof_real << std::endl;
	}
	return 0;
    }

    inline int HEPEVT_Wrapper::byte_num_to_int( unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_int == sizeof(short int) ) {
	    short int* myshortint = (short int*)&hepevt.data[b];
	    return (int)(*myshortint);
	} else if ( s_sizeof_int == sizeof(long int) ) {
	    long int* mylongint = (long int*)&hepevt.data[b];
	    return (*mylongint);
       // on some 64 bit machines, int, short, and long are all different
	} else if ( s_sizeof_int == sizeof(int) ) {
	    int* myint = (int*)&hepevt.data[b];
	    return (*myint);
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal integer number length." 
		<< s_sizeof_int << std::endl;
	}
	return 0;
    }

    inline void HEPEVT_Wrapper::write_byte_num( double in, unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_real == sizeof(float) ) {
	    float* myfloat = (float*)&hepevt.data[b];
	    (*myfloat) = (float)in;
	} else if ( s_sizeof_real == sizeof(double) ) {
	    double* mydouble = (double*)&hepevt.data[b];
	    (*mydouble) = (double)in;
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal floating point number length." 
		<< s_sizeof_real << std::endl;
	}
    }

    inline void HEPEVT_Wrapper::write_byte_num( int in, unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_int == sizeof(short int) ) {
	    short int* myshortint = (short int*)&hepevt.data[b];
	    (*myshortint) = (short int)in;
	} else if ( s_sizeof_int == sizeof(long int) ) {
	    long int* mylongint = (long int*)&hepevt.data[b];
	    (*mylongint) = (int)in;
       // on some 64 bit machines, int, short, and long are all different
	} else if ( s_sizeof_int == sizeof(int) ) {
	    int* myint = (int*)&hepevt.data[b];
	    (*myint) = (int)in;
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal integer number length." 
		<< s_sizeof_int << std::endl;
	}
    }

    //////////////
    // INLINES  //
    //////////////

    inline bool HEPEVT_Wrapper::is_double_precision() 
    { 
	// true if 8byte floating point numbers are used in the HepEVT common.
	return ( sizeof(double) == sizeof_real() );
    }

    inline int HEPEVT_Wrapper::event_number()
    { return byte_num_to_int(0); }

    inline int HEPEVT_Wrapper::number_entries() 
    { 
	int nhep = byte_num_to_int( 1*sizeof_int() );
	return ( nhep <= max_number_entries() ?
		 nhep : max_number_entries() );
    }

    inline int HEPEVT_Wrapper::status( int index )   
    { return byte_num_to_int( (2+index-1) * sizeof_int() ); }

    inline int HEPEVT_Wrapper::id( int index )
    { 
	return byte_num_to_int( (2+max_number_entries()+index-1) 
				* sizeof_int() ); 
    }

    inline int HEPEVT_Wrapper::first_parent( int index )
    { 
	int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1)) 
				      * sizeof_int() ); 
	return ( parent > 0 && parent <= number_entries() ) ?
					 parent : 0; 
    }

    inline int HEPEVT_Wrapper::last_parent( int index )
    { 
	// Returns the Index of the LAST parent in the HEPEVT record
	// for particle with Index index.
	// If there is only one parent, the last parent is forced to 
	// be the same as the first parent.
	// If there are no parents for this particle, both the first_parent
	// and the last_parent with return 0.
	// Error checking is done to ensure the parent is always
	// within range ( 0 <= parent <= nhep )
	//
	int firstparent = first_parent(index);
	int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1)+1) 
				      * sizeof_int() ); 
	return ( parent > firstparent && parent <= number_entries() ) 
						   ? parent : firstparent; 
    }

    inline int HEPEVT_Wrapper::number_parents( int index ) {
	int firstparent = first_parent(index);
	return ( firstparent>0 ) ? 
	    ( 1+last_parent(index)-firstparent ) : 0;
    }

    inline int HEPEVT_Wrapper::first_child( int index )
    { 
	int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1)) 
				     * sizeof_int() ); 
	return ( child > 0 && child <= number_entries() ) ?
				       child : 0; 
    }

    inline int HEPEVT_Wrapper::last_child( int index )
    { 
	// Returns the Index of the LAST child in the HEPEVT record
	// for particle with Index index.
	// If there is only one child, the last child is forced to 
	// be the same as the first child.
	// If there are no children for this particle, both the first_child
	// and the last_child with return 0.
	// Error checking is done to ensure the child is always
	// within range ( 0 <= parent <= nhep )
	//
	int firstchild = first_child(index);
	int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1)+1) 
				     * sizeof_int() ); 
	return ( child > firstchild && child <= number_entries() ) 
						? child : firstchild;
    }

    inline int HEPEVT_Wrapper::number_children( int index ) 
    {
	int firstchild = first_child(index);
	return ( firstchild>0 ) ? 
	    ( 1+last_child(index)-firstchild ) : 0;
    }

    inline double HEPEVT_Wrapper::px( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+0) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::py( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+1) *sizeof_real() );
    }


    inline double HEPEVT_Wrapper::pz( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+2) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::e( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+3) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::m( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+4) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::x( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+0) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::y( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+1) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::z( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+2) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::t( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+3) ) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_event_number( int evtno ) 
    { write_byte_num( evtno, 0 ); }

    inline void HEPEVT_Wrapper::set_number_entries( int noentries ) 
    { write_byte_num( noentries, 1*sizeof_int() ); }

    inline void HEPEVT_Wrapper::set_status( int index, int status ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( status, (2+index-1) * sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_id( int index, int id ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( id, (2+max_number_entries()+index-1) *sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_parents( int index, int firstparent, 
					     int lastparent ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( firstparent, (2+2*max_number_entries()+2*(index-1)) 
			             *sizeof_int() );
	write_byte_num( lastparent, (2+2*max_number_entries()+2*(index-1)+1) 
				    * sizeof_int() );
    }
    
    inline void HEPEVT_Wrapper::set_children( int index, int firstchild, 
					      int lastchild ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( firstchild, (2+4*max_number_entries()+2*(index-1)) 
				     *sizeof_int() );
	write_byte_num( lastchild, (2+4*max_number_entries()+2*(index-1)+1) 
				    *sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_momentum( int index, double px, 
					      double py, double pz, double e ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( px, (2+6*max_number_entries()) *sizeof_int()
			    + (5*(index-1)+0) *sizeof_real() );
	write_byte_num( py, (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+1) *sizeof_real() );
	write_byte_num( pz, (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+2) *sizeof_real() );
	write_byte_num( e,  (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+3) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_mass( int index, double mass ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( mass, (2+6*max_number_entries())*sizeof_int()
			      + (5*(index-1)+4) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_position( int index, double x, double y,
					      double z, double t ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( x, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+0) ) *sizeof_real() );
	write_byte_num( y, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+1) ) *sizeof_real() );
	write_byte_num( z, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+2) ) *sizeof_real() );
	write_byte_num( t, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+3) ) *sizeof_real() );
    }

} // HepMC

#endif  // HEPMC_HEPEVT_WRAPPER_H
//--------------------------------------------------------------------------
Commit	Line	Data
0ca57c2f	1	//--------------------------------------------------------------------------
	2
	3	#ifndef HEPEVT_EntriesAllocation
cb3eff76	4	#define HEPEVT_EntriesAllocation 20000
0ca57c2f	5	#endif // HEPEVT_EntriesAllocation
	6
	7	//--------------------------------------------------------------------------
	8	#ifndef HEPMC_HEPEVT_COMMON_H
	9	#define HEPMC_HEPEVT_COMMON_H
	10	//////////////////////////////////////////////////////////////////////////
	11	//
	12	// PARAMETER (NMXHEP=2000)
	13	// COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
	14	// & JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
	15	/**********************************************************/
	16	/* D E S C R I P T I O N : */
	17	/--------------------------------------------------------/
	18	/* NEVHEP - event number (or some special meaning*/
	19	/* (see documentation for details) */
	20	/* NHEP - actual number of entries in current */
	21	/* event. */
	22	/* ISTHEP[IHEP] - status code for IHEP'th entry - see */
	23	/* documentation for details */
	24	/* IDHEP [IHEP] - IHEP'th particle identifier according*/
	25	/* to PDG. */
	26	/* JMOHEP[IHEP][0] - pointer to position of 1st mother */
	27	/* JMOHEP[IHEP][1] - pointer to position of 2nd mother */
	28	/* JDAHEP[IHEP][0] - pointer to position of 1st daughter */
	29	/* JDAHEP[IHEP][1] - pointer to position of 2nd daughter */
	30	/* PHEP [IHEP][0] - X momentum */
	31	/* PHEP [IHEP][1] - Y momentum */
	32	/* PHEP [IHEP][2] - Z momentum */
	33	/* PHEP [IHEP][3] - Energy */
	34	/* PHEP [IHEP][4] - Mass */
	35	/* VHEP [IHEP][0] - X vertex */
	36	/* VHEP [IHEP][1] - Y vertex */
	37	/* VHEP [IHEP][2] - Z vertex */
	38	/* VHEP [IHEP][3] - production time */
	39	/========================================================/
	40	// Remember, array(1) is the first entry in a fortran array, array[0] is the
	41	// first entry in a C array.
	42	//
	43	// This interface to HEPEVT common block treats the block as
	44	// an array of bytes --- the precision and number of entries
	45	// is determined "on the fly" by the wrapper and used to decode
	46	// each entry.
	47	//
	48	// HEPEVT_EntriesAllocation is the maximum size of the HEPEVT common block
	49	// that can be interfaced.
	50	// It is NOT the actual size of the HEPEVT common used in each
	51	// individual application. The actual size can be changed on
	52	// the fly using HEPEVT_Wrapper::set_max_number_entries().
	53	// Thus HEPEVT_EntriesAllocation should typically be set
	54	// to the maximum possible number of entries --- 10000 is a good choice
	55	// (and is the number used by ATLAS versions of Pythia).
	56	//
	57	// Note: a statement like ( (int)&hepevt.data[0] )
	58	// takes the memory address of the first byte in HEPEVT,
	59	// interprets it as an integer pointer,
	60	// and dereferences the pointer.
	61	// i.e. it returns an integer corresponding to nevhep
	62	//
	63
	64	#include <ctype.h>
	65
	66	const unsigned int hepevt_bytes_allocation =
	67	sizeof(long int) * ( 2 + 6 * HEPEVT_EntriesAllocation )
	68	+ sizeof(double) * ( 9 * HEPEVT_EntriesAllocation );
69
70
71	#ifdef _WIN32 // Platform: Windows MS Visual C++
72	struct HEPEVT_DEF{
73	char data[hepevt_bytes_allocation];
74	};
75	extern "C" HEPEVT_DEF HEPEVT;
76	#define hepevt HEPEVT
77
78	#else
79	extern "C" {
80	extern struct {
81	char data[hepevt_bytes_allocation];
82	} hepevt_;
83	}
84	#define hepevt hepevt_
85
86	#endif // Platform
87
88	#endif // HEPMC_HEPEVT_COMMON_H
89
90	//--------------------------------------------------------------------------
91	#ifndef HEPMC_HEPEVT_WRAPPER_H
92	#define HEPMC_HEPEVT_WRAPPER_H
93
94	//////////////////////////////////////////////////////////////////////////
95	// Matt.Dobbs@Cern.CH, April 24, 2000, refer to:
96	// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
97	// High Energy Physics", Computer Physics Communications (to be published).
98	//
99	// Generic Wrapper for the fortran HEPEVT common block
100	// This class is intended for static use only - it makes no sense to
101	// instantiate it.
102	// Updated: June 30, 2000 (static initialization moved to separate .cxx file)
103	//////////////////////////////////////////////////////////////////////////
104	//
105	// The index refers to the fortran style index:
106	// i.e. index=1 refers to the first entry in the HEPEVT common block.
107	// all indices must be >0
108	// number_entries --> integer between 0 and max_number_entries() giving total
109	// number of sequential particle indices
110	// first_parent/child --> index of first mother/child if there is one,
111	// zero otherwise
112	// last_parent/child --> if number children is >1, address of last parent/child
113	// if number of children is 1, same as first_parent/child
114	// if there are no children, returns zero.
115	// is_double_precision --> T or F depending if floating point variables
116	// are 8 or 4 bytes
117	//
118
119	#include <iostream>
120	#include <cstdio> // needed for formatted output using sprintf
121
122	namespace HepMC {
123
124	//! Generic Wrapper for the fortran HEPEVT common block
125
126	/// \class HEPEVT_Wrapper
127	/// This class is intended for static use only - it makes no sense to
128	/// instantiate it.
129	///
130	class HEPEVT_Wrapper {
131	public:
132
133	/// write information from HEPEVT common block
134	static void print_hepevt( std::ostream& ostr = std::cout );
135	/// write particle information to ostr
136	static void print_hepevt_particle( int index,
137	std::ostream& ostr = std::cout );
138	static bool is_double_precision(); //!< True if common block uses double
139
140	/// check for problems with HEPEVT common block
141	static bool check_hepevt_consistency( std::ostream& ostr = std::cout );
142
143	/// set all entries in HEPEVT to zero
144	static void zero_everything();
145
146	////////////////////
147	// Access Methods //
148	////////////////////
149	static int event_number(); //!< event number
150	static int number_entries(); //!< num entries in current evt
151	static int status( int index ); //!< status code
152	static int id( int index ); //!< PDG particle id
153	static int first_parent( int index ); //!< index of 1st mother
154	static int last_parent( int index ); //!< index of last mother
155	static int number_parents( int index ); //!< number of parents
156	static int first_child( int index ); //!< index of 1st daughter
157	static int last_child( int index ); //!< index of last daughter
158	static int number_children( int index ); //!< number of children
159	static double px( int index ); //!< X momentum
160	static double py( int index ); //!< Y momentum
161	static double pz( int index ); //!< Z momentum
162	static double e( int index ); //!< Energy
163	static double m( int index ); //!< generated mass
164	static double x( int index ); //!< X Production vertex
165	static double y( int index ); //!< Y Production vertex
166	static double z( int index ); //!< Z Production vertex
167	static double t( int index ); //!< production time
168
169	////////////////////
170	// Set Methods //
171	////////////////////
172
173	/// set event number
174	static void set_event_number( int evtno );
175	/// set number of entries in HEPEVT
176	static void set_number_entries( int noentries );
177	/// set particle status
178	static void set_status( int index, int status );
179	/// set particle ID
180	static void set_id( int index, int id );
181	/// define parents of a particle
182	static void set_parents( int index, int firstparent, int lastparent );
183	/// define children of a particle
184	static void set_children( int index, int firstchild, int lastchild );
185	/// set particle momentum
186	static void set_momentum( int index, double px, double py,
187	double pz, double e );
188	/// set particle mass
189	static void set_mass( int index, double mass );
190	/// set particle production vertex
191	static void set_position( int index, double x, double y, double z,
192	double t );
193	//////////////////////
194	// HEPEVT Floorplan //
195	//////////////////////
196	static unsigned int sizeof_int(); //!< size of integer in bytes
197	static unsigned int sizeof_real(); //!< size of real in bytes
198	static int max_number_entries(); //!< size of common block
199	static void set_sizeof_int(unsigned int); //!< define size of integer
200	static void set_sizeof_real(unsigned int); //!< define size of real
201	static void set_max_number_entries(unsigned int); //!< define size of common block
202
203	protected:
204	/// navigate a byte array
205	static double byte_num_to_double( unsigned int );
206	/// navigate a byte array
207	static int byte_num_to_int( unsigned int );
208	/// pretend common block is an array of bytes
209	static void write_byte_num( double, unsigned int );
210	/// pretend common block is an array of bytes
211	static void write_byte_num( int, unsigned int );
212	/// print output legend
213	static void print_legend( std::ostream& ostr = std::cout );
214
215	private:
216	static unsigned int s_sizeof_int;
217	static unsigned int s_sizeof_real;
218	static unsigned int s_max_number_entries;
219
220	};
221
222	//////////////////////////////
223	// HEPEVT Floorplan Inlines //
224	//////////////////////////////
225	inline unsigned int HEPEVT_Wrapper::sizeof_int(){ return s_sizeof_int; }
226
227	inline unsigned int HEPEVT_Wrapper::sizeof_real(){ return s_sizeof_real; }
228
229	inline int HEPEVT_Wrapper::max_number_entries()
230	{ return (int)s_max_number_entries; }
231
232	inline void HEPEVT_Wrapper::set_sizeof_int( unsigned int size )
233	{
234	if ( size != sizeof(short int) && size != sizeof(long int) && size != sizeof(int) ) {
235	std::cerr << "HepMC is not able to handle integers "
236	<< " of size other than 2 or 4."
237	<< " You requested: " << size << std::endl;
238	}
239	s_sizeof_int = size;
240	}
241
242	inline void HEPEVT_Wrapper::set_sizeof_real( unsigned int size ) {
243	if ( size != sizeof(float) && size != sizeof(double) ) {
244	std::cerr << "HepMC is not able to handle floating point numbers"
245	<< " of size other than 4 or 8."
246	<< " You requested: " << size << std::endl;
247	}
248	s_sizeof_real = size;
249	}
250
251	inline void HEPEVT_Wrapper::set_max_number_entries( unsigned int size ) {
252	s_max_number_entries = size;
253	}
254
255	inline double HEPEVT_Wrapper::byte_num_to_double( unsigned int b ) {
256	if ( b >= hepevt_bytes_allocation ) std::cerr
257	<< "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
258	<< std::endl;
259	if ( s_sizeof_real == sizeof(float) ) {
260	float* myfloat = (float*)&hepevt.data[b];
261	return (double)(*myfloat);
262	} else if ( s_sizeof_real == sizeof(double) ) {
263	double* mydouble = (double*)&hepevt.data[b];
264	return (*mydouble);
265	} else {
266	std::cerr
267	<< "HEPEVT_Wrapper: illegal floating point number length."
268	<< s_sizeof_real << std::endl;
269	}
270	return 0;
271	}
272
273	inline int HEPEVT_Wrapper::byte_num_to_int( unsigned int b ) {
274	if ( b >= hepevt_bytes_allocation ) std::cerr
275	<< "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
276	<< std::endl;
277	if ( s_sizeof_int == sizeof(short int) ) {
278	short int* myshortint = (short int*)&hepevt.data[b];
279	return (int)(*myshortint);
280	} else if ( s_sizeof_int == sizeof(long int) ) {
281	long int* mylongint = (long int*)&hepevt.data[b];
282	return (*mylongint);
283	// on some 64 bit machines, int, short, and long are all different
284	} else if ( s_sizeof_int == sizeof(int) ) {
285	int* myint = (int*)&hepevt.data[b];
286	return (*myint);
287	} else {
288	std::cerr
289	<< "HEPEVT_Wrapper: illegal integer number length."
290	<< s_sizeof_int << std::endl;
291	}
292	return 0;
293	}
294
295	inline void HEPEVT_Wrapper::write_byte_num( double in, unsigned int b ) {
296	if ( b >= hepevt_bytes_allocation ) std::cerr
297	<< "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
298	<< std::endl;
299	if ( s_sizeof_real == sizeof(float) ) {
300	float* myfloat = (float*)&hepevt.data[b];
301	(*myfloat) = (float)in;
302	} else if ( s_sizeof_real == sizeof(double) ) {
303	double* mydouble = (double*)&hepevt.data[b];
304	(*mydouble) = (double)in;
305	} else {
306	std::cerr
307	<< "HEPEVT_Wrapper: illegal floating point number length."
308	<< s_sizeof_real << std::endl;
309	}
310	}
311
312	inline void HEPEVT_Wrapper::write_byte_num( int in, unsigned int b ) {
313	if ( b >= hepevt_bytes_allocation ) std::cerr
314	<< "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
315	<< std::endl;
316	if ( s_sizeof_int == sizeof(short int) ) {
317	short int* myshortint = (short int*)&hepevt.data[b];
318	(*myshortint) = (short int)in;
319	} else if ( s_sizeof_int == sizeof(long int) ) {
320	long int* mylongint = (long int*)&hepevt.data[b];
321	(*mylongint) = (int)in;
322	// on some 64 bit machines, int, short, and long are all different
323	} else if ( s_sizeof_int == sizeof(int) ) {
324	int* myint = (int*)&hepevt.data[b];
325	(*myint) = (int)in;
326	} else {
327	std::cerr
328	<< "HEPEVT_Wrapper: illegal integer number length."
329	<< s_sizeof_int << std::endl;
330	}
331	}
332
333	//////////////
334	// INLINES //
335	//////////////
336
337	inline bool HEPEVT_Wrapper::is_double_precision()
338	{
339	// true if 8byte floating point numbers are used in the HepEVT common.
340	return ( sizeof(double) == sizeof_real() );
341	}
342
343	inline int HEPEVT_Wrapper::event_number()
344	{ return byte_num_to_int(0); }
345
346	inline int HEPEVT_Wrapper::number_entries()
347	{
348	int nhep = byte_num_to_int( 1*sizeof_int() );
349	return ( nhep <= max_number_entries() ?
350	nhep : max_number_entries() );
351	}
352
353	inline int HEPEVT_Wrapper::status( int index )
354	{ return byte_num_to_int( (2+index-1) * sizeof_int() ); }
355
356	inline int HEPEVT_Wrapper::id( int index )
357	{
358	return byte_num_to_int( (2+max_number_entries()+index-1)
359	* sizeof_int() );
360	}
361
362	inline int HEPEVT_Wrapper::first_parent( int index )
363	{
364	int parent = byte_num_to_int( (2+2max_number_entries()+2(index-1))
365	* sizeof_int() );
366	return ( parent > 0 && parent <= number_entries() ) ?
367	parent : 0;
368	}
369
370	inline int HEPEVT_Wrapper::last_parent( int index )
371	{
372	// Returns the Index of the LAST parent in the HEPEVT record
373	// for particle with Index index.
374	// If there is only one parent, the last parent is forced to
375	// be the same as the first parent.
376	// If there are no parents for this particle, both the first_parent
377	// and the last_parent with return 0.
378	// Error checking is done to ensure the parent is always
379	// within range ( 0 <= parent <= nhep )
380	//
381	int firstparent = first_parent(index);
382	int parent = byte_num_to_int( (2+2max_number_entries()+2(index-1)+1)
383	* sizeof_int() );
384	return ( parent > firstparent && parent <= number_entries() )
385	? parent : firstparent;
386	}
387
388	inline int HEPEVT_Wrapper::number_parents( int index ) {
389	int firstparent = first_parent(index);
390	return ( firstparent>0 ) ?
391	( 1+last_parent(index)-firstparent ) : 0;
392	}
393
394	inline int HEPEVT_Wrapper::first_child( int index )
395	{
396	int child = byte_num_to_int( (2+4max_number_entries()+2(index-1))
397	* sizeof_int() );
398	return ( child > 0 && child <= number_entries() ) ?
399	child : 0;
400	}
401
402	inline int HEPEVT_Wrapper::last_child( int index )
403	{
404	// Returns the Index of the LAST child in the HEPEVT record
405	// for particle with Index index.
406	// If there is only one child, the last child is forced to
407	// be the same as the first child.
408	// If there are no children for this particle, both the first_child
409	// and the last_child with return 0.
410	// Error checking is done to ensure the child is always
411	// within range ( 0 <= parent <= nhep )
412	//
413	int firstchild = first_child(index);
414	int child = byte_num_to_int( (2+4max_number_entries()+2(index-1)+1)
415	* sizeof_int() );
416	return ( child > firstchild && child <= number_entries() )
417	? child : firstchild;
418	}
419
420	inline int HEPEVT_Wrapper::number_children( int index )
421	{
422	int firstchild = first_child(index);
423	return ( firstchild>0 ) ?
424	( 1+last_child(index)-firstchild ) : 0;
425	}
426
427	inline double HEPEVT_Wrapper::px( int index )
428	{
429	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
430	+ (5(index-1)+0) sizeof_real() );
431	}
432
433	inline double HEPEVT_Wrapper::py( int index )
434	{
435	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
436	+ (5(index-1)+1) sizeof_real() );
437	}
438
439
440	inline double HEPEVT_Wrapper::pz( int index )
441	{
442	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
443	+ (5(index-1)+2) sizeof_real() );
444	}
445
446	inline double HEPEVT_Wrapper::e( int index )
447	{
448	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
449	+ (5(index-1)+3) sizeof_real() );
450	}
451
452	inline double HEPEVT_Wrapper::m( int index )
453	{
454	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
455	+ (5(index-1)+4) sizeof_real() );
456	}
457
458	inline double HEPEVT_Wrapper::x( int index )
459	{
460	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
461	+ ( 5*max_number_entries()
462	+ (4(index-1)+0) ) sizeof_real() );
463	}
464
465	inline double HEPEVT_Wrapper::y( int index )
466	{
467	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
468	+ ( 5*max_number_entries()
469	+ (4(index-1)+1) ) sizeof_real() );
470	}
471
472	inline double HEPEVT_Wrapper::z( int index )
473	{
474	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
475	+ ( 5*max_number_entries()
476	+ (4(index-1)+2) ) sizeof_real() );
477	}
478
479	inline double HEPEVT_Wrapper::t( int index )
480	{
481	return byte_num_to_double( (2+6max_number_entries())sizeof_int()
482	+ ( 5*max_number_entries()
483	+ (4(index-1)+3) ) sizeof_real() );
484	}
485
486	inline void HEPEVT_Wrapper::set_event_number( int evtno )
487	{ write_byte_num( evtno, 0 ); }
488
489	inline void HEPEVT_Wrapper::set_number_entries( int noentries )
490	{ write_byte_num( noentries, 1*sizeof_int() ); }
491
492	inline void HEPEVT_Wrapper::set_status( int index, int status )
493	{
494	if ( index <= 0 \|\| index > max_number_entries() ) return;
495	write_byte_num( status, (2+index-1) * sizeof_int() );
496	}
497
498	inline void HEPEVT_Wrapper::set_id( int index, int id )
499	{
500	if ( index <= 0 \|\| index > max_number_entries() ) return;
501	write_byte_num( id, (2+max_number_entries()+index-1) *sizeof_int() );
502	}
503
504	inline void HEPEVT_Wrapper::set_parents( int index, int firstparent,
505	int lastparent )
506	{
507	if ( index <= 0 \|\| index > max_number_entries() ) return;
508	write_byte_num( firstparent, (2+2max_number_entries()+2(index-1))
509	*sizeof_int() );
510	write_byte_num( lastparent, (2+2max_number_entries()+2(index-1)+1)
511	* sizeof_int() );
512	}
513
514	inline void HEPEVT_Wrapper::set_children( int index, int firstchild,
515	int lastchild )
516	{
517	if ( index <= 0 \|\| index > max_number_entries() ) return;
518	write_byte_num( firstchild, (2+4max_number_entries()+2(index-1))
519	*sizeof_int() );
520	write_byte_num( lastchild, (2+4max_number_entries()+2(index-1)+1)
521	*sizeof_int() );
522	}
523
524	inline void HEPEVT_Wrapper::set_momentum( int index, double px,
525	double py, double pz, double e )
526	{
527	if ( index <= 0 \|\| index > max_number_entries() ) return;
528	write_byte_num( px, (2+6max_number_entries()) sizeof_int()
529	+ (5(index-1)+0) sizeof_real() );
530	write_byte_num( py, (2+6max_number_entries())sizeof_int()
531	+ (5(index-1)+1) sizeof_real() );
532	write_byte_num( pz, (2+6max_number_entries())sizeof_int()
533	+ (5(index-1)+2) sizeof_real() );
534	write_byte_num( e, (2+6max_number_entries())sizeof_int()
535	+ (5(index-1)+3) sizeof_real() );
536	}
537
538	inline void HEPEVT_Wrapper::set_mass( int index, double mass )
539	{
540	if ( index <= 0 \|\| index > max_number_entries() ) return;
541	write_byte_num( mass, (2+6max_number_entries())sizeof_int()
542	+ (5(index-1)+4) sizeof_real() );
543	}
544
545	inline void HEPEVT_Wrapper::set_position( int index, double x, double y,
546	double z, double t )
547	{
548	if ( index <= 0 \|\| index > max_number_entries() ) return;
549	write_byte_num( x, (2+6max_number_entries())sizeof_int()
550	+ ( 5*max_number_entries()
551	+ (4(index-1)+0) ) sizeof_real() );
552	write_byte_num( y, (2+6max_number_entries())sizeof_int()
553	+ ( 5*max_number_entries()
554	+ (4(index-1)+1) ) sizeof_real() );
555	write_byte_num( z, (2+6max_number_entries())sizeof_int()
556	+ ( 5*max_number_entries()
557	+ (4(index-1)+2) ) sizeof_real() );
558	write_byte_num( t, (2+6max_number_entries())sizeof_int()
559	+ ( 5*max_number_entries()
560	+ (4(index-1)+3) ) sizeof_real() );
561	}
562
563	} // HepMC
564
565	#endif // HEPMC_HEPEVT_WRAPPER_H
566	//--------------------------------------------------------------------------
567