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