DDG4/src/HitTupleAction.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     : M.Frank
0011 //
0012 //==========================================================================
0013 #ifndef DD4HEP_DDG4_HITTUPLEACTION_H
0014 #define DD4HEP_DDG4_HITTUPLEACTION_H
0015
0016 // Framework include files
0017 #include "DDG4/Geant4EventAction.h"
0018
0019 // Forward declarations
0020 class G4VHitsCollection;
0021 class TFile;
0022 class TTree;
0023 class TBranch;
0024
0025 /// Namespace example name of the user
0026 namespace myanalysis {
0027
0028   // Forward declarations
0029   class Geant4ParticleMap;
0030
0031   /// Class to measure the energy of escaping tracks
0032   /** Class to measure the energy of escaping tracks of a detector using Geant 4
0033    * Measure escaping energy....
0034    *
0035    *  \author  M.Frank
0036    *  \version 1.0
0037    *  \ingroup DD4HEP_SIMULATION
0038    */
0039   class HitTupleAction : public dd4hep::sim::Geant4EventAction {
0040   public:
0041     typedef std::vector<std::string> CollectionNames;
0042     /// Property: collection names to be dumped
0043     CollectionNames m_containers;
0044     /// Property: Output file name
0045     std::string     m_outFileName;
0046     /// ROOT TFile instance
0047     TFile*          m_outFile = 0;
0048     /// Tree object within the file
0049     TTree*          m_outTree = 0;
0050     /// We want to write a separate branch for all deposits of one container(sub-detector)
0051     typedef std::pair<double, std::vector<double> > Data;
0052     /// The intermediate storage of the hit deposits to be written to ROOT
0053     std::map<std::string, std::pair<TBranch*, Data> > m_deposits;
0054
0055   public:
0056     /// Standard constructor
0057     HitTupleAction(dd4hep::sim::Geant4Context* context, const std::string& nam);
0058     /// Default destructor
0059     virtual ~HitTupleAction();
0060     /// Begin-of-run callback: Open output file
0061     void beginRun(const G4Run* run);
0062     /// End-of-run callback: Close output file
0063     void endRun(const G4Run* run);
0064     /// Geant4EventAction interface: Begin-of-event callback
0065     virtual void begin(const G4Event* event)  override;
0066     /// Geant4EventAction interface: End-of-event callback
0067     virtual void end(const G4Event* event)  override;
0068   };
0069 }      // End namespace dd4hep
0070 #endif /* DD4HEP_DDG4_HITTUPLEACTION_H */
0071
0072 //====================================================================
0073 //  AIDA Detector description implementation
0074 //--------------------------------------------------------------------
0075 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0076 // All rights reserved.
0077 //
0078 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0079 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0080 //
0081 // Author     : M.Frank
0082 //
0083 //====================================================================
0084
0085 // Framework include files
0086 #include "DDG4/Geant4DataDump.h"
0087 #include "DDG4/Geant4RunAction.h"
0088 #include "DDG4/Geant4HitCollection.h"
0089
0090 // Geant 4 includes
0091 #include "G4HCofThisEvent.hh"
0092 #include "G4Event.hh"
0093
0094 // ROOT include files
0095 #include "TFile.h"
0096 #include "TTree.h"
0097 #include "TROOT.h"
0098 #include "TClass.h"
0099 #include "TBranch.h"
0100
0101 using namespace std;
0102 using namespace dd4hep;
0103 using namespace dd4hep::sim;
0104
0105 /// Standard constructor
0106 myanalysis::HitTupleAction::HitTupleAction(Geant4Context* ctxt, const string& nam)
0107   : Geant4EventAction(ctxt, nam), m_containers{"*"}
0108 {
0109   m_needsControl = true;
0110   m_containers.push_back("*");
0111   declareProperty("OutputFile",   m_outFileName);
0112   declareProperty("Collections",  m_containers);
0113   runAction().callAtBegin(this,  &HitTupleAction::beginRun);
0114   runAction().callAtEnd(this,    &HitTupleAction::endRun);
0115 }
0116
0117 /// Default destructor
0118 myanalysis::HitTupleAction::~HitTupleAction() {
0119 }
0120
0121 /// Begin-of-run callback: Open output file
0122 void myanalysis::HitTupleAction::beginRun(const G4Run* /* run */)   {
0123 }
0124
0125 /// End-of-run callback: Close output file
0126 void myanalysis::HitTupleAction::endRun(const G4Run* /* run */)    {
0127   if ( m_outFile )   {
0128     if ( m_outTree )   {
0129       m_outTree->Print();
0130     }
0131     m_outFile->Write();
0132     m_outFile->Close();
0133     delete m_outFile;
0134     m_outFile = 0;
0135   }
0136 }
0137
0138 /// Geant4EventAction interface: Begin-of-event callback
0139 void myanalysis::HitTupleAction::begin(const G4Event* /* event */)   {
0140 }
0141
0142 /// Geant4EventAction interface: End-of-event callback
0143 void myanalysis::HitTupleAction::end(const G4Event* event)    {
0144   G4HCofThisEvent* hce = event->GetHCofThisEvent();
0145   if ( hce )  {
0146     int nCol = hce->GetNumberOfCollections();
0147     if ( nCol <= 0 )   {
0148       return;
0149     }
0150     else if ( !m_outFile )   {
0151       //                      Name                   Option      Title
0152       m_outFile = TFile::Open(m_outFileName.c_str(), "RECREATE", "DDG4 User file");
0153       if ( m_outFile && !m_outFile->IsZombie() )   {
0154         m_outTree = new TTree("DDG4 User Test","DDG4 data");
0155         printout(ALWAYS,"HitTupleAction","+++ Successfully opened ROOT file %s and created TTree:%s",
0156                  m_outFile->GetName(), "DDG4 User Test");
0157         if ( m_containers.size() == 1 && (m_containers[0] == "*" || m_containers[0] == "ALL") )  {
0158           m_containers.clear();
0159           for (int i = 0; i < nCol; ++i)
0160             m_containers.push_back(hce->GetHC(i)->GetName());
0161         }
0162         // Seperate loop. We need fixed addresses when creating the branches
0163         for(const auto& c : m_containers)   {
0164           m_deposits[c].first = 0;
0165           m_deposits[c].second.first = 0e0;
0166           m_deposits[c].second.second.clear();
0167         }
0168         for(const auto& c : m_containers)   {
0169           std::pair<TBranch*, Data>& e = m_deposits[c];
0170           TClass* cl = gROOT->GetClass(typeid(Data));
0171           e.first = m_outTree->Branch(c.c_str(), cl->GetName(), (void*)0);
0172           e.first->SetAutoDelete(false);
0173           printout(ALWAYS,"HitTupleAction","+++ Prepare hit branch %s in root file.",c.c_str());
0174         }
0175       }
0176       else   {
0177         except("HitTupleAction","Failed to open ROOT N-tuple file: %s",m_outFileName.c_str());
0178       }
0179     }
0180     for ( auto& e : m_deposits )  {
0181       e.second.first->SetAddress(0);
0182     }
0183     for (int i = 0; i < nCol; ++i) {
0184       G4VHitsCollection* hc = hce->GetHC(i);
0185       const string& nam = hc->GetName();
0186       if ( find(m_containers.begin(),m_containers.end(),nam) != m_containers.end() )   {
0187         Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc);
0188         if ( coll )    {
0189           std::pair<TBranch*, Data>& e = m_deposits[nam];
0190           size_t nhits = coll->GetSize();
0191           Data* d = &e.second;
0192
0193           e.second.first = 0e0;
0194           e.second.second.clear();
0195           e.first->SetAddress(&d);
0196           for ( size_t j=0; j<nhits; ++j )   {
0197             double dep = 0e0;
0198             Geant4HitData* h = coll->hit(j);
0199             Geant4Tracker::Hit* trk_hit = dynamic_cast<Geant4Tracker::Hit*>(h);
0200             if ( 0 != trk_hit )   {
0201               dep = trk_hit->truth.deposit;
0202             }
0203             else  {
0204               Geant4Calorimeter::Hit* cal_hit = dynamic_cast<Geant4Calorimeter::Hit*>(h);
0205               if ( 0 != cal_hit )
0206                 dep = cal_hit->energyDeposit;
0207               else
0208                 continue;
0209             }
0210             if ( dep > 0 )  {
0211               e.second.first += dep;
0212               e.second.second.push_back(dep);
0213             }
0214           }
0215         }
0216       }
0217     }
0218     m_outTree->Fill();
0219     return;
0220   }
0221   warning("+++ [Event:%d] The value of G4HCofThisEvent is NULL.",event->GetEventID());
0222 }
0223
0224 #include "DDG4/Factories.h"
0225 DECLARE_GEANT4ACTION_NS(myanalysis,HitTupleAction)