DDG4/plugins/Geant4EventReaderHepEvt.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 /** \addtogroup Geant4EventReader
0015  *
0016  @{
0017   \package Geant4EventReaderHepEvt
0018  * \brief Plugin to read HepEvt ASCII files
0019  *
0020  *
0021 @}
0022  */
0023
0024
0025 // Framework include files
0026 #include <DDG4/Geant4InputAction.h>
0027
0028 // C/C++ include files
0029 #include <fstream>
0030
0031 /// Namespace for the AIDA detector description toolkit
0032 namespace dd4hep {
0033
0034   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0035   namespace sim {
0036
0037     /// Class to populate Geant4 primaries from StdHep files.
0038     /**
0039      * Class to populate Geant4 primary particles and vertices from a
0040      * file in HEPEvt format (ASCII)
0041      *
0042      *  \author  P.Kostka (main author)
0043      *  \author  M.Frank  (code reshuffeling into new DDG4 scheme)
0044      *  \version 1.0
0045      *  \ingroup DD4HEP_SIMULATION
0046      */
0047     class Geant4EventReaderHepEvt : public Geant4EventReader  {
0048
0049     protected:
0050       std::ifstream m_input;
0051       int           m_format;
0052
0053     public:
0054       /// Initializing constructor
0055       explicit Geant4EventReaderHepEvt(const std::string& nam, int format);
0056       /// Default destructor
0057       virtual ~Geant4EventReaderHepEvt();
0058       /// Read an event and fill a vector of MCParticles.
0059       virtual EventReaderStatus readParticles(int event_number,
0060                                               Vertices& vertices,
0061                                               std::vector<Particle*>& particles);
0062       virtual EventReaderStatus moveToEvent(int event_number);
0063       virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
0064     };
0065   }     /* End namespace sim   */
0066 }       /* End namespace dd4hep       */
0067
0068 //====================================================================
0069 //  AIDA Detector description implementation
0070 //--------------------------------------------------------------------
0071 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0072 // All rights reserved.
0073 //
0074 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0075 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0076 //
0077 // Author     : M.Frank
0078 //
0079 //====================================================================
0080 // #include "DDG4/Geant4EventReaderHepEvt"
0081
0082 // Framework include files
0083 #include <DDG4/Factories.h>
0084 #include <DD4hep/Printout.h>
0085 #include <CLHEP/Units/SystemOfUnits.h>
0086 #include <CLHEP/Units/PhysicalConstants.h>
0087
0088 // C/C++ include files
0089 #include <cerrno>
0090
0091 using namespace dd4hep::sim;
0092 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0093
0094 #define HEPEvtShort 1
0095 #define HEPEvtLong  2
0096
0097 // Local declarations in anaonymous namespace
0098 namespace {
0099   class Geant4EventReaderHepEvtShort : public Geant4EventReaderHepEvt  {
0100   public:
0101     /// Initializing constructor
0102     explicit Geant4EventReaderHepEvtShort(const std::string& nam) : Geant4EventReaderHepEvt(nam,HEPEvtShort) {}
0103     /// Default destructor
0104     virtual ~Geant4EventReaderHepEvtShort() {}
0105   };
0106   class Geant4EventReaderHepEvtLong : public Geant4EventReaderHepEvt  {
0107   public:
0108     /// Initializing constructor
0109     explicit Geant4EventReaderHepEvtLong(const std::string& nam) : Geant4EventReaderHepEvt(nam,HEPEvtLong) {}
0110     /// Default destructor
0111     virtual ~Geant4EventReaderHepEvtLong() {}
0112   };
0113 }
0114
0115 // Factory entry
0116 DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepEvtShort)
0117 // Factory entry
0118 DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepEvtLong)
0119
0120
0121 /// Initializing constructor
0122 Geant4EventReaderHepEvt::Geant4EventReaderHepEvt(const std::string& nam, int format)
0123 : Geant4EventReader(nam), m_input(), m_format(format)
0124 {
0125   // Now open the input file:
0126   m_input.open(nam.c_str(), std::ifstream::in);
0127   if ( !m_input.good() )   {
0128     except("Geant4EventReaderHepEvt","+++ Failed to open input stream: %s Error:%s",
0129            nam.c_str(), ::strerror(errno));
0130   }
0131 }
0132
0133 /// Default destructor
0134 Geant4EventReaderHepEvt::~Geant4EventReaderHepEvt()    {
0135   m_input.close();
0136 }
0137
0138 /// skipEvents if required
0139 Geant4EventReader::EventReaderStatus
0140 Geant4EventReaderHepEvt::moveToEvent(int event_number) {
0141   if( m_currEvent == 0 && event_number != 0 ) {
0142     printout(INFO,"EventReaderHepEvt::moveToEvent","Skipping the first %d events ", event_number );
0143     printout(INFO,"EventReaderHepEvt::moveToEvent","Event number before skipping: %d", m_currEvent );
0144     while ( m_currEvent < event_number ) {
0145       std::vector<Particle*> particles;
0146       Vertices vertices ;
0147       EventReaderStatus sc = readParticles(m_currEvent,vertices,particles);
0148       for_each(vertices.begin(), vertices.end(), detail::deleteObject<Vertex>);
0149       for_each(particles.begin(), particles.end(), detail::deleteObject<Particle>);
0150       if ( sc != EVENT_READER_OK ) return sc;
0151       //Current event is increased in readParticles already!
0152       // ++m_currEvent;
0153     }
0154   }
0155   printout(INFO,"EventReaderHepEvt::moveToEvent","Event number after skipping: %d", m_currEvent );
0156   return EVENT_READER_OK;
0157 }
0158
0159 /// Read an event and fill a vector of MCParticles.
0160 Geant4EventReader::EventReaderStatus
0161 Geant4EventReaderHepEvt::readParticles(int /* event_number */,
0162                                        Vertices& vertices,
0163                                        std::vector<Particle*>& particles)   {
0164
0165
0166   // First check the input file status
0167   if ( m_input.eof() )   {
0168     return EVENT_READER_EOF;
0169   }
0170   else if ( !m_input.good() )   {
0171     return EVENT_READER_IO_ERROR;
0172   }
0173   //static const double c_light = 299.792;// mm/ns
0174   //
0175   //  Read the event, check for errors
0176   //
0177   unsigned NHEP(0);  // number of entries
0178   m_input >> NHEP;
0179
0180   if( m_input.eof() )  { return EVENT_READER_EOF; }
0181
0182
0183   //check loop variable read from input file and chack that is reasonable
0184   // should fix coverity issue: "Using tainted variable NHEP as a loop boundary."
0185
0186   if( NHEP > 5e7 ){
0187     printout(ERROR,"EventReaderHepEvt::readParticles","Cannot read in too many particles, %d requested but an arbitrary limit has been set to 50 M", NHEP );
0188     return EVENT_READER_EOF;
0189   }
0190
0191   //
0192   //  Loop over particles
0193   int ISTHEP(0);   // status code
0194   int IDHEP(0);    // PDG code
0195   int JMOHEP1(0);  // first mother
0196   int JMOHEP2(0);  // last mother
0197   int JDAHEP1(0);  // first daughter
0198   int JDAHEP2(0);  // last daughter
0199   double PHEP1(0); // px in GeV/c
0200   double PHEP2(0); // py in GeV/c
0201   double PHEP3(0); // pz in GeV/c
0202   double PHEP4(0); // energy in GeV
0203   double PHEP5(0); // mass in GeV/c**2
0204   double VHEP1(0); // x vertex position in mm
0205   double VHEP2(0); // y vertex position in mm
0206   double VHEP3(0); // z vertex position in mm
0207   double VHEP4(0); // production time in mm/c
0208
0209   std::vector<int> daughter1;
0210   std::vector<int> daughter2;
0211
0212   for( unsigned IHEP=0; IHEP<NHEP; IHEP++ )    {
0213     if ( m_format == HEPEvtShort )
0214       m_input >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
0215               >> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
0216     else
0217       m_input >> ISTHEP >> IDHEP
0218               >> JMOHEP1 >> JMOHEP2
0219               >> JDAHEP1 >> JDAHEP2
0220               >> PHEP1 >> PHEP2 >> PHEP3
0221               >> PHEP4 >> PHEP5
0222               >> VHEP1 >> VHEP2 >> VHEP3
0223               >> VHEP4;
0224
0225     if(m_input.eof())
0226       return EVENT_READER_EOF;
0227
0228     if(! m_input.good())
0229       return EVENT_READER_IO_ERROR;
0230
0231     //
0232     //  create a MCParticle and fill it from stdhep info
0233     Particle* p = new Particle(IHEP);
0234     PropertyMask status(p->status);
0235     //
0236     //  PDGID
0237     p->pdgID = IDHEP;
0238     p->charge = 0;
0239     //
0240     //  Momentum vector
0241     p->pex = p->psx = PHEP1*CLHEP::GeV;
0242     p->pey = p->psy = PHEP2*CLHEP::GeV;
0243     p->pez = p->psz = PHEP3*CLHEP::GeV;
0244     //
0245     //  Mass
0246     p->mass = PHEP5*CLHEP::GeV;
0247     //
0248     //  Vertex
0249     p->vsx = VHEP1*CLHEP::mm;
0250     p->vsy = VHEP2*CLHEP::mm;
0251     p->vsz = VHEP3*CLHEP::mm;
0252     // endpoint (missing information in HEPEvt files)
0253     p->vex = 0.0;
0254     p->vey = 0.0;
0255     p->vez = 0.0;
0256     //
0257     //  Creation time (note the units [1/c_light])
0258     p->time       = VHEP4 * CLHEP::mm / CLHEP::c_light;
0259     p->properTime = VHEP4 * CLHEP::mm / CLHEP::c_light;
0260     //
0261     //  Generator status
0262     //  Simulator status 0 until simulator acts on it
0263     p->status = 0;
0264     if ( ISTHEP == 0 )      status.set(G4PARTICLE_GEN_EMPTY);
0265     else if ( ISTHEP == 1 ) status.set(G4PARTICLE_GEN_STABLE);
0266     else if ( ISTHEP == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
0267     else if ( ISTHEP == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0268     else                    status.set(G4PARTICLE_GEN_DOCUMENTATION);
0269     //  RAW Generator status
0270     p->genStatus = ISTHEP&G4PARTICLE_GEN_STATUS_MASK;
0271
0272     //
0273     // Keep daughters information for later
0274     daughter1.emplace_back(JDAHEP1);
0275     daughter2.emplace_back(JDAHEP2);
0276     //
0277     //  Add the particle to the collection vector
0278     particles.emplace_back(p);
0279     //
0280   }// End loop over particles
0281
0282   //
0283   //  Now make a second loop over the particles, checking the daughter
0284   //  information. This is not always consistent with parent
0285   //  information, and this utility assumes all parents listed are
0286   //  parents and all daughters listed are daughters
0287   //
0288   for( unsigned IHEP=0; IHEP<NHEP; IHEP++ )    {
0289     struct ParticleHandler  {
0290       Particle* m_part;
0291       ParticleHandler(Particle* p) : m_part(p) {}
0292       void addParent(const Particle* p)  {
0293         m_part->parents.insert(p->id);
0294       }
0295       void addDaughter(const Particle* p)  {
0296         if(m_part->daughters.find(p->id) == m_part->daughters.end()) {
0297           m_part->daughters.insert(p->id);
0298         }
0299       }
0300       Particle* findParent(const Particle* p)  {
0301         return m_part->parents.find(p->id)==m_part->parents.end() ? 0 : m_part;
0302       }
0303     };
0304
0305     //
0306     //  Get the MCParticle
0307     //
0308     Particle* mcp = particles[IHEP];
0309     ParticleHandler theParticle(mcp);
0310     //
0311     //  Get the daughter information, discarding extra information
0312     //  sometimes stored in daughter variables.
0313     //
0314     int fd = daughter1[IHEP] - 1;
0315     int ld = daughter2[IHEP] - 1;
0316
0317     //
0318     //  As with the parents, look for range, 2 discreet or 1 discreet
0319     //  daughter.
0320     if( (fd > -1) && (fd < int(particles.size())) && (ld > -1) && (ld < int(particles.size())) )  {
0321       if(ld >= fd)   {
0322         for(int id=fd;id<ld+1;id++)   {
0323           //
0324           //  Get the daughter, and see if it already lists this particle as
0325           //  a parent. If not, add this particle as a parent
0326           //
0327           ParticleHandler part(particles[id]);
0328           if ( !part.findParent(mcp) ) part.addParent(mcp);
0329           theParticle.addDaughter(particles[id]);
0330         }
0331       }
0332       else  {   //  Same logic, discrete cases
0333         ParticleHandler part_fd(particles[fd]);
0334         if ( !part_fd.findParent(mcp) ) part_fd.addParent(mcp);
0335         theParticle.addDaughter(particles[fd]);
0336
0337         ParticleHandler part_ld(particles[ld]);
0338         if ( !part_ld.findParent(mcp) ) part_ld.addParent(mcp);
0339         theParticle.addDaughter(particles[ld]);
0340       }
0341     }
0342     else if(fd > -1 && fd < int(particles.size()))  {
0343       ParticleHandler part(particles[fd]);
0344       if ( !part.findParent(mcp) ) part.addParent(mcp);
0345     }
0346     else if(ld > -1 && ld < int(particles.size()))  {
0347       ParticleHandler part(particles[ld]);
0348       if ( !part.findParent(mcp) ) part.addParent(mcp);
0349     }
0350   }  // End second loop over particles
0351
0352
0353   //fg: we simply add all particles without parents as with their own vertex.
0354   //    This might include the incoming beam particles, e.g. in
0355   //    the case of lcio files written with Whizard2, which is slightly odd,
0356   //    however should be treated correctly in Geant4InputHandling.cpp.
0357   //    We no longer make an attempt to identify the incoming particles
0358   //    based on the generator status, as this varies widely with different
0359   //    generators.
0360
0361   for( std::size_t i=0; i < particles.size(); ++i )   {
0362     Geant4ParticleHandle p(particles[i]);
0363     if ( p->parents.size() == 0 )  {
0364       Geant4Vertex* vtx = new Geant4Vertex ;
0365       vertices.emplace_back( vtx );
0366       vtx->x = p->vsx;
0367       vtx->y = p->vsy;
0368       vtx->z = p->vsz;
0369       vtx->time = p->time;
0370
0371       vtx->out.insert(p->id) ;
0372     }
0373   }
0374
0375   ++m_currEvent;
0376   return EVENT_READER_OK;
0377 }
0378