| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 ////////////////////////////////////////////// 26 ///////////////////////////////////////////////////////////////////////////////
27 // File: CCalEndOfEventAction.cc 27 // File: CCalEndOfEventAction.cc
28 // Description: CCalEndOfEventAction provides 28 // Description: CCalEndOfEventAction provides User actions at end of event
29 ////////////////////////////////////////////// 29 ///////////////////////////////////////////////////////////////////////////////
30 #include "CCalEventAction.hh" 30 #include "CCalEventAction.hh"
31 #include "CCaloSD.hh" 31 #include "CCaloSD.hh"
32 #include "CCalPrimaryGeneratorAction.hh" 32 #include "CCalPrimaryGeneratorAction.hh"
33 #include "CCalSteppingAction.hh" 33 #include "CCalSteppingAction.hh"
34 #include "CCalG4HitCollection.hh" 34 #include "CCalG4HitCollection.hh"
35 #include "CCalG4Hit.hh" 35 #include "CCalG4Hit.hh"
36 #include "CCaloOrganization.hh" 36 #include "CCaloOrganization.hh"
37 #include "CCalSDList.hh" 37 #include "CCalSDList.hh"
38 38
39 #include "G4AnalysisManager.hh" 39 #include "G4AnalysisManager.hh"
40 #include "G4ios.hh" 40 #include "G4ios.hh"
41 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
42 #include "G4Event.hh" 42 #include "G4Event.hh"
43 #include "G4SDManager.hh" 43 #include "G4SDManager.hh"
44 #include "G4HCofThisEvent.hh" 44 #include "G4HCofThisEvent.hh"
45 #include "G4RunManager.hh" 45 #include "G4RunManager.hh"
46 #include "G4UserSteppingAction.hh" 46 #include "G4UserSteppingAction.hh"
47 #include "G4UImanager.hh" 47 #include "G4UImanager.hh"
48 48
49 #include <iostream> 49 #include <iostream>
50 #include <vector> 50 #include <vector>
51 #include <map> 51 #include <map>
52 52
53 //#define debug 53 //#define debug
54 //#define ddebug 54 //#define ddebug
55 55
56 56
57 CCalEventAction::CCalEventAction(CCalPrimaryGe 57 CCalEventAction::CCalEventAction(CCalPrimaryGeneratorAction* pg,
58 CCalSteppingA 58 CCalSteppingAction* sa):
59 isInitialized(false),fPrimaryGenerator(pg), 59 isInitialized(false),fPrimaryGenerator(pg),
60 fSteppingAction(sa),SDnames(nullptr),numberO 60 fSteppingAction(sa),SDnames(nullptr),numberOfSD(0)
61 { 61 {
62 #ifdef debug 62 #ifdef debug
63 G4cout << "Instantiate CCalEndOfEventAction" 63 G4cout << "Instantiate CCalEndOfEventAction" << G4endl;
64 #endif 64 #endif
65 65
66 G4cout << "Get Calorimter organisation" << G 66 G4cout << "Get Calorimter organisation" << G4endl;
67 theOrg = new CCaloOrganization; 67 theOrg = new CCaloOrganization;
68 } 68 }
69 69
70 CCalEventAction::~CCalEventAction() { 70 CCalEventAction::~CCalEventAction() {
71 delete theOrg; 71 delete theOrg;
72 delete[] SDnames; 72 delete[] SDnames;
73 } 73 }
74 74
75 void CCalEventAction::initialize() { 75 void CCalEventAction::initialize() {
76 76
77 isInitialized = true; 77 isInitialized = true;
78 numberOfSD = CCalSDList::getInstance()->getN 78 numberOfSD = CCalSDList::getInstance()->getNumberOfCaloSD();
79 #ifdef debug 79 #ifdef debug
80 G4cout << "CCalEndOfEventAction look for " < 80 G4cout << "CCalEndOfEventAction look for " << numberOfSD
81 << " calorimeter-like SD" << G4endl; 81 << " calorimeter-like SD" << G4endl;
82 #endif 82 #endif
83 if (numberOfSD > 0) { 83 if (numberOfSD > 0) {
84 G4int n = numberOfSD; 84 G4int n = numberOfSD;
85 n = std::min(n, 2); 85 n = std::min(n, 2);
86 SDnames = new nameType[n]; 86 SDnames = new nameType[n];
87 } 87 }
88 for (G4int i=0; i<numberOfSD; ++i) { 88 for (G4int i=0; i<numberOfSD; ++i) {
89 SDnames[i] = G4String(CCalSDList::getInsta 89 SDnames[i] = G4String(CCalSDList::getInstance()->getCaloSDName(i));
90 #ifdef debug 90 #ifdef debug
91 G4cout << "CCalEndOfEventAction: found SD 91 G4cout << "CCalEndOfEventAction: found SD " << i << " name "
92 << SDnames[i] << G4endl; 92 << SDnames[i] << G4endl;
93 #endif 93 #endif
94 } 94 }
95 } 95 }
96 96
97 void CCalEventAction::BeginOfEventAction(const 97 void CCalEventAction::BeginOfEventAction(const G4Event* evt) {
98 98
99 if (!isInitialized) { initialize(); } 99 if (!isInitialized) { initialize(); }
100 G4cout << " --- Begin of event: " << evt->Ge 100 G4cout << " --- Begin of event: " << evt->GetEventID() << G4endl;
101 /* 101 /*
102 if(15 == evt->GetEventID()) { 102 if(15 == evt->GetEventID()) {
103 G4UImanager * UImanager = G4UImanager::Get 103 G4UImanager * UImanager = G4UImanager::GetUIpointer();
104 UImanager->ApplyCommand("/tracking/verbose 104 UImanager->ApplyCommand("/tracking/verbose 2");
105 } 105 }
106 */ 106 */
107 } 107 }
108 108
109 void CCalEventAction::EndOfEventAction(const G 109 void CCalEventAction::EndOfEventAction(const G4Event* evt){
110 110
111 fSteppingAction->endOfEvent(); 111 fSteppingAction->endOfEvent();
112 // 112 //
113 // Look for the Hit Collection 113 // Look for the Hit Collection
114 // 114 //
115 G4HCofThisEvent* allHC = evt->GetHCofThisEve 115 G4HCofThisEvent* allHC = evt->GetHCofThisEvent();
116 if (allHC == 0) { 116 if (allHC == 0) {
117 #ifdef debug 117 #ifdef debug
118 G4cout << "CCalEndOfEventAction: No Hit Co 118 G4cout << "CCalEndOfEventAction: No Hit Collection in this event"
119 << G4endl; 119 << G4endl;
120 #endif 120 #endif
121 return; 121 return;
122 } 122 }
123 123
124 // 124 //
125 // hits info 125 // hits info
126 // 126 //
127 127
128 //Now make summary 128 //Now make summary
129 G4float hcalE[28], ecalE[49], fullE=0., edec 129 G4float hcalE[28], ecalE[49], fullE=0., edec=0, edhc=0;
130 G4int i = 0; 130 G4int i = 0;
131 for (i = 0; i < 28; i++) {hcalE[i]=0.;} 131 for (i = 0; i < 28; i++) {hcalE[i]=0.;}
132 for (i = 0; i < 49; i++) {ecalE[i]=0.;} 132 for (i = 0; i < 49; i++) {ecalE[i]=0.;}
133 133
134 G4float* edep = new G4float[numberOfSD]; 134 G4float* edep = new G4float[numberOfSD];
135 for (i = 0; i < numberOfSD; ++i){ 135 for (i = 0; i < numberOfSD; ++i){
136 136
137 // 137 //
138 // Look for the Hit Collection 138 // Look for the Hit Collection
139 // 139 //
140 edep[i] = 0; 140 edep[i] = 0;
141 G4int caloHCid = G4SDManager::GetSDMpointe 141 G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
142 142
143 if (caloHCid >= 0) { 143 if (caloHCid >= 0) {
144 CCalG4HitCollection* theHC = 144 CCalG4HitCollection* theHC =
145 (CCalG4HitCollection*) allHC->GetHC(ca 145 (CCalG4HitCollection*) allHC->GetHC(caloHCid);
146 146
147 if (theHC != 0) { 147 if (theHC != 0) {
148 148
149 G4int nentries = theHC->entries(); 149 G4int nentries = theHC->entries();
150 #ifdef debug 150 #ifdef debug
151 G4cout << " There are " << nentries << 151 G4cout << " There are " << nentries << " hits in " << SDnames[i]
152 << " :" << G4endl; 152 << " :" << G4endl;
153 #endif 153 #endif
154 154
155 if (nentries > 0) { 155 if (nentries > 0) {
156 156
157 G4int j; 157 G4int j;
158 for (j=0; j<nentries; j++){ 158 for (j=0; j<nentries; j++){
159 #ifdef ddebug 159 #ifdef ddebug
160 G4cout << "Hit " << j; 160 G4cout << "Hit " << j;
161 #endif 161 #endif
162 CCalG4Hit* aHit = (*theHC)[j]; 162 CCalG4Hit* aHit = (*theHC)[j];
163 G4float En = aHit->getEnergyDeposi 163 G4float En = aHit->getEnergyDeposit();
164 G4int unitID = aHit->getUnitID(); 164 G4int unitID = aHit->getUnitID();
165 G4int id=-1; 165 G4int id=-1;
166 if (unitID > 0 && unitID < 29) { 166 if (unitID > 0 && unitID < 29) {
167 id = unitID - 1; // HCal 167 id = unitID - 1; // HCal
168 hcalE[id] += En/GeV; 168 hcalE[id] += En/GeV;
169 } else { 169 } else {
170 G4int i0 = unitID/4096; 170 G4int i0 = unitID/4096;
171 G4int i1 = (unitID/64)%64; 171 G4int i1 = (unitID/64)%64;
172 G4int i2 = unitID%64; 172 G4int i2 = unitID%64;
173 if (i0 == 1 && i1 < 8 && i2 < 8) 173 if (i0 == 1 && i1 < 8 && i2 < 8) {
174 id = i1*7 + i2; // ECal 174 id = i1*7 + i2; // ECal
175 ecalE[id] += En/GeV; 175 ecalE[id] += En/GeV;
176 } 176 }
177 } 177 }
178 #ifdef ddebug 178 #ifdef ddebug
179 G4cout << " with Energy = " << En/ 179 G4cout << " with Energy = " << En/MeV << " MeV in Unit " << unitID
180 << " " << id << G4endl; 180 << " " << id << G4endl;
181 #endif 181 #endif
182 fullE += En/GeV; 182 fullE += En/GeV;
183 edep[i] += En/GeV; 183 edep[i] += En/GeV;
184 } 184 }
185 #ifdef ddebug 185 #ifdef ddebug
186 G4cout << " ===> Total Energy Deposi 186 G4cout << " ===> Total Energy Deposit in this Calorimeter = "
187 << edep[i]*1000.0 << " MeV " 187 << edep[i]*1000.0 << " MeV " << G4endl;
188 #endif 188 #endif
189 } 189 }
190 } 190 }
191 } 191 }
192 if (SDnames[i] == "HadronCalorimeter") { 192 if (SDnames[i] == "HadronCalorimeter") {
193 edhc = edep[i]; 193 edhc = edep[i];
194 } else if (SDnames[i] == "CrystalMatrix") 194 } else if (SDnames[i] == "CrystalMatrix") {
195 edec = edep[i]; 195 edec = edep[i];
196 } 196 }
197 } 197 }
198 198
199 delete[] edep; 199 delete[] edep;
200 200
201 G4ThreeVector pos = fPrimaryGenerator->GetPa 201 G4ThreeVector pos = fPrimaryGenerator->GetParticlePosition();
202 G4float ener = fPrimaryGenerator->GetParticl 202 G4float ener = fPrimaryGenerator->GetParticleEnergy()/GeV;
203 G4float x = pos.x()/mm; 203 G4float x = pos.x()/mm;
204 G4float y = pos.y()/mm; 204 G4float y = pos.y()/mm;
205 G4float z = pos.z()/mm; 205 G4float z = pos.z()/mm;
206 206
207 //Save results 207 //Save results
208 G4AnalysisManager* man = G4AnalysisManager:: 208 G4AnalysisManager* man = G4AnalysisManager::Instance();
209 //1) 209 //1)
210 static G4int IDenergy = -1; 210 static G4int IDenergy = -1;
211 if (IDenergy < 0) 211 if (IDenergy < 0)
212 IDenergy = man->GetH1Id("h4000"); 212 IDenergy = man->GetH1Id("h4000");
213 man->FillH1(IDenergy,fullE); 213 man->FillH1(IDenergy,fullE);
214 //2) 214 //2)
215 #ifdef debug 215 #ifdef debug
216 G4double totalFilledEnergyHcal = 0.0; 216 G4double totalFilledEnergyHcal = 0.0;
217 #endif 217 #endif
218 static G4int IDhcalE = -1; 218 static G4int IDhcalE = -1;
219 if (IDhcalE < 0) 219 if (IDhcalE < 0)
220 IDhcalE = man->GetH1Id("h100"); 220 IDhcalE = man->GetH1Id("h100");
221 for (G4int j=0; j<28; j++) { 221 for (G4int j=0; j<28; j++) {
222 man->FillH1(IDhcalE+j,hcalE[j]); 222 man->FillH1(IDhcalE+j,hcalE[j]);
223 #ifdef debug 223 #ifdef debug
224 G4cout << "Fill Hcal histo " << j << " wit 224 G4cout << "Fill Hcal histo " << j << " with " << hcalE[j] << G4endl;
225 totalFilledEnergyHcal += hcalE[j]; 225 totalFilledEnergyHcal += hcalE[j];
226 #endif 226 #endif
227 } 227 }
228 #ifdef debug 228 #ifdef debug
229 G4cout << "CCalAnalysis::InsertEnergyHcal: T 229 G4cout << "CCalAnalysis::InsertEnergyHcal: Total filled Energy Hcal histo "
230 << totalFilledEnergyHcal << G4endl; 230 << totalFilledEnergyHcal << G4endl;
231 totalFilledEnergyEcal = 0.0; 231 totalFilledEnergyEcal = 0.0;
232 #endif 232 #endif
233 233
234 //3) 234 //3)
235 static G4int IDecalE = -1; 235 static G4int IDecalE = -1;
236 if (IDecalE < 0) 236 if (IDecalE < 0)
237 IDecalE = man->GetH1Id("h200"); 237 IDecalE = man->GetH1Id("h200");
238 for (G4int j=0; j<49; j++) { 238 for (G4int j=0; j<49; j++) {
239 man->FillH1(IDecalE+j,ecalE[j]); 239 man->FillH1(IDecalE+j,ecalE[j]);
240 #ifdef debug 240 #ifdef debug
241 G4cout << "Fill Ecal histo " << j << " wit 241 G4cout << "Fill Ecal histo " << j << " with " << ecalE[j] << G4endl;
242 totalFilledEnergyEcal += ecalE[j]; 242 totalFilledEnergyEcal += ecalE[j];
243 #endif 243 #endif
244 } 244 }
245 #ifdef debug 245 #ifdef debug
246 G4cout << "CCalAnalysis::InsertEnergyEal: To 246 G4cout << "CCalAnalysis::InsertEnergyEal: Total filled Energy Ecal histo "
247 << totalFilledEnergyEcal << G4endl; 247 << totalFilledEnergyEcal << G4endl;
248 #endif 248 #endif
249 // 4) 249 // 4)
250 G4int counter=0; 250 G4int counter=0;
251 for (G4int j=0; j<28; j++) 251 for (G4int j=0; j<28; j++)
252 { 252 {
253 man->FillNtupleFColumn(counter,hcalE[j]) 253 man->FillNtupleFColumn(counter,hcalE[j]);
254 counter++; 254 counter++;
255 } 255 }
256 for (G4int j=0; j<49; j++) 256 for (G4int j=0; j<49; j++)
257 { 257 {
258 man->FillNtupleFColumn(counter,ecalE[j]) 258 man->FillNtupleFColumn(counter,ecalE[j]);
259 counter++; 259 counter++;
260 } 260 }
261 man->FillNtupleFColumn(counter,ener); 261 man->FillNtupleFColumn(counter,ener);
262 man->FillNtupleFColumn(counter+1,x); 262 man->FillNtupleFColumn(counter+1,x);
263 man->FillNtupleFColumn(counter+2,y); 263 man->FillNtupleFColumn(counter+2,y);
264 man->FillNtupleFColumn(counter+3,z); 264 man->FillNtupleFColumn(counter+3,z);
265 man->FillNtupleFColumn(counter+4,fullE); 265 man->FillNtupleFColumn(counter+4,fullE);
266 man->FillNtupleFColumn(counter+5,edec); 266 man->FillNtupleFColumn(counter+5,edec);
267 man->FillNtupleFColumn(counter+6,edhc); 267 man->FillNtupleFColumn(counter+6,edhc);
268 man->AddNtupleRow(); 268 man->AddNtupleRow();
269 #ifdef debug 269 #ifdef debug
270 G4cout << "CCalAnalysis:: Fill Ntuple " << G 270 G4cout << "CCalAnalysis:: Fill Ntuple " << G4endl;
271 #endif 271 #endif
272 272
273 // 5) 273 // 5)
274 static G4int IDtimeProfile = -1; 274 static G4int IDtimeProfile = -1;
275 if (IDtimeProfile < 0) 275 if (IDtimeProfile < 0)
276 IDtimeProfile = man->GetH1Id("h901"); 276 IDtimeProfile = man->GetH1Id("h901");
277 for (i = 0; i < numberOfSD; i++){ 277 for (i = 0; i < numberOfSD; i++){
278 G4int caloHCid = G4SDManager::GetSDMpointe 278 G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
279 if (caloHCid >= 0) { 279 if (caloHCid >= 0) {
280 CCalG4HitCollection* theHC = 280 CCalG4HitCollection* theHC =
281 (CCalG4HitCollection*) allHC->GetHC(ca 281 (CCalG4HitCollection*) allHC->GetHC(caloHCid);
282 if (theHC != 0) { 282 if (theHC != 0) {
283 G4int nentries = theHC->entries(); 283 G4int nentries = theHC->entries();
284 if (nentries > 0) { 284 if (nentries > 0) {
285 for (G4int k=0; k<nentries; k++) { 285 for (G4int k=0; k<nentries; k++) {
286 CCalG4Hit* aHit = (*theHC)[k]; 286 CCalG4Hit* aHit = (*theHC)[k];
287 man->FillH1(IDtimeProfile,aHit->ge 287 man->FillH1(IDtimeProfile,aHit->getTimeSlice(),aHit->getEnergyDeposit()/GeV);
288 288
289 #ifdef debug 289 #ifdef debug
290 G4cout << "CCalAnalysis:: Fill Tim 290 G4cout << "CCalAnalysis:: Fill Time Profile with Hit " << k
291 << " Edeposit " << edep << 291 << " Edeposit " << edep << " Gev" << G4endl;
292 #endif 292 #endif
293 } 293 }
294 } 294 }
295 } 295 }
296 } 296 }
297 } 297 }
298 } 298 }
299 299
300 300
301 301