| 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 /// \file Run.cc 26 /// \file Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // << 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "Run.hh" 33 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "HistoManager.hh"
38 37
39 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh"
40 #include "G4SystemOfUnits.hh" <<
41 #include "G4Track.hh" 39 #include "G4Track.hh"
>> 40
42 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh"
>> 42 #include "G4SystemOfUnits.hh"
43 43
44 #include <iomanip> 44 #include <iomanip>
45 45
46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 47
48 Run::Run(DetectorConstruction* det) : fDetecto << 48 Run::Run(DetectorConstruction* det)
>> 49 : fDetector(det)
49 { 50 {
50 // initialize energy deposited per absorber << 51 //initialize energy deposited per absorber
51 // 52 //
52 for (G4int k = 0; k < kMaxAbsor; k++) { << 53 for (G4int k=0; k<kMaxAbsor; k++) {
53 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = << 54 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
54 } 55 }
55 << 56
56 // initialize total energy deposited 57 // initialize total energy deposited
57 // 58 //
58 fEdepTot = fEdepTot2 = 0.; 59 fEdepTot = fEdepTot2 = 0.;
59 << 60
60 // initialize leakage 61 // initialize leakage
61 // 62 //
62 fEnergyLeak[0] = fEnergyLeak[1] = 0.; 63 fEnergyLeak[0] = fEnergyLeak[1] = 0.;
63 fEleakTot = fEleakTot2 = 0.; 64 fEleakTot = fEleakTot2 = 0.;
64 << 65
65 // initialize total energy released 66 // initialize total energy released
66 // 67 //
67 fEtotal = fEtotal2 = 0.; 68 fEtotal = fEtotal2 = 0.;
68 << 69
69 // initialize Eflow << 70 //initialize Eflow
70 // 71 //
71 G4int nbPlanes = (fDetector->GetNbOfLayers() << 72 G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
72 fEnergyFlow.resize(nbPlanes); 73 fEnergyFlow.resize(nbPlanes);
73 for (G4int k = 0; k < nbPlanes; k++) { << 74 for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = 0.; }
74 fEnergyFlow[k] = 0.; <<
75 } <<
76 } 75 }
77 76
78 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 78
80 void Run::SetPrimary(G4ParticleDefinition* par 79 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
81 { << 80 {
82 fParticle = particle; 81 fParticle = particle;
83 fEkin = energy; 82 fEkin = energy;
84 } 83 }
85 << 84
86 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 86
88 void Run::CountProcesses(const G4VProcess* pro << 87 void Run::CountProcesses(const G4VProcess* process)
89 { 88 {
90 if (process == nullptr) return; 89 if (process == nullptr) return;
91 G4String procName = process->GetProcessName( 90 G4String procName = process->GetProcessName();
92 std::map<G4String, G4int>::iterator it = fPr << 91 std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
93 if (it == fProcCounter.end()) { << 92 if ( it == fProcCounter.end()) {
94 fProcCounter[procName] = 1; 93 fProcCounter[procName] = 1;
95 } 94 }
96 else { 95 else {
97 fProcCounter[procName]++; << 96 fProcCounter[procName]++;
98 } 97 }
99 } 98 }
100 << 99
101 //....oooOO0OOooo........oooOO0OOooo........oo 100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 101
103 void Run::SumEdepPerAbsorber(G4int kAbs, G4dou 102 void Run::SumEdepPerAbsorber(G4int kAbs, G4double EAbs, G4double LAbs)
104 { 103 {
105 // accumulate statistic with restriction << 104 //accumulate statistic with restriction
106 // 105 //
107 fSumEAbs[kAbs] += EAbs; << 106 fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs;
108 fSum2EAbs[kAbs] += EAbs * EAbs; << 107 fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs;
109 fSumLAbs[kAbs] += LAbs; <<
110 fSum2LAbs[kAbs] += LAbs * LAbs; <<
111 } 108 }
112 109
113 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 111
115 void Run::SumEnergies(G4double edeptot, G4doub 112 void Run::SumEnergies(G4double edeptot, G4double eleak0, G4double eleak1)
116 { 113 {
117 fEdepTot += edeptot; << 114 fEdepTot += edeptot; fEdepTot2 += edeptot*edeptot;
118 fEdepTot2 += edeptot * edeptot; << 115
119 << 116 fEnergyLeak[0] += eleak0; fEnergyLeak[1] += eleak1;
120 fEnergyLeak[0] += eleak0; <<
121 fEnergyLeak[1] += eleak1; <<
122 G4double eleaktot = eleak0 + eleak1; 117 G4double eleaktot = eleak0 + eleak1;
123 fEleakTot += eleaktot; << 118 fEleakTot += eleaktot; fEleakTot2 += eleaktot*eleaktot;
124 fEleakTot2 += eleaktot * eleaktot; << 119
125 <<
126 G4double etotal = edeptot + eleaktot; 120 G4double etotal = edeptot + eleaktot;
127 fEtotal += etotal; << 121 fEtotal += etotal; fEtotal2 += etotal*etotal;
128 fEtotal2 += etotal * etotal; <<
129 } 122 }
130 123
131 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 125
133 void Run::SumEnergyFlow(G4int plane, G4double 126 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
134 { 127 {
135 fEnergyFlow[plane] += Eflow; 128 fEnergyFlow[plane] += Eflow;
136 } 129 }
137 130
138 //....oooOO0OOooo........oooOO0OOooo........oo 131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
139 132
140 void Run::Merge(const G4Run* run) 133 void Run::Merge(const G4Run* run)
141 { 134 {
142 const Run* localRun = static_cast<const Run* 135 const Run* localRun = static_cast<const Run*>(run);
143 136
144 // pass information about primary particle 137 // pass information about primary particle
145 fParticle = localRun->fParticle; 138 fParticle = localRun->fParticle;
146 fEkin = localRun->fEkin; << 139 fEkin = localRun->fEkin;
147 140
148 // accumulate sums 141 // accumulate sums
149 // 142 //
150 for (G4int k = 0; k < kMaxAbsor; k++) { << 143 for (G4int k=0; k<kMaxAbsor; k++) {
151 fSumEAbs[k] += localRun->fSumEAbs[k]; << 144 fSumEAbs[k] += localRun->fSumEAbs[k];
152 fSum2EAbs[k] += localRun->fSum2EAbs[k]; << 145 fSum2EAbs[k] += localRun->fSum2EAbs[k];
153 fSumLAbs[k] += localRun->fSumLAbs[k]; << 146 fSumLAbs[k] += localRun->fSumLAbs[k];
154 fSum2LAbs[k] += localRun->fSum2LAbs[k]; 147 fSum2LAbs[k] += localRun->fSum2LAbs[k];
155 } 148 }
156 << 149
157 fEdepTot += localRun->fEdepTot; << 150 fEdepTot += localRun->fEdepTot;
158 fEdepTot2 += localRun->fEdepTot2; 151 fEdepTot2 += localRun->fEdepTot2;
159 << 152
160 fEnergyLeak[0] += localRun->fEnergyLeak[0]; << 153 fEnergyLeak[0] += localRun->fEnergyLeak[0];
161 fEnergyLeak[1] += localRun->fEnergyLeak[1]; << 154 fEnergyLeak[1] += localRun->fEnergyLeak[1];
162 << 155
163 fEleakTot += localRun->fEleakTot; << 156 fEleakTot += localRun->fEleakTot;
164 fEleakTot2 += localRun->fEleakTot2; 157 fEleakTot2 += localRun->fEleakTot2;
>> 158
>> 159 fEtotal += localRun->fEtotal;
>> 160 fEtotal2 += localRun->fEtotal2;
>> 161
>> 162 G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
>> 163 for (G4int k=0; k<nbPlanes; k++) {
>> 164 fEnergyFlow[k] += localRun->fEnergyFlow[k];
>> 165 }
>> 166
>> 167 //map: processes count
>> 168 std::map<G4String,G4int>::const_iterator itp;
>> 169 for ( itp = localRun->fProcCounter.begin();
>> 170 itp != localRun->fProcCounter.end(); ++itp ) {
165 171
166 fEtotal += localRun->fEtotal; <<
167 fEtotal2 += localRun->fEtotal2; <<
168 <<
169 G4int nbPlanes = (fDetector->GetNbOfLayers() <<
170 for (G4int k = 0; k < nbPlanes; k++) { <<
171 fEnergyFlow[k] += localRun->fEnergyFlow[k] <<
172 } <<
173 <<
174 // map: processes count <<
175 std::map<G4String, G4int>::const_iterator it <<
176 for (itp = localRun->fProcCounter.begin(); i <<
177 G4String procName = itp->first; 172 G4String procName = itp->first;
178 G4int localCount = itp->second; 173 G4int localCount = itp->second;
179 if (fProcCounter.find(procName) == fProcCo << 174 if ( fProcCounter.find(procName) == fProcCounter.end()) {
180 fProcCounter[procName] = localCount; 175 fProcCounter[procName] = localCount;
181 } 176 }
182 else { 177 else {
183 fProcCounter[procName] += localCount; 178 fProcCounter[procName] += localCount;
184 } << 179 }
185 } 180 }
186 << 181
187 G4Run::Merge(run); << 182 G4Run::Merge(run);
188 } << 183 }
189 184
190 //....oooOO0OOooo........oooOO0OOooo........oo 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 186
192 void Run::EndOfRun() 187 void Run::EndOfRun()
193 { 188 {
194 // run condition << 189 //run condition
195 // << 190 //
196 G4String Particle = fParticle->GetParticleNa << 191 G4String Particle = fParticle->GetParticleName();
197 G4cout << "\n ---> The run is " << numberOfE << 192 G4cout << "\n ---> The run is " << numberOfEvent << " "<< Particle << " of "
198 << G4BestUnit(fEkin, "Energy") << " t << 193 << G4BestUnit(fEkin,"Energy") << " through calorimeter" << G4endl;
199 << 194
200 // frequency of processes << 195 //frequency of processes
201 // 196 //
202 G4cout << "\n Process calls frequency :" << 197 G4cout << "\n Process calls frequency :" << G4endl;
203 G4int index = 0; 198 G4int index = 0;
204 std::map<G4String, G4int>::iterator it; << 199 std::map<G4String,G4int>::iterator it;
205 for (it = fProcCounter.begin(); it != fProcC 200 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
206 G4String procName = it->first; << 201 G4String procName = it->first;
207 G4int count = it->second; << 202 G4int count = it->second;
208 G4String space = " "; << 203 G4String space = " "; if (++index%3 == 0) space = "\n";
209 if (++index % 3 == 0) space = "\n"; << 204 G4cout << " " << std::setw(22) << procName << "="<< std::setw(10) << count
210 G4cout << " " << std::setw(22) << procName << 205 << space;
211 } 206 }
212 << 207
213 G4cout << G4endl; 208 G4cout << G4endl;
214 G4int nEvt = numberOfEvent; 209 G4int nEvt = numberOfEvent;
215 G4double norm = G4double(nEvt); << 210 G4double norm = G4double(nEvt);
216 if (norm > 0) norm = 1. / norm; << 211 if(norm > 0) norm = 1./norm;
217 G4double qnorm = std::sqrt(norm); 212 G4double qnorm = std::sqrt(norm);
218 213
219 // energy deposit per absorber << 214 //energy deposit per absorber
220 // 215 //
221 G4double beamEnergy = fEkin; 216 G4double beamEnergy = fEkin;
222 G4double sqbeam = std::sqrt(beamEnergy / GeV << 217 G4double sqbeam = std::sqrt(beamEnergy/GeV);
223 218
224 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resol << 219 G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
225 G4double MeanLAbs, MeanLAbs2, rmsLAbs; << 220 G4double MeanLAbs,MeanLAbs2,rmsLAbs;
226 221
227 std::ios::fmtflags mode = G4cout.flags(); 222 std::ios::fmtflags mode = G4cout.flags();
228 G4int prec = G4cout.precision(2); << 223 G4int prec = G4cout.precision(2);
229 G4cout << "\n------------------------------- 224 G4cout << "\n------------------------------------------------------------\n";
230 G4cout << std::setw(16) << "material" << std << 225 G4cout << std::setw(16) << "material"
231 << "sqrt(E0(GeV))*rmsE/Edep" << std:: << 226 << std::setw(22) << "Edep rmsE"
232 << 227 << std::setw(31) << "sqrt(E0(GeV))*rmsE/Edep"
233 for (G4int k = 1; k <= fDetector->GetNbOfAbs << 228 << std::setw(23) << "total tracklen \n \n";
234 MeanEAbs = fSumEAbs[k] * norm; << 229
235 MeanEAbs2 = fSum2EAbs[k] * norm; << 230 for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
236 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - M << 231 {
237 << 232 MeanEAbs = fSumEAbs[k]*norm;
238 resolution = 100. * sqbeam * rmsEAbs / Mea << 233 MeanEAbs2 = fSum2EAbs[k]*norm;
239 rmsres = resolution * qnorm; << 234 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
240 << 235
241 MeanLAbs = fSumLAbs[k] * norm; << 236 resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
242 MeanLAbs2 = fSum2LAbs[k] * norm; << 237 rmsres = resolution*qnorm;
243 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - M << 238
244 << 239 MeanLAbs = fSumLAbs[k]*norm;
245 // print << 240 MeanLAbs2 = fSum2LAbs[k]*norm;
246 // << 241 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
247 G4cout << std::setw(2) << k << std::setw(1 << 242
248 << std::setprecision(5) << std::set << 243 //print
249 << std::setprecision(4) << std::set << 244 //
250 << resolution << " +- " << std::set << 245 G4cout
251 << std::setprecision(4) << std::set << 246 << std::setw(2) << k
252 << std::setprecision(3) << std::set << 247 << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName()
253 } << 248 << std::setprecision(5)
254 << 249 << std::setw(10) << G4BestUnit(MeanEAbs,"Energy")
255 // total energy deposited << 250 << std::setprecision(4)
>> 251 << std::setw(8) << G4BestUnit( rmsEAbs,"Energy")
>> 252 << std::setw(10) << resolution << " +- "
>> 253 << std::setprecision(3)
>> 254 << std::setw(5) << rmsres << " %"
>> 255 << std::setprecision(4)
>> 256 << std::setw(12) << G4BestUnit(MeanLAbs,"Length") << " +- "
>> 257 << std::setprecision(3)
>> 258 << std::setw(5) << G4BestUnit( rmsLAbs,"Length")
>> 259 << G4endl;
>> 260 }
>> 261
>> 262 //total energy deposited
256 // 263 //
257 fEdepTot /= nEvt; << 264 fEdepTot /= nEvt;
258 fEdepTot2 /= nEvt; << 265 fEdepTot2 /= nEvt;
259 G4double rmsEdep = std::sqrt(std::abs(fEdepT << 266 G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot*fEdepTot));
260 << 267
261 G4cout << "\n Total energy deposited = " << << 268 G4cout << "\n Total energy deposited = " << std::setprecision(4)
262 << " +- " << G4BestUnit(rmsEdep, "Ene << 269 << G4BestUnit(fEdepTot,"Energy")
263 << 270 << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
264 // Energy leakage << 271
>> 272 //Energy leakage
265 // 273 //
266 fEnergyLeak[0] /= nEvt; << 274 fEnergyLeak[0] /= nEvt;
267 fEnergyLeak[1] /= nEvt; 275 fEnergyLeak[1] /= nEvt;
268 fEleakTot /= nEvt; << 276 fEleakTot /= nEvt;
269 fEleakTot2 /= nEvt; << 277 fEleakTot2 /= nEvt;
270 G4double rmsEleak = std::sqrt(std::abs(fElea << 278 G4double rmsEleak = std::sqrt(std::abs(fEleakTot2 - fEleakTot*fEleakTot));
271 << 279
272 G4cout << " Leakage : primary = " << G4Best << 280 G4cout << " Leakage : primary = "
273 << " secondaries = " << G4BestUnit( << 281 << G4BestUnit(fEnergyLeak[0],"Energy")
274 << " ---> total = " << G4BestUnit(fE << 282 << " secondaries = "
275 << G4BestUnit(rmsEleak, "Energy") << << 283 << G4BestUnit(fEnergyLeak[1],"Energy")
276 << 284 << " ---> total = " << G4BestUnit(fEleakTot, "Energy")
277 // total energy released << 285 << " +- " << G4BestUnit(rmsEleak, "Energy") << G4endl;
278 // << 286
279 fEtotal /= nEvt; << 287 //total energy released
280 fEtotal2 /= nEvt; << 288 //
281 G4double rmsEtotal = std::sqrt(std::abs(fEto << 289 fEtotal /= nEvt;
282 << 290 fEtotal2 /= nEvt;
283 G4cout << " Total energy released : Edep + << 291 G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal*fEtotal));
284 << G4BestUnit(rmsEtotal, "Energy") << << 292
>> 293 G4cout << " Total energy released : Edep + Eleak = "
>> 294 << G4BestUnit(fEtotal,"Energy")
>> 295 << " +- " << G4BestUnit(rmsEtotal, "Energy") << G4endl;
285 G4cout << "--------------------------------- 296 G4cout << "------------------------------------------------------------\n";
286 << 297
287 // Energy flow << 298 //Energy flow
288 // 299 //
289 G4AnalysisManager* analysis = G4AnalysisMana 300 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
290 G4int Idmax = (fDetector->GetNbOfLayers()) * << 301 G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
291 for (G4int Id = 1; Id <= Idmax + 1; Id++) { << 302 for (G4int Id=1; Id<=Idmax+1; Id++) {
292 analysis->FillH1(2 * kMaxAbsor + 1, (G4dou << 303 analysis->FillH1(2*kMaxAbsor+1, (G4double)Id, fEnergyFlow[Id]/nEvt);
293 } 304 }
294 << 305
295 // normalize histograms << 306 //normalize histograms
296 // 307 //
297 for (G4int ih = kMaxAbsor + 1; ih < 2 * kMax << 308 for (G4int ih = kMaxAbsor+1; ih < 2*kMaxAbsor+1; ih++) {
298 analysis->ScaleH1(ih, norm / MeV); << 309 analysis->ScaleH1(ih,norm/MeV);
299 } 310 }
300 << 311
301 // remove all contents in fProcCounter << 312 //remove all contents in fProcCounter
302 fProcCounter.clear(); 313 fProcCounter.clear();
303 << 314
304 G4cout.setf(mode, std::ios::floatfield); << 315 G4cout.setf(mode,std::ios::floatfield);
305 G4cout.precision(prec); 316 G4cout.precision(prec);
306 } 317 }
307 318
308 //....oooOO0OOooo........oooOO0OOooo........oo 319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 320