| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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........oooOO0OOooo........oooOO0OOooo......
46
47 Run::Run(DetectorConstruction* det) : fDetector(det) {}
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
52 {
53 fParticle = particle;
54 fEkin = energy;
55 }
56
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59 void Run::CountProcesses(const G4VProcess* process, G4int iVol)
60 {
61 if (process == nullptr) return;
62 G4String procName = process->GetProcessName();
63 if (iVol == 1) {
64 std::map<G4String, G4int>::iterator it1 = fProcCounter1.find(procName);
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 = fProcCounter2.find(procName);
75 if (it2 == fProcCounter2.end()) {
76 fProcCounter2[procName] = 1;
77 }
78 else {
79 fProcCounter2[procName]++;
80 }
81 }
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85
86 void Run::ParticleCount(G4String name, G4double Ekin, G4int iVol)
87 {
88 if (iVol == 1) {
89 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
90 if (it == fParticleDataMap1.end()) {
91 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
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 it = fParticleDataMap2.find(name);
107 if (it == fParticleDataMap2.end()) {
108 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
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........oooOO0OOooo........oooOO0OOooo......
124
125 void Run::AddEdep(G4double edep1, G4double edep2)
126 {
127 fEdepTarget += edep1;
128 fEdepTarget2 += edep1 * edep1;
129 fEdepDetect += edep2;
130 fEdepDetect2 += edep2 * edep2;
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134
135 void Run::Merge(const G4Run* run)
136 {
137 const Run* localRun = static_cast<const Run*>(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 itp1;
154 for (itp1 = localRun->fProcCounter1.begin(); itp1 != localRun->fProcCounter1.end(); ++itp1) {
155 G4String procName = itp1->first;
156 G4int localCount = itp1->second;
157 if (fProcCounter1.find(procName) == fProcCounter1.end()) {
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 itp2;
168 for (itp2 = localRun->fProcCounter2.begin(); itp2 != localRun->fProcCounter2.end(); ++itp2) {
169 G4String procName = itp2->first;
170 G4int localCount = itp2->second;
171 if (fProcCounter2.find(procName) == fProcCounter2.end()) {
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_iterator itc;
181 for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
182 G4String name = itc->first;
183 const ParticleData& localData = itc->second;
184 if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
185 fParticleDataMap1[name] =
186 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
187 }
188 else {
189 ParticleData& data = fParticleDataMap1[name];
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_iterator itn;
201 for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
202 G4String name = itn->first;
203 const ParticleData& localData = itn->second;
204 if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
205 fParticleDataMap2[name] =
206 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
207 }
208 else {
209 ParticleData& data = fParticleDataMap2[name];
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........oooOO0OOooo........oooOO0OOooo......
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->GetParticleName();
232 G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
233 << G4BestUnit(fEkin, "Energy") << " through : ";
234
235 G4cout << "\n Target : Length = " << G4BestUnit(fDetector->GetTargetLength(), "Length")
236 << " Radius = " << G4BestUnit(fDetector->GetTargetRadius(), "Length")
237 << " Material = " << fDetector->GetTargetMaterial()->GetName();
238 G4cout << "\n Detector : Length = " << G4BestUnit(fDetector->GetDetectorLength(), "Length")
239 << " Thickness = " << G4BestUnit(fDetector->GetDetectorThickness(), "Length")
240 << " Material = " << fDetector->GetDetectorMaterial()->GetName() << G4endl;
241
242 if (numberOfEvent == 0) {
243 G4cout.precision(dfprec);
244 return;
245 }
246
247 // compute mean Energy deposited and rms in target
248 //
249 G4int TotNbofEvents = numberOfEvent;
250 fEdepTarget /= TotNbofEvents;
251 fEdepTarget2 /= TotNbofEvents;
252 G4double rmsEdep = fEdepTarget2 - fEdepTarget * fEdepTarget;
253 if (rmsEdep > 0.)
254 rmsEdep = std::sqrt(rmsEdep);
255 else
256 rmsEdep = 0.;
257
258 G4cout << "\n Mean energy deposit in target, in time window = "
259 << G4BestUnit(fEdepTarget, "Energy") << "; rms = " << G4BestUnit(rmsEdep, "Energy")
260 << G4endl;
261
262 // compute mean Energy deposited and rms in detector
263 //
264 fEdepDetect /= TotNbofEvents;
265 fEdepDetect2 /= TotNbofEvents;
266 rmsEdep = fEdepDetect2 - fEdepDetect * fEdepDetect;
267 if (rmsEdep > 0.)
268 rmsEdep = std::sqrt(rmsEdep);
269 else
270 rmsEdep = 0.;
271
272 G4cout << " Mean energy deposit in detector, in time window = "
273 << G4BestUnit(fEdepDetect, "Energy") << "; rms = " << G4BestUnit(rmsEdep, "Energy")
274 << G4endl;
275
276 // frequency of processes in target
277 //
278 G4cout << "\n Process calls frequency in target :" << G4endl;
279 G4int index = 0;
280 std::map<G4String, G4int>::iterator it1;
281 for (it1 = fProcCounter1.begin(); it1 != fProcCounter1.end(); it1++) {
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 << "=" << std::setw(7) << count << space;
287 }
288 G4cout << G4endl;
289
290 // frequency of processes in detector
291 //
292 G4cout << "\n Process calls frequency in detector:" << G4endl;
293 index = 0;
294 std::map<G4String, G4int>::iterator it2;
295 for (it2 = fProcCounter2.begin(); it2 != fProcCounter2.end(); it2++) {
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 << "=" << std::setw(7) << count << space;
301 }
302 G4cout << G4endl;
303
304 // particles count in target
305 //
306 G4cout << "\n List of generated particles in target:" << G4endl;
307
308 std::map<G4String, ParticleData>::iterator itc;
309 for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); 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 << ": " << std::setw(7) << count
318 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
319 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
320 }
321
322 // particles count in detector
323 //
324 G4cout << "\n List of generated particles in detector:" << G4endl;
325
326 std::map<G4String, ParticleData>::iterator itn;
327 for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); 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 << ": " << std::setw(7) << count
336 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
337 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
338 }
339 G4cout << G4endl;
340
341 // activities in VR mode
342 //
343 WriteActivity(numberOfEvent);
344
345 // remove all contents in fProcCounter, fCount
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........oooOO0OOooo........oooOO0OOooo......
356
357 void Run::WriteActivity(G4int nevent)
358 {
359 G4ProcessTable* pTable = G4ProcessTable::GetProcessTable();
360 G4Radioactivation* rDecay =
361 (G4Radioactivation*)pTable->FindProcess("Radioactivation", "GenericIon");
362
363 // output the induced radioactivities (in VR mode only)
364 //
365 if ((rDecay == 0) || (rDecay->IsAnalogueMonteCarlo())) return;
366
367 G4String fileName = G4AnalysisManager::Instance()->GetFileName() + ".activity";
368 std::ofstream outfile(fileName, std::ios::out);
369
370 std::vector<G4RadioactivityTable*> theTables = rDecay->GetTheRadioactivityTables();
371
372 for (size_t i = 0; i < theTables.size(); i++) {
373 G4double rate, error;
374 outfile << "Radioactivities in decay window no. " << i << G4endl;
375 outfile << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) " << G4endl;
376
377 map<G4ThreeVector, G4TwoVector>* aMap = theTables[i]->GetTheMap();
378 map<G4ThreeVector, G4TwoVector>::iterator iter;
379 for (iter = aMap->begin(); iter != aMap->end(); iter++) {
380 rate = iter->second.x() / nevent;
381 error = std::sqrt(iter->second.y()) / nevent;
382 if (rate < 0.) rate = 0.; // statically it can be < 0.
383 outfile << iter->first.x() << "\t" << iter->first.y() << "\t" << iter->first.z() << "\t"
384 << rate << "\t" << error << G4endl;
385 }
386 outfile << G4endl;
387 }
388 outfile.close();
389 }
390
391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392