| 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 electromagnetic/TestEm18/src/RunActi 26 /// \file electromagnetic/TestEm18/src/RunAction.cc
27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class
28 // 28 //
>> 29 // $Id$
29 // 30 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "RunAction.hh" 34 #include "RunAction.hh"
34 <<
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh" <<
37 #include "PrimaryGeneratorAction.hh" 36 #include "PrimaryGeneratorAction.hh"
>> 37 #include "HistoManager.hh"
38 38
39 #include "G4EmCalculator.hh" <<
40 #include "G4Run.hh" 39 #include "G4Run.hh"
>> 40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh"
42 #include "Randomize.hh" << 42 #include "G4EmCalculator.hh"
43 43
>> 44 #include "Randomize.hh"
44 #include <iomanip> 45 #include <iomanip>
45 46
46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 48
48 RunAction::RunAction(DetectorConstruction* det 49 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
49 : fDetector(det), fPrimary(kin) << 50 :fDetector(det), fPrimary(kin)
50 { << 51 {
51 fHistoManager = new HistoManager(); << 52 fHistoManager = new HistoManager();
52 } 53 }
53 54
54 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 56
56 RunAction::~RunAction() 57 RunAction::~RunAction()
57 { << 58 {
58 delete fHistoManager; << 59 delete fHistoManager;
59 } 60 }
60 61
61 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 63
63 void RunAction::BeginOfRunAction(const G4Run*) << 64 void RunAction::BeginOfRunAction(const G4Run* run)
64 { 65 {
65 // initialisation << 66 G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
>> 67
>> 68 //initialisation
66 // 69 //
>> 70 fEnergyDeposit = 0.;
>> 71
>> 72 fNbCharged = fNbNeutral = 0;
>> 73 fEnergyCharged = fEnergyNeutral = 0.;
>> 74 fEmin[0] = fEmin[1] = DBL_MAX;
>> 75 fEmax[0] = fEmax[1] = 0.;
>> 76
67 fNbSteps = 0; 77 fNbSteps = 0;
68 fTrackLength = 0.; 78 fTrackLength = 0.;
69 fStepMin = DBL_MAX; << 79
70 fStepMax = 0.; << 80 //histograms
71 <<
72 fEdepPrimary = fEdepSecondary = fEdepTotal = <<
73 fEdepPrimMin = fEdepSecMin = fEdepTotMin = D <<
74 fEdepPrimMax = fEdepSecMax = fEdepTotMax = 0 <<
75 <<
76 fEnergyTransfered = 0.; <<
77 fEtransfMin = DBL_MAX; <<
78 fEtransfMax = 0.; <<
79 <<
80 fEnergyLost = 0.; <<
81 fElostMin = DBL_MAX; <<
82 fElostMax = 0.; <<
83 <<
84 fEnergyBalance = 0.; <<
85 fEbalMin = DBL_MAX; <<
86 fEbalMax = 0.; <<
87 <<
88 // histograms <<
89 // 81 //
90 G4AnalysisManager* analysisManager = G4Analy 82 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
91 if (analysisManager->IsActive()) { << 83 if ( analysisManager->IsActive() ) {
92 analysisManager->OpenFile(); 84 analysisManager->OpenFile();
93 } << 85 }
94 86
95 // show Rndm status << 87 // do not save Rndm status
>> 88 G4RunManager::GetRunManager()->SetRandomNumberStore(false);
96 CLHEP::HepRandom::showEngineStatus(); 89 CLHEP::HepRandom::showEngineStatus();
97 } 90 }
98 91
99 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 93
101 void RunAction::CountProcesses(G4String procNa <<
102 { <<
103 std::map<G4String, G4int>::iterator it = fPr <<
104 if (it == fProcCounter.end()) { <<
105 fProcCounter[procName] = 1; <<
106 } <<
107 else { <<
108 fProcCounter[procName]++; <<
109 } <<
110 } <<
111 <<
112 //....oooOO0OOooo........oooOO0OOooo........oo <<
113 <<
114 void RunAction::TrackLength(G4double step) <<
115 { <<
116 fTrackLength += step; <<
117 fNbSteps++; <<
118 if (step < fStepMin) fStepMin = step; <<
119 if (step > fStepMax) fStepMax = step; <<
120 } <<
121 <<
122 //....oooOO0OOooo........oooOO0OOooo........oo <<
123 <<
124 void RunAction::EnergyDeposited(G4double edepP <<
125 { <<
126 fEdepPrimary += edepPrim; <<
127 if (edepPrim < fEdepPrimMin) fEdepPrimMin = <<
128 if (edepPrim > fEdepPrimMax) fEdepPrimMax = <<
129 <<
130 fEdepSecondary += edepSecond; <<
131 if (edepSecond < fEdepSecMin) fEdepSecMin = <<
132 if (edepSecond > fEdepSecMax) fEdepSecMax = <<
133 } <<
134 <<
135 //....oooOO0OOooo........oooOO0OOooo........oo <<
136 <<
137 void RunAction::EnergyTransferedByProcess(G4St <<
138 { <<
139 std::map<G4String, MinMaxData>::iterator it <<
140 if (it == fEtransfByProcess.end()) { <<
141 fEtransfByProcess[process] = MinMaxData(1, <<
142 } <<
143 else { <<
144 MinMaxData& data = it->second; <<
145 data.fCount++; <<
146 data.fVsum += energy; <<
147 // update min max <<
148 G4double emin = data.fVmin; <<
149 if (energy < emin) data.fVmin = energy; <<
150 G4double emax = data.fVmax; <<
151 if (energy > emax) data.fVmax = energy; <<
152 } <<
153 } <<
154 <<
155 //....oooOO0OOooo........oooOO0OOooo........oo <<
156 <<
157 void RunAction::EnergyTransfered(G4double ener <<
158 { <<
159 fEnergyTransfered += energy; <<
160 if (energy < fEtransfMin) fEtransfMin = ener <<
161 if (energy > fEtransfMax) fEtransfMax = ener <<
162 } <<
163 <<
164 //....oooOO0OOooo........oooOO0OOooo........oo <<
165 <<
166 void RunAction::TotalEnergyLost(G4double energ <<
167 { <<
168 fEnergyLost += energy; <<
169 if (energy < fElostMin) fElostMin = energy; <<
170 if (energy > fElostMax) fElostMax = energy; <<
171 } <<
172 <<
173 //....oooOO0OOooo........oooOO0OOooo........oo <<
174 <<
175 void RunAction::EnergyBalance(G4double energy) <<
176 { <<
177 fEnergyBalance += energy; <<
178 if (energy < fEbalMin) fEbalMin = energy; <<
179 if (energy > fEbalMax) fEbalMax = energy; <<
180 } <<
181 <<
182 //....oooOO0OOooo........oooOO0OOooo........oo <<
183 <<
184 void RunAction::TotalEnergyDeposit(G4double en <<
185 { <<
186 fEdepTotal += energy; <<
187 if (energy < fEdepTotMin) fEdepTotMin = ener <<
188 if (energy > fEdepTotMax) fEdepTotMax = ener <<
189 } <<
190 <<
191 //....oooOO0OOooo........oooOO0OOooo........oo <<
192 <<
193 void RunAction::EnergySpectrumOfSecondaries(G4 <<
194 { <<
195 std::map<G4String, MinMaxData>::iterator it <<
196 if (it == fEkinOfSecondaries.end()) { <<
197 fEkinOfSecondaries[particle] = MinMaxData( <<
198 } <<
199 else { <<
200 MinMaxData& data = it->second; <<
201 data.fCount++; <<
202 data.fVsum += energy; <<
203 // update min max <<
204 G4double emin = data.fVmin; <<
205 if (energy < emin) data.fVmin = energy; <<
206 G4double emax = data.fVmax; <<
207 if (energy > emax) data.fVmax = energy; <<
208 } <<
209 } <<
210 <<
211 //....oooOO0OOooo........oooOO0OOooo........oo <<
212 <<
213 void RunAction::EndOfRunAction(const G4Run* aR 94 void RunAction::EndOfRunAction(const G4Run* aRun)
214 { 95 {
215 G4int nbEvents = aRun->GetNumberOfEvent(); 96 G4int nbEvents = aRun->GetNumberOfEvent();
216 if (nbEvents == 0) return; 97 if (nbEvents == 0) return;
217 << 98
218 G4Material* material = fDetector->GetMateria 99 G4Material* material = fDetector->GetMaterial();
219 G4double length = fDetector->GetSize(); << 100 G4double length = fDetector->GetSize();
220 G4double density = material->GetDensity(); 101 G4double density = material->GetDensity();
221 << 102
222 G4ParticleDefinition* particle = fPrimary->G << 103 G4ParticleDefinition* particle = fPrimary->GetParticleGun()
>> 104 ->GetParticleDefinition();
223 G4String partName = particle->GetParticleNam 105 G4String partName = particle->GetParticleName();
224 G4double ePrimary = fPrimary->GetParticleGun 106 G4double ePrimary = fPrimary->GetParticleGun()->GetParticleEnergy();
225 << 107
226 G4int prec = G4cout.precision(3); 108 G4int prec = G4cout.precision(3);
227 G4cout << "\n ======================== run s << 109 G4cout << "\n ======================== run summary ======================\n";
228 G4cout << "\n The run was " << nbEvents << " 110 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
229 << G4BestUnit(ePrimary, "Energy") << << 111 << G4BestUnit(ePrimary,"Energy") << " through "
230 << material->GetName() << " (density: << 112 << G4BestUnit(length,"Length") << " of "
>> 113 << material->GetName() << " (density: "
>> 114 << G4BestUnit(density,"Volumic Mass") << ")";
>> 115 G4cout << "\n ===========================================================\n";
231 G4cout << G4endl; 116 G4cout << G4endl;
232 << 117
>> 118 //save histograms
>> 119 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
>> 120 if ( analysisManager->IsActive() ) {
>> 121 analysisManager->Write();
>> 122 analysisManager->CloseFile();
>> 123 }
>> 124
233 if (particle->GetPDGCharge() == 0.) return; 125 if (particle->GetPDGCharge() == 0.) return;
>> 126
>> 127 G4cout.precision(5);
>> 128
>> 129 //track length
>> 130 //
>> 131 G4double trackLPerEvent = fTrackLength/nbEvents;
>> 132 G4double nbStepPerEvent = double(fNbSteps)/nbEvents;
>> 133 G4double stepSize = fTrackLength/fNbSteps;
>> 134
>> 135 G4cout
>> 136 << "\n TrackLength= "
>> 137 << G4BestUnit(trackLPerEvent, "Length")
>> 138 << "\t nb of steps= " << nbStepPerEvent
>> 139 << " stepSize= " << G4BestUnit(stepSize, "Length")
>> 140 << G4endl;
>> 141
>> 142 //charged secondaries (ionization, direct pair production)
>> 143 //
>> 144 G4double energyPerEvent = fEnergyCharged/nbEvents;
>> 145 G4double nbPerEvent = double(fNbCharged)/nbEvents;
>> 146 G4double meanEkin = 0.;
>> 147 if (fNbCharged) meanEkin = fEnergyCharged/fNbCharged;
>> 148
>> 149 G4cout
>> 150 << "\n d-rays : eLoss/primary= "
>> 151 << G4BestUnit(energyPerEvent, "Energy")
>> 152 << "\t nb of d-rays= " << nbPerEvent
>> 153 << " <Tkin>= " << G4BestUnit(meanEkin, "Energy")
>> 154 << " Tmin= " << G4BestUnit(fEmin[0], "Energy")
>> 155 << " Tmax= " << G4BestUnit(fEmax[0], "Energy")
>> 156 << G4endl;
>> 157
>> 158 //neutral secondaries (bremsstrahlung, pixe)
>> 159 //
>> 160 energyPerEvent = fEnergyNeutral/nbEvents;
>> 161 nbPerEvent = double(fNbNeutral)/nbEvents;
>> 162 meanEkin = 0.;
>> 163 if (fNbNeutral) meanEkin = fEnergyNeutral/fNbNeutral;
>> 164
>> 165 G4cout
>> 166 << "\n gamma : eLoss/primary= "
>> 167 << G4BestUnit(energyPerEvent, "Energy")
>> 168 << "\t nb of gammas= " << nbPerEvent
>> 169 << " <Tkin>= " << G4BestUnit(meanEkin, "Energy")
>> 170 << " Tmin= " << G4BestUnit(fEmin[1], "Energy")
>> 171 << " Tmax= " << G4BestUnit(fEmax[1], "Energy")
>> 172 << G4endl;
>> 173
234 174
235 G4cout.precision(4); << 175 G4EmCalculator emCal;
236 << 176
237 // frequency of processes << 177 //local energy deposit
238 // <<
239 G4cout << "\n Process defining step :" << G4 <<
240 G4int index = 0; <<
241 for (const auto& procCounter : fProcCounter) <<
242 G4String procName = procCounter.first; <<
243 G4int count = procCounter.second; <<
244 G4String space = " "; <<
245 if (++index % 4 == 0) space = "\n"; <<
246 G4cout << " " << std::setw(15) << procName <<
247 } <<
248 G4cout << G4endl; <<
249 <<
250 // track length <<
251 // <<
252 G4double trackLPerEvent = fTrackLength / nbE <<
253 G4double nbStepPerEvent = double(fNbSteps) / <<
254 G4double stepSize = fTrackLength / fNbSteps; <<
255 <<
256 G4cout << "\n TrackLength = " << G4BestUnit( <<
257 << " nb of steps = " << nbStepPerEve <<
258 << " stepSize = " << G4BestUnit(step <<
259 << G4BestUnit(fStepMin, "Length") << <<
260 << G4endl; <<
261 <<
262 // continuous energy deposited by primary tr <<
263 // 178 //
264 G4double energyPerEvent = fEdepPrimary / nbE << 179 energyPerEvent = fEnergyDeposit/nbEvents;
265 <<
266 G4cout << "\n Energy continuously deposited <<
267 << " (restricted dE/dx) dE1 = " << G <<
268 << G4BestUnit(fEdepPrimMin, "Energy") <<
269 << ")" << G4endl; <<
270 <<
271 // eveluation of dE1 from reading restricted <<
272 // 180 //
273 G4EmCalculator emCal; << 181 G4double r0 = emCal.GetRangeFromRestricteDEDX(ePrimary,particle,material);
274 <<
275 G4double r0 = emCal.GetRangeFromRestricteDED <<
276 G4double r1 = r0 - trackLPerEvent; 182 G4double r1 = r0 - trackLPerEvent;
277 G4double etry = ePrimary - energyPerEvent; << 183 G4double etry = ePrimary - energyPerEvent;
278 G4double efinal = 0.; 184 G4double efinal = 0.;
279 if (r1 > 0.) efinal = GetEnergyFromRestricte << 185 if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry);
280 G4double dEtable = ePrimary - efinal; 186 G4double dEtable = ePrimary - efinal;
281 G4double ratio = 0.; 187 G4double ratio = 0.;
282 if (dEtable > 0.) ratio = energyPerEvent / d << 188 if (dEtable > 0.) ratio = energyPerEvent/dEtable;
283 << 189
284 G4cout << "\n Evaluation of dE1 from reading << 190 G4cout
285 << G4BestUnit(dEtable, "Energy") << " << 191 << "\n deposit : eLoss/primary= "
286 << 192 << G4BestUnit(energyPerEvent, "Energy")
287 // energy transfered to secondary particles << 193 << "\t <dEcut > table= "
>> 194 << G4BestUnit(dEtable, "Energy")
>> 195 << " ---> simul/reference= " << ratio
>> 196 << G4endl;
>> 197
>> 198 //total energy transferred
288 // 199 //
289 G4cout << "\n Energy transfered to secondary << 200 G4double energyTotal = fEnergyDeposit + fEnergyCharged + fEnergyNeutral;
290 std::map<G4String, MinMaxData>::iterator it1 << 201 energyPerEvent = energyTotal/nbEvents;
291 for (it1 = fEtransfByProcess.begin(); it1 != <<
292 G4String name = it1->first; <<
293 MinMaxData data = it1->second; <<
294 energyPerEvent = data.fVsum / nbEvents; <<
295 G4double eMin = data.fVmin; <<
296 G4double eMax = data.fVmax; <<
297 <<
298 G4cout << " " << std::setw(17) << "due to <<
299 << G4BestUnit(energyPerEvent, "Ener <<
300 << G4BestUnit(eMax, "Energy") << ") <<
301 } <<
302 <<
303 // total energy tranfered : dE3 = sum of dE2 <<
304 // 202 //
305 energyPerEvent = fEnergyTransfered / nbEvent << 203 r0 = emCal.GetCSDARange(ePrimary,particle,material);
306 <<
307 G4cout << "\n Total energy transfered to sec <<
308 << G4BestUnit(energyPerEvent, "Energy <<
309 << " --> " << G4BestUnit(fEtransfMax, <<
310 <<
311 // total energy lost by incident particle : <<
312 // <<
313 energyPerEvent = fEnergyLost / nbEvents; <<
314 <<
315 G4cout << "\n Total energy lost by incident <<
316 << G4BestUnit(energyPerEvent, "Energy <<
317 << " --> " << G4BestUnit(fElostMax, " <<
318 <<
319 // calcul of energy lost from energy balance <<
320 // <<
321 energyPerEvent = fEnergyBalance / nbEvents; <<
322 <<
323 G4cout << "\n calcul of dE4 from energy bala <<
324 << G4BestUnit(energyPerEvent, "Energy <<
325 << " --> " << G4BestUnit(fEbalMax, "E <<
326 <<
327 // eveluation of dE4 from reading full Range <<
328 // <<
329 r0 = emCal.GetCSDARange(ePrimary, particle, <<
330 r1 = r0 - trackLPerEvent; 204 r1 = r0 - trackLPerEvent;
331 etry = ePrimary - energyPerEvent; 205 etry = ePrimary - energyPerEvent;
332 efinal = 0.; 206 efinal = 0.;
333 if (r1 > 0.) efinal = GetEnergyFromCSDARange << 207 if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry);
334 dEtable = ePrimary - efinal; 208 dEtable = ePrimary - efinal;
335 ratio = 0.; 209 ratio = 0.;
336 if (dEtable > 0.) ratio = energyPerEvent / d << 210 if (dEtable > 0.) ratio = energyPerEvent/dEtable;
337 << 211
338 G4cout << "\n Evaluation of dE4 from reading << 212 G4cout
339 << G4BestUnit(dEtable, "Energy") << " << 213 << "\n total : eLoss/primary= "
340 << 214 << G4BestUnit(energyPerEvent, "Energy")
341 // energy spectrum of secondary particles << 215 << "\t <dEfull> table= "
342 // << 216 << G4BestUnit(dEtable, "Energy")
343 G4cout << "\n Energy spectrum of secondary p << 217 << " ---> simul/reference= " << ratio
344 std::map<G4String, MinMaxData>::iterator it2 << 218 << G4endl;
345 for (it2 = fEkinOfSecondaries.begin(); it2 ! <<
346 G4String name = it2->first; <<
347 MinMaxData data = it2->second; <<
348 G4int count = data.fCount; <<
349 G4double eMean = data.fVsum / count; <<
350 G4double eMin = data.fVmin; <<
351 G4double eMax = data.fVmax; <<
352 <<
353 G4cout << " " << std::setw(13) << name << <<
354 << " Emean = " << std::setw(6) << <<
355 << G4BestUnit(eMin, "Energy") << " <<
356 } <<
357 G4cout << G4endl; <<
358 <<
359 // continuous energy deposited by secondary <<
360 // (only if secondary particles are tracked <<
361 // <<
362 if (fEdepSecondary > 0.) { <<
363 energyPerEvent = fEdepSecondary / nbEvents <<
364 <<
365 G4cout << "\n Energy continuously deposite <<
366 << " (restricted dE/dx) dE5 = " << <<
367 << G4BestUnit(fEdepSecMin, "Energy" <<
368 << ")" << G4endl; <<
369 <<
370 // total energy deposited : dE6 = dE1 + dE <<
371 // <<
372 energyPerEvent = fEdepTotal / nbEvents; <<
373 <<
374 G4cout << "\n Total energy deposited : dE6 <<
375 << G4BestUnit(energyPerEvent, "Ener <<
376 << " --> " << G4BestUnit(fEdepTotMa <<
377 << G4endl; <<
378 } <<
379 219
380 G4cout.precision(prec); 220 G4cout.precision(prec);
381 221
382 // clear maps <<
383 // <<
384 fProcCounter.clear(); <<
385 fEtransfByProcess.clear(); <<
386 fEkinOfSecondaries.clear(); <<
387 <<
388 // save histograms <<
389 G4AnalysisManager* analysisManager = G4Analy <<
390 if (analysisManager->IsActive()) { <<
391 analysisManager->Write(); <<
392 analysisManager->CloseFile(); <<
393 } <<
394 <<
395 // show Rndm status 222 // show Rndm status
396 CLHEP::HepRandom::showEngineStatus(); 223 CLHEP::HepRandom::showEngineStatus();
397 } 224 }
398 225
399 //....oooOO0OOooo........oooOO0OOooo........oo 226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400 227
401 G4double RunAction::GetEnergyFromRestrictedRan << 228 G4double RunAction::GetEnergyFromRestrictedRange(G4double range,
402 << 229 G4ParticleDefinition* particle, G4Material* material, G4double Etry)
403 { 230 {
404 G4EmCalculator emCal; 231 G4EmCalculator emCal;
405 << 232
406 G4double Energy = Etry, dE = 0., dEdx; 233 G4double Energy = Etry, dE = 0., dEdx;
407 G4double r, dr; 234 G4double r, dr;
408 G4double err = 1., errmax = 0.00001; << 235 G4double err = 1., errmax = 0.00001;
409 G4int iter = 0, itermax = 10; << 236 G4int iter = 0 , itermax = 10;
410 while (err > errmax && iter < itermax) { 237 while (err > errmax && iter < itermax) {
411 iter++; 238 iter++;
412 Energy -= dE; 239 Energy -= dE;
413 r = emCal.GetRangeFromRestricteDEDX(Energy << 240 r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material);
414 dr = r - range; << 241 dr = r - range;
415 dEdx = emCal.GetDEDX(Energy, particle, mat << 242 dEdx = emCal.GetDEDX(Energy,particle,material);
416 dE = dEdx * dr; << 243 dE = dEdx*dr;
417 err = std::abs(dE) / Energy; << 244 err = std::abs(dE)/Energy;
418 } 245 }
419 if (iter == itermax) { 246 if (iter == itermax) {
420 G4cout << "\n ---> warning: RunAction::Ge << 247 G4cout
421 << " Etry = " << G4BestUnit(Etry, << 248 << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
422 << " Energy = " << G4BestUnit(Ene << 249 << " Etry = " << G4BestUnit(Etry,"Energy")
423 << " iter = " << iter << G4endl; << 250 << " Energy = " << G4BestUnit(Energy,"Energy")
424 } << 251 << " err = " << err
425 << 252 << " iter = " << iter << G4endl;
426 return Energy; << 253 }
>> 254
>> 255 return Energy;
427 } 256 }
428 257
429 //....oooOO0OOooo........oooOO0OOooo........oo 258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430 259
431 G4double RunAction::GetEnergyFromCSDARange(G4d << 260 G4double RunAction::GetEnergyFromCSDARange(G4double range,
432 G4M << 261 G4ParticleDefinition* particle, G4Material* material, G4double Etry)
433 { 262 {
434 G4EmCalculator emCal; 263 G4EmCalculator emCal;
435 << 264
436 G4double Energy = Etry, dE = 0., dEdx; 265 G4double Energy = Etry, dE = 0., dEdx;
437 G4double r, dr; 266 G4double r, dr;
438 G4double err = 1., errmax = 0.00001; << 267 G4double err = 1., errmax = 0.00001;
439 G4int iter = 0, itermax = 10; << 268 G4int iter = 0 , itermax = 10;
440 while (err > errmax && iter < itermax) { 269 while (err > errmax && iter < itermax) {
441 iter++; 270 iter++;
442 Energy -= dE; 271 Energy -= dE;
443 r = emCal.GetCSDARange(Energy, particle, m << 272 r = emCal.GetCSDARange(Energy,particle,material);
444 dr = r - range; << 273 dr = r - range;
445 dEdx = emCal.ComputeTotalDEDX(Energy, part << 274 dEdx = emCal.ComputeTotalDEDX(Energy,particle,material);
446 dE = dEdx * dr; << 275 dE = dEdx*dr;
447 err = std::abs(dE) / Energy; << 276 err = std::abs(dE)/Energy;
448 } 277 }
449 if (iter == itermax) { 278 if (iter == itermax) {
450 G4cout << "\n ---> warning: RunAction::Ge << 279 G4cout
451 << " Etry = " << G4BestUnit(Etry, << 280 << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
452 << " Energy = " << G4BestUnit(Ene << 281 << " Etry = " << G4BestUnit(Etry,"Energy")
453 << " iter = " << iter << G4endl; << 282 << " Energy = " << G4BestUnit(Energy,"Energy")
454 } << 283 << " err = " << err
455 << 284 << " iter = " << iter << G4endl;
456 return Energy; << 285 }
>> 286
>> 287 return Energy;
457 } 288 }
458 289
459 //....oooOO0OOooo........oooOO0OOooo........oo 290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460 291