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