| Geant4 Cross Reference |
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 //
27 /// \file B5/src/EventAction.cc
28 /// \brief Implementation of the B5::EventAction class
29
30 #include "EventAction.hh"
31
32 #include "Constants.hh"
33 #include "DriftChamberHit.hh"
34 #include "EmCalorimeterHit.hh"
35 #include "HadCalorimeterHit.hh"
36 #include "HodoscopeHit.hh"
37
38 #include "G4AnalysisManager.hh"
39 #include "G4Event.hh"
40 #include "G4HCofThisEvent.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4PrimaryParticle.hh"
43 #include "G4PrimaryVertex.hh"
44 #include "G4RunManager.hh"
45 #include "G4SDManager.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4VHitsCollection.hh"
48
49 using std::array;
50 using std::vector;
51
52 namespace
53 {
54
55 // Utility function which finds a hit collection with the given Id
56 // and print warnings if not found
57 G4VHitsCollection* GetHC(const G4Event* event, G4int collId)
58 {
59 auto hce = event->GetHCofThisEvent();
60 if (!hce) {
61 G4ExceptionDescription msg;
62 msg << "No hits collection of this event found." << G4endl;
63 G4Exception("EventAction::EndOfEventAction()", "Code001", JustWarning, msg);
64 return nullptr;
65 }
66
67 auto hc = hce->GetHC(collId);
68 if (!hc) {
69 G4ExceptionDescription msg;
70 msg << "Hits collection " << collId << " of this event not found." << G4endl;
71 G4Exception("EventAction::EndOfEventAction()", "Code001", JustWarning, msg);
72 }
73 return hc;
74 }
75
76 } // namespace
77
78 namespace B5
79 {
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 EventAction::EventAction()
84 {
85 // set printing per each event
86 G4RunManager::GetRunManager()->SetPrintProgress(1);
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90
91 void EventAction::BeginOfEventAction(const G4Event*)
92 {
93 // Find hit collections and histogram Ids by names (just once)
94 // and save them in the data members of this class
95
96 if (fHodHCID[0] == -1) {
97 auto sdManager = G4SDManager::GetSDMpointer();
98 auto analysisManager = G4AnalysisManager::Instance();
99
100 // hits collections names
101 array<G4String, kDim> hHCName = {{"hodoscope1/hodoscopeColl", "hodoscope2/hodoscopeColl"}};
102 array<G4String, kDim> dHCName = {{"chamber1/driftChamberColl", "chamber2/driftChamberColl"}};
103 array<G4String, kDim> cHCName = {
104 {"EMcalorimeter/EMcalorimeterColl", "HadCalorimeter/HadCalorimeterColl"}};
105
106 // histograms names
107 array<array<G4String, kDim>, kDim> histoName = {
108 {{{"Chamber1", "Chamber2"}}, {{"Chamber1 XY", "Chamber2 XY"}}}};
109
110 for (G4int iDet = 0; iDet < kDim; ++iDet) {
111 // hit collections IDs
112 fHodHCID[iDet] = sdManager->GetCollectionID(hHCName[iDet]);
113 fDriftHCID[iDet] = sdManager->GetCollectionID(dHCName[iDet]);
114 fCalHCID[iDet] = sdManager->GetCollectionID(cHCName[iDet]);
115 // histograms IDs
116 fDriftHistoID[kH1][iDet] = analysisManager->GetH1Id(histoName[kH1][iDet]);
117 fDriftHistoID[kH2][iDet] = analysisManager->GetH2Id(histoName[kH2][iDet]);
118 }
119 }
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123
124 void EventAction::EndOfEventAction(const G4Event* event)
125 {
126 //
127 // Fill histograms & ntuple
128 //
129
130 // Get analysis manager
131 auto analysisManager = G4AnalysisManager::Instance();
132
133 // Drift chambers hits
134 for (G4int iDet = 0; iDet < kDim; ++iDet) {
135 auto hc = GetHC(event, fDriftHCID[iDet]);
136 if (!hc) return;
137
138 auto nhit = hc->GetSize();
139 analysisManager->FillH1(fDriftHistoID[kH1][iDet], nhit);
140 // columns 0, 1
141 analysisManager->FillNtupleIColumn(iDet, nhit);
142
143 for (unsigned long i = 0; i < nhit; ++i) {
144 auto hit = static_cast<DriftChamberHit*>(hc->GetHit(i));
145 auto localPos = hit->GetLocalPos();
146 analysisManager->FillH2(fDriftHistoID[kH2][iDet], localPos.x(), localPos.y());
147 }
148 }
149
150 // Em/Had Calorimeters hits
151 array<G4int, kDim> totalCalHit = {{0, 0}};
152 array<G4double, kDim> totalCalEdep = {{0., 0.}};
153
154 for (G4int iDet = 0; iDet < kDim; ++iDet) {
155 auto hc = GetHC(event, fCalHCID[iDet]);
156 if (!hc) return;
157
158 totalCalHit[iDet] = 0;
159 totalCalEdep[iDet] = 0.;
160 for (unsigned long i = 0; i < hc->GetSize(); ++i) {
161 G4double edep = 0.;
162 // The EM and Had calorimeter hits are of different types
163 if (iDet == 0) {
164 auto hit = static_cast<EmCalorimeterHit*>(hc->GetHit(i));
165 edep = hit->GetEdep();
166 }
167 else {
168 auto hit = static_cast<HadCalorimeterHit*>(hc->GetHit(i));
169 edep = hit->GetEdep();
170 }
171 if (edep > 0.) {
172 totalCalHit[iDet]++;
173 totalCalEdep[iDet] += edep;
174 }
175 fCalEdep[iDet][i] = edep;
176 }
177 // columns 2, 3
178 analysisManager->FillNtupleDColumn(iDet + 2, totalCalEdep[iDet]);
179 }
180
181 // Hodoscopes hits
182 for (G4int iDet = 0; iDet < kDim; ++iDet) {
183 auto hc = GetHC(event, fHodHCID[iDet]);
184 if (!hc) return;
185
186 for (unsigned int i = 0; i < hc->GetSize(); ++i) {
187 auto hit = static_cast<HodoscopeHit*>(hc->GetHit(i));
188 // columns 4, 5
189 analysisManager->FillNtupleDColumn(iDet + 4, hit->GetTime());
190 }
191 }
192 analysisManager->AddNtupleRow();
193
194 //
195 // Print diagnostics
196 //
197
198 auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
199 if (printModulo == 0 || event->GetEventID() % printModulo != 0) return;
200
201 auto primary = event->GetPrimaryVertex(0)->GetPrimary(0);
202 G4cout << G4endl << ">>> Event " << event->GetEventID()
203 << " >>> Simulation truth : " << primary->GetG4code()->GetParticleName() << " "
204 << primary->GetMomentum() << G4endl;
205
206 // Hodoscopes
207 for (G4int iDet = 0; iDet < kDim; ++iDet) {
208 auto hc = GetHC(event, fHodHCID[iDet]);
209 if (!hc) return;
210 G4cout << "Hodoscope " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
211 for (unsigned int i = 0; i < hc->GetSize(); ++i) {
212 hc->GetHit(i)->Print();
213 }
214 }
215
216 // Drift chambers
217 for (G4int iDet = 0; iDet < kDim; ++iDet) {
218 auto hc = GetHC(event, fDriftHCID[iDet]);
219 if (!hc) return;
220 G4cout << "Drift Chamber " << iDet + 1 << " has " << hc->GetSize() << " hits." << G4endl;
221 for (auto layer = 0; layer < kNofChambers; ++layer) {
222 for (unsigned int i = 0; i < hc->GetSize(); i++) {
223 auto hit = static_cast<DriftChamberHit*>(hc->GetHit(i));
224 if (hit->GetLayerID() == layer) hit->Print();
225 }
226 }
227 }
228
229 // Calorimeters
230 array<G4String, kDim> calName = {{"EM", "Hadron"}};
231 for (G4int iDet = 0; iDet < kDim; ++iDet) {
232 G4cout << calName[iDet] << " Calorimeter has " << totalCalHit[iDet] << " hits."
233 << " Total Edep is " << totalCalEdep[iDet] / MeV << " (MeV)" << G4endl;
234 }
235 }
236
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
238
239 } // namespace B5
240