Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 /// \file hadronic/Hadr00/src/HistoManager.cc 27 /// \brief Implementation of the HistoManager 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.Ivanchen 39 // 40 //-------------------------------------------- 41 // 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 46 #include "HistoManager.hh" 47 48 #include "HistoManagerMessenger.hh" 49 50 #include "G4HadronicProcessStore.hh" 51 #include "G4Neutron.hh" 52 #include "G4NistManager.hh" 53 #include "G4NucleiProperties.hh" 54 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleTable.hh" 56 #include "G4StableIsotopes.hh" 57 #include "G4SystemOfUnits.hh" 58 #include "G4UnitsTable.hh" 59 #include "G4ios.hh" 60 #include "globals.hh" 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 64 HistoManager::HistoManager() 65 { 66 fAnalysisManager = 0; 67 fHistoName = "hadr00"; 68 69 fNeutron = G4Neutron::Neutron(); 70 fMessenger = new HistoManagerMessenger(this) 71 fVerbose = 1; 72 73 fParticleName = "proton"; 74 fElementName = "Al"; 75 76 fTargetMaterial = 0; 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........oo 88 89 HistoManager::~HistoManager() 90 { 91 delete fMessenger; 92 } 93 94 //....oooOO0OOooo........oooOO0OOooo........oo 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= "< 104 105 fAnalysisManager = G4AnalysisManager::Instan 106 fAnalysisManager->OpenFile(fHistoName + ".ro 107 fAnalysisManager->SetFirstHistoId(1); 108 109 fAnalysisManager->CreateH1("h1", "Elastic cr 110 fBinsP, p1, p2); 111 fAnalysisManager->CreateH1("h2", "Elastic cr 112 fBinsE, e1, e2); 113 fAnalysisManager->CreateH1("h3", "Inelastic 114 fBinsP, p1, p2); 115 fAnalysisManager->CreateH1("h4", "Inelastic 116 fBinsE, e1, e2); 117 fAnalysisManager->CreateH1("h5", "Capture cr 118 fBinsE, e1, e2); 119 fAnalysisManager->CreateH1("h6", "Fission cr 120 fBinsE, e1, e2); 121 fAnalysisManager->CreateH1( 122 "h7", "Charge exchange cross section (barn 123 fAnalysisManager->CreateH1("h8", "Total cros 124 fBinsE, e1, e2); 125 fAnalysisManager->CreateH1( 126 "h9", "Inelastic cross section per volume 127 fAnalysisManager->CreateH1( 128 "h10", "Elastic cross section per volume a 129 } 130 131 //....oooOO0OOooo........oooOO0OOooo........oo 132 133 void HistoManager::EndOfRun() 134 { 135 if (fVerbose > 0) { 136 G4cout << "HistoManager: End of run action 137 } 138 139 const G4Element* elm = G4NistManager::Instan 140 const G4Material* mat = G4NistManager::Insta 141 const G4ParticleDefinition* particle = 142 G4ParticleTable::GetParticleTable()->FindP 143 144 G4cout << "### Fill Cross Sections for " << 145 if (fVerbose > 0) { 146 G4cout << "------------------------------- 147 G4cout << " N E(MeV) Elastic(b) 148 if (particle == fNeutron) { 149 G4cout << " Capture(b) Fission(b)"; 150 } 151 G4cout << " Total(b)" << G4endl; 152 G4cout << "------------------------------- 153 } 154 if (!particle || !elm) { 155 G4cout << "HistoManager WARNING Particle o 156 return; 157 } 158 159 G4int prec = G4cout.precision(); 160 G4cout.precision(4); 161 162 G4HadronicProcessStore* store = G4HadronicPr 163 G4double mass = particle->GetPDGMass(); 164 165 // Build histograms 166 167 G4double p1 = std::log10(fMinMomentum / GeV) 168 G4double p2 = std::log10(fMaxMomentum / GeV) 169 G4double e1 = std::log10(fMinKinEnergy / MeV 170 G4double e2 = std::log10(fMaxKinEnergy / MeV 171 G4double de = (e2 - e1) / G4double(fBinsE); 172 G4double dp = (p2 - p1) / G4double(fBinsP); 173 174 G4double x = e1 - de * 0.5; 175 G4double e, p, xs, xtot; 176 G4int i; 177 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; 185 e = std::pow(10., x) * MeV; 186 if (fVerbose > 0) G4cout << std::setw(5) < 187 xs = store->GetElasticCrossSectionPerAtom( 188 xtot = xs; 189 if (fVerbose > 0) G4cout << std::setw(12) 190 fAnalysisManager->FillH1(2, x, xs / barn); 191 xs = store->GetInelasticCrossSectionPerAto 192 xtot += xs; 193 if (fVerbose > 0) G4cout << " " << std::se 194 fAnalysisManager->FillH1(4, x, xs / barn); 195 if (fTargetMaterial) { 196 xs = store->GetInelasticCrossSectionPerV 197 fAnalysisManager->FillH1(9, x, xs / coef 198 xs = store->GetElasticCrossSectionPerVol 199 fAnalysisManager->FillH1(10, x, xs / coe 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 } 211 xs = store->GetChargeExchangeCrossSectionP 212 if (fVerbose > 0) G4cout << " " << std::se 213 fAnalysisManager->FillH1(7, x, xs / barn); 214 fAnalysisManager->FillH1(8, x, xtot / barn 215 } 216 217 x = p1 - dp * 0.5; 218 for (i = 0; i < fBinsP; i++) { 219 x += dp; 220 p = std::pow(10., x) * GeV; 221 e = std::sqrt(p * p + mass * mass) - mass; 222 xs = store->GetElasticCrossSectionPerAtom( 223 fAnalysisManager->FillH1(1, x, xs / barn); 224 xs = store->GetInelasticCrossSectionPerAto 225 fAnalysisManager->FillH1(3, x, xs / barn); 226 } 227 if (fVerbose > 0) { 228 G4cout << "------------------------------- 229 } 230 G4cout.precision(prec); 231 fAnalysisManager->Write(); 232 fAnalysisManager->CloseFile(); 233 fAnalysisManager->Clear(); 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 } 247 248 //....oooOO0OOooo........oooOO0OOooo........oo 249 250 void HistoManager::SetVerbose(G4int val) 251 { 252 fVerbose = val; 253 } 254 255 //....oooOO0OOooo........oooOO0OOooo........oo 256