From 9ce4923ee8407a3c43572148497cf463438fb021 Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Thu, 7 Aug 2014 12:43:35 +0000 Subject: [PATCH] Improved version of MC particle handler, fix calorimeter hit aggregation with multiple contributions --- DDEve/CMakeLists.txt | 10 +- DDEve/DDEve/DDEve.C | 56 ++++---- DDEve/DDEve/{DDG4Support.C => DDG4IO.C} | 71 +++++++--- DDEve/src/EventControl.cpp | 5 + DDG4/CMakeLists.txt | 25 ++-- DDG4/examples/dictionaries.C | 6 +- DDG4/examples/run.C | 6 +- DDG4/include/DDG4/Geant4Data.h | 14 +- DDG4/include/DDG4/Geant4DataDump.h | 74 ++++++++++ DDG4/include/DDG4/Geant4HitCollection.h | 3 +- DDG4/include/DDG4/Geant4MonteCarloTruth.h | 6 +- DDG4/include/DDG4/Geant4Output2ROOT.h | 3 + DDG4/include/DDG4/Geant4ParticleHandler.h | 5 +- DDG4/lcio/LCIOConversions.cpp | 65 ++++----- DDG4/plugins/Geant4EscapeCounter.cpp | 3 +- DDG4/plugins/Geant4SDActions.cpp | 23 +++- DDG4/python/DDG4Dict.C | 7 + DDG4/src/Geant4Data.cpp | 10 +- DDG4/src/Geant4DataDump.cpp | 97 +++++++++++++ DDG4/src/Geant4Output2ROOT.cpp | 54 +++++++- DDG4/src/Geant4ParticleHandler.cpp | 50 ++++--- DDG4/src/Geant4ROOTDump.cpp | 159 ++++++++++++++++++++++ UtilityApps/src/plugin_runner.cpp | 16 ++- UtilityApps/src/run_plugin.h | 17 ++- doc/release.notes | 39 ++---- examples/CLICSiD/eve/DDEve.xml | 1 + examples/ClientTests/scripts/FCC_Hcal.py | 64 --------- 27 files changed, 655 insertions(+), 234 deletions(-) rename DDEve/DDEve/{DDG4Support.C => DDG4IO.C} (67%) create mode 100644 DDG4/include/DDG4/Geant4DataDump.h create mode 100644 DDG4/src/Geant4DataDump.cpp create mode 100644 DDG4/src/Geant4ROOTDump.cpp diff --git a/DDEve/CMakeLists.txt b/DDEve/CMakeLists.txt index 7419daf1f..e624a4aa3 100644 --- a/DDEve/CMakeLists.txt +++ b/DDEve/CMakeLists.txt @@ -20,12 +20,18 @@ include_directories(${CMAKE_SOURCE_DIR}/DDCore/include #---DD4hepEve library -------------------------------------------------------------- file(GLOB headers include/DDEve/*.h) file(GLOB sources src/*.cpp) +#--------------------------- Support for the LCIO data I/O ------------------------ if(DD4HEP_USE_LCIO) find_package(LCIO REQUIRED) include_directories( ${LCIO_INCLUDE_DIRS} ) list(APPEND sources lcio/LCIOEventHandler.cpp) endif() - +#--------------------------- Support for the DDG4 data I/O ------------------------ +if(DD4HEP_WITH_GEANT4) + root_generate_dictionary( G__DDG4IO ${CMAKE_CURRENT_SOURCE_DIR}/DDEve/DDG4IO.C LINKDEF ${CMAKE_SOURCE_DIR}/DDCore/include/ROOT/LinkDef.h) + list(APPEND sources G__DDG4IO.cxx) +endif() +#--------------------------- Add DDEve dictionary for interactive ROOT sessions --- root_generate_dictionary( G__DDEve ${headers} LINKDEF ${CMAKE_SOURCE_DIR}/DDCore/include/ROOT/LinkDef.h) list(APPEND sources G__DDEve.cxx) @@ -36,7 +42,7 @@ SET( CMAKE_CXX_FLAGS "-Wall -Wextra -pedantic -Wno-long-long") SET_TARGET_PROPERTIES( DD4hepEve PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION}) #---DD4hepEve rootmap -------------------------------------------------------------- dd4hep_generate_rootmap(DD4hepEve) - +#---Package installation procedure(s) ---------------------------------------------- install(DIRECTORY include/DDEve DESTINATION include PATTERN ".svn" EXCLUDE ) diff --git a/DDEve/DDEve/DDEve.C b/DDEve/DDEve/DDEve.C index 4e52f4b32..5df2542ed 100644 --- a/DDEve/DDEve/DDEve.C +++ b/DDEve/DDEve/DDEve.C @@ -4,28 +4,36 @@ #include "TSystem.h" void DDEve(const char* xmlConfig=0) { - char text[1024]; - const char* dd4hep = gSystem->Getenv("DD4hepINSTALL"); - if ( 0 == dd4hep ) { - Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); - Error("DDEve","+++ Your DD4hep installation seems incomplete. The environment DD4hepINSTALL is not defined! +++"); - Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); - gSystem->Exit(EINVAL); - } - Info("DDEve","Has to be run in compiled mode to support DDG4 input ... doing this for you now...."); - ::snprintf(text,sizeof(text)," -I%s/include -D__DD4HEP_DDEVE_EXCLUSIVE__ -Wno-shadow -g -O0",dd4hep); - gSystem->AddIncludePath(text); - gSystem->Load("libDD4hepEve"); - ::snprintf(text,sizeof(text),".L %s/examples/DDEve/DDG4Support.C+",dd4hep); - gSystem->ResetErrno(); - Long_t result = gInterpreter->ProcessLine(text); - result = gSystem->Load("DDG4Support_C"); - if ( result != 1 ) { - Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); - Error("DDEve","+++ Result:%ld %s",result,text); - Error("DDEve","+++ errno=%d %s",gSystem->GetErrno(),gSystem->GetError()); - Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); - gSystem->Exit(gSystem->GetErrno()); - } - DD4hep::EveDisplay(xmlConfig); + char text[1024]; + const char* dd4hep = gSystem->Getenv("DD4hepINSTALL"); + if ( 0 == dd4hep ) { + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + Error("DDEve","+++ Your DD4hep installation seems incomplete. The environment DD4hepINSTALL is not defined! +++"); + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + gSystem->Exit(EINVAL); + } + ::snprintf(text,sizeof(text)," -I%s/include -D__DD4HEP_DDEVE_EXCLUSIVE__ -Wno-shadow -g -O0",dd4hep); + gSystem->AddIncludePath(text); + Long_t result = = gSystem->Load("libDD4hepEve"); + if ( 0 != result ) { + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + Error("DDEve","+++ Your DD4hep installation seems incomplete. FAILED to load the library 'libDD4hepEve'! +++"); + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + gSystem->Exit(EINVAL); + } +#if 0 + Info("DDEve","Has to be run in compiled mode to support DDG4 input ... doing this for you now...."); + ::snprintf(text,sizeof(text),".L %s/examples/DDEve/DDG4IO.C+",dd4hep); + gSystem->ResetErrno(); + result = gInterpreter->ProcessLine(text); + result = gSystem->Load("DDG4Support_C"); + if ( result != 1 ) { + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + Error("DDEve","+++ Result:%ld %s",result,text); + Error("DDEve","+++ errno=%d %s",gSystem->GetErrno(),gSystem->GetError()); + Error("DDEve","++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); + gSystem->Exit(gSystem->GetErrno()); + } +#endif + DD4hep::EveDisplay(xmlConfig); } diff --git a/DDEve/DDEve/DDG4Support.C b/DDEve/DDEve/DDG4IO.C similarity index 67% rename from DDEve/DDEve/DDG4Support.C rename to DDEve/DDEve/DDG4IO.C index a9292b6c5..ed1d9f3d8 100644 --- a/DDEve/DDEve/DDG4Support.C +++ b/DDEve/DDEve/DDG4IO.C @@ -8,6 +8,8 @@ //==================================================================== // Framework include files +#define __DD4HEP_DDEVE_EXCLUSIVE__ + // C/C++ include files #include <vector> #include "Math/Vector3D.h" @@ -42,6 +44,16 @@ namespace DD4hep { inline SimpleEvent::SimpleEvent() { } /// Default destructor inline SimpleEvent::~SimpleEvent() { } + + /// Default constructor + inline Particle::Particle() { } + /// Copy constructor + inline Particle::Particle(const Particle&) { NO_CALL } + /// Default destructor + inline Particle::~Particle() { } + /// Remove daughter from set + inline void Particle::removeDaughter(int) { NO_CALL } + /// Default constructor inline SimpleHit::SimpleHit() { } /// Default destructor @@ -74,8 +86,13 @@ namespace DD4hep { //#pragma link C++ class Position+; //#pragma link C++ class Direction+; +//#pragma link C++ class std::set<int>+; #pragma link C++ class DD4hep::Simulation::SimpleRun+; #pragma link C++ class DD4hep::Simulation::SimpleEvent+; + +#pragma link C++ class DD4hep::Simulation::Particle+; +#pragma link C++ class std::vector<DD4hep::Simulation::Particle*>+; + #pragma link C++ class DD4hep::Simulation::SimpleHit+; #pragma link C++ class std::vector<DD4hep::Simulation::SimpleHit*>+; #pragma link C++ class DD4hep::Simulation::SimpleHit::Contribution+; @@ -87,32 +104,52 @@ namespace DD4hep { #pragma link C++ class DD4hep::Simulation::SimpleCalorimeter::Hit+; #pragma link C++ class std::vector<DD4hep::Simulation::SimpleCalorimeter::Hit*>+; #else +#include <typeinfo> +#include "TROOT.h" +#include "TClass.h" using namespace std; using namespace DD4hep; using namespace DD4hep::Simulation; -template <typename T> static T* _fill(SimpleHit* ptr, DDEveHit* target) { - T* s = dynamic_cast<T*>(ptr); - if ( s ) { - Geometry::Position* p = &s->position; - target->x = p->X(); - target->y = p->Y(); - target->z = p->Z(); - target->deposit = s->energyDeposit; + +namespace { + template <typename T> T* _fill(SimpleHit* ptr, DDEveHit* target) { + T* s = dynamic_cast<T*>(ptr); + if ( s ) { + Geometry::Position* p = &s->position; + target->x = p->X(); + target->y = p->Y(); + target->z = p->Z(); + target->deposit = s->energyDeposit; + return s; + } + return 0; } - return s; -} -static void* _convertHitFunc(void* source, DDEveHit* target) { - void* result; - SimpleHit* hit = (SimpleHit*)source; - if ((result=_fill<SimpleTracker::Hit>(hit,target))) return result; - if ((result=_fill<SimpleCalorimeter::Hit>(hit,target))) return result; - return 0; + void* _convertHitFunc(void* source, DDEveHit* target) { + if (source ) { + static TClass* cl_calo = gROOT->GetClass(typeid(SimpleCalorimeter::Hit)); + static TClass* cl_tracker = gROOT->GetClass(typeid(SimpleTracker::Hit)); + //static TClass* cl_particles = gROOT->GetClass(typeid(Particle)); + + void* result = 0; + SimpleHit* hit = (SimpleHit*)source; + const std::type_info& type = typeid(*hit); + TClass* cl = gROOT->GetClass(type); + if ( cl == cl_tracker && (result=_fill<SimpleTracker::Hit>(hit,target)) ) return result; + if ( cl == cl_calo && (result=_fill<SimpleCalorimeter::Hit>(hit,target)) ) return result; + } + return 0; + } + union FCN { + void* v; + void* (*f)(void*, DDEveHit*); + FCN(void* (*ff)(void*, DDEveHit*)) { f=ff; } + }; } static void* _convertHit(const char*) { - return (void*)_convertHitFunc; + return FCN(_convertHitFunc).v; } #include "DD4hep/Factories.h" diff --git a/DDEve/src/EventControl.cpp b/DDEve/src/EventControl.cpp index 5931b712e..ee6fc318b 100644 --- a/DDEve/src/EventControl.cpp +++ b/DDEve/src/EventControl.cpp @@ -117,6 +117,11 @@ void EventControl::OnNewEvent(EventHandler* handler) { cl = cl.substr(idx); cl = cl.substr(0,cl.find('*')); } + else if ( (idx=cl.rfind("::")) != string::npos ) { + cl = cl.substr(idx+2); + if ( (idx=cl.rfind('*')) != string::npos ) cl = cl.substr(0,idx); + if ( (idx=cl.rfind('>')) != string::npos ) cl = cl.substr(0,idx); + } line.second.first->SetTextColor(kRed); line.second.second->SetTextColor(kRed); line.second.first->SetText(("Coll.Type: "+cl).c_str()); diff --git a/DDG4/CMakeLists.txt b/DDG4/CMakeLists.txt index cbe16a848..32bf3da16 100644 --- a/DDG4/CMakeLists.txt +++ b/DDG4/CMakeLists.txt @@ -1,8 +1,20 @@ +#---Find Geant4------------------------------------------------------------------- +find_package(Geant4 REQUIRED gdml ui_all vis_all) +INCLUDE(${Geant4_USE_FILE}) # this also takes care of geant 4 definitions and include dirs #---Includedirs------------------------------------------------------------------- include_directories(${CMAKE_SOURCE_DIR}/DDCore/include ${CMAKE_SOURCE_DIR}/DDSegmentation/include - ${CMAKE_CURRENT_SOURCE_DIR}/include ) + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${ROOT_INCLUDE_DIR} + ${CLHEP_INCLUDE_DIR} + ${Geant4_INCLUDE_DIRS} ) +#---Add Library------------------------------------------------------------------- +if(DD4HEP_USE_BOOST) + FIND_PACKAGE( Boost REQUIRED) + include_directories( ${Boost_INCLUDE_DIR}) + add_definitions(-DDD4HEP_USE_BOOST) +endif() #---Add Library------------------------------------------------------------------- file(GLOB sources src/*.cpp) @@ -18,8 +30,8 @@ list(APPEND sources G__DDG4.cxx) add_library(DD4hepG4 SHARED ${sources}) target_link_libraries(DD4hepG4 DD4hepCore ${ROOT_LIBRARIES} Reflex ${Geant4_LIBRARIES}) SET_TARGET_PROPERTIES( DD4hepG4 PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION}) +#----------------------------------------------------------------------------------- # No rootmap for link libraries! - #--------------------------- Legacy libraries (for Frank) ------------------------- file(GLOB legacy_sources legacy/*.cpp) add_library(DD4hepG4Legacy SHARED ${legacy_sources}) @@ -37,19 +49,14 @@ SET_TARGET_PROPERTIES(DD4hepG4Plugins PROPERTIES VERSION ${DD4hep_VERSION} SOVER dd4hep_generate_rootmap(DD4hepG4Plugins) #--------------------------- LCIO Plugins for new simulation framework ----------- -#if(LCIO_DIR) -# list(APPEND include_directories ${LCIO_DIR}/include) - if(DD4HEP_USE_LCIO) - + find_package(LCIO REQUIRED) + include_directories( ${LCIO_INCLUDE_DIRS} ) file(GLOB lcio_sources lcio/*.cpp) add_library(DD4hepG4LCIO SHARED ${lcio_sources}) - target_link_libraries(DD4hepG4LCIO DD4hepCore DD4hepG4 ${Geant4_LIBRARIES} ${LCIO_LIBRARIES} ${ROOT_LIBRARIES} Reflex) - SET_TARGET_PROPERTIES(DD4hepG4LCIO PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION}) dd4hep_generate_rootmap(DD4hepG4LCIO) - endif() #----------------------------------------------------------------------------------- diff --git a/DDG4/examples/dictionaries.C b/DDG4/examples/dictionaries.C index 6a2afdd13..e559d6247 100644 --- a/DDG4/examples/dictionaries.C +++ b/DDG4/examples/dictionaries.C @@ -12,6 +12,7 @@ // Framework include files #include "DDG4/Geant4Data.h" #include <vector> +#include <set> using namespace std; using namespace DD4hep; @@ -19,7 +20,7 @@ using namespace DD4hep::Simulation; // CINT configuration #if defined(__MAKECINT__) -//#pragma link C++ class Position+; +#pragma link C++ class Position+; //#pragma link C++ class Direction+; #pragma link C++ class SimpleRun+; #pragma link C++ class SimpleEvent+; @@ -34,6 +35,9 @@ using namespace DD4hep::Simulation; #pragma link C++ class SimpleCalorimeter+; #pragma link C++ class SimpleCalorimeter::Hit+; #pragma link C++ class std::vector<SimpleCalorimeter::Hit*>+; +#pragma link C++ class Particle+; +#pragma link C++ class std::vector<Particle*>+; +#pragma link C++ class std::set<int>+; //#pragma link C++ class ; #endif diff --git a/DDG4/examples/run.C b/DDG4/examples/run.C index ecf14ab74..20dde4f02 100644 --- a/DDG4/examples/run.C +++ b/DDG4/examples/run.C @@ -3,7 +3,11 @@ void run() { gInterpreter->ProcessLine(".X initAClick.C"); gInterpreter->ProcessLine(".L dictionaries.C+"); //gInterpreter->ProcessLine(".L exampleAClick.C+"); - gInterpreter->ProcessLine(".L xmlAClick.C+"); + //gInterpreter->ProcessLine(".L xmlAClick.C+"); //gInterpreter->ProcessLine(".L TEve.C+"); + gSystem->Load("libDD4hepCore"); + gSystem->Load("libDD4hepG4"); + gInterpreter->ProcessLine(".X DDG4Dump.C+"); + //gInterpreter->ProcessLine(".X a.C++"); } diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h index 3c30784e5..ef847ec0f 100644 --- a/DDG4/include/DDG4/Geant4Data.h +++ b/DDG4/include/DDG4/Geant4Data.h @@ -13,7 +13,11 @@ #ifndef __DD4HEP_DDEVE_EXCLUSIVE__ #include "DD4hep/Objects.h" #include "G4Step.hh" +#else +typedef void G4VProcess; +typedef void G4ParticleDefinition; #endif +#include <set> /* * DD4hep namespace declaration @@ -117,12 +121,12 @@ namespace DD4hep { */ class Particle { public: - int id, g4Parent, parent, reason, steps, secondaries; + int id, g4Parent, parent, reason, steps, secondaries, pdgID; double vsx, vsy, vsz; double vex, vey, vez; double psx, psy, psz, pex, pey, pez, energy, time; - const G4VProcess *process; - const G4ParticleDefinition *definition; + const G4VProcess *process; //! + const G4ParticleDefinition *definition; //! std::set<int> daughters; /// Default constructor Particle(); @@ -177,9 +181,9 @@ namespace DD4hep { MonteCarloContrib& operator=(const MonteCarloContrib& c) { if ( this != &c ) { trackID = c.trackID; - pdgID = c.pdgID; + pdgID = c.pdgID; deposit = c.deposit; - time = c.time; + time = c.time; } return *this; } diff --git a/DDG4/include/DDG4/Geant4DataDump.h b/DDG4/include/DDG4/Geant4DataDump.h new file mode 100644 index 000000000..e754e65b4 --- /dev/null +++ b/DDG4/include/DDG4/Geant4DataDump.h @@ -0,0 +1,74 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4DATADUMP_H +#define DD4HEP_DDG4_GEANT4DATADUMP_H + +// Framework include files +#include "DD4hep/Printout.h" +#include "DDG4/Geant4Data.h" + +// C/C++ include files +#include <vector> + + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + /** Class to dump the records of the intrinsic Geant4 event model. + * + * @author M.Frank + * @version 1.0 + */ + class Geant4DataDump { + public: + typedef std::vector<Particle*> Particles; + + typedef SimpleTracker::Hit TrackerHit; + typedef std::vector<SimpleTracker::Hit*> TrackerHits; + + typedef SimpleCalorimeter::Hit CalorimeterHit; + typedef std::vector<SimpleCalorimeter::Hit*> CalorimeterHits; + + protected: + /// Tag variable + std::string m_tag; + + public: + /// Default constructor + Geant4DataDump(const std::string& tag); + /// Standard destructor + virtual ~Geant4DataDump(); + + /// Print a single particle to the output logging using the specified print level + void print(PrintLevel level, const Particle* p) const; + /// Print the particle container to the output logging using the specified print level + void print(PrintLevel level, const std::string& container, const Particles* parts); + + /// Print a single tracker hit to the output logging using the specified print level + void print(PrintLevel level, const TrackerHit* h); + /// Print the tracker hits container to the output logging using the specified print level + void print(PrintLevel level, const std::string& container, const TrackerHits* hits); + + /// Print a calorimeter tracker hit to the output logging using the specified print level + void print(PrintLevel level, const CalorimeterHit* h); + /// Print the calorimeter hits container to the output logging using the specified print level + void print(PrintLevel level, const std::string& container, const CalorimeterHits* hits); + + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif /* DD4HEP_DDG4_GEANT4DATADUMP_H */ diff --git a/DDG4/include/DDG4/Geant4HitCollection.h b/DDG4/include/DDG4/Geant4HitCollection.h index 69fb9255e..198677724 100644 --- a/DDG4/include/DDG4/Geant4HitCollection.h +++ b/DDG4/include/DDG4/Geant4HitCollection.h @@ -150,7 +150,8 @@ namespace DD4hep { /// Automatic conversion to the desired type template <typename TYPE> operator TYPE*() const { return (TYPE*) m_data.second-> - cast.apply_downCast(ComponentCast::instance<TYPE>(),m_data.first); + cast.apply_dynCast(ComponentCast::instance<TYPE>(),m_data.first); + //cast.apply_downCast(ComponentCast::instance<TYPE>(),m_data.first); } }; diff --git a/DDG4/include/DDG4/Geant4MonteCarloTruth.h b/DDG4/include/DDG4/Geant4MonteCarloTruth.h index 3db9049fa..f19f64f54 100644 --- a/DDG4/include/DDG4/Geant4MonteCarloTruth.h +++ b/DDG4/include/DDG4/Geant4MonteCarloTruth.h @@ -54,6 +54,8 @@ namespace DD4hep { virtual const ParticleMap& particles() const = 0; /// Access the map of track equivalents virtual const TrackEquivalents& equivalents() const = 0; + /// Access the equivalent track id (shortcut to the usage of TrackEquivalents) + virtual int particleID(int track, bool throw_if_not_found=true) const = 0; /// Mark a Geant4 track to be kept for later MC truth analysis virtual void mark(const G4Track* track) = 0; /// Store a track, with a flag @@ -87,7 +89,9 @@ namespace DD4hep { /// Access the particle map virtual const ParticleMap& particles() const { return m_particleMap; } /// Access the map of track equivalents - virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; } + virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; } + /// Access the equivalent track id (shortcut to the usage of TrackEquivalents) + virtual int particleID(int track, bool) const { return track; } /// Mark a Geant4 track to be kept for later MC truth analysis. Default flag: CREATED_HIT virtual void mark(const G4Track* track); /// Store a track, with a flag diff --git a/DDG4/include/DDG4/Geant4Output2ROOT.h b/DDG4/include/DDG4/Geant4Output2ROOT.h index de2683d34..dd86c7731 100644 --- a/DDG4/include/DDG4/Geant4Output2ROOT.h +++ b/DDG4/include/DDG4/Geant4Output2ROOT.h @@ -63,6 +63,9 @@ namespace DD4hep { virtual void beginRun(const G4Run* run); /// Callback to store each Geant4 hit collection virtual void saveCollection(OutputContext<G4Event>& ctxt, G4VHitsCollection* collection); + /// Callback to store the Geant4 event + virtual void saveEvent(OutputContext<G4Event>& ctxt); + /// Commit data at end of filling procedure virtual void commit(OutputContext<G4Event>& ctxt); }; diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h index a31632779..3c2cacc18 100644 --- a/DDG4/include/DDG4/Geant4ParticleHandler.h +++ b/DDG4/include/DDG4/Geant4ParticleHandler.h @@ -90,7 +90,10 @@ namespace DD4hep { /// Access the particle map virtual const ParticleMap& particles() const { return m_particleMap; } /// Access the map of track equivalents - virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; } + virtual const TrackEquivalents& equivalents() const { return m_equivalentTracks; } + /// Access the equivalent track id (shortcut to the usage of TrackEquivalents) + virtual int particleID(int track, bool throw_if_not_found=true) const; + /// Mark a Geant4 track to be kept for later MC truth analysis. Default flag: CREATED_HIT virtual void mark(const G4Track* track); /// Store a track diff --git a/DDG4/lcio/LCIOConversions.cpp b/DDG4/lcio/LCIOConversions.cpp index ff82dae18..0cd2b1b1e 100644 --- a/DDG4/lcio/LCIOConversions.cpp +++ b/DDG4/lcio/LCIOConversions.cpp @@ -73,27 +73,29 @@ namespace DD4hep { template <> lcio::LCCollectionVec* Geant4DataConversion<lcio::LCCollectionVec, pair<VolMgr,Geant4HitCollection*>, - SimpleTracker::Hit>::operator()(const arg_t& args) const - { + SimpleTracker::Hit>::operator()(const arg_t& args) const { + Geant4HitCollection* coll = args.second; size_t nhits = coll->GetSize(); - SimpleHit* hit = coll->hit(0); - string dsc = encoding(args.first, hit->cellID); lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::SIMTRACKERHIT); - UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc_coll); lc_coll->reserve(nhits); - for(size_t i=0; i<nhits; ++i) { - const SimpleTracker::Hit* g4_hit = coll->hit(i); - double pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()}; - lcio::SimTrackerHitImpl* lc_hit = new lcio::SimTrackerHitImpl; - lc_hit->setCellID0( (g4_hit->cellID >> 0 ) & 0xFFFFFFFF); - lc_hit->setCellID1( (g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF); - lc_hit->setPosition(pos); - lc_hit->setEDep(g4_hit->energyDeposit); - lc_hit->setTime(g4_hit->truth.time); - lc_hit->setMomentum( g4_hit->momentum.x(), g4_hit->momentum.y() , g4_hit->momentum.z() ); - lc_hit->setPathLength( g4_hit->length); - lc_coll->addElement(lc_hit); + if ( nhits > 0 ) { + SimpleHit* hit = coll->hit(0); + string dsc = encoding(args.first, hit->cellID); + UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc_coll); + for(size_t i=0; i<nhits; ++i) { + const SimpleTracker::Hit* g4_hit = coll->hit(i); + double pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()}; + lcio::SimTrackerHitImpl* lc_hit = new lcio::SimTrackerHitImpl; + lc_hit->setCellID0( (g4_hit->cellID >> 0 ) & 0xFFFFFFFF); + lc_hit->setCellID1( (g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF); + lc_hit->setPosition(pos); + lc_hit->setEDep(g4_hit->energyDeposit); + lc_hit->setTime(g4_hit->truth.time); + lc_hit->setMomentum( g4_hit->momentum.x(), g4_hit->momentum.y() , g4_hit->momentum.z() ); + lc_hit->setPathLength( g4_hit->length); + lc_coll->addElement(lc_hit); + } } return lc_coll; } @@ -111,25 +113,26 @@ namespace DD4hep { template <> lcio::LCCollectionVec* Geant4DataConversion<lcio::LCCollectionVec, pair<VolMgr,Geant4HitCollection*>, - SimpleCalorimeter::Hit>::operator()(const arg_t& args) const - { + SimpleCalorimeter::Hit>::operator()(const arg_t& args) const { Geant4HitCollection* coll = args.second; size_t nhits = coll->GetSize(); - SimpleHit* hit = coll->hit(0); - string dsc = encoding(args.first, hit->cellID); lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::SIMCALORIMETERHIT); - UTIL::CellIDEncoder<SimCalorimeterHit> decoder(dsc,lc_coll); lc_coll->setFlag(UTIL::make_bitset32(LCIO::CHBIT_LONG,LCIO::CHBIT_STEP,LCIO::CHBIT_ID1)); lc_coll->reserve(nhits); - for(size_t i=0; i<nhits; ++i) { - const SimpleCalorimeter::Hit* g4_hit = coll->hit(i); - float pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()}; - lcio::SimCalorimeterHitImpl* lc_hit = new lcio::SimCalorimeterHitImpl; - lc_hit->setCellID0( ( g4_hit->cellID >> 0 ) & 0xFFFFFFFF ); - lc_hit->setCellID1( ( g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF ); - lc_hit->setPosition(pos); - lc_hit->setEnergy( g4_hit->energyDeposit ); - lc_coll->addElement(lc_hit); + if ( nhits > 0 ) { + SimpleHit* hit = coll->hit(0); + string dsc = encoding(args.first, hit->cellID); + UTIL::CellIDEncoder<SimCalorimeterHit> decoder(dsc,lc_coll); + for(size_t i=0; i<nhits; ++i) { + const SimpleCalorimeter::Hit* g4_hit = coll->hit(i); + float pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()}; + lcio::SimCalorimeterHitImpl* lc_hit = new lcio::SimCalorimeterHitImpl; + lc_hit->setCellID0( ( g4_hit->cellID >> 0 ) & 0xFFFFFFFF ); + lc_hit->setCellID1( ( g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF ); + lc_hit->setPosition(pos); + lc_hit->setEnergy( g4_hit->energyDeposit ); + lc_coll->addElement(lc_hit); + } } return lc_coll; } diff --git a/DDG4/plugins/Geant4EscapeCounter.cpp b/DDG4/plugins/Geant4EscapeCounter.cpp index 9ec56ac4b..4b9a46ee3 100644 --- a/DDG4/plugins/Geant4EscapeCounter.cpp +++ b/DDG4/plugins/Geant4EscapeCounter.cpp @@ -9,6 +9,7 @@ #ifndef DD4HEP_DDG4_GEANT4ESCAPECOUNTER_H #define DD4HEP_DDG4_GEANT4ESCAPECOUNTER_H +// Framework include files #include "DD4hep/Detector.h" #include "DDG4/Geant4SensDetAction.h" #include "DDG4/Geant4SteppingAction.h" @@ -23,7 +24,7 @@ namespace DD4hep { */ namespace Simulation { - /** @class Geant4EscapeCounter Geant4EscapeCounter.h DDG4/Geant4EscapeCounter.h + /** Class to measure the energy of escaping tracks of a detector using Geant 4 * * Measure escaping energy.... * diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp index fad59055b..0a4dfe0a3 100644 --- a/DDG4/plugins/Geant4SDActions.cpp +++ b/DDG4/plugins/Geant4SDActions.cpp @@ -113,6 +113,19 @@ namespace DD4hep { return pos == h->position ? h : 0; } }; + template<typename TYPE> struct CellIDCompare : public Geant4HitCollection::Compare { + long long int id; + Geant4HitWrapper::HitManipulator* manip; + /// Constructor + CellIDCompare(long long int i) : id(i), manip(Geant4HitWrapper::HitManipulator::instance<TYPE>()) { } + /// Comparison function + virtual void* operator()(const Geant4HitWrapper& w) const { + TYPE* h = w;//(const TYPE*)manip->castHit<TYPE>(w); + if ( id == h->cellID ) + return h; + return 0; + } + }; /// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ /// Geant4SensitiveAction<SimpleTracker> @@ -179,13 +192,13 @@ namespace DD4hep { Position pos = 0.5 * (h.prePos() + h.postPos()); HitContribution contrib = Hit::extractContribution(step); HitCollection* coll = collection(m_collectionID); - Hit* hit = 0;//coll->find<Hit>(PositionCompare<Hit>(pos)); + long long int cell = cellID(step); + Hit* hit = coll->find<Hit>(CellIDCompare<Hit>(cell)); if ( step->GetTotalEnergyDeposit() < 1e-5 ) return true; if ( !hit ) { Geant4TouchableHandler handler(step); - //hit = new Hit(pos); hit = new Hit(h.prePos()); - hit->cellID = cellID(step); + hit->cellID = cell; coll->add(hit); print("SimpleCalorimeter","%s> CREATE hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f %s", c_name(),contrib.deposit,pos.X(),pos.Y(),pos.Z(),handler.path().c_str()); @@ -195,8 +208,8 @@ namespace DD4hep { } } else { - print("SimpleCalorimeter","%s> UPDATE hit with deposit:%7.3f MeV Pos:%8.2f %8.2f %8.2f", - c_name(),contrib.deposit,pos.X(),pos.Y(),pos.Z()); + print("SimpleCalorimeter","%s> UPDATE hit with deposit:%7.3f MeV Pos:%8.2f %8.2f %8.2f id:%016llX cell:%016llX", + c_name(),contrib.deposit,pos.X(),pos.Y(),pos.Z(),cell,hit->cellID); } hit->truth.push_back(contrib); hit->energyDeposit += contrib.deposit; diff --git a/DDG4/python/DDG4Dict.C b/DDG4/python/DDG4Dict.C index fff56b40b..632eeeb75 100644 --- a/DDG4/python/DDG4Dict.C +++ b/DDG4/python/DDG4Dict.C @@ -17,6 +17,7 @@ using namespace std; using namespace DD4hep; using namespace DD4hep::Simulation; + // CINT configuration #if defined(__MAKECINT__) //#pragma link C++ class Position+; @@ -25,6 +26,7 @@ using namespace DD4hep::Simulation; #pragma link C++ class SimpleEvent+; //#pragma link C++ class SimpleEvent::Seeds+; #pragma link C++ class SimpleHit+; +#pragma link C++ class Particle+; #pragma link C++ class std::vector<SimpleHit*>+; #pragma link C++ class SimpleHit::Contribution+; #pragma link C++ class SimpleHit::Contributions+; @@ -34,10 +36,13 @@ using namespace DD4hep::Simulation; #pragma link C++ class SimpleCalorimeter+; #pragma link C++ class SimpleCalorimeter::Hit+; #pragma link C++ class std::vector<SimpleCalorimeter::Hit*>+; +#pragma link C++ class std::vector<Particle*>+; + //#pragma link C++ class ; #endif #include "DDG4/Geant4Config.h" +#include "DDG4/Geant4DataDump.h" #include <iostream> namespace DD4hep { @@ -206,6 +211,8 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation; #pragma link C++ class PhysicsListActionSequenceHandle; #pragma link C++ class SensDetActionSequenceHandle; +#pragma link C++ class Geant4DataDump; + #pragma link C++ class Geant4ActionCreation; #pragma link C++ class Geant4Action; #pragma link C++ class Geant4Kernel; diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp index 684d79567..92459d447 100644 --- a/DDG4/src/Geant4Data.cpp +++ b/DDG4/src/Geant4Data.cpp @@ -43,7 +43,8 @@ SimpleEvent::~SimpleEvent() { /// Copy constructor Particle::Particle(const Particle& c) - : id(c.id), g4Parent(c.g4Parent), parent(c.parent), reason(c.reason), steps(c.steps), + : id(c.id), g4Parent(c.g4Parent), parent(c.parent), reason(c.reason), + steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID), vsx(c.vsx), vsy(c.vsy), vsz(c.vsz), vex(c.vex), vey(c.vey), vez(c.vez), psx(c.psx), psy(c.psy), psz(c.psz), @@ -57,7 +58,8 @@ Particle::Particle(const Particle& c) /// Default constructor Particle::Particle() - : id(0), g4Parent(0), parent(0), reason(0), steps(0), + : id(0), g4Parent(0), parent(0), reason(0), + steps(0), secondaries(0), pdgID(0), vsx(0.0), vsy(0.0), vsz(0.0), vex(0.0), vey(0.0), vez(0.0), psx(0.0), psy(0.0), psz(0.0), @@ -145,9 +147,9 @@ SimpleTracker::Hit& SimpleTracker::Hit::storePoint(G4Step* step, G4StepPoint* pn G4ThreeVector mom = pnt->GetMomentum(); truth.trackID = trk->GetTrackID(); - truth.pdgID = trk->GetDefinition()->GetPDGEncoding(); + truth.pdgID = trk->GetDefinition()->GetPDGEncoding(); truth.deposit = step->GetTotalEnergyDeposit(); - truth.time = trk->GetGlobalTime(); + truth.time = trk->GetGlobalTime(); position.SetXYZ(pos.x(), pos.y(), pos.z()); momentum.SetXYZ(mom.x(), mom.y(), mom.z()); length = 0; diff --git a/DDG4/src/Geant4DataDump.cpp b/DDG4/src/Geant4DataDump.cpp new file mode 100644 index 000000000..176dcdae3 --- /dev/null +++ b/DDG4/src/Geant4DataDump.cpp @@ -0,0 +1,97 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +// Framework include files +#include "DD4hep/Primitives.h" +#include "DDG4/Geant4DataDump.h" + +using namespace std; +using namespace DD4hep; +using namespace DD4hep::Simulation; + +typedef ReferenceBitMask<const int> PropertyMask; + +/// Default constructor +Geant4DataDump::Geant4DataDump(const std::string& tag) : m_tag(tag) { +} + +/// Standard destructor +Geant4DataDump::~Geant4DataDump() { +} + +/// Print a single particle to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const Particle* p) const { + PropertyMask mask(p->reason); + printout(level, m_tag, " +++ TrackID: %6d %12d %6d %-7s %3s %5d %6s %8.3g %-4s %-7s %-7s %-3s", + p->id, + p->pdgID, + p->parent, + yes_no(mask.isSet(G4PARTICLE_PRIMARY)), + yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)), + int(p->daughters.size()), + yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)), + p->energy, + yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)), + yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)), + yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)), + mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "" + ); +} + +/// Print the particle container to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const std::string& container, const Particles* parts) { + if ( parts ) { + printout(level,m_tag,"+++ Track container: %-21s --------------- Track KEEP reasoning ---------------",container.c_str()); + printout(level,m_tag,"+++ # of Tracks:%6d PDG Parent Primary Secondary Energy %-8s Calo Tracker Process/Par", + int(parts->size()),"in [MeV]"); + for(Particles::const_iterator i=parts->begin(); i!= parts->end(); ++i) + print(PrintLevel(level-1), *i); + return; + } +} + +/// Print a single tracker hit to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const TrackerHit* h) { + const SimpleHit::Contribution& t = h->truth; + printout(level,m_tag," +++ Hit: Cell: %016llX Pos:(%9.3g,%9.3g,%9.3g) Len:%9.3g [mm] E:%9.3g MeV TrackID:%6d PDG:%12d dep:%9.3g time:%9.3g [ns]", + h->cellID,h->position.x(),h->position.y(),h->position.z(),h->length,h->energyDeposit,t.trackID,t.pdgID,t.deposit,t.time); +} + +/// Print the tracker hits container to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const std::string& container, const TrackerHits* hits) { + if ( hits ) { + printout(level,m_tag,"+++ %s: # Tracker hits %d",container.c_str(),int(hits->size())); + for(TrackerHits::const_iterator i=hits->begin(); i!= hits->end(); ++i) + print(PrintLevel(level-1), *i); + return; + } +} + +/// Print a calorimeter tracker hit to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const CalorimeterHit* h) { + printout(level,m_tag," +++ Hit: Cell: %016llX Pos:(%9.3g,%9.3g,%9.3g) [mm] E:%9.3g MeV #Contributions:%3d", + h->cellID,h->position.x(),h->position.y(),h->position.z(),h->energyDeposit,h->truth.size()); + const SimpleHit::Contributions& t = h->truth; + int cnt=0; + for(SimpleHit::Contributions::const_iterator i=t.begin(); i!=t.end(); ++i,++cnt) { + const SimpleHit::Contribution& c = *i; + printout(PrintLevel(level-1),m_tag," Contribution #%3d TrackID:%6d PDG:%12d %9.3g MeV %9.3g ns", + cnt,c.trackID,c.pdgID,c.deposit,c.time); + } +} + +/// Print the calorimeter hits container to the output logging using the specified print level +void Geant4DataDump::print(PrintLevel level, const std::string& container, const CalorimeterHits* hits) { + if ( hits ) { + printout(level,m_tag,"+++ %s: # Calorimeter hits %d",container.c_str(),int(hits->size())); + for(CalorimeterHits::const_iterator i=hits->begin(); i!= hits->end(); ++i) + print(PrintLevel(level-1), *i); + return; + } +} diff --git a/DDG4/src/Geant4Output2ROOT.cpp b/DDG4/src/Geant4Output2ROOT.cpp index e5bc2e6a5..118465c28 100644 --- a/DDG4/src/Geant4Output2ROOT.cpp +++ b/DDG4/src/Geant4Output2ROOT.cpp @@ -8,10 +8,13 @@ //==================================================================== // Framework include files +#include "DD4hep/Printout.h" #include "DD4hep/Primitives.h" #include "DD4hep/InstanceCount.h" #include "DDG4/Geant4HitCollection.h" #include "DDG4/Geant4Output2ROOT.h" +#include "DDG4/Geant4MonteCarloTruth.h" +#include "DDG4/Geant4Data.h" // Geant4 include files #include "G4HCofThisEvent.hh" @@ -44,7 +47,7 @@ Geant4Output2ROOT::~Geant4Output2ROOT() { } /// Create/access tree by name -TTree* Geant4Output2ROOT::section(const std::string& nam) { +TTree* Geant4Output2ROOT::section(const string& nam) { Sections::const_iterator i = m_sections.find(nam); if (i == m_sections.end()) { TDirectory::TContext context(m_file); @@ -70,7 +73,7 @@ void Geant4Output2ROOT::beginRun(const G4Run* run) { } /// Fill single EVENT branch entry (Geant4 collection data) -int Geant4Output2ROOT::fill(const std::string& nam, const ComponentCast& type, void* ptr) { +int Geant4Output2ROOT::fill(const string& nam, const ComponentCast& type, void* ptr) { if (m_file) { TBranch* b = 0; Branches::const_iterator i = m_branches.find(nam); @@ -130,13 +133,56 @@ void Geant4Output2ROOT::commit(OutputContext<G4Event>& ctxt) { Geant4OutputAction::commit(ctxt); } +/// Callback to store the Geant4 event +void Geant4Output2ROOT::saveEvent(OutputContext<G4Event>& /* ctxt */) { + Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false); + if ( truth ) { + typedef Geant4HitWrapper::HitManipulator Manip; + typedef Geant4MonteCarloTruth::ParticleMap ParticleMap; + Manip* manipulator = Geant4HitWrapper::manipulator<Particle>(); + const ParticleMap& pm = truth->particles(); + vector<void*> particles; + for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end(); ++i) { + particles.push_back((ParticleMap::mapped_type*)(*i).second); + } + fill("MCParticles",manipulator->vec_type,&particles); + } +} + /// Callback to store each Geant4 hit collection void Geant4Output2ROOT::saveCollection(OutputContext<G4Event>& /* ctxt */, G4VHitsCollection* collection) { Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection); - std::string hc_nam = collection->GetName(); - std::vector<void*> hits; + string hc_nam = collection->GetName(); + vector<void*> hits; if (coll) { coll->getHitsUnchecked(hits); + size_t nhits = coll->GetSize(); + Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false); + if ( truth && nhits > 0 ) { + try { + for(size_t i=0; i<nhits; ++i) { + SimpleHit* h = coll->hit(i); + SimpleTracker::Hit* trk_hit = dynamic_cast<SimpleTracker::Hit*>(h); + if ( 0 != trk_hit ) { + SimpleHit::Contribution& t = trk_hit->truth; + int trackID = t.trackID; + t.trackID = truth->particleID(trackID); + } + SimpleCalorimeter::Hit* cal_hit = dynamic_cast<SimpleCalorimeter::Hit*>(h); + if ( 0 != cal_hit ) { + SimpleHit::Contributions& c = cal_hit->truth; + for(SimpleHit::Contributions::iterator j=c.begin(); j!=c.end(); ++j) { + SimpleHit::Contribution& t = *j; + int trackID = t.trackID; + t.trackID = truth->particleID(trackID); + } + } + } + } + catch(...) { + printout(ERROR,name(),"+++ Exception while saving collection %s.",hc_nam.c_str()); + } + } fill(hc_nam, coll->vector_type(), &hits); } } diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp index e268e2d88..7f345fd78 100644 --- a/DDG4/src/Geant4ParticleHandler.cpp +++ b/DDG4/src/Geant4ParticleHandler.cpp @@ -162,26 +162,29 @@ void Geant4ParticleHandler::begin(const G4Track* track) { // Extract here all the information from the start of the tracking action // which we will need later to create a proper MC particle. - m_currTrack.id = h.id(); - m_currTrack.reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0; - m_currTrack.energy = kine; - m_currTrack.g4Parent = h.parent(); - m_currTrack.parent = h.parent(); - m_currTrack.definition = h.trackDef(); - m_currTrack.process = h.creatorProcess(); - m_currTrack.time = h.globalTime(); - m_currTrack.vsx = v.x(); - m_currTrack.vsy = v.y(); - m_currTrack.vsz = v.z(); - m_currTrack.vex = 0.0; - m_currTrack.vey = 0.0; - m_currTrack.vez = 0.0; - m_currTrack.psx = m.x(); - m_currTrack.psy = m.y(); - m_currTrack.psz = m.z(); - m_currTrack.pex = 0.0; - m_currTrack.pey = 0.0; - m_currTrack.pez = 0.0; + m_currTrack.id = h.id(); + m_currTrack.reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0; + m_currTrack.steps = 0; + m_currTrack.secondaries = 0; + m_currTrack.pdgID = h.trackDef()->GetPDGEncoding(); + m_currTrack.energy = kine; + m_currTrack.g4Parent = h.parent(); + m_currTrack.parent = h.parent(); + m_currTrack.definition = h.trackDef(); + m_currTrack.process = h.creatorProcess(); + m_currTrack.time = h.globalTime(); + m_currTrack.vsx = v.x(); + m_currTrack.vsy = v.y(); + m_currTrack.vsz = v.z(); + m_currTrack.vex = 0.0; + m_currTrack.vey = 0.0; + m_currTrack.vez = 0.0; + m_currTrack.psx = m.x(); + m_currTrack.psy = m.y(); + m_currTrack.psz = m.z(); + m_currTrack.pex = 0.0; + m_currTrack.pey = 0.0; + m_currTrack.pez = 0.0; m_currTrack.daughters.clear(); // If the creator process of the track is in the list of process products to be kept, set the proper flag if ( m_currTrack.process ) { @@ -362,6 +365,8 @@ int Geant4ParticleHandler::recombineParents() { PropertyMask(parent_part->reason).set(mask.value()); // Remove track from parent's list of daughters parent_part->removeDaughter(id); + parent_part->steps += par->steps; + parent_part->secondaries += par->secondaries; } } } @@ -439,3 +444,8 @@ int Geant4ParticleHandler::equivalentTrack(int g4_id) const { } return equiv_id; } + +/// Access the equivalent track id (shortcut to the usage of TrackEquivalents) +int Geant4ParticleHandler::particleID(int track, bool) const { + return equivalentTrack(track); +} diff --git a/DDG4/src/Geant4ROOTDump.cpp b/DDG4/src/Geant4ROOTDump.cpp new file mode 100644 index 000000000..3d48b9358 --- /dev/null +++ b/DDG4/src/Geant4ROOTDump.cpp @@ -0,0 +1,159 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +// Framework include files +#include "DD4hep/LCDD.h" +#include "DD4hep/Factories.h" +#include "DD4hep/Primitives.h" +#include "DDG4/Geant4DataDump.h" + +#include "TInterpreter.h" +#include "TSystem.h" +#include "TFile.h" +#include "TTree.h" +#include "TROOT.h" + +using namespace std; +using namespace DD4hep; +using namespace DD4hep::Simulation; + +typedef Geant4DataDump::Particles Particles; +typedef Geant4DataDump::TrackerHits TrackerHits; +typedef Geant4DataDump::CalorimeterHits CalorimeterHits; + +static long usage() { + printout(FATAL,"Geant4ROOTDump","usage: Geant4ROOTDump -opt (app-opts) --opt=value (plugin-opts)"); + printout(FATAL,"Geant4ROOTDump"," app-opts: "); + printout(FATAL,"Geant4ROOTDump"," -print INFO Print container summaries only."); + printout(FATAL,"Geant4ROOTDump"," -print DEBUG Print container content as well."); + printout(FATAL,"Geant4ROOTDump"," -print VERBOSE Print full data."); + printout(FATAL,"Geant4ROOTDump"," plugin opts: "); + printout(FATAL,"Geant4ROOTDump"," --usage: Print this output."); + printout(FATAL,"Geant4ROOTDump"," --input=<file> Print content of this file."); + printout(FATAL,"Geant4ROOTDump"," --entry=<number> Print specified event in the file. (starts with 0!)"); + return 'H'; +} + +static pair<TClass*,void*> load(TBranch* branch, int entry) { + TClass* cl = gROOT->GetClass(branch->GetClassName(),kTRUE); + if ( !cl ) { + return pair<TClass*,void*>(0,0); + } + void *obj = cl->New(); + branch->SetAddress(&obj); + branch->SetAutoDelete(kFALSE); + int nb = branch->GetEntry(entry); + if ( nb < 0 ) { + cl->Destructor(obj); + return pair<TClass*,void*>(0,0); + } + return pair<TClass*,void*>(cl,obj); +} + +static long dump_root(DD4hep::Geometry::LCDD&, int argc, char** argv) { + std::string input = "", tag="Geant4ROOTDump"; + int entry = -1; + + for(int j=0; j<argc; ++j) { + std::string a = argv[j]; + if ( a[0]=='-' ) a = argv[j]+1; + if ( a[0]=='-' ) a = argv[j]+2; + size_t idx = a.find('='); + if ( strncmp(a.c_str(),"usage",5)==0 ) { + usage(); + return 1; + } + if ( idx > 0 ) { + string p1 = a.substr(0,idx); + string p2 = a.substr(idx+1); + if ( strncmp(p1.c_str(),"input",3)==0 ) { + input = p2; + } + if ( strncmp(p1.c_str(),"entry",5)==0 ) { + if ( 1 != ::sscanf(p2.c_str(),"%d",&entry) ) { + printout(FATAL,tag,"+++ Argument %s is not properly formatted. must be --entry=<number>.",argv[j]); + return usage(); + } + } + continue; + } + printout(FATAL,tag,"+++ Argument %s is IGNORED! No value is found (string has no '=')",argv[j]); + return usage(); + } + const char* line = "+----------------------------------------------------------------------------------+"; + printout(INFO,tag,line); + printout(INFO,tag,"| Input:%-74s|",input.c_str()); + printout(INFO,tag,line); + + if ( input.empty() ) { + return usage(); + } + + + Geant4DataDump dump("Geant4Data"); + TFile* f = TFile::Open(input.c_str()); + + if ( f && !f->IsZombie() ) { + TTree* tree = (TTree*)f->Get("EVENT"); + TClass* cl_calo = gROOT->GetClass(typeid(CalorimeterHits)); + TClass* cl_tracker = gROOT->GetClass(typeid(TrackerHits)); + TClass* cl_particles = gROOT->GetClass(typeid(Particles)); + TObjArray* branches = tree->GetListOfBranches(); + Int_t nbranches = branches->GetEntriesFast(); + typedef pair<TClass*,void*> ENTRY; + typedef map<string,ENTRY> ENTRIES; + + for(Int_t ievt=entry<0 ? 0 : entry, nevt=entry<0 ? tree->GetEntries() : entry+1; ievt<nevt; ++ievt) { + ENTRIES event; + printout(INFO,tag,line); + printout(INFO,tag,"| Event: %6d |",ievt); + printout(INFO,tag,line); + + // First suck in all data + for (Int_t i=0;i<nbranches;i++) { + TBranch* branch = (TBranch*)branches->UncheckedAt(i); + pair<TClass*,void*> data = load(branch,ievt); + if ( data.first ) event[branch->GetName()] = data; + } + // Now dump the stuff + for(ENTRIES::const_iterator i=event.begin(); i!=event.end(); ++i) { + pair<TClass*,void*> data = (*i).second; + if ( data.first == cl_particles ) { + Particles* parts = (Particles*)data.second; + dump.print(INFO, (*i).first, parts); + for_each(parts->begin(), parts->end(), DestroyObject<Particle*>()); + } + } + for(ENTRIES::const_iterator i=event.begin(); i!=event.end(); ++i) { + pair<TClass*,void*> data = (*i).second; + if ( data.first == cl_particles ) { + } + else if ( data.first == cl_tracker ) { + TrackerHits* hits = (TrackerHits*)data.second; + dump.print(INFO, (*i).first, hits); + for_each(hits->begin(), hits->end(), DestroyObject<SimpleTracker::Hit*>()); + } + else if ( data.first == cl_calo ) { + CalorimeterHits* hits = (CalorimeterHits*)data.second; + dump.print(INFO, (*i).first, hits); + for_each(hits->begin(), hits->end(), DestroyObject<SimpleCalorimeter::Hit*>()); + } + if ( data.first ) data.first->Destructor(data.second); + } + } + delete f; + } + else { + printout(FATAL,tag,"+++ FAILED to open input file:%s",input.c_str()); + return usage(); + } + return 1; +} + +DECLARE_APPLY(Geant4ROOTDump,dump_root) diff --git a/UtilityApps/src/plugin_runner.cpp b/UtilityApps/src/plugin_runner.cpp index 4931d5871..6faae2f4c 100644 --- a/UtilityApps/src/plugin_runner.cpp +++ b/UtilityApps/src/plugin_runner.cpp @@ -13,9 +13,9 @@ //______________________________________________________________________________ namespace { void usage() { - cout << "geoPlugin -opt [-opt] \n" + cout << "geoPluginRun -opt [-opt] \n" " -plugin <name> [REQUIRED] Plugin to be executed and applied. \n" - " -input <file> [REQUIRED] Specify input file. \n"; + " -input <file> [OPTIONAL] Specify geometry input file. \n"; print_default_args() << endl; exit(EINVAL); } @@ -41,12 +41,18 @@ int main(int argc,char** argv) { usage(); } } - if ( arguments.geo_files.empty() || plugin.empty() ) + if ( plugin.empty() ) usage(); + options.push_back(0); LCDD& lcdd = dd4hep_instance(); - // Load compact files - load_compact(lcdd, arguments); + // Load compact files if required by plugin + if ( !arguments.geo_files.empty() ) { + load_compact(lcdd, arguments); + } + else { + cout << "geoPluginRun: No geometry input supplied. No geometry will be loaded." << endl; + } // Create volume manager and populate it required if ( arguments.volmgr ) run_plugin(lcdd,"DD4hepVolumeManager",0,0); // Execute plugin diff --git a/UtilityApps/src/run_plugin.h b/UtilityApps/src/run_plugin.h index b77871c5b..48a40eeec 100644 --- a/UtilityApps/src/run_plugin.h +++ b/UtilityApps/src/run_plugin.h @@ -51,16 +51,18 @@ namespace { return LCDD::getInstance(); } - int run_plugin(LCDD& lcdd, const char* name, int argc, char** argv) { + long run_plugin(LCDD& lcdd, const char* name, int argc, char** argv) { try { lcdd.apply(name,argc,argv); return 0; } catch(const exception& e) { cout << e.what() << endl; + return EINVAL; } catch(...) { cout << "UNKNOWN Exception" << endl; + return EINVAL; } ::exit(EINVAL); return EINVAL; @@ -106,7 +108,7 @@ namespace { print = DD4hep::INFO; } int handle(int& i, int argc, char** argv) { - if ( strncmp(argv[i],"-compact",2)==0 || strncmp(argv[i],"-input",2)==0 ) { + if ( strncmp(argv[i],"-compact",5)==0 || strncmp(argv[i],"-input",4)==0 ) { geo_files.push_back(argv[++i]); if ( argc>i+2 && strncmp(argv[i+1],"-build_type",6)==0 ) { build_types.push_back(argv[i+2]); @@ -116,13 +118,13 @@ namespace { build_types.push_back("BUILD_DEFAULT"); } } - else if ( strncmp(argv[i],"-load_only",2)==0 ) + else if ( strncmp(argv[i],"-load_only",5)==0 ) dry_run = true; - else if ( strncmp(argv[i],"-print",2)==0 ) + else if ( strncmp(argv[i],"-print",4)==0 ) DD4hep::setPrintLevel(DD4hep::PrintLevel(print = decodePrintLevel(argv[++i]))); - else if ( strncmp(argv[i],"-destroy",2)==0 ) + else if ( strncmp(argv[i],"-destroy",5)==0 ) destroy = true; - else if ( strncmp(argv[i],"-volmgr",2)==0 ) + else if ( strncmp(argv[i],"-volmgr",4)==0 ) volmgr = true; else return 0; @@ -193,7 +195,8 @@ namespace { if ( !args.dry_run ) { pair<int, char**> a(0,0); TRint app(name, &a.first, a.second); - run_plugin(lcdd,name,a.first,a.second); + long result = run_plugin(lcdd,name,a.first,a.second); + if ( result == EINVAL ) usage_default(name); app.Run(); } else { diff --git a/doc/release.notes b/doc/release.notes index 188ecb24d..38f473cd2 100644 --- a/doc/release.notes +++ b/doc/release.notes @@ -1,38 +1,15 @@ DD4hep ---- Release Notes ================================= - -------- -| v00-07 | seventh beta release ... - -------- - - - Some minor fixes: - - made compatible with older geant4 versions (9.5) - - add Bitflag to store CellID1 in SimCalorimeterHit collections - - fix position conversion from Geant4 to ROOT - - add cellID determination to SensitiveAction - - ... - +2014/08/07 Markus Frank +----------------------- + - DDG4: First version to support MC truth in DDG4 including + particle filtering to optimize the size of the MC record. + - DDG4 fix SimpleCalorimter sensitive action and properly support + hit aggregations. + - DDEve smaller modifications to support DDG4IO if DD4hep was + built with the Geant4 option ON. -Andre Sailer, 2014-07-17 ------------------------- - Unify cmake option variables, small cmake corrections - Change options _WITH_ to _USE_ - Print Warning that variables with _WITH_ are deprecated - Change Defintions to _USE_ as well - Add REQUIRED to find_package geant4 and xercesc if they are turned on - Updated documentation - Updated ILDExDet example - -Christian.Grefe, 2014-07-15 ---------------------------- - made DDSegmentation optionally a stand-alone package - create DDSgementationConfig.cmake when build as part of DD4hep - - - Markus Frank, 2014-07-02 - ------------------------ - - add LCIO conversions from DDSim - sensitive detectors - can now simply instantiate LCIO Sim hits 2014/06/30 Markus Frank ----------------------- diff --git a/examples/CLICSiD/eve/DDEve.xml b/examples/CLICSiD/eve/DDEve.xml index 496ebbdb7..f004bb6e5 100644 --- a/examples/CLICSiD/eve/DDEve.xml +++ b/examples/CLICSiD/eve/DDEve.xml @@ -71,6 +71,7 @@ <view name="Calo 3D" type="Calo3DProjection"> <calodata name="Ecal"/> + <calodata name="Hcal"/> <calodata name="EcalEndcap"/> <calodata name="HcalEndcap"/> </view> diff --git a/examples/ClientTests/scripts/FCC_Hcal.py b/examples/ClientTests/scripts/FCC_Hcal.py index 063aaf893..4be027936 100644 --- a/examples/ClientTests/scripts/FCC_Hcal.py +++ b/examples/ClientTests/scripts/FCC_Hcal.py @@ -19,7 +19,6 @@ def run(): kernel = DDG4.Kernel() kernel.setOutputLevel('Geant4Converter',Output.DEBUG) kernel.setOutputLevel('RootOutput',Output.INFO) - kernel.setOutputLevel('LcioOutput',Output.INFO) kernel.setOutputLevel('ShellHandler',Output.DEBUG) kernel.setOutputLevel('Gun',Output.INFO) kernel.UI = "UI" @@ -42,28 +41,6 @@ def run(): ui.SessionType = 'csh' kernel.registerGlobalAction(ui) - - # Configure Run actions: example only! - """ - run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') - run1.Property_int = 12345 - run1.Property_double = -5e15*keV - run1.Property_string = 'Startrun: Hello_2' - print run1.Property_string, run1.Property_double, run1.Property_int - run1.enableUI() - kernel.registerGlobalAction(run1) - kernel.runAction().add(run1) - """ - # Configure Event actions - """ - trk = DDG4.Action(kernel,"Geant4TrackPersistency/MonteCarloTruthHandler") - mc = DDG4.Action(kernel,"Geant4MonteCarloRecordManager/MonteCarloRecordManager") - kernel.registerGlobalAction(trk) - kernel.registerGlobalAction(mc) - trk.release() - mc.release() - """ - # Configure I/O evt_root = DDG4.EventAction(kernel,'Geant4Output2ROOT/RootOutput') evt_root.Control = True @@ -71,11 +48,6 @@ def run(): evt_root.enableUI() kernel.eventAction().add(evt_root) - #evt_lcio = DDG4.EventAction(kernel,'Geant4Output2LCIO/LcioOutput') - #evt_lcio.Output = "simple_lcio" - #evt_lcio.enableUI() - #kernel.eventAction().add(evt_lcio) - # Setup particle gun gun = DDG4.GeneratorAction(kernel,"Geant4ParticleGun/Gun") gun.energy = 100*GeV @@ -85,29 +57,15 @@ def run(): gun.isotrop = True gun.enableUI() kernel.generatorAction().add(gun) - """ - rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader") - rdr.zSpread = 0.0 - rdr.lorentzAngle = 0.0 - rdr.OutputLevel = DDG4.OutputLevel.INFO - rdr.Input = "LcioEventReader|test.data" - rdr.enableUI() - kernel.generatorAction().add(rdr) - """ # Setup global filters fur use in sensntive detectors f1 = DDG4.Filter(kernel,'GeantinoRejectFilter/GeantinoRejector') - f2 = DDG4.Filter(kernel,'EnergyDepositMinimumCut') - f2.Cut = 1e-14*keV - f2.enableUI() kernel.registerGlobalFilter(f1) - kernel.registerGlobalFilter(f2) # Now the calorimeters seq = DDG4.SensitiveSequence(kernel,'Geant4SensDetActionSequence/HcalBarrel') act = DDG4.SensitiveAction(kernel,'Geant4SimpleCalorimeterAction/HcalBarrelHandler','HcalBarrel') seq.add(f1) # ignore geantinos - #seq.add(f2) # ignore particles below threshold seq.add(act) seq = DDG4.SensitiveSequence(kernel,'Geant4SensDetActionSequence/ContainmentShell') @@ -117,29 +75,7 @@ def run(): # Now build the physics list: phys = kernel.physicsList() phys.extends = 'QGSP_BERT' - #phys.transportation = True - #phys.decays = True phys.enableUI() - """ - ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics') - ph.addParticleConstructor('G4BosonConstructor') - ph.addParticleConstructor('G4LeptonConstructor') - ph.addParticleConstructor('G4MesonConstructor') - ph.addParticleConstructor('G4BaryonConstructor') - ph.addParticleProcess('e[+-]','G4eMultipleScattering',-1,1,1) - #ph.addParticleProcess('e[+-]','G4eIonisation',-1,2,2) - ph.addParticleProcess('mu[+-]','G4MuMultipleScattering',-1,1,1) - #ph.addParticleProcess('mu[+-]','G4MuIonisation',-1,2,2) - ph.addParticleProcess('pi[+-]','G4hMultipleScattering',-1,2,2) - #ph.addParticleProcess('pi[+-]','G4hIonisation',-1,2,2) - ph.addParticleProcess('pi[+-]','G4hBremsstrahlung',-1,3,3) - ph.addParticleProcess('proton','G4hMultipleScattering',-1,2,2) - #ph.addParticleProcess('proton','G4hIonisation',-1,2,2) - ph.addParticleProcess('proton','G4hBremsstrahlung',-1,3,3) - ph.enableUI() - phys.add(ph) - """ - phys.dump() kernel.configure() -- GitLab