DDG4/edm4hep/EDM4hepFileReader.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     : A.Sailer
0011 //
0012 //==========================================================================
0013
0014 /** \addtogroup Geant4EventReader
0015  *
0016  @{
0017  \package EDM4hepFileReader
0018  * \brief Plugin to read EDM4hep files
0019  *
0020  *
0021  @}
0022 */
0023
0024 #include <CLHEP/Units/SystemOfUnits.h>
0025
0026 #include <DDG4/EventParameters.h>
0027 #include <DDG4/Factories.h>
0028 #include <DDG4/FileParameters.h>
0029 #include <DDG4/Geant4InputAction.h>
0030 #include <DDG4/RunParameters.h>
0031
0032 #include <edm4hep/EventHeaderCollection.h>
0033 #include <edm4hep/MCParticleCollection.h>
0034
0035 #include <podio/Frame.h>
0036 #include <podio/Reader.h>
0037
0038 typedef dd4hep::detail::ReferenceBitMask<int> PropertyMask;
0039
0040
0041 namespace dd4hep::sim {
0042
0043   // we use the index of the objectID to identify particles
0044   // we will not support MCParticles from different collections
0045   using MCPARTICLE_MAP=std::map<int, int>;
0046
0047   /// get the parameters from the GenericParameters of the input EDM4hep frame and store them in the EventParameters
0048   /// extension
0049   template <class T=podio::GenericParameters> void EventParameters::ingestParameters(T const& source) {
0050     for(auto const& key: source.template getKeys<int>()) {
0051       m_intValues[key] = source.template get<std::vector<int>>(key).value();
0052     }
0053     for(auto const& key: source.template getKeys<float>()) {
0054       m_fltValues[key] = source.template get<std::vector<float>>(key).value();
0055     }
0056     for(auto const& key: source.template getKeys<double>()) {
0057       m_dblValues[key] = source.template get<std::vector<double>>(key).value();
0058     }
0059     for(auto const& key: source.template getKeys<std::string>()) {
0060       m_strValues[key] = source.template get<std::vector<std::string>>(key).value();
0061     }
0062   }
0063
0064   /// get the parameters from the GenericParameters of the input EDM4hep run frame and store them in the RunParameters
0065   /// extension
0066   template <class T=podio::GenericParameters> void RunParameters::ingestParameters(T const& source) {
0067     for(auto const& key: source.template getKeys<int>()) {
0068       m_intValues[key] = source.template get<std::vector<int>>(key).value();
0069     }
0070     for(auto const& key: source.template getKeys<float>()) {
0071       m_fltValues[key] = source.template get<std::vector<float>>(key).value();
0072     }
0073     for(auto const& key: source.template getKeys<double>()) {
0074       m_dblValues[key] = source.template get<std::vector<double>>(key).value();
0075     }
0076     for(auto const& key: source.template getKeys<std::string>()) {
0077       m_strValues[key] = source.template get<std::vector<std::string>>(key).value();
0078     }
0079   }
0080
0081   template <class T=podio::GenericParameters> void FileParameters::ingestParameters(T const& source) {
0082     for(auto const& key: source.template getKeys<int>()) {
0083       m_intValues[key] = source.template get<std::vector<int>>(key).value();
0084     }
0085     for(auto const& key: source.template getKeys<float>()) {
0086       m_fltValues[key] = source.template get<std::vector<float>>(key).value();
0087     }
0088     for(auto const& key: source.template getKeys<double>()) {
0089       m_dblValues[key] = source.template get<std::vector<double>>(key).value();
0090     }
0091     for(auto const& key: source.template getKeys<std::string>()) {
0092       m_strValues[key] = source.template get<std::vector<std::string>>(key).value();
0093     }
0094   }
0095
0096   /// Class to read EDM4hep files
0097   /**
0098    *  \version 1.0
0099    *  \ingroup DD4HEP_SIMULATION
0100    */
0101   class EDM4hepFileReader : public Geant4EventReader {
0102   protected:
0103     /// Reference to reader object
0104     podio::Reader m_reader;
0105     /// Name of the MCParticle collection to read
0106     std::string m_collectionName;
0107     /// Name of the EventHeader collection to read
0108     std::string m_eventHeaderCollectionName;
0109
0110   public:
0111     /// Initializing constructor
0112     EDM4hepFileReader(const std::string& nam);
0113
0114     /// Read parameters set for this action
0115     virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters);
0116
0117     /// Read an event and fill a vector of MCParticles.
0118     virtual EventReaderStatus readParticles(int event_number, Vertices& vertices, std::vector<Particle*>& particles);
0119     /// Call to register the run parameters
0120     void registerRunParameters();
0121
0122   };
0123
0124   /// Initializing constructor
0125   dd4hep::sim::EDM4hepFileReader::EDM4hepFileReader(const std::string& nam)
0126     : Geant4EventReader(nam)
0127     , m_reader(podio::makeReader(nam))
0128     , m_collectionName("MCParticles")
0129     , m_eventHeaderCollectionName("EventHeader")
0130   {
0131     printout(INFO,"EDM4hepFileReader","Created file reader. Try to open input %s",nam.c_str());
0132     m_directAccess = true;
0133   }
0134
0135   void EDM4hepFileReader::registerRunParameters() {
0136     try {
0137       auto *parameters = new RunParameters();
0138       try {
0139         podio::Frame runFrame = m_reader.readFrame("runs", 0);
0140         parameters->ingestParameters(runFrame.getParameters());
0141       } catch (std::runtime_error& e) {
0142         // we ignore if we do not have runs information
0143       } catch(std::invalid_argument&) {
0144         // we ignore if we do not have runs information
0145       }
0146       context()->run().addExtension<RunParameters>(parameters);
0147     } catch(std::exception &e) {
0148       printout(ERROR,"EDM4hepFileReader::registerRunParameters","Failed to register run parameters: %s", e.what());
0149     }
0150
0151     try {
0152       auto *fileParameters = new FileParameters();
0153       try {
0154         podio::Frame metaFrame = m_reader.readFrame("metadata", 0);
0155         fileParameters->ingestParameters(metaFrame.getParameters());
0156       } catch (std::runtime_error& e) {
0157         // we ignore if we do not have metadata information
0158       } catch(std::invalid_argument&) {
0159         // we ignore if we do not have metadata information
0160       }
0161       context()->run().addExtension<FileParameters>(fileParameters);
0162     } catch(std::exception &e) {
0163       printout(ERROR,"EDM4hepFileReader::registerRunParameters","Failed to register file parameters: %s", e.what());
0164     }
0165   }
0166
0167   namespace {
0168     /// Helper function to look up MCParticles from mapping
0169     inline int GET_ENTRY(MCPARTICLE_MAP const& mcparts, int partID)  {
0170       MCPARTICLE_MAP::const_iterator ip = mcparts.find(partID);
0171       if ( ip == mcparts.end() )  {
0172         throw std::runtime_error("Unknown particle identifier look-up!");
0173       }
0174       return (*ip).second;
0175     }
0176   }
0177
0178   /// Read an event and fill a vector of MCParticles
0179   EDM4hepFileReader::EventReaderStatus
0180   EDM4hepFileReader::readParticles(int event_number, Vertices& vertices, std::vector<Particle*>& particles) {
0181     m_currEvent = event_number;
0182     podio::Frame frame = m_reader.readFrame("events", event_number);
0183     const auto& primaries = frame.get<edm4hep::MCParticleCollection>(m_collectionName);
0184     int eventNumber = event_number, runNumber = 0;
0185     if (primaries.isValid()) {
0186       //Read the event header collection and add it to the context as an extension
0187       const auto& eventHeaderCollection = frame.get<edm4hep::EventHeaderCollection>(m_eventHeaderCollectionName);
0188       if(eventHeaderCollection.isValid() && eventHeaderCollection.size() == 1){
0189         const auto& eh = eventHeaderCollection.at(0);
0190         eventNumber = eh.getEventNumber();
0191         runNumber = eh.getRunNumber();
0192         try {
0193           Geant4Context* ctx = context();
0194           ctx->event().addExtension<edm4hep::MutableEventHeader>(new edm4hep::MutableEventHeader(eh.clone()));
0195         } catch(std::exception& ) {}
0196         // Create input event parameters context
0197         try {
0198           Geant4Context* ctx = context();
0199           EventParameters *parameters = new EventParameters();
0200           parameters->setRunNumber(runNumber);
0201           parameters->setEventNumber(eventNumber);
0202           parameters->ingestParameters(frame.getParameters());
0203           ctx->event().addExtension<EventParameters>(parameters);
0204         } catch(std::exception &) {}
0205       } else {
0206         printout(WARNING,"EDM4hepFileReader","No EventHeader collection found or too many event headers!");
0207         try {
0208           Geant4Context* ctx = context();
0209           EventParameters *parameters = new EventParameters();
0210           parameters->setRunNumber(0);
0211           parameters->setEventNumber(event_number);
0212           parameters->ingestParameters(frame.getParameters());
0213           ctx->event().addExtension<EventParameters>(parameters);
0214         } catch(std::exception &) {}
0215       }
0216       printout(INFO,"EDM4hepFileReader","read collection %s from event %d in run %d ",
0217                m_collectionName.c_str(), eventNumber, runNumber);
0218
0219     } else {
0220       return EVENT_READER_EOF;
0221     }
0222
0223     printout(INFO,"EDM4hepFileReader", "We read the particle collection");
0224     int NHEP = primaries.size();
0225     // check if there is at least one particle
0226     if ( NHEP == 0 ) {
0227       printout(WARNING,"EDM4hepFileReader", "We have no particles");
0228       return EVENT_READER_NO_PRIMARIES;
0229     }
0230
0231     MCPARTICLE_MAP mcparts;
0232     std::vector<int>  mcpcoll;
0233     mcpcoll.resize(NHEP);
0234     for(int i=0; i<NHEP; ++i ) {
0235       edm4hep::MCParticle p = primaries.at(i);
0236       mcparts[p.getObjectID().index] = i;
0237       mcpcoll[i] = p.getObjectID().index;
0238     }
0239
0240     // build collection of MCParticles
0241     for(int i=0; i<NHEP; ++i )   {
0242       const auto& mcp = primaries.at(i);
0243       Geant4ParticleHandle p(new Particle(i));
0244       const auto mom   = mcp.getMomentum();
0245       const auto vsx   = mcp.getVertex();
0246       const auto vex   = mcp.getEndpoint();
0247       const auto spin  = mcp.getSpin();
0248       const int  pdg   = mcp.getPDG();
0249       p->pdgID        = pdg;
0250       p->charge       = int(mcp.getCharge()*3.0);
0251       p->psx          = mom[0]*CLHEP::GeV;
0252       p->psy          = mom[1]*CLHEP::GeV;
0253       p->psz          = mom[2]*CLHEP::GeV;
0254       p->time         = mcp.getTime()*CLHEP::ns;
0255       p->properTime   = mcp.getTime()*CLHEP::ns;
0256       p->vsx          = vsx[0]*CLHEP::mm;
0257       p->vsy          = vsx[1]*CLHEP::mm;
0258       p->vsz          = vsx[2]*CLHEP::mm;
0259       p->vex          = vex[0]*CLHEP::mm;
0260       p->vey          = vex[1]*CLHEP::mm;
0261       p->vez          = vex[2]*CLHEP::mm;
0262       p->process      = 0;
0263       p->spin[0]      = spin[0];
0264       p->spin[1]      = spin[1];
0265       p->spin[2]      = spin[2];
0266       p->colorFlow[0] = 0;
0267       p->colorFlow[1] = 0;
0268       p->mass         = mcp.getMass()*CLHEP::GeV;
0269       const auto par = mcp.getParents(), &dau=mcp.getDaughters();
0270       for(int num=dau.size(),k=0; k<num; ++k) {
0271         edm4hep::MCParticle dau_k = dau[k];
0272         p->daughters.insert(GET_ENTRY(mcparts,dau_k.getObjectID().index));
0273       }
0274       for(int num=par.size(),k=0; k<num; ++k) {
0275         auto par_k = par[k];
0276         p->parents.insert(GET_ENTRY(mcparts, par_k.getObjectID().index));
0277       }
0278
0279       PropertyMask status(p->status);
0280       int genStatus = mcp.getGeneratorStatus();
0281       // Copy raw generator status
0282       p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
0283       if(m_inputAction) {
0284         // in some tests we do not set up the inputAction
0285         m_inputAction->setGeneratorStatus(genStatus, status);
0286       }
0287
0288       //fg: we simply add all particles without parents as with their own vertex.
0289       //    This might include the incoming beam particles, e.g. in
0290       //    the case of lcio files written with Whizard2, which is slightly odd,
0291       //    however should be treated correctly in Geant4InputHandling.cpp.
0292       //    We no longer make an attempt to identify the incoming particles
0293       //    based on the generator status, as this varies widely with different
0294       //    generators.
0295
0296       if ( p->parents.size() == 0 )  {
0297
0298         Geant4Vertex* vtx = new Geant4Vertex ;
0299         vertices.emplace_back( vtx );
0300         vtx->x = p->vsx;
0301         vtx->y = p->vsy;
0302         vtx->z = p->vsz;
0303         vtx->time = p->time;
0304
0305         vtx->out.insert(p->id) ;
0306       }
0307
0308       if ( mcp.isCreatedInSimulation() )       status.set(G4PARTICLE_SIM_CREATED);
0309       if ( mcp.isBackscatter() )               status.set(G4PARTICLE_SIM_BACKSCATTER);
0310       if ( mcp.vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED);
0311       if ( mcp.isDecayedInTracker() )          status.set(G4PARTICLE_SIM_DECAY_TRACKER);
0312       if ( mcp.isDecayedInCalorimeter() )      status.set(G4PARTICLE_SIM_DECAY_CALO);
0313       if ( mcp.hasLeftDetector() )             status.set(G4PARTICLE_SIM_LEFT_DETECTOR);
0314       if ( mcp.isStopped() )                   status.set(G4PARTICLE_SIM_STOPPED);
0315       if ( mcp.isOverlay() )                   status.set(G4PARTICLE_SIM_OVERLAY);
0316       particles.emplace_back(p);
0317     }
0318     return EVENT_READER_OK;
0319   }
0320
0321   /// Set the parameters for the class, in particular the name of the MCParticle
0322   /// list
0323   Geant4EventReader::EventReaderStatus
0324   dd4hep::sim::EDM4hepFileReader::setParameters( std::map< std::string, std::string > & parameters ) {
0325     _getParameterValue(parameters, "MCParticleCollectionName", m_collectionName, m_collectionName);
0326     _getParameterValue(parameters, "EventHeaderCollectionName", m_eventHeaderCollectionName, m_eventHeaderCollectionName);
0327     return EVENT_READER_OK;
0328   }
0329
0330 } //end dd4hep::sim
0331
0332 DECLARE_GEANT4_EVENT_READER_NS(dd4hep::sim,EDM4hepFileReader)