| 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 electromagnetic/TestEm11/src/Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // << 29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
>> 30 //
>> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 33
31 #include "Run.hh" 34 #include "Run.hh"
32 <<
33 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
34 #include "HistoManager.hh" <<
35 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
>> 37 #include "HistoManager.hh"
36 38
37 #include "G4ProcessTable.hh" 39 #include "G4ProcessTable.hh"
38 #include "G4Radioactivation.hh" << 40 #include "G4RadioactiveDecay.hh"
39 #include "G4SystemOfUnits.hh" <<
40 #include "G4TwoVector.hh" 41 #include "G4TwoVector.hh"
41 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh"
>> 43 #include "G4SystemOfUnits.hh"
42 44
43 #include <fstream> 45 #include <fstream>
44 46
45 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 48
47 Run::Run(DetectorConstruction* det) : fDetecto << 49 Run::Run(DetectorConstruction* det)
>> 50 : G4Run(),
>> 51 fDetector(det), fParticle(0), fEkin(0.)
>> 52 {
>> 53 fEdepTarget = fEdepTarget2 = 0.;
>> 54 fEdepDetect = fEdepDetect2 = 0.;
>> 55 }
>> 56
>> 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 58
>> 59 Run::~Run()
>> 60 { }
48 61
49 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 63
51 void Run::SetPrimary(G4ParticleDefinition* par 64 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
52 { << 65 {
53 fParticle = particle; 66 fParticle = particle;
54 fEkin = energy; 67 fEkin = energy;
55 } 68 }
56 << 69
57 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 71
59 void Run::CountProcesses(const G4VProcess* pro << 72 void Run::CountProcesses(const G4VProcess* process, G4int iVol)
60 { 73 {
61 if (process == nullptr) return; <<
62 G4String procName = process->GetProcessName( 74 G4String procName = process->GetProcessName();
>> 75
63 if (iVol == 1) { 76 if (iVol == 1) {
64 std::map<G4String, G4int>::iterator it1 = << 77 std::map<G4String,G4int>::iterator it1 = fProcCounter1.find(procName);
65 if (it1 == fProcCounter1.end()) { << 78 if ( it1 == fProcCounter1.end()) {
66 fProcCounter1[procName] = 1; 79 fProcCounter1[procName] = 1;
67 } 80 }
68 else { 81 else {
69 fProcCounter1[procName]++; << 82 fProcCounter1[procName]++;
70 } 83 }
71 } 84 }
72 << 85
73 if (iVol == 2) { 86 if (iVol == 2) {
74 std::map<G4String, G4int>::iterator it2 = << 87 std::map<G4String,G4int>::iterator it2 = fProcCounter2.find(procName);
75 if (it2 == fProcCounter2.end()) { << 88 if ( it2 == fProcCounter2.end()) {
76 fProcCounter2[procName] = 1; 89 fProcCounter2[procName] = 1;
77 } 90 }
78 else { 91 else {
79 fProcCounter2[procName]++; << 92 fProcCounter2[procName]++;
80 } << 93 }
81 } 94 }
82 } 95 }
83 << 96
84 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 98
86 void Run::ParticleCount(G4String name, G4doubl 99 void Run::ParticleCount(G4String name, G4double Ekin, G4int iVol)
87 { 100 {
88 if (iVol == 1) { 101 if (iVol == 1) {
89 std::map<G4String, ParticleData>::iterator << 102 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
90 if (it == fParticleDataMap1.end()) { << 103 if ( it == fParticleDataMap1.end()) {
91 fParticleDataMap1[name] = ParticleData(1 << 104 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
92 } << 105 }
93 else { << 106 else {
94 ParticleData& data = it->second; << 107 ParticleData& data = it->second;
95 data.fCount++; << 108 data.fCount++;
96 data.fEmean += Ekin; << 109 data.fEmean += Ekin;
97 // update min max << 110 //update min max
98 G4double emin = data.fEmin; << 111 G4double emin = data.fEmin;
99 if (Ekin < emin) data.fEmin = Ekin; << 112 if (Ekin < emin) data.fEmin = Ekin;
100 G4double emax = data.fEmax; << 113 G4double emax = data.fEmax;
101 if (Ekin > emax) data.fEmax = Ekin; << 114 if (Ekin > emax) data.fEmax = Ekin;
102 } << 115 }
103 } 116 }
104 << 117
105 if (iVol == 2) { 118 if (iVol == 2) {
106 std::map<G4String, ParticleData>::iterator << 119 std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
107 if (it == fParticleDataMap2.end()) { << 120 if ( it == fParticleDataMap2.end()) {
108 fParticleDataMap2[name] = ParticleData(1 << 121 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
109 } << 122 }
110 else { << 123 else {
111 ParticleData& data = it->second; << 124 ParticleData& data = it->second;
112 data.fCount++; << 125 data.fCount++;
113 data.fEmean += Ekin; << 126 data.fEmean += Ekin;
114 // update min max << 127 //update min max
115 G4double emin = data.fEmin; << 128 G4double emin = data.fEmin;
116 if (Ekin < emin) data.fEmin = Ekin; << 129 if (Ekin < emin) data.fEmin = Ekin;
117 G4double emax = data.fEmax; << 130 G4double emax = data.fEmax;
118 if (Ekin > emax) data.fEmax = Ekin; << 131 if (Ekin > emax) data.fEmax = Ekin;
119 } << 132 }
120 } << 133 }
121 } 134 }
122 << 135
123 //....oooOO0OOooo........oooOO0OOooo........oo 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 137
125 void Run::AddEdep(G4double edep1, G4double ede 138 void Run::AddEdep(G4double edep1, G4double edep2)
126 { << 139 {
127 fEdepTarget += edep1; << 140 fEdepTarget += edep1;
128 fEdepTarget2 += edep1 * edep1; << 141 fEdepTarget2 += edep1*edep1;
129 fEdepDetect += edep2; << 142 fEdepDetect += edep2;
130 fEdepDetect2 += edep2 * edep2; << 143 fEdepDetect2 += edep2*edep2;
131 } 144 }
132 << 145
133 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 147
135 void Run::Merge(const G4Run* run) 148 void Run::Merge(const G4Run* run)
136 { 149 {
137 const Run* localRun = static_cast<const Run* 150 const Run* localRun = static_cast<const Run*>(run);
138 << 151
139 // primary particle info << 152 //primary particle info
140 // 153 //
141 fParticle = localRun->fParticle; 154 fParticle = localRun->fParticle;
142 fEkin = localRun->fEkin; << 155 fEkin = localRun->fEkin;
143 << 156
144 // accumulate sums 157 // accumulate sums
145 // 158 //
146 fEdepTarget += localRun->fEdepTarget; << 159 fEdepTarget += localRun->fEdepTarget;
147 fEdepTarget2 += localRun->fEdepTarget2; << 160 fEdepTarget2 += localRun->fEdepTarget2;
148 fEdepDetect += localRun->fEdepDetect; << 161 fEdepDetect += localRun->fEdepDetect;
149 fEdepDetect2 += localRun->fEdepDetect2; << 162 fEdepDetect2 += localRun->fEdepDetect2;
>> 163
>> 164 //map: processes count in target
>> 165
>> 166 std::map<G4String,G4int>::const_iterator itp1;
>> 167 for ( itp1 = localRun->fProcCounter1.begin();
>> 168 itp1 != localRun->fProcCounter1.end(); ++itp1 ) {
150 169
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; 170 G4String procName = itp1->first;
156 G4int localCount = itp1->second; 171 G4int localCount = itp1->second;
157 if (fProcCounter1.find(procName) == fProcC << 172 if ( fProcCounter1.find(procName) == fProcCounter1.end()) {
158 fProcCounter1[procName] = localCount; 173 fProcCounter1[procName] = localCount;
159 } 174 }
160 else { 175 else {
161 fProcCounter1[procName] += localCount; 176 fProcCounter1[procName] += localCount;
162 } << 177 }
163 } 178 }
>> 179
>> 180 //map: processes count in detector
>> 181
>> 182 std::map<G4String,G4int>::const_iterator itp2;
>> 183 for ( itp2 = localRun->fProcCounter2.begin();
>> 184 itp2 != localRun->fProcCounter2.end(); ++itp2 ) {
164 185
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; 186 G4String procName = itp2->first;
170 G4int localCount = itp2->second; 187 G4int localCount = itp2->second;
171 if (fProcCounter2.find(procName) == fProcC << 188 if ( fProcCounter2.find(procName) == fProcCounter2.end()) {
172 fProcCounter2[procName] = localCount; 189 fProcCounter2[procName] = localCount;
173 } 190 }
174 else { 191 else {
175 fProcCounter2[procName] += localCount; 192 fProcCounter2[procName] += localCount;
176 } << 193 }
177 } 194 }
178 << 195
179 // map: created particles in target << 196 //map: created particles in target
180 std::map<G4String, ParticleData>::const_iter << 197 std::map<G4String,ParticleData>::const_iterator itc;
181 for (itc = localRun->fParticleDataMap1.begin << 198 for (itc = localRun->fParticleDataMap1.begin();
>> 199 itc != localRun->fParticleDataMap1.end(); ++itc) {
>> 200
182 G4String name = itc->first; 201 G4String name = itc->first;
183 const ParticleData& localData = itc->secon << 202 const ParticleData& localData = itc->second;
184 if (fParticleDataMap1.find(name) == fParti << 203 if ( fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
185 fParticleDataMap1[name] = << 204 fParticleDataMap1[name]
186 ParticleData(localData.fCount, localDa << 205 = ParticleData(localData.fCount,
>> 206 localData.fEmean,
>> 207 localData.fEmin,
>> 208 localData.fEmax);
187 } 209 }
188 else { 210 else {
189 ParticleData& data = fParticleDataMap1[n << 211 ParticleData& data = fParticleDataMap1[name];
190 data.fCount += localData.fCount; 212 data.fCount += localData.fCount;
191 data.fEmean += localData.fEmean; 213 data.fEmean += localData.fEmean;
192 G4double emin = localData.fEmin; 214 G4double emin = localData.fEmin;
193 if (emin < data.fEmin) data.fEmin = emin 215 if (emin < data.fEmin) data.fEmin = emin;
194 G4double emax = localData.fEmax; 216 G4double emax = localData.fEmax;
195 if (emax > data.fEmax) data.fEmax = emax << 217 if (emax > data.fEmax) data.fEmax = emax;
196 } << 218 }
197 } 219 }
198 << 220
199 // map: created particle in detector << 221 //map: created particle in detector
200 std::map<G4String, ParticleData>::const_iter << 222 std::map<G4String,ParticleData>::const_iterator itn;
201 for (itn = localRun->fParticleDataMap2.begin << 223 for (itn = localRun->fParticleDataMap2.begin();
>> 224 itn != localRun->fParticleDataMap2.end(); ++itn) {
>> 225
202 G4String name = itn->first; 226 G4String name = itn->first;
203 const ParticleData& localData = itn->secon << 227 const ParticleData& localData = itn->second;
204 if (fParticleDataMap2.find(name) == fParti << 228 if ( fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
205 fParticleDataMap2[name] = << 229 fParticleDataMap2[name]
206 ParticleData(localData.fCount, localDa << 230 = ParticleData(localData.fCount,
>> 231 localData.fEmean,
>> 232 localData.fEmin,
>> 233 localData.fEmax);
207 } 234 }
208 else { 235 else {
209 ParticleData& data = fParticleDataMap2[n << 236 ParticleData& data = fParticleDataMap2[name];
210 data.fCount += localData.fCount; 237 data.fCount += localData.fCount;
211 data.fEmean += localData.fEmean; 238 data.fEmean += localData.fEmean;
212 G4double emin = localData.fEmin; 239 G4double emin = localData.fEmin;
213 if (emin < data.fEmin) data.fEmin = emin 240 if (emin < data.fEmin) data.fEmin = emin;
214 G4double emax = localData.fEmax; 241 G4double emax = localData.fEmax;
215 if (emax > data.fEmax) data.fEmax = emax << 242 if (emax > data.fEmax) data.fEmax = emax;
216 } << 243 }
217 } 244 }
218 245
219 G4Run::Merge(run); << 246 G4Run::Merge(run);
220 } << 247 }
221 248
222 //....oooOO0OOooo........oooOO0OOooo........oo 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 250
224 void Run::EndOfRun() << 251 void Run::EndOfRun()
225 { 252 {
226 G4int prec = 5, wid = prec + 2; << 253 G4int prec = 5, wid = prec + 2;
227 G4int dfprec = G4cout.precision(prec); 254 G4int dfprec = G4cout.precision(prec);
228 << 255
229 // run condition 256 // run condition
230 // << 257 //
231 G4String Particle = fParticle->GetParticleNa << 258 G4String Particle = fParticle->GetParticleName();
232 G4cout << "\n The run is " << numberOfEvent << 259 G4cout << "\n The run is " << numberOfEvent << " "<< Particle << " of "
233 << G4BestUnit(fEkin, "Energy") << " t << 260 << G4BestUnit(fEkin,"Energy") << " through : ";
234 << 261
235 G4cout << "\n Target : Length = " << G4Bes << 262 G4cout << "\n Target : Length = "
236 << " Radius = " << G4BestUnit(fDet << 263 << G4BestUnit(fDetector->GetTargetLength(),"Length")
237 << " Material = " << fDetector->GetTa << 264 << " Radius = "
238 G4cout << "\n Detector : Length = " << G4Bes << 265 << G4BestUnit(fDetector->GetTargetRadius(),"Length")
239 << " Thickness = " << G4BestUnit(fDet << 266 << " Material = "
240 << " Material = " << fDetector->GetDe << 267 << fDetector->GetTargetMaterial()->GetName();
241 << 268 G4cout << "\n Detector : Length = "
242 if (numberOfEvent == 0) { << 269 << G4BestUnit(fDetector->GetDetectorLength(),"Length")
243 G4cout.precision(dfprec); << 270 << " Thickness = "
244 return; << 271 << G4BestUnit(fDetector->GetDetectorThickness(),"Length")
245 } << 272 << " Material = "
>> 273 << fDetector->GetDetectorMaterial()->GetName() << G4endl;
246 274
>> 275 if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
>> 276
247 // compute mean Energy deposited and rms in 277 // compute mean Energy deposited and rms in target
248 // 278 //
249 G4int TotNbofEvents = numberOfEvent; 279 G4int TotNbofEvents = numberOfEvent;
250 fEdepTarget /= TotNbofEvents; << 280 fEdepTarget /= TotNbofEvents; fEdepTarget2 /= TotNbofEvents;
251 fEdepTarget2 /= TotNbofEvents; << 281 G4double rmsEdep = fEdepTarget2 - fEdepTarget*fEdepTarget;
252 G4double rmsEdep = fEdepTarget2 - fEdepTarge << 282 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
253 if (rmsEdep > 0.) << 283 else rmsEdep = 0.;
254 rmsEdep = std::sqrt(rmsEdep); << 284
255 else <<
256 rmsEdep = 0.; <<
257 <<
258 G4cout << "\n Mean energy deposit in target, 285 G4cout << "\n Mean energy deposit in target, in time window = "
259 << G4BestUnit(fEdepTarget, "Energy") << 286 << G4BestUnit(fEdepTarget,"Energy") << "; rms = "
>> 287 << G4BestUnit(rmsEdep, "Energy")
260 << G4endl; 288 << G4endl;
261 289
262 // compute mean Energy deposited and rms in 290 // compute mean Energy deposited and rms in detector
263 // 291 //
264 fEdepDetect /= TotNbofEvents; << 292 fEdepDetect /= TotNbofEvents; fEdepDetect2 /= TotNbofEvents;
265 fEdepDetect2 /= TotNbofEvents; << 293 rmsEdep = fEdepDetect2 - fEdepDetect*fEdepDetect;
266 rmsEdep = fEdepDetect2 - fEdepDetect * fEdep << 294 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
267 if (rmsEdep > 0.) << 295 else rmsEdep = 0.;
268 rmsEdep = std::sqrt(rmsEdep); << 296
269 else <<
270 rmsEdep = 0.; <<
271 <<
272 G4cout << " Mean energy deposit in detector, 297 G4cout << " Mean energy deposit in detector, in time window = "
273 << G4BestUnit(fEdepDetect, "Energy") << 298 << G4BestUnit(fEdepDetect,"Energy") << "; rms = "
>> 299 << G4BestUnit(rmsEdep, "Energy")
274 << G4endl; 300 << G4endl;
275 301
276 // frequency of processes in target 302 // frequency of processes in target
277 // 303 //
278 G4cout << "\n Process calls frequency in tar 304 G4cout << "\n Process calls frequency in target :" << G4endl;
279 G4int index = 0; 305 G4int index = 0;
280 std::map<G4String, G4int>::iterator it1; << 306 std::map<G4String,G4int>::iterator it1;
281 for (it1 = fProcCounter1.begin(); it1 != fPr 307 for (it1 = fProcCounter1.begin(); it1 != fProcCounter1.end(); it1++) {
282 G4String procName = it1->first; << 308 G4String procName = it1->first;
283 G4int count = it1->second; << 309 G4int count = it1->second;
284 G4String space = " "; << 310 G4String space = " "; if (++index%3 == 0) space = "\n";
285 if (++index % 3 == 0) space = "\n"; << 311 G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
286 G4cout << " " << std::setw(20) << procName << 312 << space;
287 } 313 }
288 G4cout << G4endl; 314 G4cout << G4endl;
289 << 315
290 // frequency of processes in detector 316 // frequency of processes in detector
291 // 317 //
292 G4cout << "\n Process calls frequency in det 318 G4cout << "\n Process calls frequency in detector:" << G4endl;
293 index = 0; 319 index = 0;
294 std::map<G4String, G4int>::iterator it2; << 320 std::map<G4String,G4int>::iterator it2;
295 for (it2 = fProcCounter2.begin(); it2 != fPr 321 for (it2 = fProcCounter2.begin(); it2 != fProcCounter2.end(); it2++) {
296 G4String procName = it2->first; << 322 G4String procName = it2->first;
297 G4int count = it2->second; << 323 G4int count = it2->second;
298 G4String space = " "; << 324 G4String space = " "; if (++index%3 == 0) space = "\n";
299 if (++index % 3 == 0) space = "\n"; << 325 G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
300 G4cout << " " << std::setw(20) << procName << 326 << space;
301 } 327 }
302 G4cout << G4endl; 328 G4cout << G4endl;
303 << 329
304 // particles count in target 330 // particles count in target
305 // 331 //
306 G4cout << "\n List of generated particles in 332 G4cout << "\n List of generated particles in target:" << G4endl;
307 << 333
308 std::map<G4String, ParticleData>::iterator i << 334 std::map<G4String,ParticleData>::iterator itc;
309 for (itc = fParticleDataMap1.begin(); itc != 335 for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
310 G4String name = itc->first; << 336 G4String name = itc->first;
311 ParticleData data = itc->second; << 337 ParticleData data = itc->second;
312 G4int count = data.fCount; << 338 G4int count = data.fCount;
313 G4double eMean = data.fEmean / count; << 339 G4double eMean = data.fEmean/count;
314 G4double eMin = data.fEmin; << 340 G4double eMin = data.fEmin;
315 G4double eMax = data.fEmax; << 341 G4double eMax = data.fEmax;
316 << 342
317 G4cout << " " << std::setw(13) << name << 343 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
318 << " Emean = " << std::setw(wid) < << 344 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
319 << G4BestUnit(eMin, "Energy") << " << 345 << "\t( " << G4BestUnit(eMin, "Energy")
320 } << 346 << " --> " << G4BestUnit(eMax, "Energy")
321 << 347 << ")" << G4endl;
322 // particles count in detector << 348 }
323 // << 349
324 G4cout << "\n List of generated particles in << 350 // particles count in detector
325 << 351 //
326 std::map<G4String, ParticleData>::iterator i << 352 G4cout << "\n List of generated particles in detector:" << G4endl;
327 for (itn = fParticleDataMap2.begin(); itn != << 353
>> 354 std::map<G4String,ParticleData>::iterator itn;
>> 355 for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
328 G4String name = itn->first; 356 G4String name = itn->first;
329 ParticleData data = itn->second; 357 ParticleData data = itn->second;
330 G4int count = data.fCount; 358 G4int count = data.fCount;
331 G4double eMean = data.fEmean / count; << 359 G4double eMean = data.fEmean/count;
332 G4double eMin = data.fEmin; 360 G4double eMin = data.fEmin;
333 G4double eMax = data.fEmax; 361 G4double eMax = data.fEmax;
334 << 362
335 G4cout << " " << std::setw(13) << name << 363 G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
336 << " Emean = " << std::setw(wid) < << 364 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
337 << G4BestUnit(eMin, "Energy") << " << 365 << "\t( " << G4BestUnit(eMin, "Energy")
>> 366 << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
338 } 367 }
339 G4cout << G4endl; 368 G4cout << G4endl;
340 << 369
341 // activities in VR mode 370 // activities in VR mode
342 // 371 //
343 WriteActivity(numberOfEvent); 372 WriteActivity(numberOfEvent);
344 << 373
345 // remove all contents in fProcCounter, fCou << 374 //remove all contents in fProcCounter, fCount
346 fProcCounter1.clear(); 375 fProcCounter1.clear();
347 fProcCounter2.clear(); 376 fProcCounter2.clear();
348 fParticleDataMap1.clear(); << 377 fParticleDataMap1.clear();
349 fParticleDataMap2.clear(); 378 fParticleDataMap2.clear();
350 << 379
351 // restore default format << 380 //restore default format
352 G4cout.precision(dfprec); << 381 G4cout.precision(dfprec);
353 } 382 }
354 383
355 //....oooOO0OOooo........oooOO0OOooo........oo 384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
356 385
357 void Run::WriteActivity(G4int nevent) 386 void Run::WriteActivity(G4int nevent)
358 { 387 {
359 G4ProcessTable* pTable = G4ProcessTable::Get << 388 G4ProcessTable *pTable = G4ProcessTable::GetProcessTable();
360 G4Radioactivation* rDecay = << 389 G4RadioactiveDecay * rDecay = (G4RadioactiveDecay *)
361 (G4Radioactivation*)pTable->FindProcess("R << 390 pTable->FindProcess("RadioactiveDecay", "GenericIon");
362 << 391
363 // output the induced radioactivities (in VR << 392 // output the induced radioactivities (in VR mode only)
364 // << 393 //
365 if ((rDecay == 0) || (rDecay->IsAnalogueMont << 394 if ((rDecay == 0) || (rDecay->IsAnalogueMonteCarlo())) return;
366 << 395
367 G4String fileName = G4AnalysisManager::Insta << 396 G4String fileName = G4AnalysisManager::Instance()->GetFileName() + ".activity";
368 std::ofstream outfile(fileName, std::ios::ou << 397 std::ofstream outfile (fileName, std::ios::out );
>> 398
>> 399 std::vector<G4RadioactivityTable*> theTables =
>> 400 rDecay->GetTheRadioactivityTables();
369 401
370 std::vector<G4RadioactivityTable*> theTables << 402 for (size_t i = 0 ; i < theTables.size(); i++) {
371 <<
372 for (size_t i = 0; i < theTables.size(); i++ <<
373 G4double rate, error; 403 G4double rate, error;
374 outfile << "Radioactivities in decay windo 404 outfile << "Radioactivities in decay window no. " << i << G4endl;
375 outfile << "Z \tA \tE \tActivity (decays/w << 405 outfile << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) "
>> 406 << G4endl;
376 407
377 map<G4ThreeVector, G4TwoVector>* aMap = th << 408 map<G4ThreeVector,G4TwoVector> *aMap = theTables[i]->GetTheMap();
378 map<G4ThreeVector, G4TwoVector>::iterator << 409 map<G4ThreeVector,G4TwoVector>::iterator iter;
379 for (iter = aMap->begin(); iter != aMap->e << 410 for (iter=aMap->begin(); iter != aMap->end(); iter++) {
380 rate = iter->second.x() / nevent; << 411 rate = iter->second.x()/nevent;
381 error = std::sqrt(iter->second.y()) / ne << 412 error = std::sqrt(iter->second.y())/nevent;
382 if (rate < 0.) rate = 0.; // statically << 413 if (rate < 0.) rate = 0.; // statically it can be < 0.
383 outfile << iter->first.x() << "\t" << it << 414 outfile << iter->first.x() <<"\t"<< iter->first.y() <<"\t"
384 << rate << "\t" << error << G4en << 415 << iter->first.z() << "\t" << rate <<"\t" << error << G4endl;
385 } 416 }
386 outfile << G4endl; 417 outfile << G4endl;
387 } << 418 }
388 outfile.close(); << 419 outfile.close();
389 } 420 }
390 421
391 //....oooOO0OOooo........oooOO0OOooo........oo 422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392 423