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