| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file electromagnetic/TestEm12/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4Material.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4UnitsTable.hh"
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 Run::Run(DetectorConstruction* detector) : fDe
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 void Run::SetPrimary(G4ParticleDefinition* par
50 {
51 fParticle = particle;
52 fEkin = energy;
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 void Run::AddEdep(G4double e)
58 {
59 fEdeposit += e;
60 fEdeposit2 += e * e;
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 void Run::AddTrackLength(G4double t)
66 {
67 fTrackLen += t;
68 fTrackLen2 += t * t;
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 void Run::AddProjRange(G4double x)
74 {
75 fProjRange += x;
76 fProjRange2 += x * x;
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oo
80
81 void Run::AddStepSize(G4int nb, G4double st)
82 {
83 fNbOfSteps += nb;
84 fNbOfSteps2 += nb * nb;
85 fStepSize += st;
86 fStepSize2 += st * st;
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 void Run::SetCsdaRange(G4double value)
92 {
93 fCsdaRange = value;
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97
98 G4double Run::GetCsdaRange()
99 {
100 return fCsdaRange;
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oo
104
105 void Run::Merge(const G4Run* run)
106 {
107 const Run* localRun = static_cast<const Run*
108
109 // pass information about primary particle
110 fParticle = localRun->fParticle;
111 fEkin = localRun->fEkin;
112
113 // accumulate sums
114 fEdeposit += localRun->fEdeposit;
115 fEdeposit2 += localRun->fEdeposit2;
116 fTrackLen += localRun->fTrackLen;
117 fTrackLen2 += localRun->fTrackLen2;
118 fProjRange += localRun->fProjRange;
119 fProjRange2 += localRun->fProjRange2;
120 fNbOfSteps += localRun->fNbOfSteps;
121 fNbOfSteps2 += localRun->fNbOfSteps2;
122 fStepSize += localRun->fStepSize;
123 fStepSize2 += localRun->fStepSize2;
124
125 fCsdaRange = localRun->fCsdaRange;
126
127 G4Run::Merge(run);
128 }
129
130 //....oooOO0OOooo........oooOO0OOooo........oo
131
132 void Run::EndOfRun()
133 {
134 std::ios::fmtflags mode = G4cout.flags();
135 G4cout.setf(std::ios::fixed, std::ios::float
136 G4int prec = G4cout.precision(2);
137
138 // run conditions
139 //
140 G4Material* material = fDetector->GetAbsorMa
141 G4double density = material->GetDensity();
142 G4String partName = fParticle->GetParticleNa
143
144 G4cout << "\n ======================== run s
145 G4cout << "\n The run is " << numberOfEvent
146 << G4BestUnit(fEkin, "Energy") << " t
147 << G4BestUnit(fDetector->GetAbsorRadi
148 << " (density: " << G4BestUnit(densit
149
150 if (numberOfEvent == 0) {
151 G4cout.setf(mode, std::ios::floatfield);
152 G4cout.precision(prec);
153 return;
154 }
155
156 fEdeposit /= numberOfEvent;
157 fEdeposit2 /= numberOfEvent;
158 G4double rms = fEdeposit2 - fEdeposit * fEde
159 if (rms > 0.)
160 rms = std::sqrt(rms);
161 else
162 rms = 0.;
163
164 G4cout.precision(3);
165 G4cout << "\n Total Energy deposited
166 << G4BestUnit(rms, "Energy") << G4end
167
168 // compute track length of primary track
169 //
170 fTrackLen /= numberOfEvent;
171 fTrackLen2 /= numberOfEvent;
172 rms = fTrackLen2 - fTrackLen * fTrackLen;
173 if (rms > 0.)
174 rms = std::sqrt(rms);
175 else
176 rms = 0.;
177
178 G4cout.precision(3);
179 G4cout << "\n Track length of primary track
180 << G4BestUnit(rms, "Length");
181
182 // compare with csda range
183 //
184 G4cout << "\n Range from EmCalculator = " <<
185 << " (from full dE/dx)" << G4endl;
186
187 // compute projected range of primary track
188 //
189 fProjRange /= numberOfEvent;
190 fProjRange2 /= numberOfEvent;
191 rms = fProjRange2 - fProjRange * fProjRange;
192 if (rms > 0.)
193 rms = std::sqrt(rms);
194 else
195 rms = 0.;
196
197 G4cout << "\n Projected range
198 << G4BestUnit(rms, "Length") << G4end
199
200 // nb of steps and step size of primary trac
201 //
202 G4double dNofEvents = double(numberOfEvent);
203 G4double fNbSteps = fNbOfSteps / dNofEvents,
204 rms = fNbSteps2 - fNbSteps * fNbSteps;
205 if (rms > 0.)
206 rms = std::sqrt(rms);
207 else
208 rms = 0.;
209
210 G4cout.precision(2);
211 G4cout << "\n Nb of steps of primary track
212
213 fStepSize /= numberOfEvent;
214 fStepSize2 /= numberOfEvent;
215 rms = fStepSize2 - fStepSize * fStepSize;
216 if (rms > 0.)
217 rms = std::sqrt(rms);
218 else
219 rms = 0.;
220
221 G4cout.precision(3);
222 G4cout << "\t Step size= " << G4BestUnit(fSt
223 << G4BestUnit(rms, "Length") << G4end
224
225 // normalize histograms of longitudinal ener
226 //
227 G4AnalysisManager* analysisManager = G4Analy
228 G4int ih = 1;
229 G4double binWidth = analysisManager->GetH1Wi
230 G4double fac = (1. / (numberOfEvent * binWid
231 analysisManager->ScaleH1(ih, fac);
232
233 // normalize histogram d(E/E0)/d(r/r0)
234 //
235 ih = 8;
236 binWidth = analysisManager->GetH1Width(ih);
237 fac = 1. / (numberOfEvent * binWidth * fEkin
238 analysisManager->ScaleH1(ih, fac);
239
240 // reset default formats
241 G4cout.setf(mode, std::ios::floatfield);
242 G4cout.precision(prec);
243 }
244
245 //....oooOO0OOooo........oooOO0OOooo........oo
246