| 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/TestEm12/src/Run.cc 26 /// \file electromagnetic/TestEm12/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 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 33
33 #include "Run.hh" 34 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 35 #include "DetectorConstruction.hh"
>> 36
36 #include "HistoManager.hh" 37 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh" 38 #include "PrimaryGeneratorAction.hh"
38 39
39 #include "G4Material.hh" 40 #include "G4Material.hh"
40 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
41 #include "G4UnitsTable.hh" 42 #include "G4UnitsTable.hh"
42 43
43 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 45
45 Run::Run(DetectorConstruction* detector) : fDe << 46 Run::Run(DetectorConstruction* detector)
>> 47 : G4Run(),
>> 48 fDetector(detector),
>> 49 fParticle(0), fEkin(0.),
>> 50 fEdeposit(0.), fEdeposit2(0.),
>> 51 fTrackLen(0.), fTrackLen2(0.),
>> 52 fProjRange(0.), fProjRange2(0.),
>> 53 fNbOfSteps(0), fNbOfSteps2(0),
>> 54 fStepSize(0.), fStepSize2(0.),
>> 55 fCsdaRange(0.)
>> 56 { }
46 57
47 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48 59
49 void Run::SetPrimary(G4ParticleDefinition* par << 60 Run::~Run()
50 { << 61 { }
>> 62
>> 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 64
>> 65 void Run::SetPrimary (G4ParticleDefinition* particle, G4double energy)
>> 66 {
51 fParticle = particle; 67 fParticle = particle;
52 fEkin = energy; << 68 fEkin = energy;
53 } 69 }
54 70
55 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 72
57 void Run::AddEdep(G4double e) << 73 void Run::AddEdep (G4double e)
58 { 74 {
59 fEdeposit += e; << 75 fEdeposit += e;
60 fEdeposit2 += e * e; << 76 fEdeposit2 += e*e;
61 } 77 }
62 78
63 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 << 80
65 void Run::AddTrackLength(G4double t) << 81 void Run::AddTrackLength (G4double t)
66 { 82 {
67 fTrackLen += t; << 83 fTrackLen += t;
68 fTrackLen2 += t * t; << 84 fTrackLen2 += t*t;
69 } 85 }
70 86
71 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 << 88
73 void Run::AddProjRange(G4double x) << 89 void Run::AddProjRange (G4double x)
74 { 90 {
75 fProjRange += x; << 91 fProjRange += x;
76 fProjRange2 += x * x; << 92 fProjRange2 += x*x;
77 } 93 }
78 94
79 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 << 96
81 void Run::AddStepSize(G4int nb, G4double st) << 97 void Run::AddStepSize (G4int nb, G4double st)
82 { 98 {
83 fNbOfSteps += nb; << 99 fNbOfSteps += nb;
84 fNbOfSteps2 += nb * nb; << 100 fNbOfSteps2 += nb*nb;
85 fStepSize += st; << 101 fStepSize += st ;
86 fStepSize2 += st * st; << 102 fStepSize2 += st*st;
87 } 103 }
88 104
89 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 << 106
91 void Run::SetCsdaRange(G4double value) << 107 void Run::SetCsdaRange(G4double value)
92 { 108 {
93 fCsdaRange = value; 109 fCsdaRange = value;
94 } 110 }
95 111
96 //....oooOO0OOooo........oooOO0OOooo........oo 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 << 113
98 G4double Run::GetCsdaRange() << 114 G4double Run::GetCsdaRange()
99 { 115 {
100 return fCsdaRange; 116 return fCsdaRange;
101 } 117 }
102 << 118
103 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 120
105 void Run::Merge(const G4Run* run) 121 void Run::Merge(const G4Run* run)
106 { 122 {
107 const Run* localRun = static_cast<const Run* 123 const Run* localRun = static_cast<const Run*>(run);
108 << 124
109 // pass information about primary particle 125 // pass information about primary particle
110 fParticle = localRun->fParticle; 126 fParticle = localRun->fParticle;
111 fEkin = localRun->fEkin; << 127 fEkin = localRun->fEkin;
112 128
113 // accumulate sums 129 // accumulate sums
114 fEdeposit += localRun->fEdeposit; << 130 fEdeposit += localRun->fEdeposit;
115 fEdeposit2 += localRun->fEdeposit2; << 131 fEdeposit2 += localRun->fEdeposit2;
116 fTrackLen += localRun->fTrackLen; << 132 fTrackLen += localRun->fTrackLen;
117 fTrackLen2 += localRun->fTrackLen2; << 133 fTrackLen2 += localRun->fTrackLen2;
118 fProjRange += localRun->fProjRange; << 134 fProjRange += localRun->fProjRange;
119 fProjRange2 += localRun->fProjRange2; 135 fProjRange2 += localRun->fProjRange2;
120 fNbOfSteps += localRun->fNbOfSteps; << 136 fNbOfSteps += localRun->fNbOfSteps ;
121 fNbOfSteps2 += localRun->fNbOfSteps2; 137 fNbOfSteps2 += localRun->fNbOfSteps2;
122 fStepSize += localRun->fStepSize; << 138 fStepSize += localRun->fStepSize;
123 fStepSize2 += localRun->fStepSize2; << 139 fStepSize2 += localRun->fStepSize2;
124 << 140
125 fCsdaRange = localRun->fCsdaRange; << 141 fCsdaRange = localRun->fCsdaRange;
126 142
127 G4Run::Merge(run); << 143 G4Run::Merge(run);
128 } << 144 }
129 145
130 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 147
132 void Run::EndOfRun() << 148 void Run::EndOfRun()
133 { 149 {
134 std::ios::fmtflags mode = G4cout.flags(); 150 std::ios::fmtflags mode = G4cout.flags();
135 G4cout.setf(std::ios::fixed, std::ios::float << 151 G4cout.setf(std::ios::fixed,std::ios::floatfield);
136 G4int prec = G4cout.precision(2); 152 G4int prec = G4cout.precision(2);
137 << 153
138 // run conditions << 154 //run conditions
139 // 155 //
140 G4Material* material = fDetector->GetAbsorMa 156 G4Material* material = fDetector->GetAbsorMaterial();
141 G4double density = material->GetDensity(); << 157 G4double density = material->GetDensity();
142 G4String partName = fParticle->GetParticleNa 158 G4String partName = fParticle->GetParticleName();
143 << 159
144 G4cout << "\n ======================== run s << 160 G4cout << "\n ======================== run summary =====================\n";
145 G4cout << "\n The run is " << numberOfEvent << 161 G4cout
146 << G4BestUnit(fEkin, "Energy") << " t << 162 << "\n The run is " << numberOfEvent << " "<< partName << " of "
147 << G4BestUnit(fDetector->GetAbsorRadi << 163 << G4BestUnit(fEkin,"Energy") << " through "
148 << " (density: " << G4BestUnit(densit << 164 << G4BestUnit(fDetector->GetAbsorRadius(),"Length") << " of "
>> 165 << material->GetName() << " (density: "
>> 166 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
149 167
150 if (numberOfEvent == 0) { 168 if (numberOfEvent == 0) {
151 G4cout.setf(mode, std::ios::floatfield); << 169 G4cout.setf(mode,std::ios::floatfield);
152 G4cout.precision(prec); << 170 G4cout.precision(prec);
153 return; 171 return;
154 } 172 }
155 << 173
156 fEdeposit /= numberOfEvent; << 174 fEdeposit /= numberOfEvent; fEdeposit2 /= numberOfEvent;
157 fEdeposit2 /= numberOfEvent; << 175 G4double rms = fEdeposit2 - fEdeposit*fEdeposit;
158 G4double rms = fEdeposit2 - fEdeposit * fEde << 176 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
159 if (rms > 0.) << 177
160 rms = std::sqrt(rms); << 178 G4cout.precision(3);
161 else << 179 G4cout
162 rms = 0.; << 180 << "\n Total Energy deposited = " << G4BestUnit(fEdeposit,"Energy")
163 << 181 << " +- " << G4BestUnit( rms,"Energy")
164 G4cout.precision(3); << 182 << G4endl;
165 G4cout << "\n Total Energy deposited << 183
166 << G4BestUnit(rms, "Energy") << G4end << 184 //compute track length of primary track
167 << 185 //
168 // compute track length of primary track << 186 fTrackLen /= numberOfEvent; fTrackLen2 /= numberOfEvent;
169 // << 187 rms = fTrackLen2 - fTrackLen*fTrackLen;
170 fTrackLen /= numberOfEvent; << 188 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
171 fTrackLen2 /= numberOfEvent; << 189
172 rms = fTrackLen2 - fTrackLen * fTrackLen; << 190 G4cout.precision(3);
173 if (rms > 0.) << 191 G4cout
174 rms = std::sqrt(rms); << 192 << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length")
175 else << 193 << " +- " << G4BestUnit( rms,"Length");
176 rms = 0.; << 194
177 << 195 //compare with csda range
178 G4cout.precision(3); << 196 //
179 G4cout << "\n Track length of primary track << 197 G4cout
180 << G4BestUnit(rms, "Length"); << 198 << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange,"Length")
181 << 199 << " (from full dE/dx)" << G4endl;
182 // compare with csda range << 200
183 // << 201
184 G4cout << "\n Range from EmCalculator = " << << 202 //compute projected range of primary track
185 << " (from full dE/dx)" << G4endl; << 203 //
186 << 204 fProjRange /= numberOfEvent; fProjRange2 /= numberOfEvent;
187 // compute projected range of primary track << 205 rms = fProjRange2 - fProjRange*fProjRange;
188 // << 206 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
189 fProjRange /= numberOfEvent; << 207
190 fProjRange2 /= numberOfEvent; << 208 G4cout
191 rms = fProjRange2 - fProjRange * fProjRange; << 209 << "\n Projected range = " << G4BestUnit(fProjRange,"Length")
192 if (rms > 0.) << 210 << " +- " << G4BestUnit( rms,"Length")
193 rms = std::sqrt(rms); << 211 << G4endl;
194 else << 212
195 rms = 0.; << 213 //nb of steps and step size of primary track
196 <<
197 G4cout << "\n Projected range <<
198 << G4BestUnit(rms, "Length") << G4end <<
199 <<
200 // nb of steps and step size of primary trac <<
201 // 214 //
202 G4double dNofEvents = double(numberOfEvent); 215 G4double dNofEvents = double(numberOfEvent);
203 G4double fNbSteps = fNbOfSteps / dNofEvents, << 216 G4double fNbSteps = fNbOfSteps/dNofEvents,
204 rms = fNbSteps2 - fNbSteps * fNbSteps; << 217 fNbSteps2 = fNbOfSteps2/dNofEvents;
205 if (rms > 0.) << 218 rms = fNbSteps2 - fNbSteps*fNbSteps;
206 rms = std::sqrt(rms); << 219 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
207 else <<
208 rms = 0.; <<
209 220
210 G4cout.precision(2); << 221 G4cout.precision(2);
211 G4cout << "\n Nb of steps of primary track 222 G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms;
212 << 223
213 fStepSize /= numberOfEvent; << 224 fStepSize /= numberOfEvent; fStepSize2 /= numberOfEvent;
214 fStepSize2 /= numberOfEvent; << 225 rms = fStepSize2 - fStepSize*fStepSize;
215 rms = fStepSize2 - fStepSize * fStepSize; << 226 if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
216 if (rms > 0.) << 227
217 rms = std::sqrt(rms); << 228 G4cout.precision(3);
218 else << 229 G4cout
219 rms = 0.; << 230 << "\t Step size= " << G4BestUnit(fStepSize,"Length")
220 << 231 << " +- " << G4BestUnit( rms,"Length")
221 G4cout.precision(3); << 232 << G4endl;
222 G4cout << "\t Step size= " << G4BestUnit(fSt <<
223 << G4BestUnit(rms, "Length") << G4end <<
224 233
225 // normalize histograms of longitudinal ener 234 // normalize histograms of longitudinal energy profile
226 // 235 //
227 G4AnalysisManager* analysisManager = G4Analy 236 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
228 G4int ih = 1; 237 G4int ih = 1;
229 G4double binWidth = analysisManager->GetH1Wi << 238 G4double binWidth = analysisManager->GetH1Width(ih)
230 G4double fac = (1. / (numberOfEvent * binWid << 239 *analysisManager->GetH1Unit(ih);
231 analysisManager->ScaleH1(ih, fac); << 240 G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
232 << 241 analysisManager->ScaleH1(ih,fac);
>> 242
233 // normalize histogram d(E/E0)/d(r/r0) 243 // normalize histogram d(E/E0)/d(r/r0)
234 // << 244 //
235 ih = 8; 245 ih = 8;
236 binWidth = analysisManager->GetH1Width(ih); 246 binWidth = analysisManager->GetH1Width(ih);
237 fac = 1. / (numberOfEvent * binWidth * fEkin << 247 fac = 1./(numberOfEvent*binWidth*fEkin);
238 analysisManager->ScaleH1(ih, fac); << 248 analysisManager->ScaleH1(ih,fac);
239 << 249
240 // reset default formats 250 // reset default formats
241 G4cout.setf(mode, std::ios::floatfield); << 251 G4cout.setf(mode,std::ios::floatfield);
242 G4cout.precision(prec); 252 G4cout.precision(prec);
243 } 253 }
244 254
245 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246 256