include/HepMC3/HEPEVT_Wrapper_Runtime_Static.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_RUNTIME_STATIC_H
0007 #define HEPMC3_HEPEVT_WRAPPER_RUNTIME_STATIC_H
0008 #include <iostream>
0009 #include <cstdio>
0010 #include <set>
0011 #include <map>
0012 #include <cstring> // memset
0013 #include <algorithm> //min max for VS2017
0014 #include "HepMC3/GenEvent.h"
0015 #include "HepMC3/GenParticle.h"
0016 #include "HepMC3/GenVertex.h"
0017 #include "HepMC3/HEPEVT_Helpers.h"
0018
0019 /**
0020  *  @file HEPEVT_Wrapper_Runtime_Static.h
0021  *  @brief Definition of \b class HEPEVT_Wrapper_Runtime_Static
0022  *
0023  *  @class HepMC3::HEPEVT_Wrapper_Runtime_Static
0024  *  @brief A static interface to HEPEVT common block implemented to deal with varying block size in runtime
0025  */
0026 namespace HepMC3
0027 {
0028
0029
0030 class HEPEVT_Wrapper_Runtime_Static
0031 {
0032
0033 //
0034 // Functions
0035 //
0036 public:
0037     /** @brief Print information from HEPEVT common block */
0038     static void print_hepevt( std::ostream& ostr = std::cout );
0039     /** @brief Print particle information */
0040     static void print_hepevt_particle( int index, std::ostream& ostr = std::cout );
0041     /** @brief Set all entries in HEPEVT to zero */
0042     static void zero_everything();
0043     /** @brief Convert GenEvent to HEPEVT*/
0044     static bool GenEvent_to_HEPEVT( const GenEvent* evt ) { return GenEvent_to_HEPEVT_static<HEPEVT_Wrapper_Runtime_Static>(evt);};
0045     /** @brief Convert HEPEVT to GenEvent*/
0046     static bool HEPEVT_to_GenEvent( GenEvent* evt ) { return HEPEVT_to_GenEvent_static<HEPEVT_Wrapper_Runtime_Static>(evt);};
0047     /** @brief Tries to fix list of daughters */
0048     static bool fix_daughters();
0049 private:
0050     /** @brief  Fortran common block HEPEVT */
0051     HEPMC3_EXPORT_API static std::shared_ptr<struct HEPEVT_Pointers<double> >  m_hepevtptr;
0052     /** @brief Block size */
0053     HEPMC3_EXPORT_API static int m_max_particles;
0054 //
0055 // Accessors
0056 //
0057 public:
0058     static void   set_max_number_entries( unsigned int size ) { m_max_particles = size; }//!< Set block size
0059     static void   set_hepevt_address(char *c); //!< Set Fortran block address
0060     static int    max_number_entries()      { return m_max_particles;                              } //!< Block size
0061     static int    event_number()            { return *(m_hepevtptr->nevhep);             } //!< Get event number
0062     static int    number_entries()          { return *(m_hepevtptr->nhep);               } //!< Get number of entries
0063     static int    status(const int index )       { return m_hepevtptr->isthep[index-1];    } //!< Get status code
0064     static int    id(const int index )           { return m_hepevtptr->idhep[index-1];     } //!< Get PDG particle id
0065     static int    first_parent(const int index ) { return m_hepevtptr->jmohep[2*(index-1)]; } //!< Get index of 1st mother
0066     static int    last_parent(const int index )  { return m_hepevtptr->jmohep[2*(index-1)+1]; } //!< Get index of last mother
0067     static int    first_child(const int index )  { return m_hepevtptr->jdahep[2*(index-1)]; } //!< Get index of 1st daughter
0068     static int    last_child(const int index )   { return m_hepevtptr->jdahep[2*(index-1)+1]; } //!< Get index of last daughter
0069     static double px(const int index )           { return m_hepevtptr->phep[5*(index-1)];   } //!< Get X momentum
0070     static double py(const int index )           { return m_hepevtptr->phep[5*(index-1)+1];   } //!< Get Y momentum
0071     static double pz(const int index )           { return m_hepevtptr->phep[5*(index-1)+2];   } //!< Get Z momentum
0072     static double e(const int index )            { return m_hepevtptr->phep[5*(index-1)+3];   } //!< Get Energy
0073     static double m(const int index )            { return m_hepevtptr->phep[5*(index-1)+4];   } //!< Get generated mass
0074     static double x(const int index )            { return m_hepevtptr->vhep[4*(index-1)];   } //!< Get X Production vertex
0075     static double y(const int index )            { return m_hepevtptr->vhep[4*(index-1)+1];   } //!< Get Y Production vertex
0076     static double z(const int index )            { return m_hepevtptr->vhep[4*(index-1)+2];   } //!< Get Z Production vertex
0077     static double t(const int index )            { return m_hepevtptr->vhep[4*(index-1)+3];   } //!< Get production time
0078     static int    number_parents(const int index );                                   //!< Get number of parents
0079     static int    number_children(const int index );                                  //!< Get number of children from the range of daughters
0080     static int    number_children_exact(const int index );                            //!< Get number of children by counting
0081     static void set_event_number( const int evtno )       { *(m_hepevtptr->nevhep) = evtno;         } //!< Set event number
0082     static void set_number_entries( const int noentries ) { *(m_hepevtptr->nhep) = noentries;       } //!< Set number of entries
0083     static void set_status( const int index, const int status ) { m_hepevtptr->isthep[index-1] = status; } //!< Set status code
0084     static void set_id(const int index, const int id )         { m_hepevtptr->idhep[index-1] = id;      } //!< Set PDG particle id
0085     static void set_parents( const int index, const int firstparent, const int lastparent );              //!< Set parents
0086     static void set_children( const int index, const int firstchild, const int lastchild );               //!< Set children
0087     static void set_momentum( const int index, const double px, const double py, const double pz, const double e );   //!< Set 4-momentum
0088     static void set_mass( const int index, double mass );                                     //!< Set mass
0089     static void set_position( const int index, const double x, const double y, const double z, const double t );      //!< Set position in time-space
0090 };
0091
0092 //
0093 // inline definitions
0094 //
0095
0096 inline void HEPEVT_Wrapper_Runtime_Static::set_hepevt_address(char *c) {
0097     m_hepevtptr = std::make_shared<struct HEPEVT_Pointers<double> >();
0098     char* x = c;
0099     m_hepevtptr->nevhep = (int*)x;
0100     x += sizeof(int);
0101     m_hepevtptr->nhep = (int*)(x);
0102     x += sizeof(int);
0103     m_hepevtptr->isthep = (int*)(x);
0104     x += sizeof(int)*m_max_particles;
0105     m_hepevtptr->idhep = (int*)(x);
0106     x += sizeof(int)*m_max_particles;
0107     m_hepevtptr->jmohep = (int*)(x);
0108     x += sizeof(int)*m_max_particles*2;
0109     m_hepevtptr->jdahep = (int*)(x);
0110     x += sizeof(int)*m_max_particles*2;
0111     m_hepevtptr->phep = (double*)(x);
0112     x += sizeof(double)*m_max_particles*5;
0113     m_hepevtptr->vhep = (double*)(x);
0114 }
0115
0116 inline void HEPEVT_Wrapper_Runtime_Static::print_hepevt( std::ostream& ostr )
0117 {
0118     ostr << " Event No.: " << *(m_hepevtptr->nevhep) << std::endl;
0119     ostr << "  Nr   Type   Parent(s)  Daughter(s)      Px       Py       Pz       E    Inv. M." << std::endl;
0120     for ( int i = 1; i <= *(m_hepevtptr->nhep); ++i )
0121     {
0122         print_hepevt_particle( i, ostr );
0123     }
0124 }
0125
0126 inline void HEPEVT_Wrapper_Runtime_Static::print_hepevt_particle( int index, std::ostream& ostr )
0127 {
0128     char buf[255];//Note: the format is fixed, so no reason for complicated treatment
0129
0130     sprintf(buf, "%5i %6i", index, m_hepevtptr->idhep[index-1]);
0131     ostr << buf;
0132     sprintf(buf, "%4i - %4i  ", m_hepevtptr->jmohep[2*(index-1)], m_hepevtptr->jmohep[2*(index-1)+1]);
0133     ostr << buf;
0134     sprintf(buf, "%4i - %4i ", m_hepevtptr->jdahep[2*(index-1)], m_hepevtptr->jdahep[2*(index-1)+1]);
0135     ostr << buf;
0136     sprintf(buf, "%8.2f %8.2f %8.2f %8.2f %8.2f", m_hepevtptr->phep[5*(index-1)], m_hepevtptr->phep[5*(index-1)+1], m_hepevtptr->phep[5*(index-1)+2],
0137             m_hepevtptr->phep[5*(index-1)+3], m_hepevtptr->phep[5*(index-1)+4]);
0138     ostr << buf << std::endl;
0139 }
0140
0141 inline void HEPEVT_Wrapper_Runtime_Static::zero_everything()
0142 {
0143     *(m_hepevtptr->nevhep) = 0;
0144     *(m_hepevtptr->nhep) = 0;
0145     memset(m_hepevtptr->isthep, 0, sizeof(int)*m_max_particles);
0146     memset(m_hepevtptr->idhep, 0, sizeof(int)*m_max_particles);
0147     memset(m_hepevtptr->jmohep, 0, sizeof(int)*m_max_particles*2);
0148     memset(m_hepevtptr->jdahep, 0, sizeof(int)*m_max_particles*2);
0149     memset(m_hepevtptr->phep, 0, sizeof(double)*m_max_particles*5);
0150     memset(m_hepevtptr->vhep, 0, sizeof(double)*m_max_particles*4);
0151 }
0152
0153 inline int HEPEVT_Wrapper_Runtime_Static::number_parents( const int index )
0154 {
0155     return (m_hepevtptr->jmohep[2*(index-1)]) ? (m_hepevtptr->jmohep[2*(index-1)+1]) ? m_hepevtptr->jmohep[2*(index-1)+1]
0156            -m_hepevtptr->jmohep[2*(index-1)] : 1 : 0;
0157 }
0158
0159 inline int HEPEVT_Wrapper_Runtime_Static::number_children( const int index )
0160 {
0161     return (m_hepevtptr->jdahep[2*(index-1)]) ? (m_hepevtptr->jdahep[2*(index-1)+1]) ? m_hepevtptr->jdahep[2*(index-1)+1]-m_hepevtptr->jdahep[2*(index-1)] : 1 : 0;
0162 }
0163
0164 inline int HEPEVT_Wrapper_Runtime_Static::number_children_exact( const int index )
0165 {
0166     int nc = 0;
0167     for ( int i = 1; i <= *(m_hepevtptr->nhep); ++i )
0168         if (((m_hepevtptr->jmohep[2*(i-1)] <= index && m_hepevtptr->jmohep[2*(i-1)+1] >= index)) || (m_hepevtptr->jmohep[2*(i-1)] == index) ||
0169                 (m_hepevtptr->jmohep[2*(index-1)+1]==index)) nc++;
0170     return nc;
0171 }
0172
0173
0174 inline void HEPEVT_Wrapper_Runtime_Static::set_parents( const int index,  const int firstparent, const int lastparent )
0175 {
0176     m_hepevtptr->jmohep[2*(index-1)] = firstparent;
0177     m_hepevtptr->jmohep[2*(index-1)+1] = lastparent;
0178 }
0179
0180 inline void HEPEVT_Wrapper_Runtime_Static::set_children(  const int index,  const int  firstchild,  const int lastchild )
0181 {
0182     m_hepevtptr->jdahep[2*(index-1)] = firstchild;
0183     m_hepevtptr->jdahep[2*(index-1)+1] = lastchild;
0184 }
0185
0186 inline void HEPEVT_Wrapper_Runtime_Static::set_momentum( const int index, const double px, const double py, const double pz, const double e )
0187 {
0188     m_hepevtptr->phep[5*(index-1)] = px;
0189     m_hepevtptr->phep[5*(index-1)+1] = py;
0190     m_hepevtptr->phep[5*(index-1)+2] = pz;
0191     m_hepevtptr->phep[5*(index-1)+3] = e;
0192 }
0193
0194 inline void HEPEVT_Wrapper_Runtime_Static::set_mass( const int index, double mass )
0195 {
0196     m_hepevtptr->phep[5*(index-1)+4] = mass;
0197 }
0198
0199 inline void HEPEVT_Wrapper_Runtime_Static::set_position( const int index, const double x, const double y, const double z, const double t )
0200 {
0201     m_hepevtptr->vhep[4*(index-1)] = x;
0202     m_hepevtptr->vhep[4*(index-1)+1] = y;
0203     m_hepevtptr->vhep[4*(index-1)+2] = z;
0204     m_hepevtptr->vhep[4*(index-1)+3] = t;
0205 }
0206
0207 inline bool HEPEVT_Wrapper_Runtime_Static::fix_daughters()
0208 {
0209     /*AV The function should be called  for a record that has correct particle ordering and mother ids.
0210     As a result it produces a record with ranges where the daughters can be found.
0211     Not every particle in the range will be a daughter. It is true only for proper events.
0212     The return tells if the record was fixed succesfully.
0213     */
0214     for ( int i = 1; i <= number_entries(); i++ )
0215         for ( int k=1; k <= number_entries(); k++ ) if (i != k)
0216                 if ((first_parent(k) <= i) && (i <= last_parent(k)))
0217                     set_children(i, (first_child(i) == 0 ? k : std::min(first_child(i), k)), (last_child(i) == 0 ? k : std::max(last_child(i), k)));
0218     bool is_fixed = true;
0219     for ( int i = 1; i <= number_entries(); i++ )
0220         is_fixed = (is_fixed && (number_children_exact(i) == number_children(i)));
0221     return is_fixed;
0222 }
0223
0224 } // namespace HepMC3
0225 #endif