Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-31 09:21:47

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // ********************************************************************
0027 //
0028 //  CaTS (Calorimetry and Tracking Simulation)
0029 //
0030 //  Authors : Hans Wenzel
0031 //            Soon Yung Jun
0032 //            (Fermi National Accelerator Laboratory)
0033 //
0034 // History
0035 //   October 18th, 2021 : first implementation
0036 //
0037 // ********************************************************************
0038 //
0039 /// \file EventAction.cc
0040 /// \brief Implementation of the CaTS::EventAction class
0041 
0042 // Geant4 headers
0043 #include <G4UserEventAction.hh>
0044 #include <G4ios.hh>
0045 #include "G4Event.hh"
0046 #include "G4HCofThisEvent.hh"
0047 #include <G4String.hh>
0048 #include <G4Types.hh>
0049 #include <G4VHit.hh>
0050 #include <G4VHitsCollection.hh>
0051 #include "G4SDManager.hh"
0052 
0053 // project headers:
0054 #include "EventAction.hh"
0055 #include "ConfigurationManager.hh"
0056 #include "Event.hh"
0057 #include "PhotonSD.hh"
0058 #include "PhotonHit.hh"
0059 #include "InteractionHit.hh"
0060 #include "lArTPCHit.hh"
0061 #include "TrackerHit.hh"
0062 #include "MscHit.hh"
0063 #include "CalorimeterHit.hh"
0064 #include "DRCalorimeterHit.hh"
0065 #ifdef WITH_ROOT
0066 #include "RootIO.hh"
0067 #endif
0068 // stl headers
0069 #include <map>
0070 #include <utility>
0071 #include <algorithm>
0072 #include <istream>
0073 #ifdef WITH_G4OPTICKS
0074 #include "OpticksFlags.hh"
0075 #include "G4Opticks.hh"
0076 #include "G4OpticksHit.hh"
0077 #endif
0078 
0079 #ifdef WITH_ROOT
0080 namespace {
0081     // Mutex to lock updating the global ion map
0082     G4Mutex ionIdMapMutex = G4MUTEX_INITIALIZER;
0083 } // namespace
0084 #endif
0085 
0086 EventAction::EventAction()
0087 : G4UserEventAction() {
0088 #ifdef WITH_ROOT
0089     RootIO::GetInstance();
0090 #endif
0091 }
0092 
0093 void EventAction::BeginOfEventAction(const G4Event*) {
0094 }
0095 
0096 void EventAction::EndOfEventAction(const G4Event* event) {
0097     G4bool verbose = ConfigurationManager::getInstance()->isEnable_verbose();
0098     if (verbose)
0099         G4cout << "EventAction::EndOfEventAction Event:   " << event->GetEventID()
0100         << G4endl;
0101     G4HCofThisEvent* HCE = event->GetHCofThisEvent();
0102     if (HCE == nullptr)
0103         return;
0104 #ifdef WITH_ROOT
0105     Event* CaTSEvt = new Event();
0106     CaTSEvt->SetEventNr(event->GetEventID());
0107     std::map<G4String, std::vector < G4VHit*>>*hcmap = CaTSEvt->GetHCMap();
0108 #endif  // end WITH_ROOT
0109 #ifdef WITH_G4OPTICKS
0110     if (ConfigurationManager::getInstance()->isEnable_opticks()) {
0111         G4Opticks* g4ok = G4Opticks::Get();
0112         G4int eventid = event->GetEventID();
0113         g4ok->propagateOpticalPhotons(eventid);
0114         unsigned num_hits = g4ok->getNumHit();
0115         G4cout << "EndOfEventAction: num_hits: " << num_hits << G4endl;
0116         if (num_hits > 0) {
0117             G4HCtable* hctable = G4SDManager::GetSDMpointer()->GetHCtable();
0118             for (G4int i = 0; i < hctable->entries(); ++i) {
0119                 std::string sdn = hctable->GetSDname(i);
0120                 std::size_t found = sdn.find("Photondetector");
0121                 if (found != std::string::npos) {
0122                     PhotonSD* aSD =
0123                             (PhotonSD*) G4SDManager::GetSDMpointer()->FindSensitiveDetector(
0124                             sdn);
0125                     aSD->AddOpticksHits();
0126                 }
0127             }
0128         }
0129         if (verbose) {
0130             G4cout << "***********************************************************"
0131                     "********************************************************"
0132                     << G4endl;
0133             G4cout << " EndOfEventAction: numphotons:   " << g4ok->getNumPhotons()
0134                     << " Gensteps: " << g4ok->getNumGensteps()
0135                     << "  Maxgensteps:  " << g4ok->getMaxGensteps() << G4endl;
0136             G4cout << " EndOfEventAction: num_hits: " << g4ok->getNumHit() << G4endl;
0137             G4cout << g4ok->dbgdesc() << G4endl;
0138         }
0139         g4ok->reset();
0140         if (verbose) {
0141             G4cout << "========================== After reset: " << G4endl;
0142             G4cout << " EndOfEventAction: numphotons:   " << g4ok->getNumPhotons()
0143                     << " Gensteps: " << g4ok->getNumGensteps()
0144                     << "  Maxgensteps:  " << g4ok->getMaxGensteps() << G4endl;
0145             G4cout << "EndOfEventAction: num_hits: " << g4ok->getNumHit() << G4endl;
0146             G4cout << g4ok->dbgdesc() << G4endl;
0147             G4cout << "***********************************************************"
0148                     "********************************************************"
0149                     << G4endl;
0150         }
0151     } // end isEnable_opticks
0152 #endif  // end   WITH_G4OPTICKS
0153     //
0154     // Now we deal with the Geant4 Hit collections.
0155     //
0156     if (verbose)
0157         G4cout << "Number of collections:  " << HCE->GetNumberOfCollections()
0158         << G4endl;
0159 #ifdef WITH_ROOT
0160     if (ConfigurationManager::getInstance()->isWriteHits()) {
0161         for (int i = 0; i < HCE->GetNumberOfCollections(); i++) {
0162             G4VHitsCollection* hc = HCE->GetHC(i);
0163             G4String hcname = hc->GetName();
0164             std::vector<std::string> y = split(hcname, '_');
0165             std::string Classname = y[1];
0166             if (verbose)
0167                 G4cout << "Classname: " << Classname << G4endl;
0168             if (Classname == "lArTPC") {
0169                 std::vector<G4VHit*> hitsVector;
0170                 G4int NbHits = hc->GetSize();
0171                 for (G4int ii = 0; ii < NbHits; ii++) {
0172                     G4VHit* hit = hc->GetHit(ii);
0173                     lArTPCHit* tpcHit = dynamic_cast<lArTPCHit*> (hit);
0174                     hitsVector.push_back(tpcHit);
0175                 }
0176                 hcmap->insert(std::make_pair(hcname, hitsVector));
0177             } else if (Classname == "Photondetector") {
0178                 std::vector<G4VHit*> hitsVector;
0179                 G4int NbHits = hc->GetSize();
0180                 if (verbose)
0181                     G4cout << "Photondetector size: " << hc->GetSize() << G4endl;
0182                 for (G4int ii = 0; ii < NbHits; ii++) {
0183                     G4VHit* hit = hc->GetHit(ii);
0184                     PhotonHit* pHit = dynamic_cast<PhotonHit*> (hit);
0185                     hitsVector.push_back(pHit);
0186                 }
0187                 hcmap->insert(std::make_pair(hcname, hitsVector));
0188             } else if (Classname == "Target") {
0189                 std::vector<G4VHit*> hitsVector;
0190                 G4int NbHits = hc->GetSize();
0191                 if (verbose)
0192                     G4cout << "Interaction size: " << hc->GetSize() << G4endl;
0193                 for (G4int ii = 0; ii < NbHits; ii++) {
0194                     G4VHit* hit = hc->GetHit(ii);
0195                     InteractionHit* iaHit = dynamic_cast<InteractionHit*> (hit);
0196                     hitsVector.push_back(iaHit);
0197                 }
0198                 hcmap->insert(std::make_pair(hcname, hitsVector));
0199             } else if (Classname == "Tracker") {
0200                 std::vector<G4VHit*> hitsVector;
0201                 G4int NbHits = hc->GetSize();
0202                 if (verbose)
0203                     G4cout << "Tracker size: " << hc->GetSize() << G4endl;
0204                 for (G4int ii = 0; ii < NbHits; ii++) {
0205                     G4VHit* hit = hc->GetHit(ii);
0206                     TrackerHit* tHit = dynamic_cast<TrackerHit*> (hit);
0207                     hitsVector.push_back(tHit);
0208                 }
0209                 hcmap->insert(std::make_pair(hcname, hitsVector));
0210             } else if (Classname == "Msc") {
0211                 std::vector<G4VHit*> hitsVector;
0212                 G4int NbHits = hc->GetSize();
0213                 if (verbose)
0214                     G4cout << "Msc size: " << hc->GetSize() << G4endl;
0215                 for (G4int ii = 0; ii < NbHits; ii++) {
0216                     G4VHit* hit = hc->GetHit(ii);
0217                     MscHit* mscHit = dynamic_cast<MscHit*> (hit);
0218                     hitsVector.push_back(mscHit);
0219                 }
0220                 hcmap->insert(std::make_pair(hcname, hitsVector));
0221             } else if (Classname == "Calorimeter") {
0222                 std::vector<G4VHit*> hitsVector;
0223                 G4int NbHits = hc->GetSize();
0224                 if (verbose)
0225                     G4cout << "Calorimeter size: " << hc->GetSize() << G4endl;
0226                 for (G4int ii = 0; ii < NbHits; ii++) {
0227                     G4VHit* hit = hc->GetHit(ii);
0228                     CalorimeterHit* cHit = dynamic_cast<CalorimeterHit*> (hit);
0229                     hitsVector.push_back(cHit);
0230                 }
0231                 hcmap->insert(std::make_pair(hcname, hitsVector));
0232             } else if (Classname == "DRCalorimeter") {
0233                 std::vector<G4VHit*> hitsVector;
0234                 G4int NbHits = hc->GetSize();
0235                 if (verbose)
0236                     G4cout << "DRCalorimeter size: " << hc->GetSize() << G4endl;
0237                 for (G4int ii = 0; ii < NbHits; ii++) {
0238                     G4VHit* hit = hc->GetHit(ii);
0239                     DRCalorimeterHit* drHit = dynamic_cast<DRCalorimeterHit*> (hit);
0240                     hitsVector.push_back(drHit);
0241                 }
0242                 hcmap->insert(std::make_pair(hcname, hitsVector));
0243             } else {
0244                 G4cout << "SD type: " << Classname << " unknown" << G4endl;
0245             }
0246         }
0247         G4AutoLock lock(&ionIdMapMutex);
0248         RootIO::GetInstance()->Write(CaTSEvt);
0249         CaTSEvt->Reset();
0250         delete CaTSEvt;
0251     } // end enableio
0252 #endif
0253 }
0254 
0255 std::vector<std::string>& EventAction::split(const std::string& s, char delim,
0256         std::vector<std::string>& elems) {
0257     std::stringstream ss(s);
0258     std::string item;
0259     while (std::getline(ss, item, delim)) {
0260         elems.push_back(item);
0261     }
0262     return elems;
0263 }
0264 
0265 std::vector<std::string> EventAction::split(const std::string& s, char delim) {
0266     std::vector<std::string> elems;
0267     return split(s, delim, elems);
0268 }