DDG4/plugins/Geant4EventReaderGuineaPig.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 Geant4EventReaderGuineaPig
0018  * \brief Reader for ascii files with e+e- pairs created from GuineaPig.
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 #include <algorithm>
0031 #include <sstream>
0032
0033 /// Namespace for the AIDA detector description toolkit
0034 namespace dd4hep {
0035
0036   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0037   namespace sim {
0038
0039     /// Reader for ascii files with e+e- pairs created from GuineaPig.
0040     /**
0041      *  Reader for ascii files with e+e- pairs created from GuineaPig.
0042      *  Will read complete the file into one event - unless skip N events is
0043      *  called, then N particles are compiled into one event.
0044      *
0045      *  \author  F.Gaede, DESY
0046      *  \author  A. Perez Perez IPHC
0047      *  \version 1.0
0048      *  \ingroup DD4HEP_SIMULATION
0049      */
0050     class Geant4EventReaderGuineaPig : public Geant4EventReader  {
0051
0052     protected:
0053       std::ifstream m_input;
0054       int m_part_num ;
0055
0056     public:
0057       /// Initializing constructor
0058       explicit Geant4EventReaderGuineaPig(const std::string& nam);
0059       /// Default destructor
0060       virtual ~Geant4EventReaderGuineaPig();
0061       /// Read an event and fill a vector of MCParticles.
0062       virtual EventReaderStatus readParticles(int event_number,
0063                                               Vertices& vertices,
0064                                               std::vector<Particle*>& particles) override ;
0065       virtual EventReaderStatus moveToEvent(int event_number) override ;
0066       virtual EventReaderStatus skipEvent() override  { return EVENT_READER_OK; }
0067       virtual EventReaderStatus setParameters( std::map< std::string, std::string > & parameters ) override ;
0068     };
0069   }     /* End namespace sim   */
0070 }       /* End namespace dd4hep       */
0071
0072
0073 // Framework include files
0074 #include <DDG4/Factories.h>
0075 #include <DD4hep/Printout.h>
0076 #include <CLHEP/Units/SystemOfUnits.h>
0077 #include <CLHEP/Units/PhysicalConstants.h>
0078
0079 // C/C++ include files
0080 #include <cerrno>
0081
0082 using namespace dd4hep::sim;
0083 typedef dd4hep::detail::ReferenceBitMask<int> PropertyMask;
0084
0085 // Factory entry
0086 DECLARE_GEANT4_EVENT_READER(Geant4EventReaderGuineaPig)
0087
0088 /// Initializing constructor
0089 Geant4EventReaderGuineaPig::Geant4EventReaderGuineaPig(const std::string& nam)
0090 : Geant4EventReader(nam), m_input(), m_part_num(-1)
0091 {
0092   // Now open the input file:
0093   m_input.open(nam.c_str(), std::ifstream::in);
0094   if ( !m_input.good() )   {
0095     std::string err = "+++ Geant4EventReaderGuineaPig: Failed to open input stream:"+nam+
0096       " Error:"+std::string(strerror(errno));
0097     throw std::runtime_error(err);
0098   }
0099 }
0100
0101 /// Default destructor
0102 Geant4EventReaderGuineaPig::~Geant4EventReaderGuineaPig()    {
0103   m_input.close();
0104 }
0105
0106
0107 Geant4EventReader::EventReaderStatus
0108 Geant4EventReaderGuineaPig::setParameters( std::map< std::string, std::string > & parameters ) {
0109
0110   _getParameterValue( parameters, "ParticlesPerEvent", m_part_num, -1);
0111
0112   if( m_part_num <  0 )
0113     printout(INFO,"EventReader","--- Will read all particles in pairs file into one event " );
0114   else
0115     printout(INFO,"EventReader","--- Will read %d particles per event from pairs file ", m_part_num );
0116
0117   return EVENT_READER_OK;
0118 }
0119
0120 Geant4EventReader::EventReaderStatus
0121 Geant4EventReaderGuineaPig::moveToEvent(int event_number) {
0122
0123   printout(DEBUG,"EventReader"," move to event_number: %d , m_currEvent %d",
0124            event_number,m_currEvent ) ;
0125
0126   if( m_currEvent == 0 && event_number > 0 ){
0127
0128     if( m_part_num <  1 ) {
0129
0130       printout(ERROR,"EventReader","--- Cannot skip to event %d in GuineaPig file without parameter 'ParticlesPerEvent' being set ! ", event_number );
0131
0132       return EVENT_READER_IO_ERROR;
0133
0134     } else {
0135
0136
0137       unsigned nSkipParticles = m_part_num * event_number ;
0138
0139       printout(INFO,"EventReader","--- Will skip first %d events, i.e. %d particles ", event_number , nSkipParticles  );
0140
0141       // First check the input file status
0142       if ( !m_input.good() || m_input.eof() )   {
0143         return EVENT_READER_IO_ERROR;
0144       }
0145
0146       for (unsigned i = 0; i<nSkipParticles; ++i){
0147         if (m_input.ignore(std::numeric_limits<std::streamsize>::max(), m_input.widen('\n'))){
0148           //just skipping the line
0149         }
0150         else
0151           return EVENT_READER_IO_ERROR ;
0152       }
0153     }
0154   }
0155   // else: nothing to do ...
0156
0157   return EVENT_READER_OK;
0158 }
0159
0160 /// Read an event and fill a vector of MCParticles.
0161 Geant4EventReader::EventReaderStatus
0162 Geant4EventReaderGuineaPig::readParticles(int /* event_number */,
0163                                           Vertices& vertices,
0164                                           std::vector<Particle*>& particles)   {
0165
0166
0167   // if no number of particles per event set, we will read the whole file
0168   if ( m_part_num < 0 )
0169     m_part_num = std::numeric_limits<int>::max() ;
0170
0171   // First check the input file status
0172   if ( m_input.eof() )   {
0173     return EVENT_READER_EOF;
0174   }
0175   else if ( !m_input.good() )   {
0176     return EVENT_READER_IO_ERROR;
0177   }
0178
0179   double Energy;
0180   double betaX;
0181   double betaY;
0182   double betaZ;
0183   double posX;
0184   double posY;
0185   double posZ;
0186
0187   //  Loop over particles
0188   for( int counter = 0; counter < m_part_num ; ++counter ){
0189
0190     // need to check for NAN as not all ifstream implementations can handle this directly
0191     std::string lineStr ;
0192     std::getline( m_input, lineStr ) ;
0193
0194     if( m_input.eof() ) {
0195       if( counter==0 ) {
0196         return EVENT_READER_IO_ERROR ;  // reading first particle of event failed
0197       } else{
0198         ++m_currEvent;
0199         return EVENT_READER_OK ; // simply EOF
0200       }
0201     }
0202
0203     std::transform(lineStr.begin(), lineStr.end(), lineStr.begin(), ::tolower);
0204     if( lineStr.find("nan") != std::string::npos){
0205
0206       printout(WARNING,"EventReader","### Read line with 'nan' entries - particle will be ignored  ! " ) ;
0207       continue ;
0208     }
0209     std::stringstream m_input_str( lineStr ) ;
0210
0211     m_input_str  >> Energy
0212                  >> betaX   >> betaY >> betaZ
0213                  >> posX    >> posY  >> posZ ;
0214
0215
0216     //    printf(" ------- %e  %e  %e  %e  %e  %e  %e \n", Energy,betaX, betaY,betaZ,posX,posY,posZ ) ;
0217
0218     //
0219     //  Create a MCParticle and fill it from stdhep info
0220     Particle* p = new Particle(counter);
0221     PropertyMask status(p->status);
0222
0223     //  PDGID: If Energy positive (negative) particle is electron (positron)
0224     //         Remember: Geant4Particle charge is in units of 1/3 elementary/positron charge
0225     p->pdgID  = 11;
0226     p->charge = -1 * 3;
0227     if(Energy < 0.0) {
0228       p->pdgID = -11;
0229       p->charge = +1 * 3;
0230     }
0231
0232     //  Momentum vector
0233     p->pex = p->psx = betaX*TMath::Abs(Energy)*CLHEP::GeV;
0234     p->pey = p->psy = betaY*TMath::Abs(Energy)*CLHEP::GeV;
0235     p->pez = p->psz = betaZ*TMath::Abs(Energy)*CLHEP::GeV;
0236
0237     //  Mass
0238     p->mass = 0.0005109989461*CLHEP::GeV;
0239     //
0240
0241
0242     //  Creation time (note the units [1/c_light])
0243     // ( not information in GuineaPig files )
0244     p->time       = 0.0;
0245     p->properTime = 0.0;
0246
0247
0248     //  Vertex
0249     p->vsx = posX*CLHEP::nm;
0250     p->vsy = posY*CLHEP::nm;
0251     p->vsz = posZ*CLHEP::nm;
0252
0253     Vertex* vtx = new Vertex ;
0254     vtx->x = p->vsx ;
0255     vtx->y = p->vsy ;
0256     vtx->z = p->vsz ;
0257     vtx->time = p->time ;
0258
0259     vtx->out.insert( p->id );
0260
0261     //
0262     //  Generator status
0263     //  Simulator status 0 until simulator acts on it
0264     p->status = 0;
0265     status.set(G4PARTICLE_GEN_STABLE);
0266
0267
0268     //  Add the particle to the collection vector
0269     particles.emplace_back(p);
0270
0271     // create a new vertex for this particle
0272     vertices.emplace_back(vtx);
0273
0274
0275   } // End loop over particles
0276
0277   ++m_currEvent;
0278
0279   return EVENT_READER_OK;
0280
0281 }
0282