Geant4 Cross Reference |
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 // 26 // ------------------------------------------------------------------- 27 // 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4BGGNucleonInelasticXS 32 // 33 // Author: Vladimir Ivanchenko 34 // 35 // Creation date: 13.03.2007 36 // Modifications: 37 // 38 // 39 // ------------------------------------------------------------------- 40 // 41 42 #include "G4BGGNucleonInelasticXS.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4ComponentGGHadronNucleusXsc.hh" 45 #include "G4NucleonNuclearCrossSection.hh" 46 #include "G4HadronNucleonXsc.hh" 47 #include "G4ComponentSAIDTotalXS.hh" 48 #include "G4Proton.hh" 49 #include "G4Neutron.hh" 50 #include "G4NistManager.hh" 51 #include "G4Material.hh" 52 #include "G4Element.hh" 53 #include "G4Isotope.hh" 54 #include "G4Log.hh" 55 #include "G4Exp.hh" 56 #include "G4NuclearRadii.hh" 57 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 59 60 namespace 61 { 62 const G4double llog10 = G4Log(10.); 63 } 64 65 G4double G4BGGNucleonInelasticXS::theGlauberFacP[93] = {0.0}; 66 G4double G4BGGNucleonInelasticXS::theCoulombFacP[93] = {0.0}; 67 G4double G4BGGNucleonInelasticXS::theGlauberFacN[93] = {0.0}; 68 G4double G4BGGNucleonInelasticXS::theCoulombFacN[93] = {0.0}; 69 G4int G4BGGNucleonInelasticXS::theA[93] = {0}; 70 71 G4BGGNucleonInelasticXS::G4BGGNucleonInelasticXS(const G4ParticleDefinition* p) 72 : G4VCrossSectionDataSet("BarashenkovGlauberGribov") 73 { 74 verboseLevel = 0; 75 fGlauberEnergy = 91.*CLHEP::GeV; 76 fLowEnergy = 14.*CLHEP::MeV; 77 78 fNucleon = new G4NucleonNuclearCrossSection(); 79 fGlauber = new G4ComponentGGHadronNucleusXsc(); 80 fHadron = new G4HadronNucleonXsc(); 81 82 theProton= G4Proton::Proton(); 83 isProton = (theProton == p); 84 SetForAllAtomsAndEnergies(true); 85 86 if (0 == theA[0]) { Initialise(); } 87 } 88 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 90 91 G4BGGNucleonInelasticXS::~G4BGGNucleonInelasticXS() 92 { 93 delete fHadron; 94 } 95 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 97 98 G4bool G4BGGNucleonInelasticXS::IsElementApplicable(const G4DynamicParticle*, 99 G4int, const G4Material*) 100 { 101 return true; 102 } 103 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 105 106 G4bool G4BGGNucleonInelasticXS::IsIsoApplicable(const G4DynamicParticle*, 107 G4int Z, G4int, 108 const G4Element*, 109 const G4Material*) 110 { 111 return (1 == Z); 112 } 113 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 115 116 G4double 117 G4BGGNucleonInelasticXS::GetElementCrossSection(const G4DynamicParticle* dp, 118 G4int ZZ, const G4Material*) 119 { 120 G4double cross = 0.0; 121 G4double ekin = dp->GetKineticEnergy(); 122 G4int Z = std::min(ZZ, 92); 123 if (1 == Z) { 124 cross = 1.0115*GetIsoCrossSection(dp,1,1); 125 } else { 126 if (ekin <= fLowEnergy) { 127 cross = CoulombFactor(ekin, Z); 128 cross *= (isProton) ? theCoulombFacP[Z] : theCoulombFacN[Z]; 129 } else if (ekin > fGlauberEnergy) { 130 cross = fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]); 131 cross *= (isProton) ? theGlauberFacP[Z] : theGlauberFacN[Z]; 132 } else { 133 cross = fNucleon->GetElementCrossSection(dp, Z); 134 } 135 } 136 137 #ifdef G4VERBOSE 138 if (verboseLevel > 1) { 139 G4cout << "G4BGGNucleonInelasticXS::GetCrossSection for " 140 << dp->GetDefinition()->GetParticleName() 141 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV 142 << " in nucleus Z= " << Z << " A= " << theA[Z] 143 << " XS(b)= " << cross/barn 144 << G4endl; 145 } 146 #endif 147 return cross; 148 } 149 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 151 152 G4double 153 G4BGGNucleonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* dp, 154 G4int, G4int A, 155 const G4Isotope*, 156 const G4Element*, 157 const G4Material*) 158 { 159 // this method should be called only for Z = 1 160 fHadron->HadronNucleonXscNS(dp->GetDefinition(), theProton, 161 dp->GetKineticEnergy()); 162 G4double cross = A*fHadron->GetInelasticHadronNucleonXsc(); 163 164 #ifdef G4VERBOSE 165 if(verboseLevel > 1) { 166 G4cout << "G4BGGNucleonInelasticXS::GetIsoCrossSection for " 167 << dp->GetDefinition()->GetParticleName() 168 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV 169 << " in nucleus Z=1 A=" << A 170 << " XS(b)= " << cross/barn 171 << G4endl; 172 } 173 #endif 174 return cross; 175 } 176 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 178 179 void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p) 180 { 181 if(&p == theProton || &p == G4Neutron::Neutron()) { 182 isProton = (theProton == &p); 183 } else { 184 G4ExceptionDescription ed; 185 ed << "This BGG cross section is applicable only to nucleons and not to " 186 << p.GetParticleName() << G4endl; 187 G4Exception("G4BGGNucleonInelasticXS::BuildPhysicsTable", "had001", 188 FatalException, ed); 189 return; 190 } 191 } 192 193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 194 195 void G4BGGNucleonInelasticXS::Initialise() 196 { 197 theA[0] = theA[1] = 1; 198 G4ThreeVector mom(0.0,0.0,1.0); 199 G4DynamicParticle dp(theProton, mom, fGlauberEnergy); 200 201 G4NistManager* nist = G4NistManager::Instance(); 202 G4double csup, csdn; 203 204 for (G4int iz=2; iz<93; ++iz) { 205 206 G4int A = G4lrint(nist->GetAtomicMassAmu(iz)); 207 theA[iz] = A; 208 209 csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A); 210 csdn = fNucleon->GetElementCrossSection(&dp, iz); 211 theGlauberFacP[iz] = csdn/csup; 212 } 213 214 dp.SetDefinition(G4Neutron::Neutron()); 215 for (G4int iz=2; iz<93; ++iz) { 216 csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, theA[iz]); 217 csdn = fNucleon->GetElementCrossSection(&dp, iz); 218 theGlauberFacN[iz] = csdn/csup; 219 220 if(verboseLevel > 1) { 221 G4cout << "G4BGGNucleonInelasticXS::Initialise Z= " << iz << " A= " << theA[iz] 222 << " GFactorP= " << theGlauberFacP[iz] 223 << " GFactorN= " << theGlauberFacN[iz] << G4endl; 224 } 225 } 226 227 theCoulombFacP[1] = theCoulombFacN[1] = 1.0; 228 dp.SetDefinition(theProton); 229 dp.SetKineticEnergy(fLowEnergy); 230 for (G4int iz=2; iz<93; ++iz) { 231 theCoulombFacP[iz] = fNucleon->GetElementCrossSection(&dp, iz) 232 /CoulombFactor(fLowEnergy, iz); 233 } 234 dp.SetDefinition(G4Neutron::Neutron()); 235 for (G4int iz=2; iz<93; ++iz) { 236 theCoulombFacN[iz] = fNucleon->GetElementCrossSection(&dp, iz) 237 /CoulombFactor(fLowEnergy, iz); 238 239 if (verboseLevel > 1) { 240 G4cout << "G4BGGNucleonInelasticXS::Initialise Z= " << iz << " A= " << theA[iz] 241 << " CFactorP= " << theCoulombFacP[iz] 242 << " CFactorN= " << theCoulombFacN[iz] << G4endl; 243 } 244 } 245 } 246 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 248 249 G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4int Z) 250 { 251 G4double res = 0.0; 252 253 if(kinEnergy <= 0.0) { return res; } 254 255 G4double elog = G4Log(kinEnergy/GeV)/llog10; 256 G4double aa = theA[Z]; 257 258 if(isProton) { 259 260 res = G4NuclearRadii::CoulombFactor(Z, aa, theProton, kinEnergy); 261 262 // from G4ProtonInelasticCrossSection 263 if(res > 0.0) { 264 G4double ff1 = 5.6 - 0.016*aa; // slope of the drop at medium energies. 265 G4double ff2 = 1.37 + 1.37/aa; // start of the slope. 266 G4double ff3 = 0.8 + 18./aa - 0.002*aa; // stephight 267 res *= (1.0 + ff3*(1.0 - (1.0/(1+G4Exp(-ff1*(elog + ff2)))))); 268 ff1 = 8. - 8./aa - 0.008*aa; // slope of the rise 269 ff2 = 2.34 - 5.4/aa - 0.0028*aa; // start of the rise 270 res /= (1.0 + G4Exp(-ff1*(elog + ff2))); 271 } 272 } else { 273 // from G4NeutronInelasticCrossSection 274 G4double p3 = 0.6 + 13./aa - 0.0005*aa; 275 G4double p4 = 7.2449 - 0.018242*aa; 276 G4double p5 = 1.36 + 1.8/aa + 0.0005*aa; 277 G4double p6 = 1. + 200./aa + 0.02*aa; 278 G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.; 279 280 G4double firstexp = G4Exp(-p4*(elog + p5)); 281 G4double secondexp = G4Exp(-p6*(elog + p7)); 282 283 res = (1. + p3*firstexp/(1. + firstexp))/(1. + secondexp); 284 } 285 return res; 286 } 287 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 289 290 void G4BGGNucleonInelasticXS::CrossSectionDescription(std::ostream& outFile) const 291 { 292 outFile << "The Barashenkov-Glauber-Gribov cross section calculates inelastic\n" 293 << "scattering of protons and neutrons from nuclei using the\n" 294 << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n" 295 << "parameterization above 91 GeV. It uses the G4HadronNucleonXsc\n" 296 << "cross section component for hydrogen targets, and the\n" 297 << "G4ComponentGGHadronNucleusXsc component for other targets.\n"; 298 } 299 300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 301