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   declareProperty("AlternativeStableStatuses", m_alternativeStableStatuses = {});
0130   m_needsControl = true;
0131
0132   runAction().callAtBegin(this, &Geant4InputAction::beginRun);
0133 }
0134
0135 /// Default destructor
0136 Geant4InputAction::~Geant4InputAction()   {
0137 }
0138
0139 ///Intialize the event reader before the run starts
0140 void Geant4InputAction::beginRun(const G4Run*) {
0141   createReader();
0142 }
0143
0144 void Geant4InputAction::createReader() {
0145   if(m_reader) {
0146     return;
0147   }
0148   if ( m_input.empty() )  {
0149     except("InputAction: No input file declared!");
0150   }
0151   std::string err;
0152   TypeName tn = TypeName::split(m_input,"|");
0153   try  {
0154     m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
0155     if ( 0 == m_reader )   {
0156       PluginDebug dbg;
0157       m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
0158       abortRun("Error creating reader plugin.",
0159                "Failed to create file reader of type %s. Cannot open dataset %s",
0160                tn.first.c_str(),tn.second.c_str());
0161     }
0162     m_reader->setParameters( m_parameters );
0163     m_reader->checkParameters( m_parameters );
0164     m_reader->setInputAction( this );
0165     m_reader->registerRunParameters();
0166   } catch(const std::exception& e)  {
0167     err = e.what();
0168   }
0169   if ( !err.empty() )  {
0170     abortRun(err,"Error when creating reader for file %s",m_input.c_str());
0171   }
0172 }
0173
0174
0175 /// helper to report Geant4 exceptions
0176 std::string Geant4InputAction::issue(int i)  const  {
0177   std::stringstream str;
0178   str << "Geant4InputAction[" << name() << "]: Event " << i << " ";
0179   return str.str();
0180 }
0181
0182 /// Read an event and return a LCCollection of MCParticles.
0183 int Geant4InputAction::readParticles(int evt_number,
0184                                      Vertices& vertices,
0185                                      std::vector<Particle*>& particles)
0186 {
0187   //in case readParticles is called directly outside of having a run, we make sure a reader exists
0188   createReader();
0189   int evid = evt_number + m_firstEvent;
0190   int status = m_reader->moveToEvent(evid);
0191   if(status == Geant4EventReader::EVENT_READER_EOF ) {
0192     long nEvents = context()->kernel().property("NumEvents").value<long>();
0193     if(nEvents < 0) {
0194       //context()->kernel().runManager().AbortRun(true);
0195       throw DD4hep_End_Of_File();
0196     }
0197   }
0198
0199   if ( Geant4EventReader::EVENT_READER_OK != status )  {
0200     std::string msg = issue(evid)+"Error when moving to event - ";
0201     if ( status == Geant4EventReader::EVENT_READER_EOF ) msg += " EOF: [end of file].";
0202     else msg += " Unknown error condition";
0203     if ( m_abort )  {
0204       abortRun(msg,"Error when reading file %s",m_input.c_str());
0205       return status;
0206     }
0207     error(msg.c_str());
0208     except("Error when reading file %s.", m_input.c_str());
0209     return status;
0210   }
0211   status = m_reader->readParticles(evid, vertices, particles);
0212   if(status == Geant4EventReader::EVENT_READER_EOF ) {
0213     long nEvents = context()->kernel().property("NumEvents").value<long>();
0214     if(nEvents < 0) {
0215       //context()->kernel().runManager().AbortRun(true);
0216       throw DD4hep_End_Of_File();
0217     }
0218   }
0219
0220   if ( Geant4EventReader::EVENT_READER_OK != status )  {
0221     std::string msg = issue(evid)+"Error when moving to event - ";
0222     if ( status == Geant4EventReader::EVENT_READER_EOF ) msg += " EOF: [end of file].";
0223     else msg += " Unknown error condition";
0224     if ( m_abort )  {
0225       abortRun(msg,"Error when reading file %s",m_input.c_str());
0226       return status;
0227     }
0228     error(msg.c_str());
0229     except("Error when reading file %s.", m_input.c_str());
0230   }
0231   return status;
0232 }
0233
0234 /// Callback to generate primary particles
0235 void Geant4InputAction::operator()(G4Event* event)   {
0236   std::vector<Particle*>    primaries;
0237   Geant4Event&              evt = context()->event();
0238   Geant4PrimaryEvent*       prim = evt.extension<Geant4PrimaryEvent>();
0239   Vertices                  vertices ;
0240   int result;
0241
0242   result = readParticles(m_currentEventNumber, vertices, primaries);
0243
0244   event->SetEventID(m_firstEvent + m_currentEventNumber);
0245   ++m_currentEventNumber;
0246
0247   if ( result != Geant4EventReader::EVENT_READER_OK )   {    // handle I/O error, but how?
0248     return;
0249   }
0250
0251   Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
0252   prim->add(m_mask, inter);
0253
0254   // check if there is at least one particle
0255   if ( primaries.empty() ) return;
0256
0257   // check if there is at least one primary vertex
0258   if ( vertices.empty() ) return;
0259
0260   print("+++ Particle interaction with %d generator particles and %d vertices ++++++++++++++++++++++++",
0261         int(primaries.size()), int(vertices.size()) );
0262
0263
0264   for(size_t i=0; i<vertices.size(); ++i )   {
0265     inter->vertices[m_mask].emplace_back( vertices[i] );
0266   }
0267
0268   // build collection of MCParticles
0269   for(auto* primPart : primaries)   {
0270     Geant4ParticleHandle p(primPart);
0271     const double mom_scale = m_momScale;
0272     PropertyMask status(p->status);
0273     p->psx  = mom_scale*p->psx;
0274     p->psy  = mom_scale*p->psy;
0275     p->psz  = mom_scale*p->psz;
0276
0277     //FIXME: this needs to be done now in the readers ...
0278     // // if ( p->parents.size() == 0 )  {
0279     // //   if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
0280     // //     vtx->in.insert(p->id);  // Beam particles and primary quarks etc.
0281     // //   else
0282     // //     vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
0283     // // }
0284
0285     inter->particles.emplace(p->id,p);
0286     p.dumpWithMomentumAndVertex(outputLevel()-1,name(),"->");
0287   }
0288 }
0289
0290 void Geant4InputAction::setGeneratorStatus(int genStatus, PropertyMask& status) {
0291   if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
0292   else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
0293   else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
0294   else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
0295   else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM);
0296   else if ( m_alternativeDecayStatuses.count(genStatus) ) status.set(G4PARTICLE_GEN_DECAYED);
0297   else if ( m_alternativeStableStatuses.count(genStatus) ) status.set(G4PARTICLE_GEN_STABLE);
0298   else
0299     status.set(G4PARTICLE_GEN_OTHER);
0300
0301   return;
0302 }