ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CCalEventAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CCalEventAction.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 //
27 // File: CCalEndOfEventAction.cc
28 // Description: CCalEndOfEventAction provides User actions at end of event
30 #include "CCalEventAction.hh"
31 #include "CCaloSD.hh"
32 #include "CCalPrimaryGeneratorAction.hh"
33 #include "CCalSteppingAction.hh"
34 #include "CCalG4HitCollection.hh"
35 #include "CCalG4Hit.hh"
36 #include "CCaloOrganization.hh"
37 #include "CCalSDList.hh"
38 #include "CCalAnalysis.hh"
39 
40 
41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Event.hh"
44 #include "G4SDManager.hh"
45 #include "G4HCofThisEvent.hh"
46 #include "G4RunManager.hh"
47 #include "G4UserSteppingAction.hh"
48 #include "G4UImanager.hh"
49 
50 #include <iostream>
51 #include <vector>
52 #include <map>
53 
54 #include "CCalAnalysis.hh"
55 
56 //#define debug
57 //#define ddebug
58 
59 
60 CCalEventAction::CCalEventAction(CCalPrimaryGeneratorAction* pg,
61  CCalSteppingAction* sa):
62  isInitialized(false),fPrimaryGenerator(pg),
63  fSteppingAction(sa),SDnames(nullptr),numberOfSD(0)
64 {
65 #ifdef debug
66  G4cout << "Instantiate CCalEndOfEventAction" << G4endl;
67 #endif
68 
69  G4cout << "Get Calorimter organisation" << G4endl;
70  theOrg = new CCaloOrganization;
71 }
72 
73 CCalEventAction::~CCalEventAction() {
74  delete theOrg;
75  delete[] SDnames;
76 }
77 
78 void CCalEventAction::initialize() {
79 
80  isInitialized = true;
81  numberOfSD = CCalSDList::getInstance()->getNumberOfCaloSD();
82 #ifdef debug
83  G4cout << "CCalEndOfEventAction look for " << numberOfSD
84  << " calorimeter-like SD" << G4endl;
85 #endif
86  if (numberOfSD > 0) {
87  G4int n = numberOfSD;
88  n = std::min(n, 2);
89  SDnames = new nameType[n];
90  }
91  for (G4int i=0; i<numberOfSD; ++i) {
92  SDnames[i] = G4String(CCalSDList::getInstance()->getCaloSDName(i));
93 #ifdef debug
94  G4cout << "CCalEndOfEventAction: found SD " << i << " name "
95  << SDnames[i] << G4endl;
96 #endif
97  }
98 }
99 
100 void CCalEventAction::BeginOfEventAction(const G4Event* evt) {
101 
102  if (!isInitialized) { initialize(); }
103  G4cout << " --- Begin of event: " << evt->GetEventID() << G4endl;
104  /*
105  if(15 == evt->GetEventID()) {
106  G4UImanager * UImanager = G4UImanager::GetUIpointer();
107  UImanager->ApplyCommand("/tracking/verbose 2");
108  }
109  */
110 }
111 
112 void CCalEventAction::EndOfEventAction(const G4Event* evt){
113 
114  fSteppingAction->endOfEvent();
115  //
116  // Look for the Hit Collection
117  //
118  G4HCofThisEvent* allHC = evt->GetHCofThisEvent();
119  if (allHC == 0) {
120 #ifdef debug
121  G4cout << "CCalEndOfEventAction: No Hit Collection in this event"
122  << G4endl;
123 #endif
124  return;
125  }
126 
127  //
128  // hits info
129  //
130 
131  //Now make summary
132  G4float hcalE[28], ecalE[49], fullE=0., edec=0, edhc=0;
133  G4int i = 0;
134  for (i = 0; i < 28; i++) {hcalE[i]=0.;}
135  for (i = 0; i < 49; i++) {ecalE[i]=0.;}
136 
137  G4float* edep = new G4float[numberOfSD];
138  G4int nhit=0;
139  for (i = 0; i < numberOfSD; ++i){
140 
141  //
142  // Look for the Hit Collection
143  //
144  edep[i] = 0;
145  G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
146 
147  if (caloHCid >= 0) {
148  CCalG4HitCollection* theHC =
149  (CCalG4HitCollection*) allHC->GetHC(caloHCid);
150 
151  if (theHC != 0) {
152 
153  G4int nentries = theHC->entries();
154 #ifdef debug
155  G4cout << " There are " << nentries << " hits in " << SDnames[i]
156  << " :" << G4endl;
157 #endif
158 
159  if (nentries > 0) {
160 
161  G4int j;
162  for (j=0; j<nentries; j++){
163 #ifdef ddebug
164  G4cout << "Hit " << j;
165 #endif
166  CCalG4Hit* aHit = (*theHC)[j];
167  G4float En = aHit->getEnergyDeposit();
168  G4int unitID = aHit->getUnitID();
169  G4int id=-1;
170  if (unitID > 0 && unitID < 29) {
171  id = unitID - 1; // HCal
172  hcalE[id] += En/GeV;
173  } else {
174  G4int i0 = unitID/4096;
175  G4int i1 = (unitID/64)%64;
176  G4int i2 = unitID%64;
177  if (i0 == 1 && i1 < 8 && i2 < 8) {
178  id = i1*7 + i2; // ECal
179  ecalE[id] += En/GeV;
180  }
181  }
182 #ifdef ddebug
183  G4cout << " with Energy = " << En/MeV << " MeV in Unit " << unitID
184  << " " << id << G4endl;
185 #endif
186  fullE += En/GeV;
187  edep[i] += En/GeV;
188  nhit++;
189  }
190 #ifdef ddebug
191  G4cout << " ===> Total Energy Deposit in this Calorimeter = "
192  << edep[i]*1000.0 << " MeV " << G4endl;
193 #endif
194  }
195  }
196  }
197  if (SDnames[i] == "HadronCalorimeter") {
198  edhc = edep[i];
199  } else if (SDnames[i] == "CrystalMatrix") {
200  edec = edep[i];
201  }
202  }
203 
204  delete[] edep;
205 
206  G4ThreeVector pos = fPrimaryGenerator->GetParticlePosition();
207  G4float ener = fPrimaryGenerator->GetParticleEnergy()/GeV;
208  G4float x = pos.x()/mm;
209  G4float y = pos.y()/mm;
210  G4float z = pos.z()/mm;
211 
212  //Save results
213  G4AnalysisManager* man = G4AnalysisManager::Instance();
214  //1)
215  static G4int IDenergy = -1;
216  if (IDenergy < 0)
217  IDenergy = man->GetH1Id("h4000");
218  man->FillH1(IDenergy,fullE);
219  //2)
220  G4double totalFilledEnergyHcal = 0.0;
221  static G4int IDhcalE = -1;
222  if (IDhcalE < 0)
223  IDhcalE = man->GetH1Id("h100");
224  for (G4int j=0; j<28; j++) {
225  man->FillH1(IDhcalE+j,hcalE[j]);
226 #ifdef debug
227  G4cout << "Fill Hcal histo " << j << " with " << hcalE[j] << G4endl;
228 #endif
229  totalFilledEnergyHcal += hcalE[j];
230  }
231 #ifdef debug
232  G4cout <<
233  "CCalAnalysis::InsertEnergyHcal: Total filled Energy Hcal histo "
234  << totalFilledEnergyHcal << G4endl;
235 #endif
236 
237  //3)
238  static G4int IDecalE = -1;
239  if (IDecalE < 0)
240  IDecalE = man->GetH1Id("h200");
241  G4double totalFilledEnergyEcal = 0.0;
242  for (G4int j=0; j<49; j++) {
243  man->FillH1(IDecalE+j,ecalE[j]);
244 #ifdef debug
245  G4cout << "Fill Ecal histo " << j << " with " << ecalE[j] << G4endl;
246 #endif
247  totalFilledEnergyEcal += ecalE[j];
248  }
249 #ifdef debug
250  G4cout <<
251  "CCalAnalysis::InsertEnergyEal: Total filled Energy Ecal histo "
252  << totalFilledEnergyEcal << G4endl;
253 #endif
254  // 4)
255  G4int counter=0;
256  for (G4int j=0; j<28; j++)
257  {
258  man->FillNtupleFColumn(counter,hcalE[j]);
259  counter++;
260  }
261  for (G4int j=0; j<49; j++)
262  {
263  man->FillNtupleFColumn(counter,ecalE[j]);
264  counter++;
265  }
266  man->FillNtupleFColumn(counter,ener);
267  man->FillNtupleFColumn(counter+1,x);
268  man->FillNtupleFColumn(counter+2,y);
269  man->FillNtupleFColumn(counter+3,z);
270  man->FillNtupleFColumn(counter+4,fullE);
271  man->FillNtupleFColumn(counter+5,edec);
272  man->FillNtupleFColumn(counter+6,edhc);
273  man->AddNtupleRow();
274 #ifdef debug
275  G4cout << "CCalAnalysis:: Fill Ntuple " << G4endl;
276 #endif
277 
278  // 5)
279  static G4int IDtimeProfile = -1;
280  if (IDtimeProfile < 0)
281  IDtimeProfile = man->GetH1Id("h901");
282  for (i = 0; i < numberOfSD; i++){
283  G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
284  if (caloHCid >= 0) {
285  CCalG4HitCollection* theHC =
286  (CCalG4HitCollection*) allHC->GetHC(caloHCid);
287  if (theHC != 0) {
288  G4int nentries = theHC->entries();
289  if (nentries > 0) {
290  for (G4int k=0; k<nentries; k++) {
291  CCalG4Hit* aHit = (*theHC)[k];
292  man->FillH1(IDtimeProfile,aHit->getTimeSlice(),aHit->getEnergyDeposit()/GeV);
293 
294 #ifdef debug
295  G4cout << "CCalAnalysis:: Fill Time Profile with Hit " << k
296  << " Edeposit " << edep << " Gev" << G4endl;
297 #endif
298  }
299  }
300  }
301  }
302  }
303 }
304 
305