| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file hadronic/Hadr00/src/HistoManager.cc << 26 // $Id: HistoManager.cc,v 1.4 2008/11/20 11:59:06 vnivanch Exp $
27 /// \brief Implementation of the HistoManager << 27 // GEANT4 tag $Name: geant4-09-02 $
28 // <<
29 // 28 //
30 //-------------------------------------------- 29 //---------------------------------------------------------------------------
31 // 30 //
32 // ClassName: HistoManager 31 // ClassName: HistoManager
33 // 32 //
34 // 33 //
35 // Author: V.Ivanchenko 30/01/01 34 // Author: V.Ivanchenko 30/01/01
36 // 35 //
37 // Modified: 36 // Modified:
38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchen 37 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko)
39 // 38 //
40 //-------------------------------------------- 39 //----------------------------------------------------------------------------
41 // 40 //
42 41
43 //....oooOO0OOooo........oooOO0OOooo........oo << 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oo << 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 44
46 #include "HistoManager.hh" 45 #include "HistoManager.hh"
47 << 46 #include "G4UnitsTable.hh"
48 #include "HistoManagerMessenger.hh" <<
49 <<
50 #include "G4HadronicProcessStore.hh" <<
51 #include "G4Neutron.hh" 47 #include "G4Neutron.hh"
52 #include "G4NistManager.hh" << 48 #include "globals.hh"
53 #include "G4NucleiProperties.hh" << 49 #include "G4ios.hh"
54 #include "G4ParticleDefinition.hh" 50 #include "G4ParticleDefinition.hh"
55 #include "G4ParticleTable.hh" 51 #include "G4ParticleTable.hh"
>> 52 #include "G4NistManager.hh"
>> 53 #include "G4HadronicProcessStore.hh"
>> 54
>> 55 #include "G4NucleiProperties.hh"
>> 56 #include "G4NistManager.hh"
56 #include "G4StableIsotopes.hh" 57 #include "G4StableIsotopes.hh"
57 #include "G4SystemOfUnits.hh" <<
58 #include "G4UnitsTable.hh" <<
59 #include "G4ios.hh" <<
60 #include "globals.hh" <<
61 58
62 //....oooOO0OOooo........oooOO0OOooo........oo << 59 #include "Histo.hh"
63 60
64 HistoManager::HistoManager() << 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 { <<
66 fAnalysisManager = 0; <<
67 fHistoName = "hadr00"; <<
68 62
69 fNeutron = G4Neutron::Neutron(); << 63 HistoManager* HistoManager::fManager = 0;
70 fMessenger = new HistoManagerMessenger(this) <<
71 fVerbose = 1; <<
72 64
73 fParticleName = "proton"; << 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 fElementName = "Al"; <<
75 66
76 fTargetMaterial = 0; << 67 HistoManager* HistoManager::GetPointer()
>> 68 {
>> 69 if(!fManager) {
>> 70 static HistoManager manager;
>> 71 fManager = &manager;
>> 72 }
>> 73 return fManager;
>> 74 }
>> 75
>> 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 77
>> 78 HistoManager::HistoManager()
>> 79 {
>> 80 histo = new Histo(verbose);
>> 81 neutron = G4Neutron::Neutron();
>> 82 verbose = 1;
>> 83
>> 84 particleName = "proton";
>> 85 elementName = "Al";
>> 86
>> 87 minKinEnergy = 0.1*MeV;
>> 88 maxKinEnergy = 1*TeV;
>> 89 minMomentum = 1*MeV;
>> 90 maxMomentum = 100*TeV;
77 91
78 fMinKinEnergy = 0.1 * MeV; << 92 nBinsE = 700;
79 fMaxKinEnergy = 10 * TeV; << 93 nBinsP = 800;
80 fMinMomentum = 1 * MeV; <<
81 fMaxMomentum = 10 * TeV; <<
82 94
83 fBinsE = 800; << 95 needsReset = false;
84 fBinsP = 700; << 96 isInitialised = false;
85 } 97 }
86 98
87 //....oooOO0OOooo........oooOO0OOooo........oo << 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 100
89 HistoManager::~HistoManager() 101 HistoManager::~HistoManager()
90 { 102 {
91 delete fMessenger; << 103 delete histo;
92 } 104 }
93 105
94 //....oooOO0OOooo........oooOO0OOooo........oo << 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95 107
96 void HistoManager::BeginOfRun() 108 void HistoManager::BeginOfRun()
97 { 109 {
98 G4double p1 = std::log10(fMinMomentum / GeV) << 110 needsReset = false;
99 G4double p2 = std::log10(fMaxMomentum / GeV) << 111
100 G4double e1 = std::log10(fMinKinEnergy / MeV << 112 G4double p1 = std::log10(minMomentum/GeV);
101 G4double e2 = std::log10(fMaxKinEnergy / MeV << 113 G4double p2 = std::log10(maxMomentum/GeV);
102 << 114 G4double e1 = std::log10(minKinEnergy/MeV);
103 // G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "< << 115 G4double e2 = std::log10(maxKinEnergy/MeV);
104 << 116
105 fAnalysisManager = G4AnalysisManager::Instan << 117 //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
106 fAnalysisManager->OpenFile(fHistoName + ".ro << 118
107 fAnalysisManager->SetFirstHistoId(1); << 119 if(isInitialised) {
108 << 120 histo->setHisto1D(0,nBinsP,p1,p2);
109 fAnalysisManager->CreateH1("h1", "Elastic cr << 121 histo->setHisto1D(1,nBinsE,e1,e2);
110 fBinsP, p1, p2); << 122 histo->setHisto1D(2,nBinsP,p1,p2);
111 fAnalysisManager->CreateH1("h2", "Elastic cr << 123 histo->setHisto1D(3,nBinsE,e1,e2);
112 fBinsE, e1, e2); << 124 histo->setHisto1D(4,nBinsE,e1,e2);
113 fAnalysisManager->CreateH1("h3", "Inelastic << 125 histo->setHisto1D(5,nBinsE,e1,e2);
114 fBinsP, p1, p2); << 126 histo->setHisto1D(6,nBinsE,e1,e2);
115 fAnalysisManager->CreateH1("h4", "Inelastic << 127
116 fBinsE, e1, e2); << 128 } else {
117 fAnalysisManager->CreateH1("h5", "Capture cr << 129 histo->add1D("h1","Elastic cross section (barn) as a functions of log10(p/GeV)",
118 fBinsE, e1, e2); << 130 nBinsP,p1,p2);
119 fAnalysisManager->CreateH1("h6", "Fission cr << 131 histo->add1D("h2","Elastic cross section (barn) as a functions of log10(E/MeV)",
120 fBinsE, e1, e2); << 132 nBinsE,e1,e2);
121 fAnalysisManager->CreateH1( << 133 histo->add1D("h3","Inelastic cross section (barn) as a functions of log10(p/GeV)",
122 "h7", "Charge exchange cross section (barn << 134 nBinsP,p1,p2);
123 fAnalysisManager->CreateH1("h8", "Total cros << 135 histo->add1D("h4","Inelastic cross section (barn) as a functions of log10(E/MeV)",
124 fBinsE, e1, e2); << 136 nBinsE,e1,e2);
125 fAnalysisManager->CreateH1( << 137 histo->add1D("h5","Capture cross section (barn) as a functions of log10(E/MeV)",
126 "h9", "Inelastic cross section per volume << 138 nBinsE,e1,e2);
127 fAnalysisManager->CreateH1( << 139 histo->add1D("h6","Fission cross section (barn) as a functions of log10(E/MeV)",
128 "h10", "Elastic cross section per volume a << 140 nBinsE,e1,e2);
>> 141 histo->add1D("h7","Charge exchange cross section (barn) as a functions of log10(E/MeV)",
>> 142 nBinsE,e1,e2);
>> 143 }
>> 144
>> 145 isInitialised = true;
>> 146 histo->book();
129 } 147 }
130 148
131 //....oooOO0OOooo........oooOO0OOooo........oo << 149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
132 150
133 void HistoManager::EndOfRun() 151 void HistoManager::EndOfRun()
134 { 152 {
135 if (fVerbose > 0) { << 153 if(verbose > 0) {
136 G4cout << "HistoManager: End of run action 154 G4cout << "HistoManager: End of run actions are started" << G4endl;
137 } 155 }
>> 156 if(needsReset) BeginOfRun();
138 157
139 const G4Element* elm = G4NistManager::Instan << 158 const G4Element* elm =
140 const G4Material* mat = G4NistManager::Insta << 159 G4NistManager::Instance()->FindOrBuildElement(elementName);
141 const G4ParticleDefinition* particle = << 160 const G4ParticleDefinition* particle =
142 G4ParticleTable::GetParticleTable()->FindP << 161 G4ParticleTable::GetParticleTable()->FindParticle(particleName);
143 << 162
144 G4cout << "### Fill Cross Sections for " << << 163 G4cout << "### Fill Cross Sections for " << particleName
145 if (fVerbose > 0) { << 164 << " off " << elementName
146 G4cout << "------------------------------- << 165 << G4endl;
147 G4cout << " N E(MeV) Elastic(b) << 166 if(verbose > 0) {
148 if (particle == fNeutron) { << 167 G4cout << "-------------------------------------------------------------"
149 G4cout << " Capture(b) Fission(b)"; << 168 << G4endl;
150 } << 169 G4cout << " N E(MeV) Elastic(barn) Inelastic(barn)";
151 G4cout << " Total(b)" << G4endl; << 170 if(particle == neutron) G4cout << " Capture(barn) Fission(barn)";
152 G4cout << "------------------------------- << 171 G4cout << " ChargeExchange(barn)" << G4endl;
>> 172 G4cout << "-------------------------------------------------------------"
>> 173 << G4endl;
153 } 174 }
154 if (!particle || !elm) { << 175 if(!particle || !elm) {
155 G4cout << "HistoManager WARNING Particle o 176 G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
156 return; 177 return;
157 } 178 }
158 179
159 G4int prec = G4cout.precision(); 180 G4int prec = G4cout.precision();
160 G4cout.precision(4); << 181 G4cout.precision(7);
161 182
162 G4HadronicProcessStore* store = G4HadronicPr 183 G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
163 G4double mass = particle->GetPDGMass(); 184 G4double mass = particle->GetPDGMass();
164 185
165 // Build histograms 186 // Build histograms
166 187
167 G4double p1 = std::log10(fMinMomentum / GeV) << 188 G4double p1 = std::log10(minMomentum/GeV);
168 G4double p2 = std::log10(fMaxMomentum / GeV) << 189 G4double p2 = std::log10(maxMomentum/GeV);
169 G4double e1 = std::log10(fMinKinEnergy / MeV << 190 G4double e1 = std::log10(minKinEnergy/MeV);
170 G4double e2 = std::log10(fMaxKinEnergy / MeV << 191 G4double e2 = std::log10(maxKinEnergy/MeV);
171 G4double de = (e2 - e1) / G4double(fBinsE); << 192 G4double de = (e2 - e1)/G4double(nBinsE);
172 G4double dp = (p2 - p1) / G4double(fBinsP); << 193 G4double dp = (p2 - p1)/G4double(nBinsP);
173 194
174 G4double x = e1 - de * 0.5; << 195 G4double x = e1 - de*0.5;
175 G4double e, p, xs, xtot; << 196 G4double e, p, xs;
176 G4int i; 197 G4int i;
177 << 198 for(i=0; i<nBinsE; i++) {
178 G4double coeff = 1.0; <<
179 if (fTargetMaterial) { <<
180 coeff = fTargetMaterial->GetDensity() * cm <<
181 } <<
182 <<
183 for (i = 0; i < fBinsE; i++) { <<
184 x += de; 199 x += de;
185 e = std::pow(10., x) * MeV; << 200 e = std::pow(10.,x)*MeV;
186 if (fVerbose > 0) G4cout << std::setw(5) < << 201 if(verbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
187 xs = store->GetElasticCrossSectionPerAtom( << 202 xs = store->GetElasticCrossSectionPerAtom(particle,e,elm);
188 xtot = xs; << 203 if(verbose>0) G4cout << std::setw(14) << xs/barn;
189 if (fVerbose > 0) G4cout << std::setw(12) << 204 histo->fill(1, x, xs/barn);
190 fAnalysisManager->FillH1(2, x, xs / barn); << 205 xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm);
191 xs = store->GetInelasticCrossSectionPerAto << 206 if(verbose>0) G4cout << " " << std::setw(17) << xs/barn;
192 xtot += xs; << 207 histo->fill(3, x, xs/barn);
193 if (fVerbose > 0) G4cout << " " << std::se << 208 if(particle == neutron) {
194 fAnalysisManager->FillH1(4, x, xs / barn); << 209 xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm);
195 if (fTargetMaterial) { << 210 if(verbose>0) G4cout << " " << std::setw(17) << xs/barn;
196 xs = store->GetInelasticCrossSectionPerV << 211 histo->fill(4, x, xs/barn);
197 fAnalysisManager->FillH1(9, x, xs / coef << 212 xs = store->GetFissionCrossSectionPerAtom(particle,e,elm);
198 xs = store->GetElasticCrossSectionPerVol << 213 if(verbose>0) G4cout << " " << std::setw(17) << xs/barn;
199 fAnalysisManager->FillH1(10, x, xs / coe << 214 histo->fill(5, x, xs/barn);
200 } <<
201 if (particle == fNeutron) { <<
202 xs = store->GetCaptureCrossSectionPerAto <<
203 xtot += xs; <<
204 if (fVerbose > 0) G4cout << " " << std:: <<
205 fAnalysisManager->FillH1(5, x, xs / barn <<
206 xs = store->GetFissionCrossSectionPerAto <<
207 xtot += xs; <<
208 if (fVerbose > 0) G4cout << " " << std:: <<
209 fAnalysisManager->FillH1(6, x, xs / barn <<
210 } 215 }
211 xs = store->GetChargeExchangeCrossSectionP << 216 xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm);
212 if (fVerbose > 0) G4cout << " " << std::se << 217 if(verbose>0) G4cout << " " << std::setw(17) << xs/barn;
213 fAnalysisManager->FillH1(7, x, xs / barn); << 218 histo->fill(6, x, xs/barn);
214 fAnalysisManager->FillH1(8, x, xtot / barn << 219 if(verbose>0) G4cout << G4endl;
215 } 220 }
216 221
217 x = p1 - dp * 0.5; << 222 x = p1 - dp*0.5;
218 for (i = 0; i < fBinsP; i++) { << 223 for(i=0; i<nBinsP; i++) {
219 x += dp; 224 x += dp;
220 p = std::pow(10., x) * GeV; << 225 p = std::pow(10.,x)*GeV;
221 e = std::sqrt(p * p + mass * mass) - mass; << 226 e = std::sqrt(p*p + mass*mass) - mass;
222 xs = store->GetElasticCrossSectionPerAtom( << 227 xs = store->GetElasticCrossSectionPerAtom(particle,e,elm);
223 fAnalysisManager->FillH1(1, x, xs / barn); << 228 histo->fill(0, x, xs/barn);
224 xs = store->GetInelasticCrossSectionPerAto << 229 xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm);
225 fAnalysisManager->FillH1(3, x, xs / barn); << 230 histo->fill(2, x, xs/barn);
226 } << 231 }
227 if (fVerbose > 0) { << 232 if(verbose > 0) {
228 G4cout << "------------------------------- << 233 G4cout << "---------------------------------------------------------"
>> 234 << G4endl;
229 } 235 }
230 G4cout.precision(prec); 236 G4cout.precision(prec);
231 fAnalysisManager->Write(); << 237 if(verbose > 1) histo->print(0);
232 fAnalysisManager->CloseFile(); << 238 histo->save();
233 fAnalysisManager->Clear(); << 239 histo->reset();
234 <<
235 G4bool extra = true; <<
236 if (fTargetMaterial && extra) { <<
237 G4double E = 5 * GeV; <<
238 G4double cross = store->GetInelasticCrossS <<
239 if (cross <= 0.0) { <<
240 cross = 1.e-100; <<
241 } <<
242 G4cout << "### " << particle->GetParticleN <<
243 << fTargetMaterial->GetName() <<
244 << " xs/X0= " << 1.0 / (cross * fTa <<
245 } <<
246 } 240 }
247 241
248 //....oooOO0OOooo........oooOO0OOooo........oo << 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249 243
250 void HistoManager::SetVerbose(G4int val) << 244 void HistoManager::SetVerbose(G4int val)
251 { 245 {
252 fVerbose = val; << 246 verbose = val;
>> 247 histo->setVerbose(val);
253 } 248 }
254 249
255 //....oooOO0OOooo........oooOO0OOooo........oo << 250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 251
256 252