From 077e664ce6a42c08831aa3ff8f47fe8369f56e14 Mon Sep 17 00:00:00 2001
From: Andre Sailer <andre.philippe.sailer@cern.ch>
Date: Tue, 25 Feb 2020 14:24:39 +0100
Subject: [PATCH] HepMC3Reader: initial prototype
---
DDG4/CMakeLists.txt | 13 +++
DDG4/hepmc/HepMC3EventReader.cpp | 135 ++++++++++++++++++++++++++
DDG4/hepmc/HepMC3EventReader.h | 48 +++++++++
DDG4/hepmc/HepMC3FileReader.cpp | 109 +++++++++++++++++++++
DDG4/python/DDSim/DD4hepSimulation.py | 5 +-
5 files changed, 309 insertions(+), 1 deletion(-)
create mode 100644 DDG4/hepmc/HepMC3EventReader.cpp
create mode 100644 DDG4/hepmc/HepMC3EventReader.h
create mode 100644 DDG4/hepmc/HepMC3FileReader.cpp
diff --git a/DDG4/CMakeLists.txt b/DDG4/CMakeLists.txt
index 3ee66b461..5d9c1f51a 100644
--- a/DDG4/CMakeLists.txt
+++ b/DDG4/CMakeLists.txt
@@ -98,6 +98,19 @@ IF(TARGET LCIO::LCIO)
set_target_properties(DDG4LCIO PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION})
ENDIF()
+
+find_package(HepMC3)
+IF(HEPMC3_VERSION)
+ dd4hep_add_plugin(DDG4HepMC3
+ SOURCES hepmc/*.cpp
+ USES DD4hep::DDG4
+ INCLUDES $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/hepmc> ${HEPMC3_INCLUDE_DIR}
+ USES ${HEPMC3_LIBRARIES}
+ )
+ install(TARGETS DDG4HepMC3 EXPORT DD4hep LIBRARY DESTINATION lib)
+ set_target_properties(DDG4HepMC3 PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION})
+ENDIF()
+
#
#--------------------------- DDRec dependent Plugins -----------------------------
#This does not compile at the moment
diff --git a/DDG4/hepmc/HepMC3EventReader.cpp b/DDG4/hepmc/HepMC3EventReader.cpp
new file mode 100644
index 000000000..63f4ceb5b
--- /dev/null
+++ b/DDG4/hepmc/HepMC3EventReader.cpp
@@ -0,0 +1,135 @@
+//==========================================================================
+// 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.
+//
+//
+//====================================================================
+
+// Framework include files
+#include "HepMC3EventReader.h"
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Factories.h"
+
+#include <G4ParticleTable.hh>
+
+#include <HepMC3/FourVector.h>
+#include <HepMC3/GenEvent.h>
+#include <HepMC3/GenParticle.h>
+#include <HepMC3/GenVertex.h>
+#include <HepMC3/Units.h>
+
+#include "G4Event.hh"
+
+using dd4hep::sim::HEPMC3EventReader;
+using SGenParticle = std::shared_ptr<HepMC3::GenParticle>;
+using PropertyMask = dd4hep::detail::ReferenceBitMask<int>;
+
+namespace {
+ template<class T> inline int GET_ENTRY(const std::map<T,int>& mcparts, T part) {
+ auto ip = mcparts.find(part);
+ if (ip == mcparts.end()) {
+ throw std::runtime_error("Unknown particle identifier look-up!");
+ }
+ return (*ip).second;
+ }
+}
+
+/// Initializing constructor
+HEPMC3EventReader::HEPMC3EventReader(const string& fileName): Geant4EventReader(fileName) {}
+
+/// Read an event and fill a vector of MCParticles.
+HEPMC3EventReader::EventReaderStatus
+HEPMC3EventReader::readParticles(int event_number, Vertices& vertices, Particles& particles) {
+ HepMC3::GenEvent genEvent;
+ std::map<SGenParticle, int> mcparts;
+ std::vector<SGenParticle> mcpcoll;
+ EventReaderStatus ret = readGenEvent(event_number, genEvent);
+
+ if ( ret != EVENT_READER_OK ) return ret;
+
+ int NHEP = genEvent.particles().size();
+ // check if there is at least one particle
+ if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
+
+ mcpcoll.resize(NHEP,0);
+ for(int i=0; i<NHEP; ++i ) {
+ auto p = genEvent.particles().at(i);
+ mcparts[p] = i;
+ mcpcoll[i] = p;
+ }
+
+ double mom_unit = (genEvent.momentum_unit() == HepMC3::Units::MomentumUnit::GEV) ? CLHEP::GeV : CLHEP::MeV;
+ double len_unit = (genEvent.length_unit() == HepMC3::Units::LengthUnit::MM) ? CLHEP::mm : CLHEP::cm;
+ // build collection of MCParticles
+ for(int i=0; i<NHEP; ++i ) {
+ auto const& mcp = mcpcoll[i];
+ Geant4ParticleHandle p(new Particle(i));
+ auto const& mom = mcp->momentum();
+ auto vsx = mcp->production_vertex()->position();
+ auto vex = (mcp->end_vertex()) ? mcp->end_vertex()->position() : HepMC3::FourVector();
+
+ const float spin[] = {0.0, 0.0, 0.0}; //mcp->getSpin(); // FIXME
+ const int color[] = {0, 0}; //mcp->getColorFlow(); // FIXME
+ const int pdg = mcp->pid();
+ PropertyMask status(p->status);
+ p->pdgID = pdg;
+ p->charge = 0; // int(mcp->getCharge()*3.0); // FIXME
+ p->psx = mom.get_component(0) * mom_unit;
+ p->psy = mom.get_component(1) * mom_unit;
+ p->psz = mom.get_component(2) * mom_unit;
+ p->time = vsx.get_component(3) * CLHEP::ns; // FIXME
+ p->properTime = vsx.get_component(3) * CLHEP::ns; // FIXME
+ p->vsx = vsx.get_component(0) * len_unit;
+ p->vsy = vsx.get_component(1) * len_unit;
+ p->vsz = vsx.get_component(2) * len_unit;
+ p->vex = vex.get_component(0) * len_unit;
+ p->vey = vex.get_component(1) * len_unit;
+ p->vez = vex.get_component(2) * len_unit;
+ p->process = 0;
+ p->spin[0] = spin[0];
+ p->spin[1] = spin[1];
+ p->spin[2] = spin[2];
+ p->colorFlow[0] = color[0];
+ p->colorFlow[0] = color[1];
+ p->mass = mcp->generated_mass()*mom_unit; // FIXME
+ auto const &par = mcp->parents(), &dau=mcp->children();
+ for(int num=dau.size(), k=0; k<num; ++k)
+ p->daughters.insert(GET_ENTRY(mcparts,dau[k]));
+ for(int num=par.size(), k=0; k<num; ++k)
+ p->parents.insert(GET_ENTRY(mcparts,par[k]));
+
+ int genStatus = mcp->status();
+ if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
+ else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
+ else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
+ else if ( genStatus == 11 ) status.set(G4PARTICLE_GEN_DECAYED); // FIXME
+ else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ else if ( genStatus == 4 ) status.set(G4PARTICLE_GEN_BEAM);
+ else
+ status.set(G4PARTICLE_GEN_OTHER);
+ // Copy raw generator status
+ p->genStatus = genStatus&G4PARTICLE_GEN_STATUS_MASK;
+
+ if ( p->parents.size() == 0 ) {
+
+ Geant4Vertex* vtx = new Geant4Vertex ;
+ vertices.emplace_back( vtx );
+ vtx->x = p->vsx;
+ vtx->y = p->vsy;
+ vtx->z = p->vsz;
+ vtx->time = p->time;
+
+ vtx->out.insert(p->id) ;
+ }
+
+ particles.emplace_back(p);
+ }
+ return EVENT_READER_OK;
+}
diff --git a/DDG4/hepmc/HepMC3EventReader.h b/DDG4/hepmc/HepMC3EventReader.h
new file mode 100644
index 000000000..cb1f04b03
--- /dev/null
+++ b/DDG4/hepmc/HepMC3EventReader.h
@@ -0,0 +1,48 @@
+//==========================================================================
+// 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 DD4HEP_DDG4_HEPMC3EVENTREADER_H
+#define DD4HEP_DDG4_HEPMC3EVENTREADER_H
+
+// Framework include files
+#include "DDG4/Geant4InputAction.h"
+
+namespace HepMC3{ struct GenEvent; }
+
+/// Namespace for the AIDA detector description toolkit
+namespace dd4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace sim {
+
+ /// Base class to read hepmc3 files.
+ /**
+ * \version 1.0
+ * \ingroup DD4HEP_SIMULATION
+ */
+ class HEPMC3EventReader : public Geant4EventReader {
+ public:
+ /// Initializing constructor
+ HEPMC3EventReader(const std::string& nam);
+ /// Default destructor
+ virtual ~HEPMC3EventReader();
+
+ /// Read an event and fill a vector of MCParticles.
+ virtual EventReaderStatus readParticles(int event_number,
+ Vertices& vertices,
+ std::vector<Particle*>& particles);
+ /// Read an event and return a MCParticles.
+ virtual EventReaderStatus readGenEvent(int event_number, HepMC3::GenEvent& genEvent) = 0;
+ };
+
+ } /* End namespace sim */
+} /* End namespace dd4hep */
+#endif /* DD4HEP_DDG4_HEPMC3EVENTREADER_H */
diff --git a/DDG4/hepmc/HepMC3FileReader.cpp b/DDG4/hepmc/HepMC3FileReader.cpp
new file mode 100644
index 000000000..81634ba74
--- /dev/null
+++ b/DDG4/hepmc/HepMC3FileReader.cpp
@@ -0,0 +1,109 @@
+//==========================================================================
+// 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.
+//
+//
+//==========================================================================
+
+/** \addtogroup Geant4EventReader
+ *
+ @{
+ \package HepMC3FileReader
+ * \brief Plugin to read HepMC3 Ascii files
+ *
+ And here we can put a longer description, parameters, examples...
+ *
+@}
+ */
+
+#include "HepMC3EventReader.h"
+
+#include <HepMC3/ReaderFactory.h>
+
+/// Namespace for the AIDA detector description toolkit
+namespace dd4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace sim {
+
+ /// Base class to read hepmc3 event files
+ /**
+ * \version 1.0
+ * \ingroup DD4HEP_SIMULATION
+ */
+ class HEPMC3FileReader : public HEPMC3EventReader {
+ protected:
+ /// Reference to reader object
+ std::shared_ptr<HepMC3::Reader> m_reader;
+ public:
+ /// Initializing constructor
+ HEPMC3FileReader(const std::string& nam);
+ /// Default destructor
+ virtual ~HEPMC3FileReader() = default;
+
+ /// Read an event and fill a vector of MCParticles.
+ virtual EventReaderStatus readGenEvent(int event_number, HepMC3::GenEvent& genEvent);
+ /// skip to event with event_number
+ virtual EventReaderStatus moveToEvent(int event_number);
+ //virtual EventReaderStatus skipEvent() { return EVENT_READER_OK; }
+ virtual EventReaderStatus setParameters(std::map< std::string, std::string >& parameters);
+ };
+ }
+}
+
+#include "DD4hep/Printout.h"
+#include "DDG4/Factories.h"
+
+using dd4hep::sim::HEPMC3FileReader;
+using dd4hep::sim::Geant4EventReader;
+
+// Factory entry
+DECLARE_GEANT4_EVENT_READER_NS(dd4hep::sim,HEPMC3FileReader)
+
+/// Initializing constructor
+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_directAccess = false;
+}
+
+/// moveToSpecifiedEvent, a.k.a. skipNEvents
+Geant4EventReader::EventReaderStatus
+HEPMC3FileReader::moveToEvent(int event_number) {
+ printout(INFO,"HEPMC3FileReader::moveToEvent","Skipping the first %d events ", event_number);
+ while( m_currEvent != event_number) {
+ printout(INFO,"HEPMC3FileReader::moveToEvent","Event number before skipping: %d", m_currEvent );
+ HepMC3::GenEvent genEvent;
+ m_reader->read_event(genEvent) ;
+ ++m_currEvent;
+ printout(INFO,"HEPMC3FileReader::moveToEvent","Event number after skipping: %d", m_currEvent );
+ }
+ return EVENT_READER_OK;
+}
+
+/// Read an event and fill a vector of MCParticles.
+Geant4EventReader::EventReaderStatus
+HEPMC3FileReader::readGenEvent(int /*event_number*/, HepMC3::GenEvent& genEvent) {
+ m_reader->read_event(genEvent);
+ ++m_currEvent;
+ if (genEvent.particles().size()) {
+ printout(INFO,"HEPMC3FileReader","Read event from file");
+ return EVENT_READER_OK;
+ }
+ return EVENT_READER_EOF;
+}
+
+/// Set the parameters for the class
+Geant4EventReader::EventReaderStatus
+HEPMC3FileReader::setParameters( std::map< std::string, std::string > & parameters ) {
+ _getParameterValue(parameters, "Flow1", m_flow1, std::string("flow1"));
+ _getParameterValue(parameters, "Flow2", m_flow2, std::string("flow2"));
+ return EVENT_READER_OK;
+}
diff --git a/DDG4/python/DDSim/DD4hepSimulation.py b/DDG4/python/DDSim/DD4hepSimulation.py
index 0e5f584ae..881d4f043 100644
--- a/DDG4/python/DDSim/DD4hepSimulation.py
+++ b/DDG4/python/DDSim/DD4hepSimulation.py
@@ -34,7 +34,7 @@ try:
except ImportError:
ARGCOMPLETEENABLED = False
-POSSIBLEINPUTFILES = (".stdhep", ".slcio", ".HEPEvt", ".hepevt", ".hepmc", ".pairs")
+POSSIBLEINPUTFILES = (".stdhep", ".slcio", ".HEPEvt", ".hepevt", ".hepmc", ".pairs", ".hepmc3")
def outputLevel(level):
@@ -397,6 +397,9 @@ class DD4hepSimulation(object):
elif inputFile.endswith(".hepevt"):
gen = DDG4.GeneratorAction(kernel, "Geant4InputAction/hepevt%d" % index)
gen.Input = "Geant4EventReaderHepEvtLong|" + inputFile
+ elif inputFile.endswith(".hepmc3"):
+ gen = DDG4.GeneratorAction(kernel, "Geant4InputAction/hepmc%d" % index)
+ gen.Input = "HEPMC3FileReader|" + inputFile
elif inputFile.endswith(".hepmc"):
gen = DDG4.GeneratorAction(kernel, "Geant4InputAction/hepmc%d" % index)
gen.Input = "Geant4EventReaderHepMC|" + inputFile
--
GitLab