Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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/TestEm9/src/HistoManager.cc
 27 /// \brief Implementation of the HistoManager class
 28 //
 29 //
 30 //---------------------------------------------------------------------------
 31 //
 32 // ClassName:   HistoManager
 33 //
 34 //
 35 // Author:      V.Ivanchenko 30/01/01
 36 //
 37 //----------------------------------------------------------------------------
 38 //
 39 
 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 42 
 43 #include "HistoManager.hh"
 44 
 45 #include "EmAcceptance.hh"
 46 #include "Histo.hh"
 47 
 48 #include "G4Electron.hh"
 49 #include "G4EmProcessSubType.hh"
 50 #include "G4Gamma.hh"
 51 #include "G4GammaGeneralProcess.hh"
 52 #include "G4MaterialCutsCouple.hh"
 53 #include "G4Positron.hh"
 54 #include "G4SystemOfUnits.hh"
 55 #include "G4UnitsTable.hh"
 56 #include "G4VEmProcess.hh"
 57 #include "G4VEnergyLossProcess.hh"
 58 #include "G4VProcess.hh"
 59 
 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61 
 62 HistoManager* HistoManager::fManager = nullptr;
 63 
 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 65 
 66 HistoManager* HistoManager::GetPointer()
 67 {
 68   if (nullptr == fManager) {
 69     fManager = new HistoManager();
 70   }
 71   return fManager;
 72 }
 73 
 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75 
 76 HistoManager::HistoManager()
 77   : fGamma(G4Gamma::Gamma()),
 78     fElectron(G4Electron::Electron()),
 79     fPositron(G4Positron::Positron()),
 80     fHisto(new Histo())
 81 {
 82   fVerbose = 1;
 83   fEvt1 = -1;
 84   fEvt2 = -1;
 85   fNmax = 3;
 86   fMaxEnergy = 50.0 * MeV;
 87   fBeamEnergy = 1. * GeV;
 88   fMaxEnergyAbs = 10. * MeV;
 89   fBinsE = 100;
 90   fBinsEA = 40;
 91   fBinsED = 100;
 92   fNHisto = 13;
 93 
 94   BookHisto();
 95 }
 96 
 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98 
 99 HistoManager::~HistoManager()
100 {
101   delete fHisto;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105 
106 void HistoManager::BookHisto()
107 {
108   fHisto->Add1D("10", "Evis/E0 in central crystal", fBinsED, 0.0, 1, 1.0);
109   fHisto->Add1D("11", "Evis/E0 in 3x3", fBinsED, 0.0, 1.0, 1.0);
110   fHisto->Add1D("12", "Evis/E0 in 5x5", fBinsED, 0.0, 1.0, 1.0);
111   fHisto->Add1D("13", "Energy (MeV) of delta-electrons", fBinsE, 0.0, fMaxEnergy, MeV);
112   fHisto->Add1D("14", "Energy (MeV) of gammas", fBinsE, 0.0, fMaxEnergy, MeV);
113   fHisto->Add1D("15", "Energy (MeV) in abs1", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
114   fHisto->Add1D("16", "Energy (MeV) in abs2", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
115   fHisto->Add1D("17", "Energy (MeV) in abs3", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
116   fHisto->Add1D("18", "Energy (MeV) in abs4", fBinsEA, 0.0, fMaxEnergyAbs, MeV);
117   fHisto->Add1D("19", "Number of vertex hits", 20, -0.5, 19.5, 1.0);
118   fHisto->Add1D("20", "E1/E9 Ratio", fBinsED, 0.0, 1, 1.0);
119   fHisto->Add1D("21", "E1/E25 Ratio", fBinsED, 0.0, 1.0, 1.0);
120   fHisto->Add1D("22", "E9/E25 Ratio", fBinsED, 0.0, 1.0, 1.0);
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
125 void HistoManager::BeginOfRun()
126 {
127   // initilise scoring
128   fEvt = 0;
129   fElec = 0;
130   fPosit = 0;
131   fGam = 0;
132   fStep = 0;
133   fLowe = 0;
134 
135   for (G4int i = 0; i < 6; i++) {
136     fStat[i] = 0;
137     fEdep[i] = 0.0;
138     fErms[i] = 0.0;
139     if (i < 3) {
140       fEdeptr[i] = 0.0;
141       fErmstr[i] = 0.0;
142     }
143   }
144 
145   // initialise counters
146   fBrem.resize(93, 0.0);
147   fPhot.resize(93, 0.0);
148   fComp.resize(93, 0.0);
149   fConv.resize(93, 0.0);
150 
151   // initialise acceptance - by default is not applied
152   for (G4int i = 0; i < fNmax; i++) {
153     fEdeptrue[i] = 1.0;
154     fRmstrue[i] = 1.0;
155     fLimittrue[i] = 10.;
156   }
157 
158   if (fHisto->IsActive()) {
159     for (G4int i = 0; i < fNHisto; ++i) {
160       fHisto->Activate(i, true);
161     }
162     fHisto->Book();
163 
164     if (fVerbose > 0) {
165       G4cout << "HistoManager: Histograms are booked and run has been started" << G4endl;
166     }
167   }
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
172 void HistoManager::EndOfRun(G4int runID)
173 {
174   G4cout << "HistoManager: End of run actions are started   RunID# " << runID << G4endl;
175   G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
176 
177   // average
178 
179   G4cout << "=================================================================" << G4endl;
180   G4double x = (G4double)fEvt;
181   if (fEvt > 0) x = 1.0 / x;
182   G4int j;
183   for (j = 0; j < fNmax; j++) {
184     // total mean
185     fEdep[j] *= x;
186     G4double y = fErms[j] * x - fEdep[j] * fEdep[j];
187     if (y < 0.0) y = 0.0;
188     fErms[j] = std::sqrt(y);
189 
190     // trancated mean
191     G4double xx = G4double(fStat[j]);
192     if (xx > 0.0) xx = 1.0 / xx;
193     fEdeptr[j] *= xx;
194     y = fErmstr[j] * xx - fEdeptr[j] * fEdeptr[j];
195     if (y < 0.0) y = 0.0;
196     fErmstr[j] = std::sqrt(y);
197   }
198   G4double xe = x * (G4double)fElec;
199   G4double xg = x * (G4double)fGam;
200   G4double xp = x * (G4double)fPosit;
201   G4double xs = x * fStep;
202 
203   G4double f = 100. * std::sqrt(fBeamEnergy / GeV);
204 
205   G4cout << "Number of events             " << fEvt << G4endl;
206   G4cout << std::setprecision(4) << "Average number of e-         " << xe << G4endl;
207   G4cout << std::setprecision(4) << "Average number of gamma      " << xg << G4endl;
208   G4cout << std::setprecision(4) << "Average number of e+         " << xp << G4endl;
209   G4cout << std::setprecision(4) << "Average number of steps      " << xs << G4endl;
210 
211   for (j = 0; j < 3; ++j) {
212     G4double ex = fEdeptr[j];
213     G4double sx = fErmstr[j];
214     G4double xx = G4double(fStat[j]);
215     if (xx > 0.0) xx = 1.0 / xx;
216     G4double r = sx * std::sqrt(xx);
217     G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << ex << " +- "
218            << r;
219     if (ex > 0.1) G4cout << "  res=  " << f * sx / ex << " %   " << fStat[j];
220     G4cout << G4endl;
221   }
222   if (fLimittrue[0] < 10. || fLimittrue[1] < 10. || fLimittrue[2] < 10.) {
223     G4cout << "===========  Mean values without trancating =====================" << G4endl;
224     for (j = 0; j < fNmax; j++) {
225       G4double ex = fEdep[j];
226       G4double sx = fErms[j];
227       G4double rx = sx * std::sqrt(x);
228       G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << ex << " +- "
229              << rx;
230       if (ex > 0.0) G4cout << "  res=  " << f * sx / ex << " %";
231       G4cout << G4endl;
232     }
233   }
234   G4cout << "===========  Ratios without trancating ===========================" << G4endl;
235   for (j = 3; j < 6; ++j) {
236     G4double e = fEdep[j];
237     G4double xx = G4double(fStat[j]);
238     if (xx > 0.0) xx = 1.0 / xx;
239     e *= xx;
240     G4double y = fErms[j] * xx - e * e;
241     G4double r = 0.0;
242     if (y > 0.0) r = std::sqrt(y * xx);
243     G4cout << "  " << nam[j] << " =                   " << e << " +- " << r;
244     G4cout << G4endl;
245   }
246   G4cout << std::setprecision(4) << "Beam Energy                  " << fBeamEnergy / GeV << " GeV"
247          << G4endl;
248   if (fLowe > 0) G4cout << "Number of events E/E0<0.8    " << fLowe << G4endl;
249   G4cout << "==================================================================" << G4endl;
250   G4cout << G4endl;
251 
252   // normalise histograms
253   if (fHisto->IsActive()) {
254     for (G4int i = 0; i < fNHisto; ++i) {
255       fHisto->ScaleH1(i, x);
256     }
257     fHisto->Save();
258   }
259   if (0 < runID) {
260     return;
261   }
262 
263   // Acceptance only for the first run
264   EmAcceptance acc;
265   G4bool isStarted = false;
266   for (j = 0; j < fNmax; j++) {
267     G4double ltrue = fLimittrue[j];
268     if (ltrue < DBL_MAX) {
269       if (!isStarted) {
270         acc.BeginOfAcceptance("Crystal Calorimeter", fEvt);
271         isStarted = true;
272       }
273       G4double etrue = fEdeptrue[j];
274       G4double rtrue = fRmstrue[j];
275       acc.EmAcceptanceGauss("Edep" + nam[j], fEvt, fEdeptr[j], etrue, rtrue, ltrue);
276       acc.EmAcceptanceGauss("Erms" + nam[j], fEvt, fErmstr[j], rtrue, rtrue, 2.0 * ltrue);
277     }
278   }
279   if (isStarted) acc.EndOfAcceptance();
280 
281   // atom frequency
282   G4cout << "   Z  bremsstrahlung photoeffect  compton    conversion" << G4endl;
283   for (j = 1; j < 93; ++j) {
284     G4int n1 = G4int(fBrem[j] * x);
285     G4int n2 = G4int(fPhot[j] * x);
286     G4int n3 = G4int(fComp[j] * x);
287     G4int n4 = G4int(fConv[j] * x);
288     if (n1 + n2 + n3 + n4 > 0) {
289       G4cout << std::setw(4) << j << std::setw(12) << n1 << std::setw(12) << n2 << std::setw(12)
290              << n3 << std::setw(12) << n4 << G4endl;
291     }
292   }
293 }
294 
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296 
297 void HistoManager::BeginOfEvent()
298 {
299   ++fEvt;
300 
301   fEabs1 = 0.0;
302   fEabs2 = 0.0;
303   fEabs3 = 0.0;
304   fEabs4 = 0.0;
305   fEvertex.clear();
306   fNvertex.clear();
307   for (G4int i = 0; i < 25; i++) {
308     fE[i] = 0.0;
309   }
310 }
311 
312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
313 
314 void HistoManager::EndOfEvent()
315 {
316   G4double e9 = 0.0;
317   G4double e25 = 0.0;
318   for (G4int i = 0; i < 25; i++) {
319     fE[i] /= fBeamEnergy;
320     e25 += fE[i];
321     if ((6 <= i && 8 >= i) || (11 <= i && 13 >= i) || (16 <= i && 18 >= i)) e9 += fE[i];
322   }
323 
324   if (1 < fVerbose && e25 < 0.8) {
325     ++fLowe;
326     G4cout << "### in the event# " << fEvt << "  E25= " << e25 << G4endl;
327   }
328 
329   // compute ratios
330   G4double e0 = fE[12];
331   G4double e19 = 0.0;
332   G4double e125 = 0.0;
333   G4double e925 = 0.0;
334   if (e9 > 0.0) {
335     e19 = e0 / e9;
336     e125 = e0 / e25;
337     e925 = e9 / e25;
338     fEdep[3] += e19;
339     fErms[3] += e19 * e19;
340     fEdep[4] += e125;
341     fErms[4] += e125 * e125;
342     fEdep[5] += e925;
343     fErms[5] += e925 * e925;
344     fStat[3] += 1;
345     fStat[4] += 1;
346     fStat[5] += 1;
347   }
348 
349   // Fill histo
350   fHisto->Fill(0, e0, 1.0);
351   fHisto->Fill(1, e9, 1.0);
352   fHisto->Fill(2, e25, 1.0);
353   fHisto->Fill(5, fEabs1, 1.0);
354   fHisto->Fill(6, fEabs2, 1.0);
355   fHisto->Fill(7, fEabs3, 1.0);
356   fHisto->Fill(8, fEabs4, 1.0);
357   fHisto->Fill(9, G4double(fNvertex.size()), 1.0);
358   fHisto->Fill(10, e19, 1.0);
359   fHisto->Fill(11, e125, 1.0);
360   fHisto->Fill(12, e925, 1.0);
361 
362   // compute sums
363   fEdep[0] += e0;
364   fErms[0] += e0 * e0;
365   fEdep[1] += e9;
366   fErms[1] += e9 * e9;
367   fEdep[2] += e25;
368   fErms[2] += e25 * e25;
369 
370   // trancated mean
371   if (std::abs(e0 - fEdeptrue[0]) < fRmstrue[0] * fLimittrue[0]) {
372     fStat[0] += 1;
373     fEdeptr[0] += e0;
374     fErmstr[0] += e0 * e0;
375   }
376   if (std::abs(e9 - fEdeptrue[1]) < fRmstrue[1] * fLimittrue[1]) {
377     fStat[1] += 1;
378     fEdeptr[1] += e9;
379     fErmstr[1] += e9 * e9;
380   }
381   if (std::abs(e25 - fEdeptrue[2]) < fRmstrue[2] * fLimittrue[2]) {
382     fStat[2] += 1;
383     fEdeptr[2] += e25;
384     fErmstr[2] += e25 * e25;
385   }
386 }
387 
388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 
390 void HistoManager::ScoreNewTrack(const G4Track* aTrack)
391 {
392   // Save primary parameters
393   ResetTrackLength();
394   const G4ParticleDefinition* particle = aTrack->GetDefinition();
395   const G4DynamicParticle* dynParticle = aTrack->GetDynamicParticle();
396 
397   G4int pid = aTrack->GetParentID();
398   G4double kinE = dynParticle->GetKineticEnergy();
399   G4ThreeVector pos = aTrack->GetVertexPosition();
400 
401   // primary
402   if (0 == pid) {
403     G4double mass = 0.0;
404     if (particle) {
405       mass = particle->GetPDGMass();
406     }
407 
408     G4ThreeVector dir = dynParticle->GetMomentumDirection();
409     if (1 < fVerbose) {
410       G4cout << "TrackingAction: Primary kinE(MeV)= " << kinE / MeV << "; m(MeV)= " << mass / MeV
411              << "; pos= " << pos << ";  dir= " << dir << G4endl;
412     }
413 
414     // secondary
415   }
416   else {
417     const G4VProcess* proc = aTrack->GetCreatorProcess();
418     G4int type = proc->GetProcessSubType();
419 
420     if (type == fBremsstrahlung) {
421       auto elm = static_cast<const G4VEnergyLossProcess*>(proc)->GetCurrentElement();
422       if (nullptr != elm) {
423         G4int Z = elm->GetZasInt();
424         if (Z > 0 && Z < 93) {
425           fBrem[Z] += 1.0;
426         }
427       }
428     }
429     else if (type == fPhotoElectricEffect) {
430       auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
431       if (nullptr != elm) {
432         G4int Z = elm->GetZasInt();
433         if (Z > 0 && Z < 93) {
434           fPhot[Z] += 1.0;
435         }
436       }
437     }
438     else if (type == fGammaConversion) {
439       auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
440       if (nullptr != elm) {
441         G4int Z = elm->GetZasInt();
442         if (Z > 0 && Z < 93) {
443           fConv[Z] += 1.0;
444         }
445       }
446     }
447     else if (type == fComptonScattering) {
448       auto elm = static_cast<const G4VEmProcess*>(proc)->GetCurrentElement();
449       if (nullptr != elm) {
450         G4int Z = elm->GetZasInt();
451         if (Z > 0 && Z < 93) {
452           fComp[Z] += 1.0;
453         }
454       }
455     }
456 
457     // delta-electron
458     if (particle == fElectron) {
459       if (1 < fVerbose) {
460         G4cout << "TrackingAction: Secondary electron " << G4endl;
461       }
462       AddDeltaElectron(dynParticle);
463     }
464     else if (particle == fPositron) {
465       if (1 < fVerbose) {
466         G4cout << "TrackingAction: Secondary positron " << G4endl;
467       }
468       AddPositron(dynParticle);
469     }
470     else if (particle == fGamma) {
471       if (1 < fVerbose) {
472         G4cout << "TrackingAction: Secondary gamma; parentID= " << pid << " E= " << kinE << G4endl;
473       }
474       AddPhoton(dynParticle);
475     }
476   }
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
480 
481 void HistoManager::AddEnergy(G4double edep, G4int volIndex, G4int copyNo)
482 {
483   if (1 < fVerbose) {
484     G4cout << "HistoManager::AddEnergy: e(keV)= " << edep / keV << "; volIdx= " << volIndex
485            << "; copyNo= " << copyNo << G4endl;
486   }
487   if (0 == volIndex) {
488     fE[copyNo] += edep;
489   }
490   else if (1 == volIndex) {
491     fEabs1 += edep;
492   }
493   else if (2 == volIndex) {
494     fEabs2 += edep;
495   }
496   else if (3 == volIndex) {
497     fEabs3 += edep;
498   }
499   else if (4 == volIndex) {
500     fEabs4 += edep;
501   }
502   else if (5 == volIndex) {
503     G4int n = fNvertex.size();
504     G4bool newPad = true;
505     if (n > 0) {
506       for (G4int i = 0; i < n; i++) {
507         if (copyNo == fNvertex[i]) {
508           newPad = false;
509           fEvertex[i] += edep;
510           break;
511         }
512       }
513     }
514     if (newPad) {
515       fNvertex.push_back(copyNo);
516       fEvertex.push_back(edep);
517     }
518   }
519 }
520 
521 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
522 
523 void HistoManager::AddDeltaElectron(const G4DynamicParticle* elec)
524 {
525   G4double e = elec->GetKineticEnergy() / MeV;
526   if (e > 0.0) {
527     ++fElec;
528     fHisto->Fill(3, e, 1.0);
529   }
530 }
531 
532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
533 
534 void HistoManager::AddPhoton(const G4DynamicParticle* ph)
535 {
536   G4double e = ph->GetKineticEnergy() / MeV;
537   if (e > 0.0) {
538     ++fGam;
539     fHisto->Fill(4, e, 1.0);
540   }
541 }
542 
543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
544 
545 void HistoManager::SetEdepAndRMS(G4int i, const G4ThreeVector& val)
546 {
547   if (i < fNmax && i >= 0) {
548     if (val[0] > 0.0) fEdeptrue[i] = val[0];
549     if (val[1] > 0.0) fRmstrue[i] = val[1];
550     if (val[2] > 0.0) fLimittrue[i] = val[2];
551   }
552 }
553 
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
555