diff --git a/sensitive/Geant4EventReaderGuineaPig.cpp b/sensitive/Geant4EventReaderGuineaPig.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9664d2b6142e2946ddb6d440916c532e42146a6d --- /dev/null +++ b/sensitive/Geant4EventReaderGuineaPig.cpp @@ -0,0 +1,213 @@ +// $Id: $ +//========================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// 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> + +/// Namespace for the AIDA detector description toolkit +namespace DD4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace Simulation { + + /// Class to populate Geant4 primaries from StdHep files. + /** + * Class to populate Geant4 primary particles and vertices from a + * file in HEPEvt format (ASCII) + * + * \author P.Kostka (main author) + * \author M.Frank (code reshuffeling into new DDG4 scheme) + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4EventReaderGuineaPig : public Geant4EventReader { + + protected: + std::ifstream m_input; + + public: + /// Initializing constructor + explicit Geant4EventReaderGuineaPig(const std::string& nam); + /// Default destructor + virtual ~Geant4EventReaderGuineaPig(); + /// Read an event and fill a vector of MCParticles. + virtual EventReaderStatus readParticles(int event_number, + Vertex& primary_vertex, + std::vector<Particle*>& particles); + virtual EventReaderStatus moveToEvent(int event_number); + virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; } + }; + } /* End namespace Simulation */ +} /* End namespace DD4hep */ + +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +//==================================================================== +// #include "DDG4/Geant4EventReaderGuineaPig" + +// 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 std; +using namespace CLHEP; +using namespace DD4hep::Simulation; +typedef DD4hep::ReferenceBitMask<int> PropertyMask; + +// Factory entry +DECLARE_GEANT4_EVENT_READER(Geant4EventReaderGuineaPig) + +/// Initializing constructor +Geant4EventReaderGuineaPig::Geant4EventReaderGuineaPig(const string& nam) +: Geant4EventReader(nam), m_input() +{ + // Now open the input file: + m_input.open(nam.c_str(),ifstream::in); + if ( !m_input.good() ) { + string err = "+++ Geant4EventReaderGuineaPig: Failed to open input stream:"+nam+ + " Error:"+string(strerror(errno)); + throw runtime_error(err); + } +} + +/// Default destructor +Geant4EventReaderGuineaPig::~Geant4EventReaderGuineaPig() { + m_input.close(); +} + +/// skipEvents if required +Geant4EventReader::EventReaderStatus +Geant4EventReaderGuineaPig::moveToEvent(int event_number) { + if( m_currEvent == 0 && event_number != 0 ) { + printout(INFO,"EventReaderGuineaPig::moveToEvent","Skipping the first %d events ", event_number ); + printout(INFO,"EventReaderGuineaPig::moveToEvent","Event number before skipping: %d", m_currEvent ); + while ( m_currEvent < event_number ) { + Geant4Vertex vertex; + std::vector<Particle*> particles; + EventReaderStatus sc = readParticles(m_currEvent,vertex,particles); + for_each(particles.begin(),particles.end(),deleteObject<Particle>); + if ( sc != EVENT_READER_OK ) return sc; + //Current event is increased in readParticles already! + // ++m_currEvent; + } + } + printout(INFO,"EventReaderGuineaPig::moveToEvent","Event number after skipping: %d", m_currEvent ); + return EVENT_READER_OK; +} + +/// Read an event and fill a vector of MCParticles. +Geant4EventReader::EventReaderStatus +Geant4EventReaderGuineaPig::readParticles(int /* event_number */, + Vertex& /* primary_vertex */, + vector<Particle*>& particles) { + + // First check the input file status + if ( !m_input.good() || m_input.eof() ) { + return EVENT_READER_IO_ERROR; + } + + //cout << "Hello! I am in the readParticles function of Geant4EventReaderGuineaPig" << endl; + + //static const double c_light = 299.792;// mm/ns + // + // Read the event, check for errors + // + + //if( m_input.eof() ) { + // + // End of File :: ??? Exception ??? + // -> FG: EOF is not an exception as it happens for every file at the end ! + //return EVENT_READER_IO_ERROR; + //} + + double Energy; + double betaX; + double betaY; + double betaZ; + double posX; + double posY; + double posZ; + + // Loop over particles + int counter = 0; + while(m_input >> Energy + >> betaX >> betaY >> betaZ + >> posX >> posY >> posZ) { + + cout << endl; + cout << "Reading line " << counter+1 + << ": (E,betaX,betaY,betaZ,posX,posY,posZ) = (" << Energy << "," << betaX << "," <<betaY << "," << betaZ << "," << posX << "," << posY << "," << posZ << ")" + << endl; + cout << endl; + + //if(m_input.eof()) return EVENT_READER_IO_ERROR; + // + // 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) + p->pdgID = 11; + p->charge = -1; + if(Energy < 0.0) { + p->pdgID = -11; + p->charge = +1; + } + + // Momentum vector + p->pex = p->psx = betaX*TMath::Abs(Energy)*GeV; + p->pey = p->psy = betaY*TMath::Abs(Energy)*GeV; + p->pez = p->psz = betaZ*TMath::Abs(Energy)*GeV; + + // Mass + p->mass = 0.0005109989461*GeV; + // + // Vertex + // (missing information in HEPEvt files) + p->vsx = posX*nm; + p->vsy = posY*nm; + p->vsz = posZ*nm; + // + // Generator status + // Simulator status 0 until simulator acts on it + p->status = 0; + status.set(G4PARTICLE_GEN_STABLE); + + // Creation time (note the units [1/c_light]) + // (No information in HEPEvt files) + p->time = 0.0; + p->properTime = 0.0; + + // Add the particle to the collection vector + particles.push_back(p); + counter++; + } // End loop over particles + + ++m_currEvent; + + return EVENT_READER_OK; + +} +