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