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 ExGflashEventAction.cc 27 /// \file ExGflashEventAction.cc 28 /// \brief Implementation of the ExGflashEvent 28 /// \brief Implementation of the ExGflashEventAction class 29 // 29 // 30 // Created by Joanna Weng 26.11.2004 30 // Created by Joanna Weng 26.11.2004 31 31 32 #include "ExGflashEventAction.hh" << 33 32 >> 33 #include "ExGflashEventAction.hh" 34 #include "ExGflashHit.hh" 34 #include "ExGflashHit.hh" 35 << 36 #include "G4Event.hh" << 37 #include "G4EventManager.hh" 35 #include "G4EventManager.hh" 38 #include "G4SDManager.hh" 36 #include "G4SDManager.hh" 39 #include "G4SystemOfUnits.hh" << 40 #include "G4TrajectoryContainer.hh" << 41 #include "G4UImanager.hh" 37 #include "G4UImanager.hh" 42 // std << 38 #include "G4TrajectoryContainer.hh" 43 #include <algorithm> << 39 #include "G4Event.hh" >> 40 #include "G4SystemOfUnits.hh" >> 41 //std 44 #include <iostream> 42 #include <iostream> 45 // Gflash << 43 #include <algorithm> >> 44 //Gflash 46 using namespace std; 45 using namespace std; 47 46 >> 47 48 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 49 49 50 ExGflashEventAction::ExGflashEventAction() = d << 50 ExGflashEventAction::ExGflashEventAction() >> 51 : G4UserEventAction(), >> 52 fNevent(0),fDtime(0.0),fCalorimeterCollectionId(-1) >> 53 { >> 54 } 51 55 52 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 53 57 54 ExGflashEventAction::~ExGflashEventAction() 58 ExGflashEventAction::~ExGflashEventAction() 55 { 59 { 56 if (fNevent > 0) { << 60 if ( fNevent > 0 ) { 57 G4cout << "Internal Real Elapsed Time /eve << 61 G4cout << "Internal Real Elapsed Time /event is: "<< fDtime /fNevent<< G4endl; 58 } 62 } 59 } 63 } 60 64 61 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 62 66 63 void ExGflashEventAction::BeginOfEventAction(c << 67 void ExGflashEventAction::BeginOfEventAction(const G4Event *evt) 64 { 68 { 65 fTimerIntern.Start(); 69 fTimerIntern.Start(); 66 G4cout << " ------ Start ExGflashEventAction << 70 G4cout<<" ------ Start ExGflashEventAction ----- "<<G4endl; 67 fNevent = evt->GetEventID(); << 71 fNevent=evt->GetEventID(); 68 G4cout << " Start generating event Nr " << f << 72 G4cout<<" Start generating event Nr "<<fNevent<<G4endl<<G4endl; 69 } 73 } 70 74 71 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 72 76 73 void ExGflashEventAction::EndOfEventAction(con << 77 void ExGflashEventAction::EndOfEventAction(const G4Event *evt) 74 { << 78 { 75 fTimerIntern.Stop(); 79 fTimerIntern.Stop(); 76 G4cout << G4endl; 80 G4cout << G4endl; 77 G4cout << "********************************* 81 G4cout << "******************************************"; 78 G4cout << G4endl; 82 G4cout << G4endl; 79 G4cout << "Internal Real Elapsed Time is: " << 83 G4cout << "Internal Real Elapsed Time is: "<< fTimerIntern.GetRealElapsed(); 80 G4cout << G4endl; 84 G4cout << G4endl; 81 G4cout << "Internal System Elapsed Time: " < 85 G4cout << "Internal System Elapsed Time: " << fTimerIntern.GetSystemElapsed(); 82 G4cout << G4endl; 86 G4cout << G4endl; 83 G4cout << "Internal GetUserElapsed Time: " < 87 G4cout << "Internal GetUserElapsed Time: " << fTimerIntern.GetUserElapsed(); 84 G4cout << G4endl; 88 G4cout << G4endl; 85 G4cout << "********************************* << 89 G4cout << "******************************************"<< G4endl; 86 fDtime += fTimerIntern.GetRealElapsed(); << 90 fDtime+=fTimerIntern.GetRealElapsed(); 87 G4cout << " ------ ExGflashEventAction::End << 91 G4cout<<" ------ ExGflashEventAction::End of event nr. "<<fNevent<<" -----"<< G4endl; 88 << 92 89 G4SDManager* SDman = G4SDManager::GetSDMpoin << 93 G4SDManager * SDman = G4SDManager::GetSDMpointer(); 90 G4String colNam; 94 G4String colNam; 91 fCalorimeterCollectionId = SDman->GetCollect << 95 fCalorimeterCollectionId=SDman->GetCollectionID(colNam="ExGflashCollection"); 92 if (fCalorimeterCollectionId < 0) return; << 96 if (fCalorimeterCollectionId<0) return; 93 G4HCofThisEvent* HCE = evt->GetHCofThisEvent << 97 G4HCofThisEvent * HCE = evt->GetHCofThisEvent(); 94 ExGflashHitsCollection* THC = nullptr; << 98 ExGflashHitsCollection* THC = 0; 95 G4double totE = 0; 99 G4double totE = 0; 96 // Read out of the crysta ECAL 100 // Read out of the crysta ECAL 97 THC = (ExGflashHitsCollection*)(HCE->GetHC(f << 101 THC=(ExGflashHitsCollection *)(HCE->GetHC(fCalorimeterCollectionId)); 98 if (THC) { << 102 if (THC) 99 /// Hits in sensitive Detector << 103 { 100 int n_hit = THC->entries(); << 104 /// Hits in sensitive Detector 101 G4cout << " " << n_hit << " hits are stor << 105 int n_hit = THC->entries(); 102 G4PrimaryVertex* pvertex = evt->GetPrimary << 106 G4cout<<" " << n_hit<< " hits are stored in ExGflashHitsCollection "<<G4endl; 103 /// Computing (x,y,z) of vertex of initial << 107 G4PrimaryVertex* pvertex=evt->GetPrimaryVertex(); 104 G4ThreeVector vtx = pvertex->GetPosition() << 108 ///Computing (x,y,z) of vertex of initial particles 105 G4PrimaryParticle* pparticle = pvertex->Ge << 109 G4ThreeVector vtx=pvertex->GetPosition(); >> 110 G4PrimaryParticle* pparticle=pvertex->GetPrimary(); 106 // direction of the Shower 111 // direction of the Shower 107 G4ThreeVector mom = pparticle->GetMomentum << 112 G4ThreeVector mom=pparticle->GetMomentum()/pparticle->GetMomentum().mag(); 108 << 113 109 //@@@ ExGflashEventAction: Magicnumber 114 //@@@ ExGflashEventAction: Magicnumber 110 G4double energyincrystal[100]; 115 G4double energyincrystal[100]; 111 G4int hitsincrystal[100]; 116 G4int hitsincrystal[100]; 112 for (int i = 0; i < 100; i++) << 117 for (int i=0;i<100;i++) energyincrystal[i]=0.; 113 energyincrystal[i] = 0.; << 118 for (int i=0;i<100;i++) hitsincrystal[i]=0.; 114 for (int i = 0; i < 100; i++) << 119 115 hitsincrystal[i] = 0.; << 116 << 117 //@@@ ExGflashEventAction: Magicnumber 120 //@@@ ExGflashEventAction: Magicnumber 118 /// For all Hits in sensitive detector 121 /// For all Hits in sensitive detector 119 for (int i = 0; i < n_hit; i++) { << 122 for (int i=0;i<n_hit;i++) 120 G4double estep = (*THC)[i]->GetEdep() / << 123 { 121 if (estep > 0.0) { << 124 G4double estep = (*THC)[i]->GetEdep()/GeV; 122 totE += (*THC)[i]->GetEdep() / GeV; << 125 if (estep >0.0) 123 G4int num = (*THC)[i]->GetCrystalNum() << 126 { 124 << 127 totE += (*THC)[i]->GetEdep()/GeV; 125 energyincrystal[num] += (*THC)[i]->Get << 128 G4int num=(*THC)[i]->GetCrystalNum(); 126 hitsincrystal[num]++; << 129 127 // G4cout << num << G4endl; << 130 energyincrystal[num]+=(*THC)[i]->GetEdep()/GeV; 128 // G4cout << " Crystal Nummer " << << 131 hitsincrystal[num]++; 129 // G4cout << (*THC)[i]->GetCrystalN << 132 //G4cout << num << G4endl; 130 // " "<<(*THC)[i]->GetCrystalNum()%1 << 133 // G4cout << " Crystal Nummer " << (*THC)[i]->GetCrystalNum() << G4endl; 131 << 134 // G4cout << (*THC)[i]->GetCrystalNum() /10 << 132 G4ThreeVector hitpos = (*THC)[i]->GetP << 135 // " "<<(*THC)[i]->GetCrystalNum()%10 << G4endl; 133 G4ThreeVector l(hitpos.x(), hitpos.y() << 136 134 // distance from shower start << 137 G4ThreeVector hitpos=(*THC)[i]->GetPos(); 135 l = l - vtx; << 138 G4ThreeVector l (hitpos.x(), hitpos.y(), hitpos.z()); 136 // projection on shower axis = longitu << 139 // distance from shower start 137 G4ThreeVector longitudinal = l; << 140 l = l - vtx; 138 // shower profiles (Radial) << 141 // projection on shower axis = longitudinal profile 139 G4ThreeVector radial = vtx.cross(l); << 142 G4ThreeVector longitudinal = l; 140 } << 143 // shower profiles (Radial) 141 } << 144 G4ThreeVector radial = vtx.cross(l); 142 G4double max = 0; << 145 } 143 G4int index = 0; << 144 // Find crystal with maximum energy << 145 for (int i = 0; i < 100; i++) { << 146 // G4cout << i <<" " << energyincrystal << 147 if (max < energyincrystal[i]) { << 148 max = energyincrystal[i]; << 149 index = i; << 150 } 146 } 151 } << 147 G4double max = 0; 152 // G4cout << index <<" NMAX " << index << << 148 G4int index = 0; 153 << 149 //Find crystal with maximum energy >> 150 for (int i=0;i<100;i++) >> 151 { >> 152 //G4cout << i <<" " << energyincrystal[i] << G4endl; >> 153 if (max <energyincrystal[i]) >> 154 { >> 155 max=energyincrystal[i]; >> 156 index = i; >> 157 } >> 158 } >> 159 //G4cout << index <<" NMAX " << index << G4endl; >> 160 154 // 3x3 matrix of crystals around the cryst 161 // 3x3 matrix of crystals around the crystal with the maximum energy deposit 155 G4double e3x3 = energyincrystal[index] + e << 162 G4double e3x3 = 156 + energyincrystal[index - << 163 energyincrystal[index]+energyincrystal[index+1]+energyincrystal[index-1]+ 157 + energyincrystal[index - << 164 energyincrystal[index-10]+energyincrystal[index-9]+energyincrystal[index-11]+ 158 + energyincrystal[index + << 165 energyincrystal[index+10]+energyincrystal[index+11]+energyincrystal[index+9]; 159 << 166 160 // 5x5 matrix of crystals around the cryst << 167 // 5x5 matrix of crystals around the crystal with the maximum energy deposit 161 G4double e5x5 = 168 G4double e5x5 = 162 energyincrystal[index] + energyincrystal << 169 energyincrystal[index]+energyincrystal[index+1]+energyincrystal[index-1]+ 163 + energyincrystal[index + 2] + energyinc << 170 energyincrystal[index+2]+energyincrystal[index-2]+ 164 + energyincrystal[index - 9] + energyinc << 171 energyincrystal[index-10]+energyincrystal[index-9]+energyincrystal[index-11]+ 165 + energyincrystal[index - 12] + energyin << 172 energyincrystal[index-8]+energyincrystal[index-12]+ 166 + energyincrystal[index + 9] + energyinc << 173 energyincrystal[index+10]+energyincrystal[index+11]+energyincrystal[index+9]+ >> 174 energyincrystal[index+12]+energyincrystal[index+8]; 167 175 168 // 3x3 matrix of crystals around the cryst 176 // 3x3 matrix of crystals around the crystal with the maximum energy deposit 169 G4int num3x3 = hitsincrystal[index] + hits << 177 G4int num3x3 = 170 + hitsincrystal[index - 10] << 178 hitsincrystal[index]+hitsincrystal[index+1]+hitsincrystal[index-1]+ 171 + hitsincrystal[index - 11] << 179 hitsincrystal[index-10]+hitsincrystal[index-9]+hitsincrystal[index-11]+ 172 + hitsincrystal[index + 11] << 180 hitsincrystal[index+10]+hitsincrystal[index+11]+hitsincrystal[index+9]; 173 << 181 174 // 5x5 matrix of crystals around the cryst << 182 // 5x5 matrix of crystals around the crystal with the maximum energy deposit 175 G4int num5x5 = hitsincrystal[index] + hits << 183 G4int num5x5 = 176 + hitsincrystal[index + 2] << 184 hitsincrystal[index]+hitsincrystal[index+1]+hitsincrystal[index-1]+ 177 + hitsincrystal[index - 9] << 185 hitsincrystal[index+2]+hitsincrystal[index-2]+ 178 + hitsincrystal[index - 12] << 186 hitsincrystal[index-10]+hitsincrystal[index-9]+hitsincrystal[index-11]+ 179 + hitsincrystal[index + 11] << 187 hitsincrystal[index-8]+hitsincrystal[index-12]+ 180 + hitsincrystal[index + 12] << 188 hitsincrystal[index+10]+hitsincrystal[index+11]+hitsincrystal[index+9]+ 181 << 189 hitsincrystal[index+12]+hitsincrystal[index+8]; 182 G4cout << " e1 " << energyincrystal[ind << 190 183 << G4endl; << 191 G4cout << " e1 " << energyincrystal[index] 184 << 192 << " e3x3 " << e3x3<< " GeV e5x5 " <<e5x5 <<G4endl; 185 G4cout << " num1 " << hitsincrystal[ind << 193 186 << num5x5 << G4endl; << 194 G4cout << " num1 " << hitsincrystal[index] >> 195 << " num3x3 " << num3x3<< " num5x5 " <<num5x5 <<G4endl; 187 } 196 } 188 << 197 189 G4cout << " Total energy deposited in the ca 198 G4cout << " Total energy deposited in the calorimeter: " << totE << " (GeV)" << G4endl; 190 G4TrajectoryContainer* trajectoryContainer = << 199 G4TrajectoryContainer * trajectoryContainer = evt->GetTrajectoryContainer(); 191 G4int n_trajectories = 0; 200 G4int n_trajectories = 0; 192 if (trajectoryContainer) { << 201 if(trajectoryContainer){ n_trajectories = trajectoryContainer->entries(); } 193 n_trajectories = trajectoryContainer->entr << 202 G4cout << " " << n_trajectories << " trajectories stored in this event." << G4endl; 194 } << 203 195 G4cout << " " << n_trajectories << " traj << 196 } 204 } 197 205 198 //....oooOO0OOooo........oooOO0OOooo........oo 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 207 >> 208 >> 209 >> 210 >> 211 >> 212 >> 213 >> 214 >> 215 >> 216 >> 217 >> 218 >> 219 199 220