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 // ClassName: HistoManager
32 //
33 //
34 // Author: V.Ivanchenko 30/01/01
35 //
36 // Modified:
37 // 04.06.2006 Adoptation of Hadr01 (V.Ivanchenko)
38 // 16.11.2006 Add beamFlag (V.Ivanchenko)
39 //
40 //----------------------------------------------------------------------------
41 //
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 
46 #include "HistoManager.hh"
47 #include "G4Track.hh"
48 #include "G4Step.hh"
49 #include "G4UnitsTable.hh"
50 #include "G4Neutron.hh"
51 #include "G4Proton.hh"
52 #include "G4AntiProton.hh"
53 #include "G4Gamma.hh"
54 #include "G4Electron.hh"
55 #include "G4Positron.hh"
56 #include "G4MuonPlus.hh"
57 #include "G4MuonMinus.hh"
58 #include "G4PionPlus.hh"
59 #include "G4PionMinus.hh"
60 #include "G4PionZero.hh"
61 #include "G4KaonPlus.hh"
62 #include "G4KaonMinus.hh"
63 #include "G4KaonZeroShort.hh"
64 #include "G4KaonZeroLong.hh"
65 #include "G4Deuteron.hh"
66 #include "G4Triton.hh"
67 #include "G4He3.hh"
68 #include "G4Alpha.hh"
69 #include "Histo.hh"
70 #include "globals.hh"
71 #include "G4SystemOfUnits.hh"
72 #include "G4PhysicalConstants.hh"
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 HistoManager* HistoManager::fManager = 0;
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
80 HistoManager* HistoManager::GetPointer()
81 {
82  if(!fManager) {
83  static HistoManager manager;
84  fManager = &manager;
85  }
86  return fManager;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
91 HistoManager::HistoManager()
92 : fPrimaryDef(nullptr),
93  fEdepMax(1.0*GeV),
94  fLength (300.*mm),
95  fPrimaryKineticEnergy(0.0),
96  fVerbose(0),
97  fNBinsE (100),
98  fNSlices(300),
99  fNHisto (28),
100  fBeamFlag(true),
101  fHistoBooked(false)
102 {
103  fHisto = new Histo();
104  fHisto->SetVerbose(fVerbose);
105  BookHisto();
106  fNeutron = G4Neutron::Neutron();
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110 
111 HistoManager::~HistoManager()
112 {
113  delete fHisto;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
118 void HistoManager::BookHisto()
119 {
120  fHistoBooked = true;
121  fHisto->Add1D("1","Energy deposition (MeV/mm/event) in the target",
122  fNSlices,0.0,fLength/mm,MeV/mm);
123  fHisto->Add1D("2","Log10 Energy (MeV) of gammas",fNBinsE,-5.,5.,1.0);
124  fHisto->Add1D("3","Log10 Energy (MeV) of electrons",fNBinsE,-5.,5.,1.0);
125  fHisto->Add1D("4","Log10 Energy (MeV) of positrons",fNBinsE,-5.,5.,1.0);
126  fHisto->Add1D("5","Log10 Energy (MeV) of protons",fNBinsE,-5.,5.,1.0);
127  fHisto->Add1D("6","Log10 Energy (MeV) of neutrons",fNBinsE,-5.,5.,1.0);
128  fHisto->Add1D("7","Log10 Energy (MeV) of charged pions",fNBinsE,-4.,6.,1.0);
129  fHisto->Add1D("8","Log10 Energy (MeV) of pi0",fNBinsE,-4.,6.,1.0);
130  fHisto->Add1D("9","Log10 Energy (MeV) of charged kaons",fNBinsE,-4.,6.,1.0);
131  fHisto->Add1D("10","Log10 Energy (MeV) of neutral kaons",fNBinsE,-4.,6.,1.0);
132  fHisto->Add1D("11","Log10 Energy (MeV) of deuterons and tritons",
133  fNBinsE,-5.,5.,1.0);
134  fHisto->Add1D("12","Log10 Energy (MeV) of He3 and alpha",fNBinsE,-5.,5.,1.0);
135  fHisto->Add1D("13","Log10 Energy (MeV) of Generic Ions",fNBinsE,-5.,5.,1.0);
136  fHisto->Add1D("14","Log10 Energy (MeV) of muons",fNBinsE,-4.,6.,1.0);
137  fHisto->Add1D("15","log10 Energy (MeV) of side-leaked neutrons",
138  fNBinsE,-5.,5.,1.0);
139  fHisto->Add1D("16","log10 Energy (MeV) of forward-leaked neutrons",
140  fNBinsE,-5.,5.,1.0);
141  fHisto->Add1D("17","log10 Energy (MeV) of backward-leaked neutrons",
142  fNBinsE,-5.,5.,1.0);
143  fHisto->Add1D("18","log10 Energy (MeV) of leaking protons",
144  fNBinsE,-4.,6.,1.0);
145  fHisto->Add1D("19","log10 Energy (MeV) of leaking charged pions",
146  fNBinsE,-4.,6.,1.0);
147  fHisto->Add1D("20","Log10 Energy (MeV) of pi+",fNBinsE,-4.,6.,1.0);
148  fHisto->Add1D("21","Log10 Energy (MeV) of pi-",fNBinsE,-4.,6.,1.0);
149  fHisto->Add1D("22",
150  "Energy deposition in the target normalized to beam energy",
151  110,0.0,1.1,1.0);
152  fHisto->Add1D("23",
153  "EM energy deposition in the target normalized to beam energy",
154  110,0.0,1.1,1.0);
155  fHisto->Add1D("24",
156  "Pion energy deposition in the target normalized to beam energy",
157  110,0.0,1.1,1.0);
158  fHisto->Add1D("25",
159  "Proton energy deposition in the target normalized to beam energy",
160  110,0.0,1.1,1.0);
161  fHisto->Add1D("26","Energy (MeV) of pi+",fNBinsE,0.,500.,1.0);
162  fHisto->Add1D("27","Energy (MeV) of pi-",fNBinsE,0.,500.,1.0);
163  fHisto->Add1D("28","Energy (MeV) of pi0",fNBinsE,0.,500.,1.0);
164 
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
169 void HistoManager::BeginOfRun()
170 {
171  fAbsZ0 = -0.5*fLength;
172  fNevt = 0;
173  fNelec = 0;
174  fNposit = 0;
175  fNgam = 0;
176  fNstep = 0;
177  fNprot_leak = 0;
178  fNpiofNleak = 0;
179  fNions = 0;
180  fNdeut = 0;
181  fNalpha = 0;
182  fNkaons = 0;
183  fNmuons = 0;
184  fNcpions = 0;
185  fNpi0 = 0;
186  fNneutron = 0;
187  fNproton = 0;
188  fNaproton = 0;
189  fNneu_forw = 0;
190  fNneu_leak = 0;
191  fNneu_back = 0;
192 
193  fEdepSum = 0.0;
194  fEdepSum2 = 0.0;
195 
196  if(!fHistoBooked) { BookHisto(); }
197  fHisto->Book();
198 
199  if(fVerbose > 0) {
200  G4cout << "HistoManager: Histograms are booked and run has been started"
201  <<G4endl;
202  }
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 
207 void HistoManager::EndOfRun()
208 {
209 
210  G4cout << "HistoManager: End of run actions are started" << G4endl;
211 
212  // Average values
213  G4cout<<"========================================================"<<G4endl;
214 
215  G4double x = (G4double)fNevt;
216  if(fNevt > 0) { x = 1.0/x; }
217 
218  G4double xe = x*(G4double)fNelec;
219  G4double xg = x*(G4double)fNgam;
220  G4double xp = x*(G4double)fNposit;
221  G4double xs = x*(G4double)fNstep;
222  G4double xn = x*(G4double)fNneutron;
223  G4double xpn = x*(G4double)fNproton;
224  G4double xap = x*(G4double)fNaproton;
225  G4double xnf = x*(G4double)fNneu_forw;
226  G4double xnb = x*(G4double)fNneu_leak;
227  G4double xnbw= x*(G4double)fNneu_back;
228  G4double xpl = x*(G4double)fNprot_leak;
229  G4double xal = x*(G4double)fNpiofNleak;
230  G4double xpc = x*(G4double)fNcpions;
231  G4double xp0 = x*(G4double)fNpi0;
232  G4double xpk = x*(G4double)fNkaons;
233  G4double xpm = x*(G4double)fNmuons;
234  G4double xid = x*(G4double)fNdeut;
235  G4double xia = x*(G4double)fNalpha;
236  G4double xio = x*(G4double)fNions;
237 
238  fEdepSum *= x;
239  fEdepSum2 *= x;
240  fEdepSum2 -= fEdepSum*fEdepSum;
241  if(fEdepSum2 > 0.0) { fEdepSum2 = std::sqrt(fEdepSum2); }
242  else { fEdepSum2 = 0.0; }
243 
244  G4cout << "Beam particle "
245  << fPrimaryDef->GetParticleName() <<G4endl;
246  G4cout << "Beam Energy(MeV) "
247  << fPrimaryKineticEnergy/MeV <<G4endl;
248  G4cout << "Number of events "
249  << fNevt <<G4endl;
250  G4cout << std::setprecision(4) << "Average energy deposit (MeV) "
251  << fEdepSum/MeV
252  << " RMS(MeV) " << fEdepSum2/MeV << G4endl;
253  G4cout << std::setprecision(4) << "Average number of steps "
254  << xs << G4endl;
255  G4cout << std::setprecision(4) << "Average number of gamma "
256  << xg << G4endl;
257  G4cout << std::setprecision(4) << "Average number of e- "
258  << xe << G4endl;
259  G4cout << std::setprecision(4) << "Average number of e+ "
260  << xp << G4endl;
261  G4cout << std::setprecision(4) << "Average number of neutrons "
262  << xn << G4endl;
263  G4cout << std::setprecision(4) << "Average number of protons "
264  << xpn << G4endl;
265  G4cout << std::setprecision(4) << "Average number of antiprotons "
266  << xap << G4endl;
267  G4cout << std::setprecision(4) << "Average number of pi+ & pi- "
268  << xpc << G4endl;
269  G4cout << std::setprecision(4) << "Average number of pi0 "
270  << xp0 << G4endl;
271  G4cout << std::setprecision(4) << "Average number of kaons "
272  << xpk << G4endl;
273  G4cout << std::setprecision(4) << "Average number of muons "
274  << xpm << G4endl;
275  G4cout << std::setprecision(4) << "Average number of deuterons+tritons "
276  << xid << G4endl;
277  G4cout << std::setprecision(4) << "Average number of He3+alpha "
278  << xia << G4endl;
279  G4cout << std::setprecision(4) << "Average number of ions "
280  << xio << G4endl;
281  G4cout << std::setprecision(4) << "Average number of forward neutrons "
282  << xnf << G4endl;
283  G4cout << std::setprecision(4) << "Average number of reflected neutrons "
284  << xnb << G4endl;
285  G4cout << std::setprecision(4) << "Average number of leaked neutrons "
286  << xnbw << G4endl;
287  G4cout << std::setprecision(4) << "Average number of proton leak "
288  << xpl << G4endl;
289  G4cout << std::setprecision(4) << "Average number of pion leak "
290  << xal << G4endl;
291  G4cout<<"========================================================"<<G4endl;
292  G4cout<<G4endl;
293 
294  // normalise histograms
295  for(G4int i=0; i<fNHisto; i++) {
296  fHisto->ScaleH1(i,x);
297  }
298 
299  fHisto->Save();
300 }
301 
302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 
304 void HistoManager::BeginOfEvent()
305 {
306  fEdepEvt = 0.0;
307  fEdepEM = 0.0;
308  fEdepPI = 0.0;
309  fEdepP = 0.0;
310 }
311 
312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
313 
314 void HistoManager::EndOfEvent()
315 {
316  fEdepSum += fEdepEvt;
317  fEdepSum2 += fEdepEvt*fEdepEvt;
318  fHisto->Fill(21,fEdepEvt/fPrimaryKineticEnergy,1.0);
319  fHisto->Fill(22,fEdepEM/fPrimaryKineticEnergy,1.0);
320  fHisto->Fill(23,fEdepPI/fPrimaryKineticEnergy,1.0);
321  fHisto->Fill(24,fEdepP/fPrimaryKineticEnergy,1.0);
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325 
326 void HistoManager::ScoreNewTrack(const G4Track* track)
327 {
328  const G4ParticleDefinition* pd = track->GetDefinition();
329  G4String name = pd->GetParticleName();
330  G4double ee = track->GetKineticEnergy();
331  G4double e = ee;
332 
333  // Primary track
334  if(0 == track->GetParentID()) {
335 
336  fNevt++;
337  fPrimaryKineticEnergy = e;
338  fPrimaryDef = pd;
339  G4ThreeVector dir = track->GetMomentumDirection();
340  if(1 < fVerbose)
341  G4cout << "### Primary " << name
342  << " kinE(MeV)= " << e/MeV
343  << "; m(MeV)= " << pd->GetPDGMass()/MeV
344  << "; pos(mm)= " << track->GetPosition()/mm
345  << "; dir= " << track->GetMomentumDirection()
346  << G4endl;
347 
348  // Secondary track
349  } else {
350  if(1 < fVerbose)
351  G4cout << "=== Secondary " << name
352  << " kinE(MeV)= " << e/MeV
353  << "; m(MeV)= " << pd->GetPDGMass()/MeV
354  << "; pos(mm)= " << track->GetPosition()/mm
355  << "; dir= " << track->GetMomentumDirection()
356  << G4endl;
357  e = std::log10(e/MeV);
358  if(pd == G4Gamma::Gamma()) {
359  fNgam++;
360  fHisto->Fill(1,e,1.0);
361  } else if ( pd == G4Electron::Electron()) {
362  fNelec++;
363  fHisto->Fill(2,e,1.0);
364  } else if ( pd == G4Positron::Positron()) {
365  fNposit++;
366  fHisto->Fill(3,e,1.0);
367  } else if ( pd == G4Proton::Proton()) {
368  fNproton++;
369  fHisto->Fill(4,e,1.0);
370  } else if ( pd == fNeutron) {
371  fNneutron++;
372  fHisto->Fill(5,e,1.0);
373  } else if ( pd == G4AntiProton::AntiProton()) {
374  fNaproton++;
375  } else if ( pd == G4PionPlus::PionPlus() ) {
376  fNcpions++;
377  fHisto->Fill(6,e,1.0);
378  fHisto->Fill(19,e,1.0);
379  fHisto->Fill(25,ee,1.0);
380 
381  } else if ( pd == G4PionMinus::PionMinus()) {
382  fNcpions++;
383  fHisto->Fill(6,e,1.0);
384  fHisto->Fill(20,e,1.0);
385  fHisto->Fill(26,ee,1.0);
386 
387  } else if ( pd == G4PionZero::PionZero()) {
388  fNpi0++;
389  fHisto->Fill(7,e,1.0);
390  fHisto->Fill(27,ee,1.0);
391 
392  } else if ( pd == G4KaonPlus::KaonPlus() ||
393  pd == G4KaonMinus::KaonMinus()) {
394  fNkaons++;
395  fHisto->Fill(8,e,1.0);
396  } else if ( pd == G4KaonZeroShort::KaonZeroShort() ||
397  pd == G4KaonZeroLong::KaonZeroLong()) {
398  fNkaons++;
399  fHisto->Fill(9,e,1.0);
400  } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
401  fNdeut++;
402  fHisto->Fill(10,e,1.0);
403  } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
404  fNalpha++;
405  fHisto->Fill(11,e,1.0);
406  } else if ( pd->GetParticleType() == "nucleus") {
407  fNions++;
408  fHisto->Fill(12,e,1.0);
409  } else if ( pd == G4MuonPlus::MuonPlus() ||
410  pd == G4MuonMinus::MuonMinus()) {
411  fNmuons++;
412  fHisto->Fill(13,e,1.0);
413  }
414  }
415 }
416 
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418 
419 void HistoManager::AddTargetStep(const G4Step* step)
420 {
421  fNstep++;
422  G4double fEdep = step->GetTotalEnergyDeposit();
423  if(1 < fVerbose) {
424  G4cout << "TargetSD::ProcessHits: beta1= "
425  << step->GetPreStepPoint()->GetVelocity()/c_light
426  << " beta2= " << step->GetPostStepPoint()->GetVelocity()/c_light
427  << " weight= " << step->GetTrack()->GetWeight()
428  << G4endl;
429  }
430  if(fEdep >= DBL_MIN) {
431  const G4Track* track = step->GetTrack();
432 
433  G4ThreeVector pos =
434  (step->GetPreStepPoint()->GetPosition() +
435  step->GetPostStepPoint()->GetPosition())*0.5;
436 
437  G4double z = pos.z() - fAbsZ0;
438 
439  // scoring
440  fEdepEvt += fEdep;
441  fHisto->Fill(0,z,fEdep);
442  const G4ParticleDefinition* pd = track->GetDefinition();
443 
444  if(pd == G4Gamma::Gamma() || pd == G4Electron::Electron()
445  || pd == G4Positron::Positron()) {
446  fEdepEM += fEdep;
447  } else if ( pd == G4PionPlus::PionPlus() ||
448  pd == G4PionMinus::PionMinus()) {
449  fEdepPI += fEdep;
450  } else if ( pd == G4Proton::Proton() ||
451  pd == G4AntiProton::AntiProton()) {
452  fEdepP += fEdep;
453  }
454 
455  if(1 < fVerbose) {
456  G4cout << "HistoManager::AddEnergy: e(keV)= " << fEdep/keV
457  << "; z(mm)= " << z/mm
458  << "; step(mm)= " << step->GetStepLength()/mm
459  << " by " << pd->GetParticleName()
460  << " E(MeV)= " << track->GetKineticEnergy()/MeV
461  << G4endl;
462  }
463  }
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467 
468 void HistoManager::AddLeakingParticle(const G4Track* track)
469 {
470  const G4ParticleDefinition* pd = track->GetDefinition();
471  const G4StepPoint* sp = track->GetStep()->GetPreStepPoint();
472  G4double e = std::log10(sp->GetKineticEnergy()/MeV);
473 
474  G4ThreeVector pos = sp->GetPosition();
475  G4ThreeVector dir = sp->GetMomentumDirection();
476  G4double x = pos.x();
477  G4double y = pos.y();
478  G4double z = pos.z();
479 
480  G4bool isLeaking = false;
481 
482  // Forward
483  if(z > -fAbsZ0 && dir.z() > 0.0) {
484  isLeaking = true;
485  if(pd == fNeutron) {
486  ++fNneu_forw;
487  fHisto->Fill(15,e,1.0);
488  } else isLeaking = true;
489 
490  // Backward
491  } else if (z < fAbsZ0 && dir.z() < 0.0) {
492  isLeaking = true;
493  if(pd == fNeutron) {
494  ++fNneu_back;
495  fHisto->Fill(16,e,1.0);
496  } else isLeaking = true;
497 
498  // Side
499  } else if (std::abs(z) <= -fAbsZ0 && x*dir.x() + y*dir.y() > 0.0) {
500  isLeaking = true;
501  if(pd == fNeutron) {
502  ++fNneu_leak;
503  fHisto->Fill(14,e,1.0);
504  } else isLeaking = true;
505  }
506 
507  // protons and pions
508  if(isLeaking) {
509  if(pd == G4Proton::Proton()) {
510  fHisto->Fill(17,e,1.0);
511  ++fNprot_leak;
512  } else if (pd == G4PionPlus::PionPlus() ||
513  pd == G4PionMinus::PionMinus()) {
514  fHisto->Fill(18,e,1.0);
515  ++fNpiofNleak;
516  }
517  }
518 }
519 
520 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
521 
522 void HistoManager::SetVerbose(G4int val)
523 {
524  fVerbose = val;
525  fHisto->SetVerbose(val);
526 }
527 
528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
529 
530 void HistoManager::Fill(G4int id, G4double x, G4double w)
531 {
532  fHisto->Fill(id, x, w);
533 }
534 
535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
536