| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30
31 #include "Run.hh"
32
33 #include "DetectorConstruction.hh"
34 #include "HistoManager.hh"
35 #include "PrimaryGeneratorAction.hh"
36
37 #include "G4ProcessTable.hh"
38 #include "G4Radioactivation.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4TwoVector.hh"
41 #include "G4UnitsTable.hh"
42
43 #include <fstream>
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 Run::Run(DetectorConstruction* det) : fDetecto
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50
51 void Run::SetPrimary(G4ParticleDefinition* par
52 {
53 fParticle = particle;
54 fEkin = energy;
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 void Run::CountProcesses(const G4VProcess* pro
60 {
61 if (process == nullptr) return;
62 G4String procName = process->GetProcessName(
63 if (iVol == 1) {
64 std::map<G4String, G4int>::iterator it1 =
65 if (it1 == fProcCounter1.end()) {
66 fProcCounter1[procName] = 1;
67 }
68 else {
69 fProcCounter1[procName]++;
70 }
71 }
72
73 if (iVol == 2) {
74 std::map<G4String, G4int>::iterator it2 =
75 if (it2 == fProcCounter2.end()) {
76 fProcCounter2[procName] = 1;
77 }
78 else {
79 fProcCounter2[procName]++;
80 }
81 }
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oo
85
86 void Run::ParticleCount(G4String name, G4doubl
87 {
88 if (iVol == 1) {
89 std::map<G4String, ParticleData>::iterator
90 if (it == fParticleDataMap1.end()) {
91 fParticleDataMap1[name] = ParticleData(1
92 }
93 else {
94 ParticleData& data = it->second;
95 data.fCount++;
96 data.fEmean += Ekin;
97 // update min max
98 G4double emin = data.fEmin;
99 if (Ekin < emin) data.fEmin = Ekin;
100 G4double emax = data.fEmax;
101 if (Ekin > emax) data.fEmax = Ekin;
102 }
103 }
104
105 if (iVol == 2) {
106 std::map<G4String, ParticleData>::iterator
107 if (it == fParticleDataMap2.end()) {
108 fParticleDataMap2[name] = ParticleData(1
109 }
110 else {
111 ParticleData& data = it->second;
112 data.fCount++;
113 data.fEmean += Ekin;
114 // update min max
115 G4double emin = data.fEmin;
116 if (Ekin < emin) data.fEmin = Ekin;
117 G4double emax = data.fEmax;
118 if (Ekin > emax) data.fEmax = Ekin;
119 }
120 }
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oo
124
125 void Run::AddEdep(G4double edep1, G4double ede
126 {
127 fEdepTarget += edep1;
128 fEdepTarget2 += edep1 * edep1;
129 fEdepDetect += edep2;
130 fEdepDetect2 += edep2 * edep2;
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oo
134
135 void Run::Merge(const G4Run* run)
136 {
137 const Run* localRun = static_cast<const Run*
138
139 // primary particle info
140 //
141 fParticle = localRun->fParticle;
142 fEkin = localRun->fEkin;
143
144 // accumulate sums
145 //
146 fEdepTarget += localRun->fEdepTarget;
147 fEdepTarget2 += localRun->fEdepTarget2;
148 fEdepDetect += localRun->fEdepDetect;
149 fEdepDetect2 += localRun->fEdepDetect2;
150
151 // map: processes count in target
152
153 std::map<G4String, G4int>::const_iterator it
154 for (itp1 = localRun->fProcCounter1.begin();
155 G4String procName = itp1->first;
156 G4int localCount = itp1->second;
157 if (fProcCounter1.find(procName) == fProcC
158 fProcCounter1[procName] = localCount;
159 }
160 else {
161 fProcCounter1[procName] += localCount;
162 }
163 }
164
165 // map: processes count in detector
166
167 std::map<G4String, G4int>::const_iterator it
168 for (itp2 = localRun->fProcCounter2.begin();
169 G4String procName = itp2->first;
170 G4int localCount = itp2->second;
171 if (fProcCounter2.find(procName) == fProcC
172 fProcCounter2[procName] = localCount;
173 }
174 else {
175 fProcCounter2[procName] += localCount;
176 }
177 }
178
179 // map: created particles in target
180 std::map<G4String, ParticleData>::const_iter
181 for (itc = localRun->fParticleDataMap1.begin
182 G4String name = itc->first;
183 const ParticleData& localData = itc->secon
184 if (fParticleDataMap1.find(name) == fParti
185 fParticleDataMap1[name] =
186 ParticleData(localData.fCount, localDa
187 }
188 else {
189 ParticleData& data = fParticleDataMap1[n
190 data.fCount += localData.fCount;
191 data.fEmean += localData.fEmean;
192 G4double emin = localData.fEmin;
193 if (emin < data.fEmin) data.fEmin = emin
194 G4double emax = localData.fEmax;
195 if (emax > data.fEmax) data.fEmax = emax
196 }
197 }
198
199 // map: created particle in detector
200 std::map<G4String, ParticleData>::const_iter
201 for (itn = localRun->fParticleDataMap2.begin
202 G4String name = itn->first;
203 const ParticleData& localData = itn->secon
204 if (fParticleDataMap2.find(name) == fParti
205 fParticleDataMap2[name] =
206 ParticleData(localData.fCount, localDa
207 }
208 else {
209 ParticleData& data = fParticleDataMap2[n
210 data.fCount += localData.fCount;
211 data.fEmean += localData.fEmean;
212 G4double emin = localData.fEmin;
213 if (emin < data.fEmin) data.fEmin = emin
214 G4double emax = localData.fEmax;
215 if (emax > data.fEmax) data.fEmax = emax
216 }
217 }
218
219 G4Run::Merge(run);
220 }
221
222 //....oooOO0OOooo........oooOO0OOooo........oo
223
224 void Run::EndOfRun()
225 {
226 G4int prec = 5, wid = prec + 2;
227 G4int dfprec = G4cout.precision(prec);
228
229 // run condition
230 //
231 G4String Particle = fParticle->GetParticleNa
232 G4cout << "\n The run is " << numberOfEvent
233 << G4BestUnit(fEkin, "Energy") << " t
234
235 G4cout << "\n Target : Length = " << G4Bes
236 << " Radius = " << G4BestUnit(fDet
237 << " Material = " << fDetector->GetTa
238 G4cout << "\n Detector : Length = " << G4Bes
239 << " Thickness = " << G4BestUnit(fDet
240 << " Material = " << fDetector->GetDe
241
242 if (numberOfEvent == 0) {
243 G4cout.precision(dfprec);
244 return;
245 }
246
247 // compute mean Energy deposited and rms in
248 //
249 G4int TotNbofEvents = numberOfEvent;
250 fEdepTarget /= TotNbofEvents;
251 fEdepTarget2 /= TotNbofEvents;
252 G4double rmsEdep = fEdepTarget2 - fEdepTarge
253 if (rmsEdep > 0.)
254 rmsEdep = std::sqrt(rmsEdep);
255 else
256 rmsEdep = 0.;
257
258 G4cout << "\n Mean energy deposit in target,
259 << G4BestUnit(fEdepTarget, "Energy")
260 << G4endl;
261
262 // compute mean Energy deposited and rms in
263 //
264 fEdepDetect /= TotNbofEvents;
265 fEdepDetect2 /= TotNbofEvents;
266 rmsEdep = fEdepDetect2 - fEdepDetect * fEdep
267 if (rmsEdep > 0.)
268 rmsEdep = std::sqrt(rmsEdep);
269 else
270 rmsEdep = 0.;
271
272 G4cout << " Mean energy deposit in detector,
273 << G4BestUnit(fEdepDetect, "Energy")
274 << G4endl;
275
276 // frequency of processes in target
277 //
278 G4cout << "\n Process calls frequency in tar
279 G4int index = 0;
280 std::map<G4String, G4int>::iterator it1;
281 for (it1 = fProcCounter1.begin(); it1 != fPr
282 G4String procName = it1->first;
283 G4int count = it1->second;
284 G4String space = " ";
285 if (++index % 3 == 0) space = "\n";
286 G4cout << " " << std::setw(20) << procName
287 }
288 G4cout << G4endl;
289
290 // frequency of processes in detector
291 //
292 G4cout << "\n Process calls frequency in det
293 index = 0;
294 std::map<G4String, G4int>::iterator it2;
295 for (it2 = fProcCounter2.begin(); it2 != fPr
296 G4String procName = it2->first;
297 G4int count = it2->second;
298 G4String space = " ";
299 if (++index % 3 == 0) space = "\n";
300 G4cout << " " << std::setw(20) << procName
301 }
302 G4cout << G4endl;
303
304 // particles count in target
305 //
306 G4cout << "\n List of generated particles in
307
308 std::map<G4String, ParticleData>::iterator i
309 for (itc = fParticleDataMap1.begin(); itc !=
310 G4String name = itc->first;
311 ParticleData data = itc->second;
312 G4int count = data.fCount;
313 G4double eMean = data.fEmean / count;
314 G4double eMin = data.fEmin;
315 G4double eMax = data.fEmax;
316
317 G4cout << " " << std::setw(13) << name <<
318 << " Emean = " << std::setw(wid) <
319 << G4BestUnit(eMin, "Energy") << "
320 }
321
322 // particles count in detector
323 //
324 G4cout << "\n List of generated particles in
325
326 std::map<G4String, ParticleData>::iterator i
327 for (itn = fParticleDataMap2.begin(); itn !=
328 G4String name = itn->first;
329 ParticleData data = itn->second;
330 G4int count = data.fCount;
331 G4double eMean = data.fEmean / count;
332 G4double eMin = data.fEmin;
333 G4double eMax = data.fEmax;
334
335 G4cout << " " << std::setw(13) << name <<
336 << " Emean = " << std::setw(wid) <
337 << G4BestUnit(eMin, "Energy") << "
338 }
339 G4cout << G4endl;
340
341 // activities in VR mode
342 //
343 WriteActivity(numberOfEvent);
344
345 // remove all contents in fProcCounter, fCou
346 fProcCounter1.clear();
347 fProcCounter2.clear();
348 fParticleDataMap1.clear();
349 fParticleDataMap2.clear();
350
351 // restore default format
352 G4cout.precision(dfprec);
353 }
354
355 //....oooOO0OOooo........oooOO0OOooo........oo
356
357 void Run::WriteActivity(G4int nevent)
358 {
359 G4ProcessTable* pTable = G4ProcessTable::Get
360 G4Radioactivation* rDecay =
361 (G4Radioactivation*)pTable->FindProcess("R
362
363 // output the induced radioactivities (in VR
364 //
365 if ((rDecay == 0) || (rDecay->IsAnalogueMont
366
367 G4String fileName = G4AnalysisManager::Insta
368 std::ofstream outfile(fileName, std::ios::ou
369
370 std::vector<G4RadioactivityTable*> theTables
371
372 for (size_t i = 0; i < theTables.size(); i++
373 G4double rate, error;
374 outfile << "Radioactivities in decay windo
375 outfile << "Z \tA \tE \tActivity (decays/w
376
377 map<G4ThreeVector, G4TwoVector>* aMap = th
378 map<G4ThreeVector, G4TwoVector>::iterator
379 for (iter = aMap->begin(); iter != aMap->e
380 rate = iter->second.x() / nevent;
381 error = std::sqrt(iter->second.y()) / ne
382 if (rate < 0.) rate = 0.; // statically
383 outfile << iter->first.x() << "\t" << it
384 << rate << "\t" << error << G4en
385 }
386 outfile << G4endl;
387 }
388 outfile.close();
389 }
390
391 //....oooOO0OOooo........oooOO0OOooo........oo
392