B5/src/EventAction.cc

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 /// \file B5/src/EventAction.cc
0028 /// \brief Implementation of the B5::EventAction class
0029
0030 #include "EventAction.hh"
0031
0032 #include "Constants.hh"
0033 #include "DriftChamberHit.hh"
0034 #include "EmCalorimeterHit.hh"
0035 #include "HadCalorimeterHit.hh"
0036 #include "HodoscopeHit.hh"
0037
0038 #include "G4AnalysisManager.hh"
0039 #include "G4Event.hh"
0040 #include "G4HCofThisEvent.hh"
0041 #include "G4ParticleDefinition.hh"
0042 #include "G4PrimaryParticle.hh"
0043 #include "G4PrimaryVertex.hh"
0044 #include "G4RunManager.hh"
0045 #include "G4SDManager.hh"
0046 #include "G4SystemOfUnits.hh"
0047 #include "G4VHitsCollection.hh"
0048
0049 using std::array;
0050 using std::vector;
0051
0052 namespace
0053 {
0054
0055 // Utility function which finds a hit collection with the given Id
0056 // and print warnings if not found
0057 G4VHitsCollection* GetHC(const G4Event* event, G4int collId)
0058 {
0059   auto hce = event->GetHCofThisEvent();
0060   if (!hce) {
0061     G4ExceptionDescription msg;
0062     msg << "No hits collection of this event found." << G4endl;
0063     G4Exception("EventAction::EndOfEventAction()", "Code001", JustWarning, msg);
0064     return nullptr;
0065   }
0066
0067   auto hc = hce->GetHC(collId);
0068   if (!hc) {
0069     G4ExceptionDescription msg;
0070     msg << "Hits collection " << collId << " of this event not found." << G4endl;
0071     G4Exception("EventAction::EndOfEventAction()", "Code001", JustWarning, msg);
0072   }
0073   return hc;
0074 }
0075
0076 }  // namespace
0077
0078 namespace B5
0079 {
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 EventAction::EventAction()
0084 {
0085   // set printing per each event
0086   G4RunManager::GetRunManager()->SetPrintProgress(1);
0087 }
0088
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0090
0091 void EventAction::BeginOfEventAction(const G4Event*)
0092 {
0093   // Find hit collections and histogram Ids by names (just once)
0094   // and save them in the data members of this class
0095
0096   if (fHodHCID[0] == -1) {
0097     auto sdManager = G4SDManager::GetSDMpointer();
0098     auto analysisManager = G4AnalysisManager::Instance();
0099
0100     // hits collections names
0101     array<G4String, kDim> hHCName = {{"hodoscope1/hodoscopeColl", "hodoscope2/hodoscopeColl"}};
0102     array<G4String, kDim> dHCName = {{"chamber1/driftChamberColl", "chamber2/driftChamberColl"}};
0103     array<G4String, kDim> cHCName = {
0104       {"EMcalorimeter/EMcalorimeterColl", "HadCalorimeter/HadCalorimeterColl"}};
0105
0106     // histograms names
0107     array<array<G4String, kDim>, kDim> histoName = {
0108       {{{"Chamber1", "Chamber2"}}, {{"Chamber1 XY", "Chamber2 XY"}}}};
0109
0110     for (G4int iDet = 0; iDet < kDim; ++iDet) {
0111       // hit collections IDs
0112       fHodHCID[iDet] = sdManager->GetCollectionID(hHCName[iDet]);
0113       fDriftHCID[iDet] = sdManager->GetCollectionID(dHCName[iDet]);
0114       fCalHCID[iDet] = sdManager->GetCollectionID(cHCName[iDet]);
0115       // histograms IDs
0116       fDriftHistoID[kH1][iDet] = analysisManager->GetH1Id(histoName[kH1][iDet]);
0117       fDriftHistoID[kH2][iDet] = analysisManager->GetH2Id(histoName[kH2][iDet]);
0118     }
0119   }
0120 }
0121
0122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0123
0124 void EventAction::EndOfEventAction(const G4Event* event)
0125 {
0126   //
0127   // Fill histograms & ntuple
0128   //
0129
0130   // Get analysis manager
0131   auto analysisManager = G4AnalysisManager::Instance();
0132
0133   // Drift chambers hits
0134   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0135     auto hc = GetHC(event, fDriftHCID[iDet]);
0136     if (!hc) return;
0137
0138     auto nhit = hc->GetSize();
0139     analysisManager->FillH1(fDriftHistoID[kH1][iDet], nhit);
0140     // columns 0, 1
0141     analysisManager->FillNtupleIColumn(iDet, nhit);
0142
0143     for (unsigned long i = 0; i < nhit; ++i) {
0144       auto hit = static_cast<DriftChamberHit*>(hc->GetHit(i));
0145       auto localPos = hit->GetLocalPos();
0146       analysisManager->FillH2(fDriftHistoID[kH2][iDet], localPos.x(), localPos.y());
0147     }
0148   }
0149
0150   // Em/Had Calorimeters hits
0151   array<G4int, kDim> totalCalHit = {{0, 0}};
0152   array<G4double, kDim> totalCalEdep = {{0., 0.}};
0153
0154   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0155     auto hc = GetHC(event, fCalHCID[iDet]);
0156     if (!hc) return;
0157
0158     totalCalHit[iDet] = 0;
0159     totalCalEdep[iDet] = 0.;
0160     for (unsigned long i = 0; i < hc->GetSize(); ++i) {
0161       G4double edep = 0.;
0162       // The EM and Had calorimeter hits are of different types
0163       if (iDet == 0) {
0164         auto hit = static_cast<EmCalorimeterHit*>(hc->GetHit(i));
0165         edep = hit->GetEdep();
0166       }
0167       else {
0168         auto hit = static_cast<HadCalorimeterHit*>(hc->GetHit(i));
0169         edep = hit->GetEdep();
0170       }
0171       if (edep > 0.) {
0172         totalCalHit[iDet]++;
0173         totalCalEdep[iDet] += edep;
0174       }
0175       fCalEdep[iDet][i] = edep;
0176     }
0177     // columns 2, 3
0178     analysisManager->FillNtupleDColumn(iDet + 2, totalCalEdep[iDet]);
0179   }
0180
0181   // Hodoscopes hits
0182   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0183     auto hc = GetHC(event, fHodHCID[iDet]);
0184     if (!hc) return;
0185
0186     for (unsigned int i = 0; i < hc->GetSize(); ++i) {
0187       auto hit = static_cast<HodoscopeHit*>(hc->GetHit(i));
0188       // columns 4, 5
0189       analysisManager->FillNtupleDColumn(iDet + 4, hit->GetTime());
0190     }
0191   }
0192   analysisManager->AddNtupleRow();
0193
0194   //
0195   // Print diagnostics
0196   //
0197
0198   auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
0199   if (printModulo == 0 || event->GetEventID() % printModulo != 0) return;
0200
0201   auto primary = event->GetPrimaryVertex(0)->GetPrimary(0);
0202   G4cout << G4endl << ">>> Event " << event->GetEventID()
0203          << " >>> Simulation truth : " << primary->GetG4code()->GetParticleName() << " "
0204          << primary->GetMomentum() << G4endl;
0205
0206   // Hodoscopes
0207   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0208     auto hc = GetHC(event, fHodHCID[iDet]);
0209     if (!hc) return;
0210     G4cout << "Hodoscope " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
0211     for (unsigned int i = 0; i < hc->GetSize(); ++i) {
0212       hc->GetHit(i)->Print();
0213     }
0214   }
0215
0216   // Drift chambers
0217   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0218     auto hc = GetHC(event, fDriftHCID[iDet]);
0219     if (!hc) return;
0220     G4cout << "Drift Chamber " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
0221     for (auto layer = 0; layer < kNofChambers; ++layer) {
0222       for (unsigned int i = 0; i < hc->GetSize(); i++) {
0223         auto hit = static_cast<DriftChamberHit*>(hc->GetHit(i));
0224         if (hit->GetLayerID() == layer) hit->Print();
0225       }
0226     }
0227   }
0228
0229   // Calorimeters
0230   array<G4String, kDim> calName = {{"EM", "Hadron"}};
0231   for (G4int iDet = 0; iDet < kDim; ++iDet) {
0232     G4cout << calName[iDet] << " Calorimeter has " << totalCalHit[iDet] << " hits."
0233            << " Total Edep is " << totalCalEdep[iDet] / MeV << " (MeV)" << G4endl;
0234   }
0235 }
0236
0237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0238
0239 }  // namespace B5