ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
B5EventAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file B5EventAction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
29 
30 #include "B5EventAction.hh"
31 #include "B5HodoscopeHit.hh"
32 #include "B5DriftChamberHit.hh"
33 #include "B5EmCalorimeterHit.hh"
34 #include "B5HadCalorimeterHit.hh"
35 #include "B5Constants.hh"
36 
37 #include "G4Event.hh"
38 #include "G4RunManager.hh"
39 #include "G4EventManager.hh"
40 #include "G4HCofThisEvent.hh"
41 #include "G4VHitsCollection.hh"
42 #include "G4SDManager.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4ios.hh"
45 #include "g4analysis.hh"
46 
47 using std::array;
48 using std::vector;
49 
50 
51 namespace {
52 
53 // Utility function which finds a hit collection with the given Id
54 // and print warnings if not found
55 G4VHitsCollection* GetHC(const G4Event* event, G4int collId) {
56  auto hce = event->GetHCofThisEvent();
57  if (!hce) {
59  msg << "No hits collection of this event found." << G4endl;
60  G4Exception("B5EventAction::EndOfEventAction()",
61  "B5Code001", JustWarning, msg);
62  return nullptr;
63  }
64 
65  auto hc = hce->GetHC(collId);
66  if ( ! hc) {
68  msg << "Hits collection " << collId << " of this event not found." << G4endl;
69  G4Exception("B5EventAction::EndOfEventAction()",
70  "B5Code001", JustWarning, msg);
71  }
72  return hc;
73 }
74 
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
81  fHodHCID {{ -1, -1 }},
82  fDriftHCID{{ -1, -1 }},
83  fCalHCID {{ -1, -1 }},
84  fDriftHistoID{{ {{ -1, -1 }}, {{ -1, -1 }} }},
85  fCalEdep{{ vector<G4double>(kNofEmCells, 0.), vector<G4double>(kNofHadCells, 0.) }}
86  // std::array<T, N> is an aggregate that contains a C array.
87  // To initialize it, we need outer braces for the class itself
88  // and inner braces for the C array
89 {
90  // set printing per each event
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
97 {}
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
102 {
103  // Find hit collections and histogram Ids by names (just once)
104  // and save them in the data members of this class
105 
106  if (fHodHCID[0] == -1) {
107  auto sdManager = G4SDManager::GetSDMpointer();
108  auto analysisManager = G4AnalysisManager::Instance();
109 
110  // hits collections names
111  array<G4String, kDim> hHCName
112  = {{ "hodoscope1/hodoscopeColl", "hodoscope2/hodoscopeColl" }};
113  array<G4String, kDim> dHCName
114  = {{ "chamber1/driftChamberColl", "chamber2/driftChamberColl" }};
115  array<G4String, kDim> cHCName
116  = {{ "EMcalorimeter/EMcalorimeterColl", "HadCalorimeter/HadCalorimeterColl" }};
117 
118  // histograms names
119  array<array<G4String, kDim>, kDim> histoName
120  = {{ {{ "Chamber1", "Chamber2" }}, {{ "Chamber1 XY", "Chamber2 XY" }} }};
121 
122  for (G4int iDet = 0; iDet < kDim; ++iDet) {
123  // hit collections IDs
124  fHodHCID[iDet] = sdManager->GetCollectionID(hHCName[iDet]);
125  fDriftHCID[iDet] = sdManager->GetCollectionID(dHCName[iDet]);
126  fCalHCID[iDet] = sdManager->GetCollectionID(cHCName[iDet]);
127  // histograms IDs
128  fDriftHistoID[kH1][iDet] = analysisManager->GetH1Id(histoName[kH1][iDet]);
129  fDriftHistoID[kH2][iDet] = analysisManager->GetH2Id(histoName[kH2][iDet]);
130  }
131  }
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
136 void B5EventAction::EndOfEventAction(const G4Event* event)
137 {
138  //
139  // Fill histograms & ntuple
140  //
141 
142  // Get analysis manager
143  auto analysisManager = G4AnalysisManager::Instance();
144 
145  // Drift chambers hits
146  for (G4int iDet = 0; iDet < kDim; ++iDet) {
147  auto hc = GetHC(event, fDriftHCID[iDet]);
148  if ( ! hc ) return;
149 
150  auto nhit = hc->GetSize();
151  analysisManager->FillH1(fDriftHistoID[kH1][iDet], nhit );
152  // columns 0, 1
153  analysisManager->FillNtupleIColumn(iDet, nhit);
154 
155  for (unsigned long i = 0; i < nhit; ++i) {
156  auto hit = static_cast<B5DriftChamberHit*>(hc->GetHit(i));
157  auto localPos = hit->GetLocalPos();
158  analysisManager->FillH2(fDriftHistoID[kH2][iDet], localPos.x(), localPos.y());
159  }
160  }
161 
162  // Em/Had Calorimeters hits
163  array<G4int, kDim> totalCalHit = {{ 0, 0 }};
164  array<G4double, kDim> totalCalEdep = {{ 0., 0. }};
165 
166  for (G4int iDet = 0; iDet < kDim; ++iDet) {
167  auto hc = GetHC(event, fCalHCID[iDet]);
168  if ( ! hc ) return;
169 
170  totalCalHit[iDet] = 0;
171  totalCalEdep[iDet] = 0.;
172  for (unsigned long i = 0; i < hc->GetSize(); ++i) {
173  G4double edep = 0.;
174  // The EM and Had calorimeter hits are of different types
175  if (iDet == 0) {
176  auto hit = static_cast<B5EmCalorimeterHit*>(hc->GetHit(i));
177  edep = hit->GetEdep();
178  } else {
179  auto hit = static_cast<B5HadCalorimeterHit*>(hc->GetHit(i));
180  edep = hit->GetEdep();
181  }
182  if ( edep > 0. ) {
183  totalCalHit[iDet]++;
184  totalCalEdep[iDet] += edep;
185  }
186  fCalEdep[iDet][i] = edep;
187  }
188  // columns 2, 3
189  analysisManager->FillNtupleDColumn(iDet + 2, totalCalEdep[iDet]);
190  }
191 
192  // Hodoscopes hits
193  for (G4int iDet = 0; iDet < kDim; ++iDet) {
194  auto hc = GetHC(event, fHodHCID[iDet]);
195  if ( ! hc ) return;
196 
197  for (unsigned int i = 0; i<hc->GetSize(); ++i) {
198  auto hit = static_cast<B5HodoscopeHit*>(hc->GetHit(i));
199  // columns 4, 5
200  analysisManager->FillNtupleDColumn(iDet + 4, hit->GetTime());
201  }
202  }
203  analysisManager->AddNtupleRow();
204 
205  //
206  // Print diagnostics
207  //
208 
209  auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
210  if ( printModulo == 0 || event->GetEventID() % printModulo != 0) return;
211 
212  auto primary = event->GetPrimaryVertex(0)->GetPrimary(0);
213  G4cout
214  << G4endl
215  << ">>> Event " << event->GetEventID() << " >>> Simulation truth : "
216  << primary->GetG4code()->GetParticleName()
217  << " " << primary->GetMomentum() << G4endl;
218 
219  // Hodoscopes
220  for (G4int iDet = 0; iDet < kDim; ++iDet) {
221  auto hc = GetHC(event, fHodHCID[iDet]);
222  if ( ! hc ) return;
223  G4cout << "Hodoscope " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
224  for (unsigned int i = 0; i<hc->GetSize(); ++i) {
225  hc->GetHit(i)->Print();
226  }
227  }
228 
229  // Drift chambers
230  for (G4int iDet = 0; iDet < kDim; ++iDet) {
231  auto hc = GetHC(event, fDriftHCID[iDet]);
232  if ( ! hc ) return;
233  G4cout << "Drift Chamber " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
234  for (auto layer = 0; layer < kNofChambers; ++layer) {
235  for (unsigned int i = 0; i < hc->GetSize(); i++) {
236  auto hit = static_cast<B5DriftChamberHit*>(hc->GetHit(i));
237  if (hit->GetLayerID() == layer) hit->Print();
238  }
239  }
240  }
241 
242  // Calorimeters
243  array<G4String, kDim> calName = {{ "EM", "Hadron" }};
244  for (G4int iDet = 0; iDet < kDim; ++iDet) {
245  G4cout << calName[iDet] << " Calorimeter has " << totalCalHit[iDet] << " hits."
246  << " Total Edep is " << totalCalEdep[iDet]/MeV << " (MeV)" << G4endl;
247  }
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......