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