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 //
40 //----------------------------------------------------------------------------
41 //
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
46 #include "HistoManager.hh"
47 #include "G4UnitsTable.hh"
48 #include "G4Neutron.hh"
49 #include "globals.hh"
50 #include "G4ios.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4NistManager.hh"
54 #include "G4HadronicProcessStore.hh"
55 
56 #include "G4NucleiProperties.hh"
57 #include "G4NistManager.hh"
58 #include "G4StableIsotopes.hh"
59 #include "G4SystemOfUnits.hh"
60 
61 #include "HistoManagerMessenger.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 HistoManager::HistoManager()
66 {
67  fAnalysisManager = 0;
68  fHistoName = "hadr00";
69 
70  fNeutron = G4Neutron::Neutron();
71  fMessenger = new HistoManagerMessenger(this);
72  fVerbose = 1;
73 
74  fParticleName = "proton";
75  fElementName = "Al";
76 
77  fTargetMaterial = 0;
78 
79  fMinKinEnergy = 0.1*MeV;
80  fMaxKinEnergy = 10*TeV;
81  fMinMomentum = 1*MeV;
82  fMaxMomentum = 10*TeV;
83 
84  fBinsE = 800;
85  fBinsP = 700;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
90 HistoManager::~HistoManager()
91 {
92  delete fMessenger;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
97 void HistoManager::BeginOfRun()
98 {
99  G4double p1 = std::log10(fMinMomentum/GeV);
100  G4double p2 = std::log10(fMaxMomentum/GeV);
101  G4double e1 = std::log10(fMinKinEnergy/MeV);
102  G4double e2 = std::log10(fMaxKinEnergy/MeV);
103 
104  //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
105 
106  fAnalysisManager = G4AnalysisManager::Instance();
107  fAnalysisManager->OpenFile(fHistoName+".root");
108  fAnalysisManager->SetFirstHistoId(1);
109 
110  fAnalysisManager->CreateH1("h1",
111  "Elastic cross section (barn) as a functions of log10(p/GeV)",
112  fBinsP,p1,p2);
113  fAnalysisManager->CreateH1("h2",
114  "Elastic cross section (barn) as a functions of log10(E/MeV)",
115  fBinsE,e1,e2);
116  fAnalysisManager->CreateH1("h3",
117  "Inelastic cross section (barn) as a functions of log10(p/GeV)",
118  fBinsP,p1,p2);
119  fAnalysisManager->CreateH1("h4",
120  "Inelastic cross section (barn) as a functions of log10(E/MeV)",
121  fBinsE,e1,e2);
122  fAnalysisManager->CreateH1("h5",
123  "Capture cross section (barn) as a functions of log10(E/MeV)",
124  fBinsE,e1,e2);
125  fAnalysisManager->CreateH1("h6",
126  "Fission cross section (barn) as a functions of log10(E/MeV)",
127  fBinsE,e1,e2);
128  fAnalysisManager->CreateH1("h7",
129  "Charge exchange cross section (barn) as a functions of log10(E/MeV)",
130  fBinsE,e1,e2);
131  fAnalysisManager->CreateH1("h8",
132  "Total cross section (barn) as a functions of log10(E/MeV)",
133  fBinsE,e1,e2);
134  fAnalysisManager->CreateH1("h9",
135  "Inelastic cross section per volume as a functions of log10(E/MeV)",
136  fBinsE,e1,e2);
137  fAnalysisManager->CreateH1("h10",
138  "Elastic cross section per volume as a functions of log10(E/MeV)",
139  fBinsE,e1,e2);
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 void HistoManager::EndOfRun()
145 {
146  if(fVerbose > 0) {
147  G4cout << "HistoManager: End of run actions are started" << G4endl;
148  }
149 
150  const G4Element* elm =
151  G4NistManager::Instance()->FindOrBuildElement(fElementName);
152  const G4Material* mat =
153  G4NistManager::Instance()->FindOrBuildMaterial("G4_"+fElementName);
154  const G4ParticleDefinition* particle =
155  G4ParticleTable::GetParticleTable()->FindParticle(fParticleName);
156 
157  G4cout << "### Fill Cross Sections for " << fParticleName
158  << " off " << fElementName
159  << G4endl;
160  if(fVerbose > 0) {
161  G4cout << "-------------------------------------------------------------"
162  << G4endl;
163  G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
164  if(particle == fNeutron) { G4cout << " Capture(b) Fission(b)"; }
165  G4cout << " Total(b)" << G4endl;
166  G4cout << "-------------------------------------------------------------"
167  << G4endl;
168  }
169  if(!particle || !elm) {
170  G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
171  return;
172  }
173 
174  G4int prec = G4cout.precision();
175  G4cout.precision(4);
176 
177  G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
178  G4double mass = particle->GetPDGMass();
179 
180  // Build histograms
181 
182  G4double p1 = std::log10(fMinMomentum/GeV);
183  G4double p2 = std::log10(fMaxMomentum/GeV);
184  G4double e1 = std::log10(fMinKinEnergy/MeV);
185  G4double e2 = std::log10(fMaxKinEnergy/MeV);
186  G4double de = (e2 - e1)/G4double(fBinsE);
187  G4double dp = (p2 - p1)/G4double(fBinsP);
188 
189  G4double x = e1 - de*0.5;
190  G4double e, p, xs, xtot;
191  G4int i;
192 
193  G4double coeff = 1.0;
194  if(fTargetMaterial) { coeff = fTargetMaterial->GetDensity()*cm2/g; }
195 
196  for(i=0; i<fBinsE; i++) {
197  x += de;
198  e = std::pow(10.,x)*MeV;
199  if(fVerbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
200  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
201  xtot = xs;
202  if(fVerbose>0) G4cout << std::setw(12) << xs/barn;
203  fAnalysisManager->FillH1(2, x, xs/barn);
204  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
205  xtot += xs;
206  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
207  fAnalysisManager->FillH1(4, x, xs/barn);
208  if(fTargetMaterial) {
209  xs =
210  store->GetInelasticCrossSectionPerVolume(particle,e,fTargetMaterial);
211  fAnalysisManager->FillH1(9, x, xs/coeff);
212  xs =
213  store->GetElasticCrossSectionPerVolume(particle,e,fTargetMaterial);
214  fAnalysisManager->FillH1(10, x, xs/coeff);
215  }
216  if(particle == fNeutron) {
217  xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm,mat);
218  xtot += xs;
219  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
220  fAnalysisManager->FillH1(5, x, xs/barn);
221  xs = store->GetFissionCrossSectionPerAtom(particle,e,elm,mat);
222  xtot += xs;
223  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
224  fAnalysisManager->FillH1(6, x, xs/barn);
225  }
226  xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm,mat);
227  if(fVerbose>0) G4cout << " " << std::setw(12) << xtot/barn << G4endl;
228  fAnalysisManager->FillH1(7, x, xs/barn);
229  fAnalysisManager->FillH1(8, x, xtot/barn);
230  }
231 
232  x = p1 - dp*0.5;
233  for(i=0; i<fBinsP; i++) {
234  x += dp;
235  p = std::pow(10.,x)*GeV;
236  e = std::sqrt(p*p + mass*mass) - mass;
237  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
238  fAnalysisManager->FillH1(1, x, xs/barn);
239  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
240  fAnalysisManager->FillH1(3, x, xs/barn);
241  }
242  if(fVerbose > 0) {
243  G4cout << "-------------------------------------------------------------"
244  << G4endl;
245  }
246  G4cout.precision(prec);
247  fAnalysisManager->Write();
248  fAnalysisManager->CloseFile();
249 
250  G4bool extra = true;
251  if(fTargetMaterial && extra) {
252  G4double E= 5*GeV;
253  G4double cross =
254  store->GetInelasticCrossSectionPerVolume(particle,E,fTargetMaterial);
255  if(cross <= 0.0) { cross = 1.e-100; }
256  G4cout << "### " << particle->GetParticleName() << " " << E/GeV
257  << " GeV on " << fTargetMaterial->GetName()
258  << " xs/X0= " << 1.0/(cross*fTargetMaterial->GetRadlen()) << G4endl;
259  }
260  delete fAnalysisManager;
261 }
262 
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264 
265 void HistoManager::SetVerbose(G4int val)
266 {
267  fVerbose = val;
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271