Geant4EventReaderGuineaPig.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : M.Frank
//
//==========================================================================
 
// Framework include files
#include <DDG4/Geant4InputAction.h>
 
// C/C++ include files
#include <fstream>
#include <algorithm>
#include <sstream>
 
namespace dd4hep {
 
  namespace sim {
 
 
    class Geant4EventReaderGuineaPig : public Geant4EventReader  {
 
    protected:
      std::ifstream m_input;
      int m_part_num ;
      
    public:
      explicit Geant4EventReaderGuineaPig(const std::string& nam);
      virtual ~Geant4EventReaderGuineaPig();
      virtual EventReaderStatus readParticles(int event_number,
                                              Vertices& vertices,
                                              std::vector<Particle*>& particles) override ;
      virtual EventReaderStatus moveToEvent(int event_number) override ;
      virtual EventReaderStatus skipEvent() override  { return EVENT_READER_OK; }
      virtual EventReaderStatus setParameters( std::map< std::string, std::string > & parameters ) override ;
    };
  }     /* End namespace sim   */
}       /* End namespace dd4hep       */
 
 
// Framework include files
#include <DDG4/Factories.h>
#include <DD4hep/Printout.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
 
// C/C++ include files
#include <cerrno>
 
using namespace dd4hep::sim;
typedef dd4hep::detail::ReferenceBitMask<int> PropertyMask;
 
// Factory entry
DECLARE_GEANT4_EVENT_READER(Geant4EventReaderGuineaPig)
 
Geant4EventReaderGuineaPig::Geant4EventReaderGuineaPig(const std::string& nam)
: Geant4EventReader(nam), m_input(), m_part_num(-1) 
{
  // Now open the input file:
  m_input.open(nam.c_str(), std::ifstream::in);
  if ( !m_input.good() )   {
    std::string err = "+++ Geant4EventReaderGuineaPig: Failed to open input stream:"+nam+
      " Error:"+std::string(strerror(errno));
    throw std::runtime_error(err);
  }
}
 
Geant4EventReaderGuineaPig::~Geant4EventReaderGuineaPig()    {
  m_input.close();
}
 
 
Geant4EventReader::EventReaderStatus
Geant4EventReaderGuineaPig::setParameters( std::map< std::string, std::string > & parameters ) {
 
  _getParameterValue( parameters, "ParticlesPerEvent", m_part_num, -1);
  
  if( m_part_num <  0 ) 
    printout(INFO,"EventReader","--- Will read all particles in pairs file into one event " );
  else
    printout(INFO,"EventReader","--- Will read %d particles per event from pairs file ", m_part_num );
    
  return EVENT_READER_OK;
}
 
Geant4EventReader::EventReaderStatus
Geant4EventReaderGuineaPig::moveToEvent(int event_number) {
  
  printout(DEBUG,"EventReader"," move to event_number: %d , m_currEvent %d",
           event_number,m_currEvent ) ;
  
  if( m_currEvent == 0 && event_number > 0 ){
 
    if( m_part_num <  1 ) {
 
      printout(ERROR,"EventReader","--- Cannot skip to event %d in GuineaPig file without parameter 'ParticlesPerEvent' being set ! ", event_number );
 
      return EVENT_READER_IO_ERROR;
 
    } else {
 
 
      unsigned nSkipParticles = m_part_num * event_number ;
 
      printout(INFO,"EventReader","--- Will skip first %d events, i.e. %d particles ", event_number , nSkipParticles  );
 
      // First check the input file status
      if ( !m_input.good() || m_input.eof() )   {
        return EVENT_READER_IO_ERROR;
      }
 
      for (unsigned i = 0; i<nSkipParticles; ++i){
        if (m_input.ignore(std::numeric_limits<std::streamsize>::max(), m_input.widen('\n'))){
          //just skipping the line
        }
        else
          return EVENT_READER_IO_ERROR ;
      }
    }
  }
  // else: nothing to do ...
 
  return EVENT_READER_OK;
}
 
Geant4EventReader::EventReaderStatus
Geant4EventReaderGuineaPig::readParticles(int /* event_number */, 
                                          Vertices& vertices,
                                          std::vector<Particle*>& particles)   {
 
 
  // if no number of particles per event set, we will read the whole file
  if ( m_part_num < 0 )
    m_part_num = std::numeric_limits<int>::max() ; 
 
  // First check the input file status
  if ( m_input.eof() )   {
    return EVENT_READER_EOF;
  }
  else if ( !m_input.good() )   {
    return EVENT_READER_IO_ERROR;
  }
 
  double Energy;
  double betaX;
  double betaY;
  double betaZ;
  double posX;
  double posY;
  double posZ;
 
  //  Loop over particles
  for( int counter = 0; counter < m_part_num ; ++counter ){      
 
    // need to check for NAN as not all ifstream implementations can handle this directly
    std::string lineStr ;
    std::getline( m_input, lineStr ) ;
 
    if( m_input.eof() ) {
      if( counter==0 ) { 
        return EVENT_READER_IO_ERROR ;  // reading first particle of event failed 
      } else{
        ++m_currEvent;
        return EVENT_READER_OK ; // simply EOF
      }
    }
 
    std::transform(lineStr.begin(), lineStr.end(), lineStr.begin(), ::tolower);
    if( lineStr.find("nan") != std::string::npos){
 
      printout(WARNING,"EventReader","### Read line with 'nan' entries - particle will be ignored  ! " ) ;
      continue ;
    }
    std::stringstream m_input_str( lineStr ) ;
 
    m_input_str  >> Energy
                 >> betaX   >> betaY >> betaZ
                 >> posX    >> posY  >> posZ ;
 
    
    //    printf(" ------- %e  %e  %e  %e  %e  %e  %e \n", Energy,betaX, betaY,betaZ,posX,posY,posZ ) ;
 
    //
    //  Create a MCParticle and fill it from stdhep info
    Particle* p = new Particle(counter);
    PropertyMask status(p->status);
 
    //  PDGID: If Energy positive (negative) particle is electron (positron)
    //         Remember: Geant4Particle charge is in units of 1/3 elementary/positron charge
    p->pdgID  = 11;
    p->charge = -1 * 3;
    if(Energy < 0.0) {
      p->pdgID = -11;
      p->charge = +1 * 3;
    }
 
    //  Momentum vector
    p->pex = p->psx = betaX*TMath::Abs(Energy)*CLHEP::GeV;
    p->pey = p->psy = betaY*TMath::Abs(Energy)*CLHEP::GeV;
    p->pez = p->psz = betaZ*TMath::Abs(Energy)*CLHEP::GeV;
 
    //  Mass
    p->mass = 0.0005109989461*CLHEP::GeV;
    //
 
 
    //  Creation time (note the units [1/c_light])
    // ( not information in GuineaPig files )
    p->time       = 0.0;
    p->properTime = 0.0;
 
 
    //  Vertex
    p->vsx = posX*CLHEP::nm;
    p->vsy = posY*CLHEP::nm;
    p->vsz = posZ*CLHEP::nm;
 
    Vertex* vtx = new Vertex ;
    vtx->x = p->vsx ;
    vtx->y = p->vsy ;
    vtx->z = p->vsz ;
    vtx->time = p->time ;
 
    vtx->out.insert( p->id ); 
 
    //
    //  Generator status
    //  Simulator status 0 until simulator acts on it
    p->status = 0;
    status.set(G4PARTICLE_GEN_STABLE);
 
 
    //  Add the particle to the collection vector
    particles.emplace_back(p);
 
    // create a new vertex for this particle
    vertices.emplace_back(vtx);
 
 
  } // End loop over particles
 
  ++m_currEvent;
 
  return EVENT_READER_OK;
 
}