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