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 // Modified:
38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko)
39 // 16.11.2006 Add beamFlag (V.Ivanchenko)
40 // 16.10.2012 Renamed G4IonFTFPBinaryCascadePhysics as G4IonPhysics (A.Ribon)
41 //
42 //----------------------------------------------------------------------------
43 //
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 #include "HistoManager.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 "G4VModularPhysicsList.hh"
72 #include "G4VPhysicsConstructor.hh"
73 #include "G4RunManager.hh"
74 #include "DetectorConstruction.hh"
75 #include "G4NistManager.hh"
76 #include "G4Material.hh"
77 #include "G4BuilderType.hh"
78 #include "G4SystemOfUnits.hh"
79 
80 #include "IonUrQMDPhysics.hh"
81 #include "IonHIJINGPhysics.hh"
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
85 HistoManager* HistoManager::fManager = 0;
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
89 HistoManager* HistoManager::GetPointer()
90 {
91  if(!fManager) {
92  static HistoManager manager;
93  fManager = &manager;
94  }
95  return fManager;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
100 HistoManager::HistoManager()
101 {
102  fVerbose= 0;
103  fNSlices = 100;
104  fNBinsE = 100;
105  fNHisto = 20;
106  fLength = 300.*mm;
107  fEdepMax = 1.0*GeV;
108 
109  fPrimaryDef= 0;
110  fPrimaryKineticEnergy = 0.0;
111  fMaterial = 0;
112  fBeamFlag = true;
113 
114  fHisto = new Histo();
115  fHisto->SetVerbose(fVerbose);
116  fNeutron = G4Neutron::Neutron();
117  fPhysList = 0;
118  fIonPhysics= 0;
119  Initialise();
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
124 void HistoManager::Initialise()
125 {
126  fAbsZ0 = -0.5*fLength;
127  fNevt = 0;
128  fNelec = 0;
129  fNposit = 0;
130  fNgam = 0;
131  fNstep = 0;
132  fNions = 0;
133  fNdeut = 0;
134  fNalpha = 0;
135  fNkaons = 0;
136  fNmuons = 0;
137  fNcpions = 0;
138  fNpi0 = 0;
139  fNneutron = 0;
140  fNproton = 0;
141  fNaproton = 0;
142 
143  fEdepEvt = 0.0;
144  fEdepSum = 0.0;
145  fEdepSum2 = 0.0;
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 
150 HistoManager::~HistoManager()
151 {
152  delete fHisto;
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
157 void HistoManager::BookHisto()
158 {
159  fHisto->Add1D("0","Energy deposition (MeV/mm/event) in the target",
160  fNSlices,0.0,fLength/mm,MeV/mm);
161  fHisto->Add1D("1","Log10 Energy (GeV) of gammas",fNBinsE,-5.,5.,1.0);
162  fHisto->Add1D("2","Log10 Energy (GeV) of electrons",fNBinsE,-5.,5.,1.0);
163  fHisto->Add1D("3","Log10 Energy (GeV) of positrons",fNBinsE,-5.,5.,1.0);
164  fHisto->Add1D("4","Log10 Energy (GeV) of protons",fNBinsE,-5.,5.,1.0);
165  fHisto->Add1D("5","Log10 Energy (GeV) of neutrons",fNBinsE,-5.,5.,1.0);
166  fHisto->Add1D("6","Log10 Energy (GeV) of charged pions",fNBinsE,-4.,6.,1.0);
167  fHisto->Add1D("7","Log10 Energy (GeV) of pi0",fNBinsE,-4.,6.,1.0);
168  fHisto->Add1D("8","Log10 Energy (GeV) of charged kaons",fNBinsE,-4.,6.,1.0);
169  fHisto->Add1D("9","Log10 Energy (GeV) of neutral kaons",fNBinsE,-4.,6.,1.0);
170  fHisto->Add1D("10","Log10 Energy (GeV) of deuterons and tritons",
171  fNBinsE,-5.,5.,1.0);
172  fHisto->Add1D("11","Log10 Energy (GeV) of He3 and alpha",fNBinsE,-5.,5.,1.0);
173  fHisto->Add1D("12","Log10 Energy (GeV) of Generic Ions",fNBinsE,-5.,5.,1.0);
174  fHisto->Add1D("13","Log10 Energy (GeV) of muons",fNBinsE,-4.,6.,1.0);
175  fHisto->Add1D("14","Log10 Energy (GeV) of pi+",fNBinsE,-4.,6.,1.0);
176  fHisto->Add1D("15","Log10 Energy (GeV) of pi-",fNBinsE,-4.,6.,1.0);
177  fHisto->Add1D("16","X Section (mb) of Secondary Fragments Z with E>1 GeV (mb)"
178  ,25,0.5,25.5,1.0);
179  fHisto->Add1D("17","Secondary Fragment A E>1 GeV",50,0.5,50.5,1.0);
180  fHisto->Add1D("18","Secondary Fragment Z E<1 GeV",25,0.5,25.5,1.0);
181  fHisto->Add1D("19","Secondary Fragment A E<1 GeV",50,0.5,50.5,1.0);
182  fHisto->Add1D("20","X Section (mb) of Secondary Fragments Z (mb) ",
183  25,0.5,25.5,1.0);
184  fHisto->Add1D("21","Secondary Fragment A ",50,0.5,50.5,1.0);
185 }
186 
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
188 
189 void HistoManager::BeginOfRun()
190 {
191  Initialise();
192  BookHisto();
193  fHisto->Book();
194 
195  if(fVerbose > 0) {
196  G4cout << "HistoManager: Histograms are booked and run has been started"
197  <<G4endl<<" BeginOfRun (After fHisto->book)"<< G4endl;
198  }
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 
203 void HistoManager::EndOfRun()
204 {
205 
206  G4cout << "HistoManager: End of run actions are started" << G4endl;
207 
208  // Average values
209  G4cout<<"========================================================"<<G4endl;
210 
211  G4double x = (G4double)fNevt;
212  if(fNevt > 0) { x = 1.0/x; }
213 
214  G4double xe = x*fNelec;
215  G4double xg = x*fNgam;
216  G4double xp = x*fNposit;
217  G4double xs = x*fNstep;
218  G4double xn = x*fNneutron;
219  G4double xpn = x*fNproton;
220  G4double xap = x*fNaproton;
221  G4double xpc = x*fNcpions;
222  G4double xp0 = x*fNpi0;
223  G4double xpk = x*fNkaons;
224  G4double xpm = x*fNmuons;
225  G4double xid = x*fNdeut;
226  G4double xia = x*fNalpha;
227  G4double xio = x*fNions;
228 
229  fEdepSum *= x;
230  fEdepSum2 *= x;
231  fEdepSum2 -= fEdepSum*fEdepSum;
232  if(fEdepSum2 > 0.0) fEdepSum2 = std::sqrt(fEdepSum2);
233  else fEdepSum2 = 0.0;
234 
235  G4cout << "Beam particle "
236  << fPrimaryDef->GetParticleName() <<G4endl;
237  G4cout << "Beam Energy(GeV) "
238  << fPrimaryKineticEnergy/GeV <<G4endl;
239  G4cout << "Number of events "
240  << fNevt <<G4endl;
241  G4cout << std::setprecision(4) << "Average energy deposit (GeV) "
242  << fEdepSum/GeV
243  << " RMS(GeV) " << fEdepSum2/GeV << G4endl;
244  G4cout << std::setprecision(4) << "Average number of steps "
245  << xs << G4endl;
246  G4cout << std::setprecision(4) << "Average number of gamma "
247  << xg << G4endl;
248  G4cout << std::setprecision(4) << "Average number of e- "
249  << xe << G4endl;
250  G4cout << std::setprecision(4) << "Average number of e+ "
251  << xp << G4endl;
252  G4cout << std::setprecision(4) << "Average number of neutrons "
253  << xn << G4endl;
254  G4cout << std::setprecision(4) << "Average number of protons "
255  << xpn << G4endl;
256  G4cout << std::setprecision(4) << "Average number of antiprotons "
257  << xap << G4endl;
258  G4cout << std::setprecision(4) << "Average number of pi+ & pi- "
259  << xpc << G4endl;
260  G4cout << std::setprecision(4) << "Average number of pi0 "
261  << xp0 << G4endl;
262  G4cout << std::setprecision(4) << "Average number of kaons "
263  << xpk << G4endl;
264  G4cout << std::setprecision(4) << "Average number of muons "
265  << xpm << G4endl;
266  G4cout << std::setprecision(4) << "Average number of deuterons+tritons "
267  << xid << G4endl;
268  G4cout << std::setprecision(4) << "Average number of He3+alpha "
269  << xia << G4endl;
270  G4cout << std::setprecision(4) << "Average number of ions "
271  << xio << G4endl;
272  G4cout<<"========================================================"<<G4endl;
273  G4cout<<G4endl;
274 
275  // normalise histograms
276  for(G4int i=0; i<fNHisto; i++) {
277  fHisto->ScaleH1(i,x);
278  }
279  // will work only for pure material - 1 element
280  // G4cout << fMaterial << G4endl;
281  G4double F= 1000/(fLength*barn*fMaterial->GetTotNbOfAtomsPerVolume());
282  if(F > 0.0) {
283  fHisto->ScaleH1(16, F);
284  fHisto->ScaleH1(20, F);
285  }
286 
287  fHisto->Save();
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 
292 void HistoManager::BeginOfEvent()
293 {
294  ++fNevt;
295  fEdepEvt = 0.0;
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 
300 void HistoManager::EndOfEvent()
301 {
302  fEdepSum += fEdepEvt;
303  fEdepSum2 += fEdepEvt*fEdepEvt;
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 
308 void HistoManager::ScoreNewTrack(const G4Track* track)
309 {
310  const G4ParticleDefinition* pd = track->GetDefinition();
311  G4String name = pd->GetParticleName();
312  G4double e = track->GetKineticEnergy();
313 
314  // Primary track
315  if(0 == track->GetParentID()) {
316 
317  fPrimaryKineticEnergy = e;
318  fPrimaryDef = pd;
319  G4ThreeVector dir = track->GetMomentumDirection();
320  if(1 < fVerbose)
321  G4cout << "### Primary " << name
322  << " kinE(GeV)= " << e/GeV
323  << "; m(GeV)= " << pd->GetPDGMass()/GeV
324  << "; pos(mm)= " << track->GetPosition()/mm
325  << "; dir= " << track->GetMomentumDirection()
326  << G4endl;
327 
328  // Secondary track
329  } else {
330  if(1 < fVerbose) {
331  G4cout << "=== Secondary " << name
332  << " kinE(GeV)= " << e/GeV
333  << "; m(GeV)= " << pd->GetPDGMass()/GeV
334  << "; pos(mm)= " << track->GetPosition()/mm
335  << "; dir= " << track->GetMomentumDirection()
336  << G4endl;
337  }
338  e = std::log10(e/GeV);
339  if(pd == G4Gamma::Gamma()) {
340  fNgam++;
341  fHisto->Fill(1,e,1.0);
342  } else if ( pd == G4Electron::Electron()) {
343  fNelec++;
344  fHisto->Fill(2,e,1.0);
345  } else if ( pd == G4Positron::Positron()) {
346  fNposit++;
347  fHisto->Fill(3,e,1.0);
348  } else if ( pd == G4Proton::Proton()) {
349  fNproton++;
350  fHisto->Fill(4,e,1.0);
351  } else if ( pd == fNeutron) {
352  fNneutron++;
353  fHisto->Fill(5,e,1.0);
354  } else if ( pd == G4AntiProton::AntiProton()) {
355  fNaproton++;
356  } else if ( pd == G4PionPlus::PionPlus() ) {
357  fNcpions++;
358  fHisto->Fill(6,e,1.0);
359  fHisto->Fill(14,e,1.0);
360 
361  } else if ( pd == G4PionMinus::PionMinus()) {
362  fNcpions++;
363  fHisto->Fill(6,e,1.0);
364  fHisto->Fill(15,e,1.0);
365 
366  } else if ( pd == G4PionZero::PionZero()) {
367  fNpi0++;
368  fHisto->Fill(7,e,1.0);
369  } else if ( pd == G4KaonPlus::KaonPlus() ||
370  pd == G4KaonMinus::KaonMinus()) {
371  fNkaons++;
372  fHisto->Fill(8,e,1.0);
373  } else if ( pd == G4KaonZeroShort::KaonZeroShort() ||
374  pd == G4KaonZeroLong::KaonZeroLong()) {
375  fNkaons++;
376  fHisto->Fill(9,e,1.0);
377  } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
378  fNdeut++;
379  fHisto->Fill(10,e,1.0);
380  } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
381  fNalpha++;
382  fHisto->Fill(11,e,1.0);
383  } else if ( pd->GetParticleType() == "nucleus") {
384  fNions++;
385  fHisto->Fill(12,e,1.0);
386  G4double Z = pd->GetPDGCharge()/eplus;
387  G4double A = (G4double)pd->GetBaryonNumber();
388  if(e > 0.0) {
389  fHisto->Fill(16,Z,1.0);
390  fHisto->Fill(17,A,1.0);
391  } else {
392  fHisto->Fill(18,Z,1.0);
393  fHisto->Fill(19,A,1.0);
394  }
395  fHisto->Fill(20,Z,1.0);
396  fHisto->Fill(21,A,1.0);
397  } else if ( pd == G4MuonPlus::MuonPlus() ||
398  pd == G4MuonMinus::MuonMinus()) {
399  fNmuons++;
400  fHisto->Fill(13,e,1.0);
401  }
402  }
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
406 
407 void HistoManager::AddTargetStep(const G4Step* step)
408 {
409  ++fNstep;
410  G4double edep = step->GetTotalEnergyDeposit();
411 
412  if(edep > 0.0) {
413  G4ThreeVector pos =
414  (step->GetPreStepPoint()->GetPosition() +
415  step->GetPostStepPoint()->GetPosition())*0.5;
416 
417  G4double z = pos.z() - fAbsZ0;
418 
419  // scoring
420  fEdepEvt += edep;
421  fHisto->Fill(0,z,edep);
422 
423  if(1 < fVerbose) {
424  const G4Track* track = step->GetTrack();
425  G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
426  << "; z(mm)= " << z/mm
427  << "; step(mm)= " << step->GetStepLength()/mm
428  << " by " << track->GetDefinition()->GetParticleName()
429  << " E(GeV)= " << track->GetKineticEnergy()/GeV
430  << G4endl;
431  }
432  }
433 }
434 
435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436 
437 void HistoManager::SetVerbose(G4int val)
438 {
439  fVerbose = val;
440  fHisto->SetVerbose(val);
441 }
442 
443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444 
445 void HistoManager::Fill(G4int id, G4double x, G4double w)
446 {
447  fHisto->Fill(id, x, w);
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
451 
452 void HistoManager::SetIonPhysics(const G4String& nam)
453 {
454  if(fIonPhysics) {
455  G4cout << "### HistoManager WARNING: Ion Physics is already defined: <"
456  << nam << "> is ignored!" << G4endl;
457  } else if(nam == "HIJING") {
458 #ifdef G4_USE_HIJING
459  fIonPhysics = new IonHIJINGPhysics();
460  fPhysList->ReplacePhysics(fIonPhysics);
461  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
462  G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added"
463  << G4endl;
464 #else
465  G4cout << "Error: Ion Physics HIJING is requested but is not available"
466  <<G4endl;
467 #endif
468  } else if(nam == "UrQMD") {
469 #ifdef G4_USE_URQMD
470  fIonPhysics = new IonUrQMDPhysics(1);
471  fPhysList->ReplacePhysics(fIonPhysics);
472  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
473  G4cout << "### SetIonPhysics: Ion Physics UrQMD is added"
474  << G4endl;
475 #else
476  G4cout << "Error: Ion Physics UrQMD is requested but is not available"
477  <<G4endl;
478 #endif
479  } else {
480  G4cout << "### HistoManager WARNING: Ion Physics <"
481  << nam << "> is unknown!" << G4endl;
482  }
483 }
484 
485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
486