| 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 electromagnetic/TestEm12/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
44
45 Run::Run(DetectorConstruction* detector) : fDetector(detector) {}
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48
49 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
50 {
51 fParticle = particle;
52 fEkin = energy;
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 void Run::AddEdep(G4double e)
58 {
59 fEdeposit += e;
60 fEdeposit2 += e * e;
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64
65 void Run::AddTrackLength(G4double t)
66 {
67 fTrackLen += t;
68 fTrackLen2 += t * t;
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72
73 void Run::AddProjRange(G4double x)
74 {
75 fProjRange += x;
76 fProjRange2 += x * x;
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
90
91 void Run::SetCsdaRange(G4double value)
92 {
93 fCsdaRange = value;
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97
98 G4double Run::GetCsdaRange()
99 {
100 return fCsdaRange;
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105 void Run::Merge(const G4Run* run)
106 {
107 const Run* localRun = static_cast<const Run*>(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........oooOO0OOooo........oooOO0OOooo......
131
132 void Run::EndOfRun()
133 {
134 std::ios::fmtflags mode = G4cout.flags();
135 G4cout.setf(std::ios::fixed, std::ios::floatfield);
136 G4int prec = G4cout.precision(2);
137
138 // run conditions
139 //
140 G4Material* material = fDetector->GetAbsorMaterial();
141 G4double density = material->GetDensity();
142 G4String partName = fParticle->GetParticleName();
143
144 G4cout << "\n ======================== run summary =====================\n";
145 G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
146 << G4BestUnit(fEkin, "Energy") << " through "
147 << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << " of " << material->GetName()
148 << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
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 * fEdeposit;
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 = " << G4BestUnit(fEdeposit, "Energy") << " +- "
166 << G4BestUnit(rms, "Energy") << G4endl;
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 = " << G4BestUnit(fTrackLen, "Length") << " +- "
180 << G4BestUnit(rms, "Length");
181
182 // compare with csda range
183 //
184 G4cout << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange, "Length")
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 = " << G4BestUnit(fProjRange, "Length") << " +- "
198 << G4BestUnit(rms, "Length") << G4endl;
199
200 // nb of steps and step size of primary track
201 //
202 G4double dNofEvents = double(numberOfEvent);
203 G4double fNbSteps = fNbOfSteps / dNofEvents, fNbSteps2 = fNbOfSteps2 / 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 = " << fNbSteps << " +- " << rms;
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(fStepSize, "Length") << " +- "
223 << G4BestUnit(rms, "Length") << G4endl;
224
225 // normalize histograms of longitudinal energy profile
226 //
227 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
228 G4int ih = 1;
229 G4double binWidth = analysisManager->GetH1Width(ih) * analysisManager->GetH1Unit(ih);
230 G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
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........oooOO0OOooo........oooOO0OOooo......
246