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