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