DDG4/src/Geant4InputAction.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  P.Kostka (main author)
0011 // @author  M.Frank  (code reshuffeling into new DDG4 scheme)
0012 //
0013 //====================================================================
0014
0015 // Framework include files
0016 #include <DD4hep/Memory.h>
0017 #include <DD4hep/Plugins.h>
0018 #include <DDG4/Geant4Primary.h>
0019 #include <DDG4/Geant4Context.h>
0020 #include <DDG4/Geant4Kernel.h>
0021 #include <DDG4/Geant4InputAction.h>
0022 #include <DDG4/Geant4RunAction.h>
0023
0024 #include <G4Event.hh>
0025
0026 using namespace dd4hep::sim;
0027 using Vertices = Geant4InputAction::Vertices;
0028 using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
0029
0030
0031 /// Initializing constructor
0032 Geant4EventReader::Geant4EventReader(const std::string& nam) : m_name(nam)
0033 {
0034 }
0035
0036 /// Default destructor
0037 Geant4EventReader::~Geant4EventReader()   {
0038 }
0039
0040 /// Get the context (from the input action)
0041 Geant4Context* Geant4EventReader::context() const {
0042   if( 0 == m_inputAction ) {
0043     printout(FATAL,"Geant4EventReader", "No input action registered!");
0044     throw std::runtime_error("Geant4EventReader: No input action registered!");
0045   }
0046   return m_inputAction->context();
0047 }
0048
0049 /// Set the input action
0050 void Geant4EventReader::setInputAction(Geant4InputAction* action) {
0051   m_inputAction = action;
0052 }
0053
0054 /// Skip event. To be implemented for sequential sources
0055 Geant4EventReader::EventReaderStatus Geant4EventReader::skipEvent()  {
0056   if ( hasDirectAccess() )   {
0057     ++m_currEvent;
0058     return EVENT_READER_OK;
0059   }
0060   std::vector<Particle*> particles;
0061   Vertices vertices ;
0062
0063   ++m_currEvent;
0064   EventReaderStatus sc = readParticles(m_currEvent,vertices,particles);
0065   for_each(particles.begin(),particles.end(),detail::deleteObject<Particle>);
0066   for_each(vertices.begin(),vertices.end(),detail::deleteObject<Vertex>);
0067   return sc;
0068 }
0069
0070 /// check if all parameters have been consumed by the reader, otherwise throws exception
0071 void Geant4EventReader::checkParameters(std::map< std::string, std::string > &parameters) {
0072
0073   if( parameters.empty() ) {
0074     return;
0075   }
0076   for (auto const& pairNV : parameters ) {
0077     printout(FATAL,"EventReader::checkParameters","Unknown parameter name: %s with value %s",
0078          pairNV.first.c_str(),
0079          pairNV.second.c_str());
0080   }
0081   throw std::runtime_error("Unknown parameter for event reader");
0082
0083 }
0084
0085 #if 0
0086 Geant4EventReader::EventReaderStatus
0087 Geant4EventReader::moveToEvent(int event_number)   {
0088   if ( event_number >= INT_MIN )   {
0089     return EVENT_READER_OK;  // Logic below does not work as expected.
0090   }                          // This shortcuts it!
0091                              // APS: would have been nice to know what exactly doesn't work...
0092   if ( m_currEvent == event_number )  {
0093     return EVENT_READER_OK;
0094   }
0095   else if ( hasDirectAccess() )   {
0096     m_currEvent = event_number;
0097     return EVENT_READER_OK;
0098   }
0099   else if ( event_number<m_currEvent )   {
0100     return EVENT_READER_ERROR;
0101   }
0102   else  {
0103     for(int i=m_currEvent; i<event_number;++i)
0104       skipEvent();
0105     m_currEvent = event_number;
0106     return EVENT_READER_OK;
0107   }
0108   return EVENT_READER_ERROR;
0109 }
0110 #else
0111 /// Move to the indicated event number.
0112 Geant4EventReader::EventReaderStatus
0113 Geant4EventReader::moveToEvent(int /* event_number */)   {
0114   return EVENT_READER_OK;
0115 }
0116 #endif
0117
0118 /// Standard constructor
0119 Geant4InputAction::Geant4InputAction(Geant4Context* ctxt, const std::string& nam)
0120   : Geant4GeneratorAction(ctxt,nam), m_reader(0), m_currentEventNumber(0)
0121 {
0122   declareProperty("Input",          m_input);
0123   declareProperty("Sync",           m_firstEvent=0);
0124   declareProperty("Mask",           m_mask = 0);
0125   declareProperty("MomentumScale",  m_momScale = 1.0);
0126   declareProperty("HaveAbort",      m_abort = true);
0127   declareProperty("Parameters",     m_parameters = {});
0128   declareProperty("AlternativeDecayStatuses", m_alternativeDecayStatuses = {});
0129   m_needsControl = true;
0130
0131   runAction().callAtBegin(this, &Geant4InputAction::beginRun);
0132 }
0133
0134 /// Default destructor
0135 Geant4InputAction::~Geant4InputAction()   {
0136 }
0137
0138 ///Intialize the event reader before the run starts
0139 void Geant4InputAction::beginRun(const G4Run*) {
0140   createReader();
0141 }
0142
0143 void Geant4InputAction::createReader() {
0144   if(m_reader) {
0145     return;
0146   }
0147   if ( m_input.empty() )  {
0148     except("InputAction: No input file declared!");
0149   }
0150   std::string err;
0151   TypeName tn = TypeName::split(m_input,"|");
0152   try  {
0153     m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
0154     if ( 0 == m_reader )   {
0155       PluginDebug dbg;
0156       m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
0157       abortRun("Error creating reader plugin.",
0158                "Failed to create file reader of type %s. Cannot open dataset %s",
0159                tn.first.c_str(),tn.second.c_str());
0160     }
0161     m_reader->setParameters( m_parameters );
0162     m_reader->checkParameters( m_parameters );
0163     m_reader->setInputAction( this );
0164     m_reader->registerRunParameters();
0165   } catch(const std::exception& e)  {
0166     err = e.what();
0167   }
0168   if ( !err.empty() )  {
0169     abortRun(err,"Error when creating reader for file %s",m_input.c_str());
0170   }
0171 }
0172
0173
0174 /// helper to report Geant4 exceptions
0175 std::string Geant4InputAction::issue(int i)  const  {
0176   std::stringstream str;
0177   str << "Geant4InputAction[" << name() << "]: Event " << i << " ";
0178   return str.str();
0179 }
0180
0181 /// Read an event and return a LCCollection of MCParticles.
0182 int Geant4InputAction::readParticles(int evt_number,
0183                                      Vertices& vertices,
0184                                      std::vector<Particle*>& particles)
0185 {
0186   //in case readParticles is called directly outside of having a run, we make sure a reader exists
0187   createReader();
0188   int evid = evt_number + m_firstEvent;
0189   int status = m_reader->moveToEvent(evid);
0190   if(status == Geant4EventReader::EVENT_READER_EOF ) {
0191     long nEvents = context()->kernel().property("NumEvents").value<long>();
0192     if(nEvents < 0) {
0193       //context()->kernel().runManager().AbortRun(true);
0194       throw DD4hep_End_Of_File();
0195     }
0196   }
0197
0198   if ( Geant4EventReader::EVENT_READER_OK != status )  {
0199     std::string msg = issue(evid)+"Error when moving to event - ";
0200     if ( status == Geant4EventReader::EVENT_READER_EOF ) msg += " EOF: [end of file].";
0201     else msg += " Unknown error condition";
0202     if ( m_abort )  {
0203       abortRun(msg,"Error when reading file %s",m_input.c_str());
0204       return status;
0205     }
0206     error(msg.c_str());
0207     except("Error when reading file %s.", m_input.c_str());
0208     return status;
0209   }
0210   status = m_reader->readParticles(evid, vertices, particles);
0211   if(status == Geant4EventReader::EVENT_READER_EOF ) {
0212     long nEvents = context()->kernel().property("NumEvents").value<long>();
0213     if(nEvents < 0) {
0214       //context()->kernel().runManager().AbortRun(true);
0215       throw DD4hep_End_Of_File();
0216     }
0217   }
0218
0219   if ( Geant4EventReader::EVENT_READER_OK != status )  {
0220     std::string msg = issue(evid)+"Error when moving to event - ";
0221     if ( status == Geant4EventReader::EVENT_READER_EOF ) msg += " EOF: [end of file].";
0222     else msg += " Unknown error condition";
0223     if ( m_abort )  {
0224       abortRun(msg,"Error when reading file %s",m_input.c_str());
0225       return status;
0226     }
0227     error(msg.c_str());
0228     except("Error when reading file %s.", m_input.c_str());
0229   }
0230   return status;
0231 }
0232
0233 /// Callback to generate primary particles
0234 void Geant4InputAction::operator()(G4Event* event)   {
0235   std::vector<Particle*>    primaries;
0236   Geant4Event&              evt = context()->event();
0237   Geant4PrimaryEvent*       prim = evt.extension<Geant4PrimaryEvent>();
0238   Vertices                  vertices ;
0239   int result;
0240
0241   result = readParticles(m_currentEventNumber, vertices, primaries);
0242
0243   event->SetEventID(m_firstEvent + m_currentEventNumber);
0244   ++m_currentEventNumber;
0245
0246   if ( result != Geant4EventReader::EVENT_READER_OK )   {    // handle I/O error, but how?
0247     return;
0248   }
0249
0250   Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
0251   prim->add(m_mask, inter);
0252
0253   // check if there is at least one particle
0254   if ( primaries.empty() ) return;
0255
0256   // check if there is at least one primary vertex
0257   if ( vertices.empty() ) return;
0258
0259   print("+++ Particle interaction with %d generator particles and %d vertices ++++++++++++++++++++++++",
0260         int(primaries.size()), int(vertices.size()) );
0261
0262
0263   for(size_t i=0; i<vertices.size(); ++i )   {
0264     inter->vertices[m_mask].emplace_back( vertices[i] );
0265   }
0266
0267   // build collection of MCParticles
0268   for(auto* primPart : primaries)   {
0269     Geant4ParticleHandle p(primPart);
0270     const double mom_scale = m_momScale;
0271     PropertyMask status(p->status);
0272     p->psx  = mom_scale*p->psx;
0273     p->psy  = mom_scale*p->psy;
0274     p->psz  = mom_scale*p->psz;
0275
0276     //FIXME: this needs to be done now in the readers ...
0277     // // if ( p->parents.size() == 0 )  {
0278     // //   if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
0279     // //     vtx->in.insert(p->id);  // Beam particles and primary quarks etc.
0280     // //   else
0281     // //     vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
0282     // // }
0283
0284     inter->particles.emplace(p->id,p);
0285     p.dumpWithMomentumAndVertex(outputLevel()-1,name(),"->");
0286   }
0287 }
0288
0289 void Geant4InputAction::setGeneratorStatus(int genStatus, PropertyMask& status) {
0290   if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
0291   else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
0292   else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
0293   else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0294   else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM);
0295   else if ( m_alternativeDecayStatuses.count(genStatus) ) status.set(G4PARTICLE_GEN_DECAYED);
0296   else
0297     status.set(G4PARTICLE_GEN_OTHER);
0298
0299   return;
0300 }