Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file.
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 //
28 //
29 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 #include "EmAcceptance.hh"
39 
40 #include "G4ParticleTable.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4Track.hh"
43 #include "G4Gamma.hh"
44 #include "G4Electron.hh"
45 #include "G4Positron.hh"
46 
47 #include "G4UnitsTable.hh"
48 #include "G4SystemOfUnits.hh"
49 
50 #include <iomanip>
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53 
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)
61 {
62  //initialize cumulative quantities
63  //
64  for (G4int k=0; k<kMaxAbsor; k++) {
65  fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
66  fEnergyDeposit[k].clear();
67  fEdeptrue[k] = fRmstrue[k] = 1.;
68  fLimittrue[k] = DBL_MAX;
69  }
70 
71  //initialize Eflow
72  //
73  G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
74  fEnergyFlow.resize(nbPlanes);
75  fLateralEleak.resize(nbPlanes);
76  for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
81 Run::~Run()
82 { }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
87 {
88  fParticle = particle;
89  fEkin = energy;
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 
95 {
96  //accumulate statistic with restriction
97  //
98  if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
99  fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs;
100  fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
106 {
107  fEnergyFlow[plane] += Eflow;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
113 {
114  fLateralEleak[cell] += Eflow;
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118 
120 {
121  fChargedStep += 1.0;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 
127 {
128  fNeutralStep += 1.0;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 void Run::AddSecondaryTrack(const G4Track* track)
134 {
135  const G4ParticleDefinition* d = track->GetDefinition();
136  if(d == G4Gamma::Gamma()) { ++fN_gamma; }
137  else if (d == G4Electron::Electron()) { ++fN_elec; }
138  else if (d == G4Positron::Positron()) { ++fN_pos; }
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 
143 void Run::Merge(const G4Run* run)
144 {
145  const Run* localRun = static_cast<const Run*>(run);
146 
147  // pass information about primary particle
148  fParticle = localRun->fParticle;
149  fEkin = localRun->fEkin;
150 
151  // accumulate sums
152  //
153  for (G4int k=0; k<kMaxAbsor; k++) {
154  fSumEAbs[k] += localRun->fSumEAbs[k];
155  fSum2EAbs[k] += localRun->fSum2EAbs[k];
156  fSumLAbs[k] += localRun->fSumLAbs[k];
157  fSum2LAbs[k] += localRun->fSum2LAbs[k];
158  }
159 
160  G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
161  for (G4int k=0; k<nbPlanes; k++) {
162  fEnergyFlow[k] += localRun->fEnergyFlow[k];
163  fLateralEleak[k] += localRun->fLateralEleak[k];
164  }
165 
166 
167  fChargedStep += localRun->fChargedStep;
168  fNeutralStep += localRun->fNeutralStep;
169 
170  fN_gamma += localRun->fN_gamma;
171  fN_elec += localRun->fN_elec;
172  fN_pos += localRun->fN_pos;
173 
174  fApplyLimit = localRun->fApplyLimit;
175 
176  for (G4int k=0; k<kMaxAbsor; k++) {
177  fEdeptrue[k] = localRun->fEdeptrue[k];
178  fRmstrue[k] = localRun->fRmstrue[k];
179  fLimittrue[k] = localRun->fLimittrue[k];
180  }
181 
182  G4Run::Merge(run);
183 }
184 
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186 
187 void Run::EndOfRun()
188 {
189  G4int nEvt = numberOfEvent;
190  G4double norm = G4double(nEvt);
191  if(norm > 0) norm = 1./norm;
192  G4double qnorm = std::sqrt(norm);
193 
194  fChargedStep *= norm;
195  fNeutralStep *= norm;
196 
197  //compute and print statistic
198  //
199  G4double beamEnergy = fEkin;
200  G4double sqbeam = std::sqrt(beamEnergy/GeV);
201 
202  G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
203  G4double MeanLAbs,MeanLAbs2,rmsLAbs;
204 
205  std::ios::fmtflags mode = G4cout.flags();
206  G4int prec = G4cout.precision(2);
207  G4cout << "\n------------------------------------------------------------\n";
208  G4cout << std::setw(14) << "material"
209  << std::setw(17) << "Edep RMS"
210  << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
211  << std::setw(23) << "total tracklen \n \n";
212 
213  for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
214  {
215  MeanEAbs = fSumEAbs[k]*norm;
216  MeanEAbs2 = fSum2EAbs[k]*norm;
217  rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
218  //G4cout << "k= " << k << " RMS= " << rmsEAbs
219  // << " fApplyLimit: " << fApplyLimit << G4endl;
220  if(fApplyLimit) {
221  G4int nn = 0;
222  G4double sume = 0.0;
223  G4double sume2 = 0.0;
224  // compute trancated means
225  G4double lim = rmsEAbs * 2.5;
226  for(G4int i=0; i<nEvt; i++) {
227  G4double e = (fEnergyDeposit[k])[i];
228  if(std::abs(e - MeanEAbs) < lim) {
229  sume += e;
230  sume2 += e*e;
231  nn++;
232  }
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));
239  }
240 
241  resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
242  rmsres = resolution*qnorm;
243 
244  // Save mean and RMS
245  fSumEAbs[k] = MeanEAbs;
246  fSum2EAbs[k] = rmsEAbs;
247 
248  MeanLAbs = fSumLAbs[k]*norm;
249  MeanLAbs2 = fSum2LAbs[k]*norm;
250  rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
251 
252  //print
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";
268 
269  G4cout << " Beam particle "
270  << fParticle->GetParticleName()
271  << " E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
272  G4cout << " Mean number of gamma " << (G4double)fN_gamma*norm << G4endl;
273  G4cout << " Mean number of e- " << (G4double)fN_elec*norm << G4endl;
274  G4cout << " Mean number of e+ " << (G4double)fN_pos*norm << G4endl;
275  G4cout << std::setprecision(6)
276  << " Mean number of charged steps " << fChargedStep << G4endl;
277  G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
278  G4cout << "------------------------------------------------------------\n";
279 
280  //Energy flow
281  //
282  G4AnalysisManager* analysis = G4AnalysisManager::Instance();
283  G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
284  for (G4int Id=1; Id<=Idmax+1; Id++) {
285  analysis->FillH1(2*kMaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
286  analysis->FillH1(2*kMaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
287  }
288 
289  //Energy deposit from energy flow balance
290  //
291  G4double EdepTot[kMaxAbsor];
292  for (G4int k=0; k<kMaxAbsor; k++) EdepTot[k] = 0.;
293 
294  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
295  for (G4int Id=1; Id<=Idmax; Id++) {
296  G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
297  EdepTot[iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]);
298  }
299 
300  G4cout << std::setprecision(3)
301  << "\n Energy deposition from Energy flow balance : \n"
302  << std::setw(10) << " material \t Total Edep \n \n";
303  G4cout.precision(6);
304 
305  for (G4int k=1; k<=nbOfAbsor; k++) {
306  EdepTot [k] *= norm;
307  G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
308  << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
309  }
310 
311  G4cout << "\n------------------------------------------------------------\n"
312  << G4endl;
313 
314  // Acceptance
315  EmAcceptance acc;
316  G4bool isStarted = false;
317  for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
318  if (fLimittrue[j] < DBL_MAX) {
319  if (!isStarted) {
320  acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
321  isStarted = true;
322  }
323  MeanEAbs = fSumEAbs[j];
324  rmsEAbs = fSum2EAbs[j];
325  G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
326  acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
327  fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
328  acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
329  fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
330  }
331  }
332  if(isStarted) acc.EndOfAcceptance();
333 
334  //normalize histograms
335  //
336  for (G4int ih = kMaxAbsor+1; ih < kMaxHisto; ih++) {
337  analysis->ScaleH1(ih,norm/MeV);
338  }
339 
340  G4cout.setf(mode,std::ios::floatfield);
341  G4cout.precision(prec);
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345 
347 {
348  if (i>=0 && i<kMaxAbsor) {
349  fEdeptrue [i] = edep;
350  fRmstrue [i] = rms;
351  fLimittrue[i] = lim;
352  }
353 }
354 
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
356 
358 {
359  fApplyLimit = val;
360 }
361 
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ParticleDefinition * GetDefinition() const
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
Run()
Definition: Run.cc:43
void EmAcceptanceGauss(const G4String &title, G4int stat, G4double avr, G4double avr0, G4double rms, G4double limit)
Definition: EmAcceptance.cc:68
const G4String & GetName() const
Definition: G4Material.hh:178
void SumEnergyFlow(G4int plane, G4double Eflow)
Definition: Run.cc:105
void BeginOfAcceptance(const G4String &title, G4int stat)
Definition: EmAcceptance.cc:49
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void SetEdepAndRMS(G4int, G4double, G4double, G4double)
Definition: Run.cc:346
int G4int
Definition: G4Types.hh:78
void AddChargedStep()
Definition: Run.cc:119
const G4String & GetParticleName() const
G4Material * GetAbsorMaterial(G4int i)
void AddNeutralStep()
Definition: Run.cc:126
static const double prec
Definition: RanecuEngine.cc:58
void SumLateralEleak(G4int cell, G4double Eflow)
Definition: Run.cc:112
G4GLOB_DLL std::ostream G4cout
void AddSecondaryTrack(const G4Track *)
Definition: Run.cc:133
void EndOfAcceptance()
Definition: EmAcceptance.cc:58
bool G4bool
Definition: G4Types.hh:79
const G4int kMaxAbsor
Definition: G4Run.hh:46
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4double energy(const ThreeVector &p, const G4double m)
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
Detector construction class to define materials and geometry.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
Definition: Run.cc:77
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
void FillPerEvent()
Definition: Run.cc:112
const G4int kMaxHisto
Definition: HistoManager.hh:46
#define DBL_MAX
Definition: templates.hh:83
virtual void Merge(const G4Run *)
Definition: Run.cc:115
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
~Run()
Definition: Run.cc:72
void SetApplyLimit(G4bool)
Definition: Run.cc:357
const G4int kAbs
Definition: B4bRunData.hh:41