| 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
>> 27 /// \brief Implementation of the HistoManager class
26 // 28 //
>> 29 // $Id: HistoManager.cc 75733 2013-11-05 17:57:53Z vnivanch $
27 // 30 //
>> 31 //---------------------------------------------------------------------------
>> 32 //
>> 33 // ClassName: HistoManager
>> 34 //
>> 35 //
>> 36 // Author: V.Ivanchenko 30/01/01
>> 37 //
>> 38 // Modified:
>> 39 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko)
>> 40 //
>> 41 //----------------------------------------------------------------------------
>> 42 //
>> 43
28 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
29 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 46
31 #include "HistoManager.hh" 47 #include "HistoManager.hh"
32 <<
33 #include "DetectorConstruction.hh" <<
34 <<
35 #include "G4UnitsTable.hh" 48 #include "G4UnitsTable.hh"
>> 49 #include "G4Neutron.hh"
>> 50 #include "globals.hh"
>> 51 #include "G4ios.hh"
>> 52 #include "G4ParticleDefinition.hh"
>> 53 #include "G4ParticleTable.hh"
>> 54 #include "G4NistManager.hh"
>> 55 #include "G4HadronicProcessStore.hh"
>> 56
>> 57 #include "G4NucleiProperties.hh"
>> 58 #include "G4NistManager.hh"
>> 59 #include "G4StableIsotopes.hh"
>> 60 #include "G4SystemOfUnits.hh"
>> 61
>> 62 #include "HistoManagerMessenger.hh"
36 63
37 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
38 65
39 HistoManager::HistoManager() 66 HistoManager::HistoManager()
40 { 67 {
41 Book(); << 68 fAnalysisManager = 0;
>> 69 fHistoName = "hadr00";
>> 70
>> 71 fNeutron = G4Neutron::Neutron();
>> 72 fMessenger = new HistoManagerMessenger(this);
>> 73 fVerbose = 1;
>> 74
>> 75 fParticleName = "proton";
>> 76 fElementName = "Al";
>> 77
>> 78 fMinKinEnergy = 0.1*MeV;
>> 79 fMaxKinEnergy = 10*TeV;
>> 80 fMinMomentum = 1*MeV;
>> 81 fMaxMomentum = 10*TeV;
>> 82
>> 83 fBinsE = 800;
>> 84 fBinsP = 700;
>> 85 }
>> 86
>> 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 88
>> 89 HistoManager::~HistoManager()
>> 90 {
>> 91 delete fMessenger;
>> 92 }
>> 93
>> 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 95
>> 96 void HistoManager::BeginOfRun()
>> 97 {
>> 98 G4double p1 = std::log10(fMinMomentum/GeV);
>> 99 G4double p2 = std::log10(fMaxMomentum/GeV);
>> 100 G4double e1 = std::log10(fMinKinEnergy/MeV);
>> 101 G4double e2 = std::log10(fMaxKinEnergy/MeV);
>> 102
>> 103 //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
>> 104
>> 105 fAnalysisManager = G4AnalysisManager::Instance();
>> 106 fAnalysisManager->OpenFile(fHistoName+".root");
>> 107 fAnalysisManager->SetFirstHistoId(1);
>> 108
>> 109 fAnalysisManager->CreateH1("h1",
>> 110 "Elastic cross section (barn) as a functions of log10(p/GeV)",
>> 111 fBinsP,p1,p2);
>> 112 fAnalysisManager->CreateH1("h2",
>> 113 "Elastic cross section (barn) as a functions of log10(E/MeV)",
>> 114 fBinsE,e1,e2);
>> 115 fAnalysisManager->CreateH1("h3",
>> 116 "Inelastic cross section (barn) as a functions of log10(p/GeV)",
>> 117 fBinsP,p1,p2);
>> 118 fAnalysisManager->CreateH1("h4",
>> 119 "Inelastic cross section (barn) as a functions of log10(E/MeV)",
>> 120 fBinsE,e1,e2);
>> 121 fAnalysisManager->CreateH1("h5",
>> 122 "Capture cross section (barn) as a functions of log10(E/MeV)",
>> 123 fBinsE,e1,e2);
>> 124 fAnalysisManager->CreateH1("h6",
>> 125 "Fission cross section (barn) as a functions of log10(E/MeV)",
>> 126 fBinsE,e1,e2);
>> 127 fAnalysisManager->CreateH1("h7",
>> 128 "Charge exchange cross section (barn) as a functions of log10(E/MeV)",
>> 129 fBinsE,e1,e2);
>> 130 fAnalysisManager->CreateH1("h8",
>> 131 "Total cross section (barn) as a functions of log10(E/MeV)",
>> 132 fBinsE,e1,e2);
42 } 133 }
43 134
44 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 136
46 void HistoManager::Book() << 137 void HistoManager::EndOfRun()
47 { 138 {
48 // Create or get analysis manager << 139 if(fVerbose > 0) {
49 // The choice of analysis technology is done << 140 G4cout << "HistoManager: End of run actions are started" << G4endl;
50 // in HistoManager.hh << 141 }
51 G4AnalysisManager* analysisManager = G4Analy << 142
52 analysisManager->SetDefaultFileType("root"); << 143 const G4Element* elm =
53 analysisManager->SetFileName(fFileName); << 144 G4NistManager::Instance()->FindOrBuildElement(fElementName);
54 analysisManager->SetVerboseLevel(1); << 145 const G4Material* mat =
55 analysisManager->SetActivation(true); // en << 146 G4NistManager::Instance()->FindOrBuildMaterial("G4_"+fElementName);
56 << 147 const G4ParticleDefinition* particle =
57 // Define histograms start values << 148 G4ParticleTable::GetParticleTable()->FindParticle(fParticleName);
58 << 149
59 const G4String id[] = {"0", "1", "2", "3" << 150 G4cout << "### Fill Cross Sections for " << fParticleName
60 "9", "10", "11", "12 << 151 << " off " << fElementName
61 "18", "19", "20", "21 << 152 << G4endl;
62 G4String title; << 153 if(fVerbose > 0) {
63 << 154 G4cout << "-------------------------------------------------------------"
64 // Default values (to be reset via /analysis << 155 << G4endl;
65 G4int nbins = 100; << 156 G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
66 G4double vmin = 0.; << 157 if(particle == fNeutron) { G4cout << " Capture(b) Fission(b)"; }
67 G4double vmax = 100.; << 158 G4cout << " Total(b)" << G4endl;
68 << 159 G4cout << "-------------------------------------------------------------"
69 // Create all histograms as inactivated << 160 << G4endl;
70 // as we have not yet set nbins, vmin, vmax << 161 }
71 for (G4int k = 0; k < kMaxHisto; k++) { << 162 if(!particle || !elm) {
72 if (k < kMaxAbsor) title = "Edep in absorb << 163 G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
73 if (k > kMaxAbsor) title = "Edep longit. p << 164 return;
74 if (k == 2 * kMaxAbsor + 1) title = "energ << 165 }
75 if (k == 2 * kMaxAbsor + 2) title = "total << 166
76 if (k == 2 * kMaxAbsor + 3) title = "total << 167 G4int prec = G4cout.precision();
77 if (k == 2 * kMaxAbsor + 4) title = "total << 168 G4cout.precision(4);
78 G4int ih = analysisManager->CreateH1(id[k] << 169
79 analysisManager->SetH1Activation(ih, false << 170 G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
>> 171 G4double mass = particle->GetPDGMass();
>> 172
>> 173 // Build histograms
>> 174
>> 175 G4double p1 = std::log10(fMinMomentum/GeV);
>> 176 G4double p2 = std::log10(fMaxMomentum/GeV);
>> 177 G4double e1 = std::log10(fMinKinEnergy/MeV);
>> 178 G4double e2 = std::log10(fMaxKinEnergy/MeV);
>> 179 G4double de = (e2 - e1)/G4double(fBinsE);
>> 180 G4double dp = (p2 - p1)/G4double(fBinsP);
>> 181
>> 182 G4double x = e1 - de*0.5;
>> 183 G4double e, p, xs, xtot;
>> 184 G4int i;
>> 185 for(i=0; i<fBinsE; i++) {
>> 186 x += de;
>> 187 e = std::pow(10.,x)*MeV;
>> 188 if(fVerbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
>> 189 xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
>> 190 xtot = xs;
>> 191 if(fVerbose>0) G4cout << std::setw(12) << xs/barn;
>> 192 fAnalysisManager->FillH1(2, x, xs/barn);
>> 193 xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
>> 194 xtot += xs;
>> 195 if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
>> 196 fAnalysisManager->FillH1(4, x, xs/barn);
>> 197 if(particle == fNeutron) {
>> 198 xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm,mat);
>> 199 xtot += xs;
>> 200 if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
>> 201 fAnalysisManager->FillH1(5, x, xs/barn);
>> 202 xs = store->GetFissionCrossSectionPerAtom(particle,e,elm,mat);
>> 203 xtot += xs;
>> 204 if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
>> 205 fAnalysisManager->FillH1(6, x, xs/barn);
>> 206 }
>> 207 xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm,mat);
>> 208 if(fVerbose>0) G4cout << " " << std::setw(12) << xtot/barn << G4endl;
>> 209 fAnalysisManager->FillH1(7, x, xs/barn);
>> 210 fAnalysisManager->FillH1(8, x, xtot/barn);
>> 211 }
>> 212
>> 213 x = p1 - dp*0.5;
>> 214 for(i=0; i<fBinsP; i++) {
>> 215 x += dp;
>> 216 p = std::pow(10.,x)*GeV;
>> 217 e = std::sqrt(p*p + mass*mass) - mass;
>> 218 xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
>> 219 fAnalysisManager->FillH1(1, x, xs/barn);
>> 220 xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
>> 221 fAnalysisManager->FillH1(3, x, xs/barn);
80 } 222 }
>> 223 if(fVerbose > 0) {
>> 224 G4cout << "-------------------------------------------------------------"
>> 225 << G4endl;
>> 226 }
>> 227 G4cout.precision(prec);
>> 228 fAnalysisManager->Write();
>> 229 fAnalysisManager->CloseFile();
>> 230
>> 231 delete fAnalysisManager;
>> 232 }
>> 233
>> 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 235
>> 236 void HistoManager::SetVerbose(G4int val)
>> 237 {
>> 238 fVerbose = val;
81 } 239 }
>> 240
>> 241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 242
82 243