| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \file B5/src/EventAction.cc
28 /// \brief Implementation of the B5::EventActi
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 collecti
56 // and print warnings if not found
57 G4VHitsCollection* GetHC(const G4Event* event,
58 {
59 auto hce = event->GetHCofThisEvent();
60 if (!hce) {
61 G4ExceptionDescription msg;
62 msg << "No hits collection of this event f
63 G4Exception("EventAction::EndOfEventAction
64 return nullptr;
65 }
66
67 auto hc = hce->GetHC(collId);
68 if (!hc) {
69 G4ExceptionDescription msg;
70 msg << "Hits collection " << collId << " o
71 G4Exception("EventAction::EndOfEventAction
72 }
73 return hc;
74 }
75
76 } // namespace
77
78 namespace B5
79 {
80
81 //....oooOO0OOooo........oooOO0OOooo........oo
82
83 EventAction::EventAction()
84 {
85 // set printing per each event
86 G4RunManager::GetRunManager()->SetPrintProgr
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 void EventAction::BeginOfEventAction(const G4E
92 {
93 // Find hit collections and histogram Ids by
94 // and save them in the data members of this
95
96 if (fHodHCID[0] == -1) {
97 auto sdManager = G4SDManager::GetSDMpointe
98 auto analysisManager = G4AnalysisManager::
99
100 // hits collections names
101 array<G4String, kDim> hHCName = {{"hodosco
102 array<G4String, kDim> dHCName = {{"chamber
103 array<G4String, kDim> cHCName = {
104 {"EMcalorimeter/EMcalorimeterColl", "Had
105
106 // histograms names
107 array<array<G4String, kDim>, kDim> histoNa
108 {{{"Chamber1", "Chamber2"}}, {{"Chamber1
109
110 for (G4int iDet = 0; iDet < kDim; ++iDet)
111 // hit collections IDs
112 fHodHCID[iDet] = sdManager->GetCollectio
113 fDriftHCID[iDet] = sdManager->GetCollect
114 fCalHCID[iDet] = sdManager->GetCollectio
115 // histograms IDs
116 fDriftHistoID[kH1][iDet] = analysisManag
117 fDriftHistoID[kH2][iDet] = analysisManag
118 }
119 }
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 void EventAction::EndOfEventAction(const G4Eve
125 {
126 //
127 // Fill histograms & ntuple
128 //
129
130 // Get analysis manager
131 auto analysisManager = G4AnalysisManager::In
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]
140 // columns 0, 1
141 analysisManager->FillNtupleIColumn(iDet, n
142
143 for (unsigned long i = 0; i < nhit; ++i) {
144 auto hit = static_cast<DriftChamberHit*>
145 auto localPos = hit->GetLocalPos();
146 analysisManager->FillH2(fDriftHistoID[kH
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(
161 G4double edep = 0.;
162 // The EM and Had calorimeter hits are o
163 if (iDet == 0) {
164 auto hit = static_cast<EmCalorimeterHi
165 edep = hit->GetEdep();
166 }
167 else {
168 auto hit = static_cast<HadCalorimeterH
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 +
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()
187 auto hit = static_cast<HodoscopeHit*>(hc
188 // columns 4, 5
189 analysisManager->FillNtupleDColumn(iDet
190 }
191 }
192 analysisManager->AddNtupleRow();
193
194 //
195 // Print diagnostics
196 //
197
198 auto printModulo = G4RunManager::GetRunManag
199 if (printModulo == 0 || event->GetEventID()
200
201 auto primary = event->GetPrimaryVertex(0)->G
202 G4cout << G4endl << ">>> Event " << event->G
203 << " >>> Simulation truth : " << prim
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 << " ha
211 for (unsigned int i = 0; i < hc->GetSize()
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 <<
221 for (auto layer = 0; layer < kNofChambers;
222 for (unsigned int i = 0; i < hc->GetSize
223 auto hit = static_cast<DriftChamberHit
224 if (hit->GetLayerID() == layer) hit->P
225 }
226 }
227 }
228
229 // Calorimeters
230 array<G4String, kDim> calName = {{"EM", "Had
231 for (G4int iDet = 0; iDet < kDim; ++iDet) {
232 G4cout << calName[iDet] << " Calorimeter h
233 << " Total Edep is " << totalCalEde
234 }
235 }
236
237 //....oooOO0OOooo........oooOO0OOooo........oo
238
239 } // namespace B5
240