| 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/TestEm3/src/Run.cc 26 /// \file electromagnetic/TestEm3/src/Run.cc
27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class
28 // 28 //
29 // << 29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 32
33 #include "Run.hh" 33 #include "Run.hh"
34 <<
35 #include "DetectorConstruction.hh" 34 #include "DetectorConstruction.hh"
36 #include "EmAcceptance.hh" <<
37 #include "HistoManager.hh" <<
38 #include "PrimaryGeneratorAction.hh" 35 #include "PrimaryGeneratorAction.hh"
>> 36 #include "HistoManager.hh"
>> 37 #include "EmAcceptance.hh"
39 38
40 #include "G4Electron.hh" <<
41 #include "G4Gamma.hh" <<
42 #include "G4ParticleDefinition.hh" <<
43 #include "G4ParticleTable.hh" 39 #include "G4ParticleTable.hh"
44 #include "G4Positron.hh" << 40 #include "G4ParticleDefinition.hh"
45 #include "G4SystemOfUnits.hh" <<
46 #include "G4Track.hh" 41 #include "G4Track.hh"
>> 42 #include "G4Gamma.hh"
>> 43 #include "G4Electron.hh"
>> 44 #include "G4Positron.hh"
>> 45
47 #include "G4UnitsTable.hh" 46 #include "G4UnitsTable.hh"
>> 47 #include "G4SystemOfUnits.hh"
48 48
49 #include <iomanip> 49 #include <iomanip>
50 50
51 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 52
53 Run::Run(DetectorConstruction* det) : fDetecto << 53 Run::Run(DetectorConstruction* det)
>> 54 : G4Run(),
>> 55 fDetector(det),
>> 56 fParticle(nullptr), fEkin(0.),
>> 57 fChargedStep(0), fNeutralStep(0),
>> 58 fN_gamma(0), fN_elec(0), fN_pos(0),
>> 59 fApplyLimit(false)
54 { 60 {
55 // initialize cumulative quantities << 61 //initialize cumulative quantities
56 // 62 //
57 for (G4int k = 0; k < kMaxAbsor; k++) { << 63 for (G4int k=0; k<kMaxAbsor; k++) {
58 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = << 64 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
59 fEnergyDeposit[k].clear(); 65 fEnergyDeposit[k].clear();
60 fEdeptrue[k] = fRmstrue[k] = 1.; 66 fEdeptrue[k] = fRmstrue[k] = 1.;
61 fLimittrue[k] = DBL_MAX; << 67 fLimittrue[k] = DBL_MAX;
62 } 68 }
63 << 69
64 fEdepTot = fEdepTot2 = 0.; << 70 //initialize Eflow
65 fEleakTot = fEleakTot2 = 0.; <<
66 fEtotal = fEtotal2 = 0.; <<
67 <<
68 // initialize Eflow <<
69 // 71 //
70 G4int nbPlanes = (fDetector->GetNbOfLayers() << 72 G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
71 fEnergyFlow.resize(nbPlanes); 73 fEnergyFlow.resize(nbPlanes);
72 fLateralEleak.resize(nbPlanes); 74 fLateralEleak.resize(nbPlanes);
73 for (G4int k = 0; k < nbPlanes; k++) { << 75 for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
74 fEnergyFlow[k] = fLateralEleak[k] = 0.; <<
75 } <<
76 } 76 }
77 77
78 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 79
>> 80 Run::~Run()
>> 81 { }
>> 82
>> 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 84
80 void Run::SetPrimary(G4ParticleDefinition* par 85 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
81 { << 86 {
82 fParticle = particle; 87 fParticle = particle;
83 fEkin = energy; 88 fEkin = energy;
84 } 89 }
85 90
86 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 92
88 void Run::FillPerEvent(G4int kAbs, G4double EA 93 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
89 { 94 {
90 // accumulate statistic with restriction << 95 //accumulate statistic with restriction
91 // 96 //
92 if (fApplyLimit) fEnergyDeposit[kAbs].push_b << 97 if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
93 fSumEAbs[kAbs] += EAbs; << 98 fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs;
94 fSum2EAbs[kAbs] += EAbs * EAbs; << 99 fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs;
95 fSumLAbs[kAbs] += LAbs; <<
96 fSum2LAbs[kAbs] += LAbs * LAbs; <<
97 } <<
98 <<
99 //....oooOO0OOooo........oooOO0OOooo........oo <<
100 <<
101 void Run::SumEnergies(G4double edeptot, G4doub <<
102 { <<
103 fEdepTot += edeptot; <<
104 fEdepTot2 += edeptot * edeptot; <<
105 fEleakTot += eleak; <<
106 fEleakTot2 += eleak * eleak; <<
107 <<
108 G4double etotal = edeptot + eleak; <<
109 fEtotal += etotal; <<
110 fEtotal2 += etotal * etotal; <<
111 } 100 }
112 101
113 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 103
115 void Run::SumEnergyFlow(G4int plane, G4double 104 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
116 { 105 {
117 fEnergyFlow[plane] += Eflow; 106 fEnergyFlow[plane] += Eflow;
118 } 107 }
119 108
120 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 << 110
122 void Run::SumLateralEleak(G4int cell, G4double 111 void Run::SumLateralEleak(G4int cell, G4double Eflow)
123 { 112 {
124 fLateralEleak[cell] += Eflow; 113 fLateralEleak[cell] += Eflow;
125 } 114 }
126 115
127 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 << 117
129 void Run::AddChargedStep() << 118 void Run::AddChargedStep()
130 { 119 {
131 fChargedStep += 1.0; << 120 fChargedStep += 1.0;
132 } 121 }
133 122
134 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 << 124
136 void Run::AddNeutralStep() << 125 void Run::AddNeutralStep()
137 { 126 {
138 fNeutralStep += 1.0; << 127 fNeutralStep += 1.0;
139 } 128 }
140 << 129
141 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 131
143 void Run::AddSecondaryTrack(const G4Track* tra 132 void Run::AddSecondaryTrack(const G4Track* track)
144 { 133 {
145 const G4ParticleDefinition* d = track->GetDe 134 const G4ParticleDefinition* d = track->GetDefinition();
146 if (d == G4Gamma::Gamma()) { << 135 if(d == G4Gamma::Gamma()) { ++fN_gamma; }
147 ++fN_gamma; << 136 else if (d == G4Electron::Electron()) { ++fN_elec; }
148 } << 137 else if (d == G4Positron::Positron()) { ++fN_pos; }
149 else if (d == G4Electron::Electron()) { <<
150 ++fN_elec; <<
151 } <<
152 else if (d == G4Positron::Positron()) { <<
153 ++fN_pos; <<
154 } <<
155 } 138 }
156 << 139
157 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
158 141
159 void Run::Merge(const G4Run* run) 142 void Run::Merge(const G4Run* run)
160 { 143 {
161 const Run* localRun = static_cast<const Run* 144 const Run* localRun = static_cast<const Run*>(run);
162 145
163 // pass information about primary particle 146 // pass information about primary particle
164 fParticle = localRun->fParticle; 147 fParticle = localRun->fParticle;
165 fEkin = localRun->fEkin; << 148 fEkin = localRun->fEkin;
166 149
167 // accumulate sums 150 // accumulate sums
168 // 151 //
169 for (G4int k = 0; k < kMaxAbsor; k++) { << 152 for (G4int k=0; k<kMaxAbsor; k++) {
170 fSumEAbs[k] += localRun->fSumEAbs[k]; << 153 fSumEAbs[k] += localRun->fSumEAbs[k];
171 fSum2EAbs[k] += localRun->fSum2EAbs[k]; << 154 fSum2EAbs[k] += localRun->fSum2EAbs[k];
172 fSumLAbs[k] += localRun->fSumLAbs[k]; << 155 fSumLAbs[k] += localRun->fSumLAbs[k];
173 fSum2LAbs[k] += localRun->fSum2LAbs[k]; 156 fSum2LAbs[k] += localRun->fSum2LAbs[k];
174 } 157 }
175 << 158
176 fEdepTot += localRun->fEdepTot; << 159 G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
177 fEdepTot2 += localRun->fEdepTot2; << 160 for (G4int k=0; k<nbPlanes; k++) {
178 << 161 fEnergyFlow[k] += localRun->fEnergyFlow[k];
179 fEleakTot += localRun->fEleakTot; <<
180 fEleakTot2 += localRun->fEleakTot2; <<
181 <<
182 fEtotal += localRun->fEtotal; <<
183 fEtotal2 += localRun->fEtotal2; <<
184 <<
185 G4int nbPlanes = (fDetector->GetNbOfLayers() <<
186 for (G4int k = 0; k < nbPlanes; k++) { <<
187 fEnergyFlow[k] += localRun->fEnergyFlow[k] <<
188 fLateralEleak[k] += localRun->fLateralElea 162 fLateralEleak[k] += localRun->fLateralEleak[k];
189 } 163 }
190 << 164
191 for (G4int k=0; k<kMaxAbsor; k++) { << 165
192 fEnergyDeposit[k].insert(fEnergyDeposit[k] << 166 fChargedStep += localRun->fChargedStep;
193 localRun->fEnergyDeposit[k].begin(), local <<
194 } <<
195 <<
196 fChargedStep += localRun->fChargedStep; <<
197 fNeutralStep += localRun->fNeutralStep; 167 fNeutralStep += localRun->fNeutralStep;
198 << 168
199 fN_gamma += localRun->fN_gamma; << 169 fN_gamma += localRun->fN_gamma;
200 fN_elec += localRun->fN_elec; << 170 fN_elec += localRun->fN_elec;
201 fN_pos += localRun->fN_pos; << 171 fN_pos += localRun->fN_pos;
202 << 172
203 fApplyLimit = localRun->fApplyLimit; 173 fApplyLimit = localRun->fApplyLimit;
204 << 174
205 for (G4int k = 0; k < kMaxAbsor; k++) { << 175 for (G4int k=0; k<kMaxAbsor; k++) {
206 fEdeptrue[k] = localRun->fEdeptrue[k]; << 176 fEdeptrue[k] = localRun->fEdeptrue[k];
207 fRmstrue[k] = localRun->fRmstrue[k]; << 177 fRmstrue[k] = localRun->fRmstrue[k];
208 fLimittrue[k] = localRun->fLimittrue[k]; << 178 fLimittrue[k] = localRun->fLimittrue[k];
209 } << 179 }
210 << 180
211 G4Run::Merge(run); << 181 G4Run::Merge(run);
212 } << 182 }
213 183
214 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215 185
216 void Run::EndOfRun() 186 void Run::EndOfRun()
217 { 187 {
218 G4int nEvt = numberOfEvent; 188 G4int nEvt = numberOfEvent;
219 G4double norm = G4double(nEvt); << 189 G4double norm = G4double(nEvt);
220 if (norm > 0) norm = 1. / norm; << 190 if(norm > 0) norm = 1./norm;
221 G4double qnorm = std::sqrt(norm); 191 G4double qnorm = std::sqrt(norm);
222 192
223 fChargedStep *= norm; 193 fChargedStep *= norm;
224 fNeutralStep *= norm; 194 fNeutralStep *= norm;
225 195
226 // compute and print statistic << 196 //compute and print statistic
227 // 197 //
228 G4double beamEnergy = fEkin; 198 G4double beamEnergy = fEkin;
229 G4double sqbeam = std::sqrt(beamEnergy / GeV << 199 G4double sqbeam = std::sqrt(beamEnergy/GeV);
230 200
231 G4double MeanEAbs, MeanEAbs2, rmsEAbs, resol << 201 G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
232 G4double MeanLAbs, MeanLAbs2, rmsLAbs; << 202 G4double MeanLAbs,MeanLAbs2,rmsLAbs;
233 203
234 std::ios::fmtflags mode = G4cout.flags(); 204 std::ios::fmtflags mode = G4cout.flags();
235 G4int prec = G4cout.precision(2); << 205 G4int prec = G4cout.precision(2);
236 G4cout << "\n------------------------------- 206 G4cout << "\n------------------------------------------------------------\n";
237 G4cout << std::setw(14) << "material" << std << 207 G4cout << std::setw(14) << "material"
238 << "sqrt(E0(GeV))*rmsE/Emean" << std: << 208 << std::setw(17) << "Edep RMS"
239 << 209 << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
240 for (G4int k = 1; k <= fDetector->GetNbOfAbs << 210 << std::setw(23) << "total tracklen \n \n";
241 MeanEAbs = fSumEAbs[k] * norm; << 211
242 MeanEAbs2 = fSum2EAbs[k] * norm; << 212 for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
243 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - M << 213 {
244 // G4cout << "k= " << k << " RMS= " << r << 214 MeanEAbs = fSumEAbs[k]*norm;
245 // << " fApplyLimit: " << fApplyLimi << 215 MeanEAbs2 = fSum2EAbs[k]*norm;
246 if (fApplyLimit) { << 216 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
247 G4int nn = 0; << 217 //G4cout << "k= " << k << " RMS= " << rmsEAbs
248 G4double sume = 0.0; << 218 // << " fApplyLimit: " << fApplyLimit << G4endl;
249 G4double sume2 = 0.0; << 219 if(fApplyLimit) {
250 // compute trancated means << 220 G4int nn = 0;
251 G4double lim = rmsEAbs * 2.5; << 221 G4double sume = 0.0;
252 for (G4int i = 0; i < nEvt; i++) { << 222 G4double sume2 = 0.0;
253 G4double e = (fEnergyDeposit[k])[i]; << 223 // compute trancated means
254 if (std::abs(e - MeanEAbs) < lim) { << 224 G4double lim = rmsEAbs * 2.5;
255 sume += e; << 225 for(G4int i=0; i<nEvt; i++) {
256 sume2 += e * e; << 226 G4double e = (fEnergyDeposit[k])[i];
257 nn++; << 227 if(std::abs(e - MeanEAbs) < lim) {
>> 228 sume += e;
>> 229 sume2 += e*e;
>> 230 nn++;
>> 231 }
258 } 232 }
>> 233 G4double norm1 = G4double(nn);
>> 234 if(norm1 > 0.0) norm1 = 1.0/norm1;
>> 235 MeanEAbs = sume*norm1;
>> 236 MeanEAbs2 = sume2*norm1;
>> 237 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
259 } 238 }
260 G4double norm1 = G4double(nn); <<
261 if (norm1 > 0.0) norm1 = 1.0 / norm1; <<
262 MeanEAbs = sume * norm1; <<
263 MeanEAbs2 = sume2 * norm1; <<
264 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - <<
265 } <<
266 <<
267 resolution = (MeanEAbs > 0.) ? 100. * sqbe <<
268 rmsres = resolution * qnorm; <<
269 <<
270 // Save mean and RMS <<
271 fSumEAbs[k] = MeanEAbs; <<
272 fSum2EAbs[k] = rmsEAbs; <<
273 <<
274 MeanLAbs = fSumLAbs[k] * norm; <<
275 MeanLAbs2 = fSum2LAbs[k] * norm; <<
276 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - M <<
277 <<
278 // print <<
279 // <<
280 G4cout << std::setw(14) << fDetector->GetA <<
281 << std::setprecision(5) << std::set <<
282 << std::setprecision(4) << std::set <<
283 << resolution << " +- " << std::set <<
284 << std::setw(10) << G4BestUnit(Mean <<
285 << G4BestUnit(rmsLAbs, "Length") << <<
286 } <<
287 239
288 // total energy deposited << 240 resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
289 // << 241 rmsres = resolution*qnorm;
290 fEdepTot *= norm; <<
291 fEdepTot2 *= norm; <<
292 G4double rmsEdep = std::sqrt(std::abs(fEdepT <<
293 <<
294 G4cout << "\n Total energy deposited = " << <<
295 << " +- " << G4BestUnit(rmsEdep, "Ene <<
296 <<
297 // Energy leakage <<
298 // <<
299 fEleakTot *= norm; <<
300 fEleakTot2 *= norm; <<
301 G4double rmsEleak = std::sqrt(std::abs(fElea <<
302 <<
303 G4cout << " Energy leakage = " << G4BestUnit <<
304 << G4BestUnit(rmsEleak, "Energy") << <<
305 242
306 // total energy << 243 // Save mean and RMS
307 // << 244 fSumEAbs[k] = MeanEAbs;
308 fEtotal *= norm; << 245 fSum2EAbs[k] = rmsEAbs;
309 fEtotal2 *= norm; << 246
310 G4double rmsEtotal = std::sqrt(std::abs(fEto << 247 MeanLAbs = fSumLAbs[k]*norm;
311 << 248 MeanLAbs2 = fSum2LAbs[k]*norm;
312 G4cout << " Total energy : Edep + Eleak = " << 249 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
313 << G4BestUnit(rmsEtotal, "Energy") << << 250
314 << 251 //print
315 G4cout << "--------------------------------- << 252 //
>> 253 G4cout
>> 254 << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
>> 255 << std::setprecision(5)
>> 256 << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " : "
>> 257 << std::setprecision(4)
>> 258 << std::setw(5) << G4BestUnit( rmsEAbs,"Energy")
>> 259 << std::setw(10) << resolution << " +- "
>> 260 << std::setw(5) << rmsres << " %"
>> 261 << std::setprecision(3)
>> 262 << std::setw(10) << G4BestUnit(MeanLAbs,"Length") << " +- "
>> 263 << std::setw(4) << G4BestUnit( rmsLAbs,"Length")
>> 264 << G4endl;
>> 265 }
>> 266 G4cout << "\n------------------------------------------------------------\n";
316 267
317 G4cout << " Beam particle " << fParticle->Ge << 268 G4cout << " Beam particle "
318 << " E = " << G4BestUnit(beamEnergy, << 269 << fParticle->GetParticleName()
319 G4cout << " Mean number of gamma " << << 270 << " E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
320 G4cout << " Mean number of e- " << << 271 G4cout << " Mean number of gamma " << (G4double)fN_gamma*norm << G4endl;
321 G4cout << " Mean number of e+ " << << 272 G4cout << " Mean number of e- " << (G4double)fN_elec*norm << G4endl;
322 G4cout << std::setprecision(6) << " Mean num << 273 G4cout << " Mean number of e+ " << (G4double)fN_pos*norm << G4endl;
>> 274 G4cout << std::setprecision(6)
>> 275 << " Mean number of charged steps " << fChargedStep << G4endl;
323 G4cout << " Mean number of neutral steps " 276 G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
324 G4cout << "--------------------------------- 277 G4cout << "------------------------------------------------------------\n";
325 << 278
326 // Energy flow << 279 //Energy flow
327 // 280 //
328 G4AnalysisManager* analysis = G4AnalysisMana 281 G4AnalysisManager* analysis = G4AnalysisManager::Instance();
329 G4int Idmax = (fDetector->GetNbOfLayers()) * << 282 G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
330 for (G4int Id = 1; Id <= Idmax + 1; Id++) { << 283 for (G4int Id=1; Id<=Idmax+1; Id++) {
331 analysis->FillH1(2 * kMaxAbsor + 1, (G4dou << 284 analysis->FillH1(2*kMaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
332 analysis->FillH1(2 * kMaxAbsor + 2, (G4dou << 285 analysis->FillH1(2*kMaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
333 } 286 }
334 << 287
335 // Energy deposit from energy flow balance << 288 //Energy deposit from energy flow balance
336 // 289 //
337 G4double EdepTot[kMaxAbsor]; 290 G4double EdepTot[kMaxAbsor];
338 for (G4int k = 0; k < kMaxAbsor; k++) << 291 for (G4int k=0; k<kMaxAbsor; k++) EdepTot[k] = 0.;
339 EdepTot[k] = 0.; << 292
340 <<
341 G4int nbOfAbsor = fDetector->GetNbOfAbsor(); 293 G4int nbOfAbsor = fDetector->GetNbOfAbsor();
342 for (G4int Id = 1; Id <= Idmax; Id++) { << 294 for (G4int Id=1; Id<=Idmax; Id++) {
343 G4int iAbsor = Id % nbOfAbsor; << 295 G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
344 if (iAbsor == 0) iAbsor = nbOfAbsor; << 296 EdepTot[iAbsor] += (fEnergyFlow[Id]-fEnergyFlow[Id+1]-fLateralEleak[Id]);
345 EdepTot[iAbsor] += (fEnergyFlow[Id] - fEne << 297 }
346 } << 298
347 << 299 G4cout << std::setprecision(3)
348 G4cout << std::setprecision(3) << "\n Energy << 300 << "\n Energy deposition from Energy flow balance : \n"
349 << std::setw(10) << " material \t To 301 << std::setw(10) << " material \t Total Edep \n \n";
350 G4cout.precision(6); 302 G4cout.precision(6);
351 << 303
352 for (G4int k = 1; k <= nbOfAbsor; k++) { << 304 for (G4int k=1; k<=nbOfAbsor; k++) {
353 EdepTot[k] *= norm; << 305 EdepTot [k] *= norm;
354 G4cout << std::setw(10) << fDetector->GetA 306 G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
355 << "\t " << G4BestUnit(EdepTot[k], << 307 << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
356 } 308 }
357 << 309
358 G4cout << "\n------------------------------- << 310 G4cout << "\n------------------------------------------------------------\n"
>> 311 << G4endl;
359 312
360 // Acceptance 313 // Acceptance
361 EmAcceptance acc; 314 EmAcceptance acc;
362 G4bool isStarted = false; 315 G4bool isStarted = false;
363 for (G4int j = 1; j <= fDetector->GetNbOfAbs << 316 for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
364 if (fLimittrue[j] < DBL_MAX) { 317 if (fLimittrue[j] < DBL_MAX) {
365 if (!isStarted) { 318 if (!isStarted) {
366 acc.BeginOfAcceptance("Sampling Calori << 319 acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
367 isStarted = true; 320 isStarted = true;
368 } 321 }
369 MeanEAbs = fSumEAbs[j]; 322 MeanEAbs = fSumEAbs[j];
370 rmsEAbs = fSum2EAbs[j]; << 323 rmsEAbs = fSum2EAbs[j];
371 G4String mat = fDetector->GetAbsorMateri 324 G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
372 acc.EmAcceptanceGauss("Edep" + mat, nEvt << 325 acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
373 acc.EmAcceptanceGauss("Erms" + mat, nEvt << 326 fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
374 2.0 * fLimittrue[j << 327 acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
>> 328 fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
375 } 329 }
376 } 330 }
377 if (isStarted) acc.EndOfAcceptance(); << 331 if(isStarted) acc.EndOfAcceptance();
378 332
379 // normalize histograms << 333 //normalize histograms
380 // 334 //
381 for (G4int ih = kMaxAbsor + 1; ih < kMaxHist << 335 for (G4int ih = kMaxAbsor+1; ih < kMaxHisto; ih++) {
382 analysis->ScaleH1(ih, norm / MeV); << 336 analysis->ScaleH1(ih,norm/MeV);
383 } 337 }
384 << 338
385 G4cout.setf(mode, std::ios::floatfield); << 339 G4cout.setf(mode,std::ios::floatfield);
386 G4cout.precision(prec); 340 G4cout.precision(prec);
387 } 341 }
388 342
389 //....oooOO0OOooo........oooOO0OOooo........oo 343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 344
391 void Run::SetEdepAndRMS(G4int i, G4double edep 345 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
392 { 346 {
393 if (i >= 0 && i < kMaxAbsor) { << 347 if (i>=0 && i<kMaxAbsor) {
394 fEdeptrue[i] = edep; << 348 fEdeptrue [i] = edep;
395 fRmstrue[i] = rms; << 349 fRmstrue [i] = rms;
396 fLimittrue[i] = lim; 350 fLimittrue[i] = lim;
397 } 351 }
398 } 352 }
399 353
400 //....oooOO0OOooo........oooOO0OOooo........oo 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
401 355
402 void Run::SetApplyLimit(G4bool val) 356 void Run::SetApplyLimit(G4bool val)
403 { 357 {
404 fApplyLimit = val; 358 fApplyLimit = val;
405 } 359 }
406 360
407 //....oooOO0OOooo........oooOO0OOooo........oo 361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 362