| 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