Propagate run Parameters from input to output

DDG4/edm4hep/Geant4Output2EDM4hep.cpp

@@ -17,6 +17,7 @@
 #include <DD4hep/Detector.h>
 #include <DDG4/EventParameters.h>
 #include <DDG4/Geant4OutputAction.h>
+#include <DDG4/RunParameters.h>
 
 /// edm4hep include files
 #include <edm4hep/MCParticleCollection.h>
@@ -125,6 +126,29 @@ namespace dd4hep {
       }
 #endif
     }
+
+    template <> void RunParameters::extractParameters(podio::Frame& frame)   {
+      for(auto const& p: this->intParameters()) {
+        printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
+        frame.putParameter(p.first, p.second);
+      }
+      for(auto const& p: this->fltParameters()) {
+        printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
+        frame.putParameter(p.first, p.second);
+      }
+      for(auto const& p: this->strParameters()) {
+        printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
+        frame.putParameter(p.first, p.second);
+      }
+#if podio_VERSION_MAJOR > 0 || podio_VERSION_MINOR > 16 || podio_VERSION_PATCH > 2
+      // This functionality is only present in podio > 0.16.2
+      for (auto const& p: this->dblParameters()) {
+        printout(DEBUG, "Geant4OutputEDM4hep", "Saving run parameter: %s", p.first.c_str());
+        frame.putParameter(p.first, p.second);
+      }
+#endif
+    }
+
   }    // End namespace sim
 }      // End namespace dd4hep
 #endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
@@ -278,6 +302,11 @@ void Geant4Output2EDM4hep::saveRun(const G4Run* run)   {
   runHeader.putParameter("DD4hepVersion", versionString());
   runHeader.putParameter("detectorName", context()->detectorDescription().header().name());
 
+  RunParameters* parameters = context()->run().extension<RunParameters>(false);
+  if ( parameters ) {
+    parameters->extractParameters(runHeader);
+  }
+
   m_file->writeFrame(runHeader, "runs");
 }
 

DDG4/hepmc/HepMC3FileReader.cpp

@@ -27,8 +27,10 @@
 #include "HepMC3EventReader.h"
 
 #include "DDG4/EventParameters.h"
+#include "DDG4/RunParameters.h"
 
 #include <HepMC3/ReaderFactory.h>
+#include <HepMC3/Version.h>
 
 /// Namespace for the AIDA detector description toolkit
 namespace dd4hep  {
@@ -68,6 +70,21 @@ namespace dd4hep  {
       }
     }
 
+    template <class T=HepMC3::GenRunInfo> void RunParameters::ingestParameters(T const& runInfo) {
+      // This attributes is not the same return type as for GenEvent!
+      for(auto const& attr: runInfo.attributes()){
+        if(auto int_attr = std::dynamic_pointer_cast<HepMC3::IntAttribute>(attr.second)) {
+          m_intValues[attr.first] = {int_attr->value()};
+        } else if(auto flt_attr = std::dynamic_pointer_cast<HepMC3::FloatAttribute>(attr.second)) {
+          m_fltValues[attr.first] = {flt_attr->value()};
+        } else if(auto dbl_attr = std::dynamic_pointer_cast<HepMC3::DoubleAttribute>(attr.second)) {
+          m_dblValues[attr.first] = {dbl_attr->value()};
+        } else { // anything else
+          m_strValues[attr.first] = {attr.second->unparsed_string()};
+        }
+      }
+    }
+
     /// Base class to read hepmc3 event files
     /**
      *  \version 1.0
@@ -89,6 +106,9 @@ namespace dd4hep  {
       virtual EventReaderStatus moveToEvent(int event_number);
       //virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
       virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters);
+      /// register the run parameters into an extension for the run context
+      virtual void registerRunParameters();
+
     };
   }
 }
@@ -106,11 +126,33 @@ DECLARE_GEANT4_EVENT_READER_NS(dd4hep::sim,HEPMC3FileReader)
 HEPMC3FileReader::HEPMC3FileReader(const std::string& nam)
 : HEPMC3EventReader(nam)
 {
-  m_reader = HepMC3::deduce_reader(nam);
   printout(INFO,"HEPMC3FileReader","Created file reader. Try to open input %s", nam.c_str());
+  m_reader = HepMC3::deduce_reader(nam);
+#if HEPMC3_VERSION_CODE >= 3002006
+  // to get the runInfo in the Ascii reader we have to force HepMC to read the first event
+  m_reader->skip(1);
+  // then we get the run info (shared pointer)
+  auto runInfo = m_reader->run_info();
+  // and close the reader
+  m_reader->close();
+  // so we can open the file again from the start
+  m_reader = HepMC3::deduce_reader(nam);
+  // and set the run info object now
+  m_reader->set_run_info(runInfo);
+#endif
   m_directAccess = false;
 }
 
+void HEPMC3FileReader::registerRunParameters() {
+  try {
+    auto *parameters = new RunParameters();
+    parameters->ingestParameters(*(m_reader->run_info()));
+    context()->run().addExtension<RunParameters>(parameters);
+  } catch(std::exception &e) {
+    printout(ERROR,"HEPMC3FileReader::registerRunParameters","Failed to register run parameters: %s", e.what());
+  }
+}
+
 /// moveToSpecifiedEvent, a.k.a. skipNEvents
 Geant4EventReader::EventReaderStatus
 HEPMC3FileReader::moveToEvent(int event_number) {

DDG4/include/DDG4/Geant4Context.h

@@ -162,7 +162,7 @@ namespace dd4hep {
      *  the user framework.
      *  - The access to the dd4hep objects is via the Geant4Context::detectorDescription() call,
      *  - the access to DDG4 as a whole is supported via Geant4Context::kernel() and
-     *  - the access to the user gframework using a specialized implementation of:
+     *  - the access to the user framework using a specialized implementation of:
      *  template <typename T> T& userFramework()  const;
      *
      *  A user defined implementations must be specialized somewhere in a compilation unit

DDG4/include/DDG4/Geant4InputAction.h

@@ -36,6 +36,7 @@
 
 // Forward declarations
 class G4Event;
+class G4Run;
 
 /// Namespace for the AIDA detector description toolkit
 namespace dd4hep  {
@@ -137,6 +138,9 @@ namespace dd4hep  {
 
       /// make sure that all parameters have been processed, otherwise throw exceptions
       virtual void checkParameters( std::map< std::string, std::string >& );
+
+      /// Register Run Parameters
+      virtual void registerRunParameters() {}
     };
 
     /// Generic input action capable of using the Geant4EventReader class.
@@ -174,6 +178,11 @@ namespace dd4hep  {
       /// Property: named parameters to configure file readers or input actions
       std::map< std::string, std::string> m_parameters;
 
+      /// Perform some actions before the run starts, like opening the event inputs
+      void beginRun(const G4Run*);
+
+      /// Create the input reader
+      void createReader();
     public:
       /// Read an event and return a LCCollectionVec of MCParticles.
       int readParticles(int event_number,

DDG4/include/DDG4/RunParameters.h

@@ -0,0 +1,68 @@
+//==========================================================================
+//  AIDA Detector description implementation
+//--------------------------------------------------------------------------
+// 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.
+//
+//
+//==========================================================================
+#ifndef DDG4_RUNPARAMETERS_H
+#define DDG4_RUNPARAMETERS_H
+
+#include <map>
+#include <string>
+#include <vector>
+
+
+/// Namespace for the AIDA detector description toolkit
+namespace dd4hep  {
+
+  /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+  namespace sim  {
+
+    /// Extension to pass input run data to output run data
+    /**
+     *  \version 1.0
+     *  \ingroup DD4HEP_SIMULATION
+     */
+    class RunParameters  {
+    protected:
+      std::map<std::string, std::vector<int>>         m_intValues {};
+      std::map<std::string, std::vector<float>>       m_fltValues {};
+      std::map<std::string, std::vector<std::string>> m_strValues {};
+      std::map<std::string, std::vector<double>>      m_dblValues {};
+      int                                             m_runNumber = -1;
+
+    public:
+      /// Initializing constructor
+      RunParameters() = default;
+      /// Default destructor
+      ~RunParameters() = default;
+
+      /// Set the Run parameters
+      void setRunNumber(int runNumber);
+      /// Get the run number
+      int runNumber() const { return m_runNumber; }
+
+      /// Copy the parameters from source
+      template <class T> void ingestParameters(T const& source);
+      /// Put parameters into destination
+      template <class T> void extractParameters(T& destination);
+
+      /// Get the int Run parameters
+      auto const& intParameters() const { return m_intValues; }
+      /// Get the float Run parameters
+      auto const& fltParameters() const { return m_fltValues; }
+      /// Get the string Run parameters
+      auto const& strParameters() const { return m_strValues; }
+      /// Get the double Run parameters
+      auto const& dblParameters() const { return m_dblValues; }
+
+    };
+
+  }     /* End namespace sim   */
+}       /* End namespace dd4hep */
+#endif // DDG4_RUNPARAMETERS_H

DDG4/lcio/Geant4Output2LCIO.cpp

@@ -19,6 +19,7 @@
 #include "DDG4/Geant4OutputAction.h"
 
 #include "DDG4/EventParameters.h"
+#include "DDG4/RunParameters.h"
 // Geant4 headers
 #include "G4Threading.hh"
 #include "G4AutoLock.hh"
@@ -62,6 +63,23 @@ namespace dd4hep {
 #endif
     }
 
+    template <class T=lcio::LCRunHeaderImpl> void RunParameters::extractParameters(T& runHeader){
+      auto& lcparameters = runHeader.parameters();
+      for(auto const& ival: this->intParameters()) {
+        lcparameters.setValues(ival.first, ival.second);
+      }
+      for(auto const& ival: this->fltParameters()) {
+        lcparameters.setValues(ival.first, ival.second);
+      }
+      for(auto const& ival: this->strParameters()) {
+        lcparameters.setValues(ival.first, ival.second);
+      }
+#if LCIO_VERSION_GE(2, 17)
+      for(auto const& ival: this->dblParameters()) {
+        lcparameters.setValues(ival.first, ival.second);
+      }
+#endif
+    }
 
     class Geant4ParticleMap;
 
@@ -250,6 +268,10 @@ void Geant4Output2LCIO::saveRun(const G4Run* run)  {
   rh->parameters().setValue("DD4HEPVersion", versionString());
   rh->setRunNumber(m_runNo);
   rh->setDetectorName(context()->detectorDescription().header().name());
+  auto* parameters = context()->run().extension<RunParameters>(false);
+  if (parameters) {
+    parameters->extractParameters(*rh);
+  }
   m_file->writeRunHeader(rh);
 }
 

DDG4/python/DDSim/Helper/OutputConfig.py

@@ -152,6 +152,7 @@ def _configureEDM4HEP(self, dds, geant4):
     logger.info("++++ Setting up EDM4hep ROOT Output ++++")
     e4Out = geant4.setupEDM4hepOutput('EDM4hepOutput', dds.outputFile)
     eventPars = dds.meta.parseEventParameters()
+    e4Out.RunHeader = dds.meta.addParametersToRunHeader(dds)
     e4Out.EventParametersString, e4Out.EventParametersInt, e4Out.EventParametersFloat = eventPars
     e4Out.RunNumberOffset = dds.meta.runNumberOffset if dds.meta.runNumberOffset > 0 else 0
     e4Out.EventNumberOffset = dds.meta.eventNumberOffset if dds.meta.eventNumberOffset > 0 else 0

DDG4/src/Geant4InputAction.cpp

@@ -19,6 +19,7 @@
 #include "DDG4/Geant4Context.h"
 #include "DDG4/Geant4Kernel.h"
 #include "DDG4/Geant4InputAction.h"
+#include "DDG4/Geant4RunAction.h"
 
 #include "G4Event.hh"
 
@@ -127,12 +128,50 @@ Geant4InputAction::Geant4InputAction(Geant4Context* ctxt, const string& nam)
   declareProperty("HaveAbort",      m_abort = true);
   declareProperty("Parameters",     m_parameters = {});
   m_needsControl = true;
+
+  runAction().callAtBegin(this, &Geant4InputAction::beginRun);
 }
 
 /// Default destructor
 Geant4InputAction::~Geant4InputAction()   {
 }
 
+///Intialize the event reader before the run starts
+void Geant4InputAction::beginRun(const G4Run*) {
+  createReader();
+}
+
+void Geant4InputAction::createReader() {
+  if(m_reader) {
+    return;
+  }
+  if ( m_input.empty() )  {
+    except("InputAction: No input file declared!");
+  }
+  string err;
+  TypeName tn = TypeName::split(m_input,"|");
+  try  {
+    m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
+    if ( 0 == m_reader )   {
+      PluginDebug dbg;
+      m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
+      abortRun("Error creating reader plugin.",
+               "Failed to create file reader of type %s. Cannot open dataset %s",
+               tn.first.c_str(),tn.second.c_str());
+    }
+    m_reader->setParameters( m_parameters );
+    m_reader->checkParameters( m_parameters );
+    m_reader->setInputAction( this );
+    m_reader->registerRunParameters();
+  } catch(const exception& e)  {
+    err = e.what();
+  }
+  if ( !err.empty() )  {
+    abortRun(err,"Error when creating reader for file %s",m_input.c_str());
+  }
+}
+
+
 /// helper to report Geant4 exceptions
 string Geant4InputAction::issue(int i)  const  {
   stringstream str;
@@ -145,35 +184,9 @@ int Geant4InputAction::readParticles(int evt_number,
                                      Vertices& vertices,
                                      std::vector<Particle*>& particles)
 {
+  //in case readParticles is called diractly outside of having a run, we make sure a reader exists
+  createReader();
   int evid = evt_number + m_firstEvent;
-  if ( 0 == m_reader )  {
-    if ( m_input.empty() )  {
-      except("InputAction: No input file declared!");
-    }
-    string err;
-    TypeName tn = TypeName::split(m_input,"|");
-    try  {
-      m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
-      if ( 0 == m_reader )   {
-        PluginDebug dbg;
-        m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
-        abortRun(issue(evid)+"Error creating reader plugin.",
-                 "Failed to create file reader of type %s. Cannot open dataset %s",
-                 tn.first.c_str(),tn.second.c_str());
-        return Geant4EventReader::EVENT_READER_NO_FACTORY;
-      }
-      m_reader->setParameters( m_parameters );
-      m_reader->checkParameters( m_parameters );
-      m_reader->setInputAction( this );
-    }
-    catch(const exception& e)  {
-      err = e.what();
-    }
-    if ( !err.empty() )  {
-      abortRun(issue(evid)+err,"Error when creating reader for file %s",m_input.c_str());
-      return Geant4EventReader::EVENT_READER_NO_FACTORY;
-    }
-  }
   int status = m_reader->moveToEvent(evid);
   if(status == Geant4EventReader::EVENT_READER_EOF ) {
     long nEvents = context()->kernel().property("NumEvents").value<long>();