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