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