include/HepMC3/HEPEVT_Wrapper.h

0001 // -*- C++ -*-
0002 //
0003 // This file is part of HepMC
0004 // Copyright (C) 2014-2023 The HepMC collaboration (see AUTHORS for details)
0005 //
0006 #ifndef HEPMC3_HEPEVT_WRAPPER_H
0007 #define HEPMC3_HEPEVT_WRAPPER_H
0008 /**
0009  *  @file HEPEVT_Wrapper.h
0010  *  @brief Definition of \b class HEPEVT_Wrapper
0011  *
0012  *  @class HepMC3::HEPEVT_Wrapper
0013  *  @brief An interface to HEPEVT common block implemented in a traditional way.
0014  *         When possible this implementation should be avoided and the templated
0015  *         version should be used instead.
0016  *
0017  *  @note This header file does not include HEPEVT definition, only declaration.
0018  *        Including this wrapper requires that HEPEVT is defined somewhere
0019  *        in the project (most likely as FORTRAN common block).
0020  *
0021  *  @note Make sure that HEPEVT definition in project matches this definition
0022  *        (NMXHEP, double precision, etc.) Change this definition if necessary.
0023  *
0024  */
0025
0026 #ifndef HEPMC3_HEPEVT_NMXHEP
0027 /** Default number of particles in the HEPEVT structure */
0028 #define HEPMC3_HEPEVT_NMXHEP 10000
0029 #endif
0030
0031 #ifndef HEPMC3_HEPEVT_PRECISION
0032 /** Default precision of the 4-momentum, time-space position and mass */
0033 #define HEPMC3_HEPEVT_PRECISION double
0034 #endif
0035
0036 /* This definition of HEPEVT corresponds to FORTRAN definition:
0037
0038       PARAMETER (NMXHEP=10000)
0039       COMMON /HEPEVT/  NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
0040      &                 JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
0041       INTEGER          NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
0042       DOUBLE PRECISION PHEP,VHEP
0043 */
0044
0045 static const int NMXHEP = HEPMC3_HEPEVT_NMXHEP; //!< Number of particles in the HEPEVT structure
0046 typedef HEPMC3_HEPEVT_PRECISION momentum_t;     //!< Precision of the 4-momentum, time-space position and mass
0047
0048 /** @struct HEPEVT
0049  *  @brief  C structure representing Fortran common block HEPEVT
0050  * T. Sjöstrand et al., "A proposed standard event record",
0051  *  in `Z physics at LEP 1', eds. G. Altarelli, R. Kleiss and C. Verzegnassi,
0052  * Geneva, Switzerland, September 4-5, 1989, CERN 89-08 (Geneva, 1989), Vol. 3, p. 327
0053  * Disk representation is given by Fortran WRITE/READ format.
0054  */
0055 struct HEPEVT
0056 {
0057     int        nevhep;             //!< Event number
0058     int        nhep;               //!< Number of entries in the event
0059     int        isthep[NMXHEP];     //!< Status code
0060     int        idhep [NMXHEP];     //!< PDG ID
0061     int        jmohep[NMXHEP][2];  //!< Pointer to position of 1st and 2nd (or last!) mother
0062     int        jdahep[NMXHEP][2];  //!< Pointer to position of 1nd and 2nd (or last!) daughter
0063     momentum_t phep  [NMXHEP][5];  //!< Momentum: px, py, pz, e, m
0064     momentum_t vhep  [NMXHEP][4];  //!< Time-space position: x, y, z, t
0065 };                               //!< Fortran common block HEPEVT
0066
0067 #include <iostream>
0068 #include <cstdio>
0069 #include <set>
0070 #include <map>
0071 #include <cstring> // memset
0072 #include <cassert>
0073 #include <algorithm> //min max for VS2017
0074 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
0075 #include "HepMC3/GenEvent.h"
0076 #include "HepMC3/GenParticle.h"
0077 #include "HepMC3/GenVertex.h"
0078 #include "HepMC3/HEPEVT_Helpers.h"
0079 #endif
0080 /** @brief Pointer to HEPEVT common block */
0081 HEPMC3_EXPORT_API  struct HEPEVT*  hepevtptr;
0082
0083 namespace HepMC3
0084 {
0085 class HEPEVT_Wrapper
0086 {
0087
0088 //
0089 // Functions
0090 //
0091 public:
0092
0093     /** @brief Print information from HEPEVT common block */
0094     static void print_hepevt( std::ostream& ostr = std::cout );
0095     /** @brief Print particle information */
0096     static void print_hepevt_particle( int index, std::ostream& ostr = std::cout );
0097     /** @brief Set all entries in HEPEVT to zero */
0098     static void zero_everything();
0099 #ifndef HEPEVT_WRAPPER_HEADER_ONLY
0100     /** @brief Convert GenEvent to HEPEVT*/
0101     static bool GenEvent_to_HEPEVT( const GenEvent* evt ) { return GenEvent_to_HEPEVT_static<HEPEVT_Wrapper>(evt);};
0102     /** @brief Convert HEPEVT to GenEvent*/
0103     static bool HEPEVT_to_GenEvent( GenEvent* evt ) { return HEPEVT_to_GenEvent_static<HEPEVT_Wrapper>(evt);};
0104 #endif
0105     /** @brief Tries to fix list of daughters */
0106     static bool fix_daughters();
0107 //
0108 // Accessors
0109 //
0110 public:
0111     static void   set_max_number_entries( unsigned int size ) { if (size != NMXHEP) printf("This implementation does not support change of the block size.\n"); assert(size == NMXHEP); }//!< Set block size
0112     static void   set_hepevt_address(char *c) { hepevtptr = (struct HEPEVT*)c;          } //!< Set Fortran block address
0113     static int    max_number_entries()      { return NMXHEP;                              } //!< Block size
0114     static int    event_number()            { return hepevtptr->nevhep;             } //!< Get event number
0115     static int    number_entries()          { return hepevtptr->nhep;               } //!< Get number of entries
0116     static int    status(const int& index )       { return hepevtptr->isthep[index-1];    } //!< Get status code
0117     static int    id(const int& index )           { return hepevtptr->idhep[index-1];     } //!< Get PDG particle id
0118     static int    first_parent(const int& index ) { return hepevtptr->jmohep[index-1][0]; } //!< Get index of 1st mother
0119     static int    last_parent(const int& index )  { return hepevtptr->jmohep[index-1][1]; } //!< Get index of last mother
0120     static int    first_child(const int& index )  { return hepevtptr->jdahep[index-1][0]; } //!< Get index of 1st daughter
0121     static int    last_child(const int& index )   { return hepevtptr->jdahep[index-1][1]; } //!< Get index of last daughter
0122     static double px(const int& index )           { return hepevtptr->phep[index-1][0];   } //!< Get X momentum
0123     static double py(const int& index )           { return hepevtptr->phep[index-1][1];   } //!< Get Y momentum
0124     static double pz(const int& index )           { return hepevtptr->phep[index-1][2];   } //!< Get Z momentum
0125     static double e(const int& index )            { return hepevtptr->phep[index-1][3];   } //!< Get Energy
0126     static double m(const int& index )            { return hepevtptr->phep[index-1][4];   } //!< Get generated mass
0127     static double x(const int& index )            { return hepevtptr->vhep[index-1][0];   } //!< Get X Production vertex
0128     static double y(const int& index )            { return hepevtptr->vhep[index-1][1];   } //!< Get Y Production vertex
0129     static double z(const int& index )            { return hepevtptr->vhep[index-1][2];   } //!< Get Z Production vertex
0130     static double t(const int& index )            { return hepevtptr->vhep[index-1][3];   } //!< Get production time
0131     static int    number_parents(const int& index );                                   //!< Get number of parents
0132     static int    number_children(const int& index );                                  //!< Get number of children from the range of daughters
0133     static int    number_children_exact(const int& index );                            //!< Get number of children by counting
0134     static void set_event_number( const int& evtno )       { hepevtptr->nevhep = evtno;         } //!< Set event number
0135     static void set_number_entries( const int& noentries ) { hepevtptr->nhep = noentries;       } //!< Set number of entries
0136     static void set_status( const int& index, const int& status ) { hepevtptr->isthep[index-1] = status; } //!< Set status code
0137     static void set_id(const int& index, const int& id )         { hepevtptr->idhep[index-1] = id;      } //!< Set PDG particle id
0138     static void set_parents( const int& index, const int& firstparent, const int& lastparent );              //!< Set parents
0139     static void set_children( const int& index, const int& firstchild, const int& lastchild );               //!< Set children
0140     static void set_momentum( const int& index, const double& px, const double& py, const double& pz, const double& e );   //!< Set 4-momentum
0141     static void set_mass( const int& index, double mass );                                     //!< Set mass
0142     static void set_position( const int& index, const double& x, const double& y, const double& z, const double& t );      //!< Set position in time-space
0143 };
0144
0145 //
0146 // inline definitions
0147 //
0148 inline void HEPEVT_Wrapper::print_hepevt( std::ostream& ostr )
0149 {
0150     ostr << " Event No.: " << hepevtptr->nevhep << std::endl;
0151     ostr<< "  Nr   Type   Parent(s)  Daughter(s)      Px       Py       Pz       E    Inv. M." << std::endl;
0152     for ( int i = 1; i <= hepevtptr->nhep; ++i )
0153     {
0154         HEPEVT_Wrapper::print_hepevt_particle( i, ostr );
0155     }
0156 }
0157
0158 inline void HEPEVT_Wrapper::print_hepevt_particle( int index, std::ostream& ostr )
0159 {
0160     char buf[255];//Note: the format is fixed, so no reason for complicated treatment
0161
0162     sprintf(buf,"%5i %6i",index,hepevtptr->idhep[index-1]);
0163     ostr << buf;
0164     sprintf(buf,"%4i - %4i  ",hepevtptr->jmohep[index-1][0],hepevtptr->jmohep[index-1][1]);
0165     ostr << buf;
0166     sprintf(buf,"%4i - %4i ",hepevtptr->jdahep[index-1][0],hepevtptr->jdahep[index-1][1]);
0167     ostr << buf;
0168     sprintf(buf,"%8.2f %8.2f %8.2f %8.2f %8.2f",hepevtptr->phep[index-1][0],hepevtptr->phep[index-1][1],hepevtptr->phep[index-1][2],hepevtptr->phep[index-1][3],hepevtptr->phep[index-1][4]);
0169     ostr << buf << std::endl;
0170 }
0171
0172 inline void HEPEVT_Wrapper::zero_everything()
0173 {
0174     memset(hepevtptr,0,sizeof(struct HEPEVT));
0175 }
0176
0177 inline int HEPEVT_Wrapper::number_parents( const int& index )
0178 {
0179     return (hepevtptr->jmohep[index-1][0]) ? (hepevtptr->jmohep[index-1][1]) ? hepevtptr->jmohep[index-1][1]-hepevtptr->jmohep[index-1][0] : 1 : 0;
0180 }
0181
0182 inline int HEPEVT_Wrapper::number_children( const int& index )
0183 {
0184     return (hepevtptr->jdahep[index-1][0]) ? (hepevtptr->jdahep[index-1][1]) ? hepevtptr->jdahep[index-1][1]-hepevtptr->jdahep[index-1][0] : 1 : 0;
0185 }
0186
0187 inline int HEPEVT_Wrapper::number_children_exact( const int& index )
0188 {
0189     int nc=0;
0190     for ( int i = 1; i <= hepevtptr->nhep; ++i )
0191         if (((hepevtptr->jmohep[i-1][0] <= index && hepevtptr->jmohep[i-1][1]>=index))||(hepevtptr->jmohep[i-1][0]==index)||(hepevtptr->jmohep[i-1][1]==index)) nc++;
0192     return nc;
0193 }
0194
0195 inline void HEPEVT_Wrapper::set_parents( const int& index,  const int& firstparent, const int&lastparent )
0196 {
0197     hepevtptr->jmohep[index-1][0] = firstparent;
0198     hepevtptr->jmohep[index-1][1] = lastparent;
0199 }
0200
0201 inline void HEPEVT_Wrapper::set_children(  const int& index,  const int&  firstchild,  const int& lastchild )
0202 {
0203     hepevtptr->jdahep[index-1][0] = firstchild;
0204     hepevtptr->jdahep[index-1][1] = lastchild;
0205 }
0206
0207 inline void HEPEVT_Wrapper::set_momentum( const int& index, const double& px,const double& py, const double& pz, const double& e )
0208 {
0209     hepevtptr->phep[index-1][0] = px;
0210     hepevtptr->phep[index-1][1] = py;
0211     hepevtptr->phep[index-1][2] = pz;
0212     hepevtptr->phep[index-1][3] = e;
0213 }
0214
0215 inline void HEPEVT_Wrapper::set_mass( const int& index, double mass )
0216 {
0217     hepevtptr->phep[index-1][4] = mass;
0218 }
0219
0220 inline void HEPEVT_Wrapper::set_position(  const int& index, const double& x, const double& y, const double& z, const double& t )
0221 {
0222     hepevtptr->vhep[index-1][0] = x;
0223     hepevtptr->vhep[index-1][1] = y;
0224     hepevtptr->vhep[index-1][2] = z;
0225     hepevtptr->vhep[index-1][3] = t;
0226 }
0227
0228 inline bool HEPEVT_Wrapper::fix_daughters()
0229 {
0230     /*AV The function should be called  for a record that has correct particle ordering and mother ids.
0231     As a result it produces a record with ranges where the daughters can be found.
0232     Not every particle in the range will be a daughter. It is true only for proper events.
0233     The return tells if the record was fixed succesfully.
0234     */
0235     for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ )
0236         for ( int k = 1; k <= HEPEVT_Wrapper::number_entries(); k++ ) if (i!=k)
0237                 if ((HEPEVT_Wrapper::first_parent(k) <= i) && (i <= HEPEVT_Wrapper::last_parent(k)))
0238                     HEPEVT_Wrapper::set_children(i,(HEPEVT_Wrapper::first_child(i)==0?k:std::min(HEPEVT_Wrapper::first_child(i),k)),(HEPEVT_Wrapper::last_child(i)==0?k:std::max(HEPEVT_Wrapper::last_child(i),k)));
0239     bool is_fixed = true;
0240     for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); i++ )
0241         is_fixed = (is_fixed && (HEPEVT_Wrapper::number_children_exact(i)==HEPEVT_Wrapper::number_children(i)));
0242     return is_fixed;
0243 }
0244
0245 } // namespace HepMC3
0246 #endif