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 // $Id: G4BGGNucleonElasticXS.cc 93682 2015-10-28 10:09:49Z gcosmo $ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class file 30 // GEANT4 Class file 29 // 31 // 30 // 32 // 31 // File name: G4BGGNucleonElasticXS 33 // File name: G4BGGNucleonElasticXS 32 // 34 // 33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 34 // 36 // 35 // Creation date: 13.03.2007 37 // Creation date: 13.03.2007 >> 38 // Modifications: >> 39 // 36 // 40 // 37 // ------------------------------------------- 41 // ------------------------------------------------------------------- 38 // 42 // 39 43 40 #include "G4BGGNucleonElasticXS.hh" 44 #include "G4BGGNucleonElasticXS.hh" 41 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh" 42 #include "G4ComponentGGHadronNucleusXsc.hh" 46 #include "G4ComponentGGHadronNucleusXsc.hh" 43 #include "G4NucleonNuclearCrossSection.hh" 47 #include "G4NucleonNuclearCrossSection.hh" 44 #include "G4HadronNucleonXsc.hh" 48 #include "G4HadronNucleonXsc.hh" 45 #include "G4NuclearRadii.hh" << 49 #include "G4ComponentSAIDTotalXS.hh" 46 #include "G4Proton.hh" 50 #include "G4Proton.hh" 47 #include "G4Neutron.hh" 51 #include "G4Neutron.hh" 48 #include "G4NistManager.hh" 52 #include "G4NistManager.hh" 49 #include "G4NuclearRadii.hh" << 53 #include "G4Log.hh" >> 54 #include "G4Exp.hh" 50 55 51 #include "G4CrossSectionDataSetRegistry.hh" 56 #include "G4CrossSectionDataSetRegistry.hh" 52 57 53 G4double G4BGGNucleonElasticXS::theGlauberFacP << 58 const G4double llog10 = G4Log(10.); 54 G4double G4BGGNucleonElasticXS::theCoulombFacP << 55 G4double G4BGGNucleonElasticXS::theGlauberFacN << 56 G4double G4BGGNucleonElasticXS::theCoulombFacN << 57 G4int G4BGGNucleonElasticXS::theA[93] = {0}; << 58 59 59 G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(c 60 G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition* p) 60 : G4VCrossSectionDataSet("BarashenkovGlauberG << 61 : G4VCrossSectionDataSet("Barashenkov-Glauber") 61 { 62 { 62 verboseLevel = 0; 63 verboseLevel = 0; 63 fGlauberEnergy = 91.*GeV; 64 fGlauberEnergy = 91.*GeV; 64 fLowEnergy = 14.0*MeV; << 65 fPDGEnergy = 5*GeV; 65 fNucleon = new G4NucleonNuclearCrossSection( << 66 fLowEnergy = 14.*MeV; 66 fGlauber = new G4ComponentGGHadronNucleusXsc << 67 fSAIDLowEnergyLimit = 1*MeV; 67 fHadron = new G4HadronNucleonXsc(); << 68 fSAIDHighEnergyLimit = 1.3*GeV; 68 << 69 fLowestXSection = millibarn; 69 theProton = G4Proton::Proton(); << 70 for (G4int i = 0; i < 93; ++i) { 70 isProton = (theProton == p); << 71 theGlauberFac[i] = 0.0; 71 SetForAllAtomsAndEnergies(true); << 72 theCoulombFac[i] = 0.0; 72 << 73 theA[i] = 1; 73 if (0 == theA[0]) { Initialise(); } << 74 } >> 75 fNucleon = 0; >> 76 fGlauber = 0; >> 77 fHadron = 0; >> 78 fSAID = 0; >> 79 particle = p; >> 80 theProton= G4Proton::Proton(); >> 81 isProton = false; >> 82 isInitialized = false; 74 } 83 } 75 84 76 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 77 86 78 G4BGGNucleonElasticXS::~G4BGGNucleonElasticXS( 87 G4BGGNucleonElasticXS::~G4BGGNucleonElasticXS() 79 { 88 { >> 89 delete fSAID; 80 delete fHadron; 90 delete fHadron; >> 91 // The cross section registry will delete fNucleon >> 92 delete fGlauber; 81 } 93 } 82 94 83 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 84 96 85 G4bool 97 G4bool 86 G4BGGNucleonElasticXS::IsElementApplicable(con 98 G4BGGNucleonElasticXS::IsElementApplicable(const G4DynamicParticle*, G4int, 87 con << 99 const G4Material*) 88 { 100 { 89 return true; 101 return true; 90 } 102 } 91 103 92 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 93 105 94 G4bool G4BGGNucleonElasticXS::IsIsoApplicable( 106 G4bool G4BGGNucleonElasticXS::IsIsoApplicable(const G4DynamicParticle*, 95 << 107 G4int Z, G4int, 96 << 108 const G4Element*, 97 << 109 const G4Material*) 98 { 110 { 99 return (1 == Z); 111 return (1 == Z); 100 } 112 } 101 113 102 //....oooOO0OOooo........oooOO0OOooo........oo 114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 103 115 104 G4double 116 G4double 105 G4BGGNucleonElasticXS::GetElementCrossSection( 117 G4BGGNucleonElasticXS::GetElementCrossSection(const G4DynamicParticle* dp, 106 << 118 G4int ZZ, const G4Material*) 107 { 119 { 108 // this method should be called only for Z > 120 // this method should be called only for Z > 1 109 121 110 G4double cross = 0.0; 122 G4double cross = 0.0; 111 G4int Z = std::min(ZZ, 92); << 112 G4double ekin = dp->GetKineticEnergy(); 123 G4double ekin = dp->GetKineticEnergy(); >> 124 G4int Z = ZZ; 113 if(1 == Z) { 125 if(1 == Z) { 114 cross = 1.0115*GetIsoCrossSection(dp,1,1); 126 cross = 1.0115*GetIsoCrossSection(dp,1,1); 115 } else { 127 } else { >> 128 if(Z > 92) { Z = 92; } >> 129 116 if(ekin <= fLowEnergy) { 130 if(ekin <= fLowEnergy) { 117 cross = (isProton) ? theCoulombFacP[Z] : << 131 cross = theCoulombFac[Z]*CoulombFactor(ekin, Z); 118 cross *= CoulombFactor(ekin, Z); << 132 119 } else if(ekin > fGlauberEnergy) { 133 } else if(ekin > fGlauberEnergy) { 120 cross = (isProton) ? theGlauberFacP[Z] : << 134 cross = theGlauberFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, theA[Z]); 121 cross *= fGlauber->GetElasticGlauberGrib << 122 } else { 135 } else { 123 cross = fNucleon->GetElasticCrossSection 136 cross = fNucleon->GetElasticCrossSection(dp, Z); 124 } 137 } 125 } 138 } 126 #ifdef G4VERBOSE << 139 127 if (verboseLevel > 1) { << 140 if(verboseLevel > 1) { 128 G4cout << "G4BGGNucleonElasticXS::GetEleme 141 G4cout << "G4BGGNucleonElasticXS::GetElementCrossSection for " 129 << dp->GetDefinition()->GetParticle << 142 << dp->GetDefinition()->GetParticleName() 130 << " Ekin(GeV)= " << dp->GetKineti << 143 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV 131 << " in nucleus Z= " << Z << " A= << 144 << " in nucleus Z= " << Z << " A= " << theA[Z] 132 << " XS(b)= " << cross/barn << 145 << " XS(b)= " << cross/barn 133 << G4endl; << 146 << G4endl; 134 } 147 } 135 #endif << 148 //AR-04Dec2013 if(cross <= fLowestXSection) { cross = 0.0; } 136 return cross; 149 return cross; 137 } 150 } 138 151 139 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 140 153 141 G4double 154 G4double 142 G4BGGNucleonElasticXS::GetIsoCrossSection(cons 155 G4BGGNucleonElasticXS::GetIsoCrossSection(const G4DynamicParticle* dp, 143 G4in << 156 G4int Z, G4int A, 144 cons << 157 const G4Isotope*, 145 cons << 158 const G4Element*, 146 cons << 159 const G4Material*) 147 { 160 { 148 // this method should be called only for Z = 161 // this method should be called only for Z = 1 149 fHadron->HadronNucleonXscNS(dp->GetDefinitio << 150 dp->GetKineticEn << 151 G4double cross = A*fHadron->GetElasticHadron << 152 162 153 #ifdef G4VERBOSE << 163 G4double cross = 0.0; 154 if (verboseLevel > 1) { << 164 G4double ekin = dp->GetKineticEnergy(); >> 165 >> 166 if(ekin <= fSAIDLowEnergyLimit) { >> 167 cross = theCoulombFac[0]*CoulombFactor(ekin, 1); >> 168 } else if(ekin <= fSAIDHighEnergyLimit) { >> 169 cross = fSAID->GetElasticIsotopeCrossSection(particle, ekin, 1, 1); >> 170 } else if(ekin <= fPDGEnergy) { >> 171 G4double cross1 = >> 172 fSAID->GetElasticIsotopeCrossSection(particle, fPDGEnergy, 1, 1); >> 173 G4double cross2 = theCoulombFac[1]; >> 174 cross = cross1 + (cross2 - cross1)*(ekin - fSAIDHighEnergyLimit) >> 175 /(fPDGEnergy - fSAIDHighEnergyLimit); >> 176 } else { >> 177 fHadron->GetHadronNucleonXscPDG(dp, theProton); >> 178 cross = fHadron->GetElasticHadronNucleonXsc(); >> 179 } >> 180 cross *= A; >> 181 >> 182 if(verboseLevel > 1) { 155 G4cout << "G4BGGNucleonElasticXS::GetIsoCr 183 G4cout << "G4BGGNucleonElasticXS::GetIsoCrossSection for " 156 << dp->GetDefinition()->GetParticle << 184 << dp->GetDefinition()->GetParticleName() 157 << " Ekin(GeV)= " << dp->GetKineti << 185 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV 158 << " in nucleus Z=1 A=" << A << 186 << " in nucleus Z= " << Z << " A= " << A 159 << " XS(b)= " << cross/barn << 187 << " XS(b)= " << cross/barn 160 << G4endl; << 188 << G4endl; 161 } 189 } 162 #endif << 190 //AR-04Dec2013 if(cross <= fLowestXSection) { cross = 0.0; } 163 return cross; 191 return cross; 164 } 192 } 165 193 166 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 167 195 168 void G4BGGNucleonElasticXS::BuildPhysicsTable( 196 void G4BGGNucleonElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p) 169 { 197 { 170 if(&p == theProton || &p == G4Neutron::Neutr 198 if(&p == theProton || &p == G4Neutron::Neutron()) { 171 isProton = (theProton == &p); << 199 particle = &p; 172 200 173 } else { 201 } else { 174 G4ExceptionDescription ed; << 202 G4cout << "### G4BGGNucleonElasticXS WARNING: is not applicable to " 175 ed << "This BGG cross section is applicabl << 203 << p.GetParticleName() 176 << p.GetParticleName() << G4endl; << 204 << G4endl; 177 G4Exception("G4BGGNucleonElasticXS::BuildP << 205 throw G4HadronicException(__FILE__, __LINE__, 178 FatalException, ed); << 206 "G4BGGNucleonElasticXS::BuildPhysicsTable is used for wrong particle"); >> 207 return; 179 } 208 } 180 } << 181 209 182 //....oooOO0OOooo........oooOO0OOooo........oo << 210 if(isInitialized) { return; } >> 211 isInitialized = true; >> 212 >> 213 fNucleon = (G4NucleonNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NucleonNuclearCrossSection::Default_Name()); >> 214 fGlauber = new G4ComponentGGHadronNucleusXsc(); >> 215 fHadron = new G4HadronNucleonXsc(); >> 216 fSAID = new G4ComponentSAIDTotalXS(); >> 217 >> 218 fNucleon->BuildPhysicsTable(*particle); >> 219 fGlauber->BuildPhysicsTable(*particle); >> 220 >> 221 if(particle == theProton) { >> 222 isProton = true; >> 223 fSAIDHighEnergyLimit = 3*GeV; >> 224 } 183 225 184 void G4BGGNucleonElasticXS::Initialise() << 185 { << 186 theA[0] = theA[1] = 1; << 187 G4ThreeVector mom(0.0,0.0,1.0); 226 G4ThreeVector mom(0.0,0.0,1.0); 188 G4DynamicParticle dp(theProton, mom, fGlaube << 227 G4DynamicParticle dp(particle, mom, fGlauberEnergy); 189 228 190 G4NistManager* nist = G4NistManager::Instanc 229 G4NistManager* nist = G4NistManager::Instance(); >> 230 191 G4double csup, csdn; 231 G4double csup, csdn; >> 232 G4int A; >> 233 >> 234 if(verboseLevel > 0) { >> 235 G4cout << "### G4BGGNucleonElasticXS::Initialise for " >> 236 << particle->GetParticleName() << G4endl; >> 237 } >> 238 >> 239 for(G4int iz=2; iz<93; iz++) { 192 240 193 for (G4int iz=2; iz<93; ++iz) { << 241 A = G4lrint(nist->GetAtomicMassAmu(iz)); 194 G4int A = G4lrint(nist->GetAtomicMassAmu(i << 195 theA[iz] = A; 242 theA[iz] = A; 196 243 197 csup = fGlauber->GetElasticGlauberGribov(& 244 csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A); 198 csdn = fNucleon->GetElasticCrossSection(&d 245 csdn = fNucleon->GetElasticCrossSection(&dp, iz); 199 theGlauberFacP[iz] = csdn/csup; << 200 } << 201 << 202 dp.SetDefinition(G4Neutron::Neutron()); << 203 for (G4int iz=2; iz<93; ++iz) { << 204 csup = fGlauber->GetElasticGlauberGribov(& << 205 csdn = fNucleon->GetElasticCrossSection(&d << 206 theGlauberFacN[iz] = csdn/csup; << 207 246 208 if (verboseLevel > 1) { << 247 theGlauberFac[iz] = csdn/csup; 209 G4cout << "G4BGGNucleonElasticXS::Init << 248 if(verboseLevel > 0) { 210 << " GFactorP=" << theGlauberFa << 249 G4cout << "Z= " << iz << " A= " << A 211 << " GFactorN=" << theGlauberFa << 250 << " factor= " << theGlauberFac[iz] << G4endl; 212 } 251 } 213 } 252 } 214 253 215 theCoulombFacP[0] = theCoulombFacP[1] = << 254 theCoulombFac[0] = 216 theCoulombFacN[0] = theCoulombFacN[1] = 1. << 255 fSAID->GetElasticIsotopeCrossSection(particle,fSAIDLowEnergyLimit,1,1) 217 dp.SetDefinition(theProton); << 256 /CoulombFactor(fSAIDLowEnergyLimit, 1); >> 257 >> 258 dp.SetKineticEnergy(fPDGEnergy); >> 259 fHadron->GetHadronNucleonXscPDG(&dp, theProton); >> 260 theCoulombFac[1] = fHadron->GetElasticHadronNucleonXsc(); >> 261 >> 262 if(verboseLevel > 0) { >> 263 G4cout << "Z=1 A=1 " << " factor0= " << theCoulombFac[0] >> 264 << " factor1= " << theCoulombFac[1] >> 265 << G4endl; >> 266 } >> 267 218 dp.SetKineticEnergy(fLowEnergy); 268 dp.SetKineticEnergy(fLowEnergy); 219 for (G4int iz=2; iz<93; ++iz) { << 269 for(G4int iz=2; iz<93; iz++) { 220 theCoulombFacP[iz] = fNucleon->GetElasticC << 270 theCoulombFac[iz] = 221 /CoulombFactor(fLowEnergy, iz); << 271 fNucleon->GetElasticCrossSection(&dp, iz)/CoulombFactor(fLowEnergy, iz); 222 } << 272 if(verboseLevel > 0) { 223 dp.SetDefinition(G4Neutron::Neutron()); << 273 G4cout << "Z= " << iz << " A= " << theA[iz] 224 for(G4int iz=2; iz<93; ++iz) { << 274 << " factor= " << theCoulombFac[iz] << G4endl; 225 theCoulombFacN[iz] = fNucleon->GetElasticC << 226 /CoulombFactor(fLowEnergy, iz); << 227 << 228 if (verboseLevel > 1) { << 229 G4cout << "G4BGGNucleonElasticXS::Initia << 230 << " CFactorP=" << theCoulombFacP[iz] << 231 << " CFactorN=" << theCoulombFacN[iz] < << 232 } 275 } 233 } 276 } 234 } 277 } 235 278 236 //....oooOO0OOooo........oooOO0OOooo........oo 279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 237 280 238 G4double G4BGGNucleonElasticXS::CoulombFactor( 281 G4double G4BGGNucleonElasticXS::CoulombFactor(G4double kinEnergy, G4int Z) 239 { 282 { 240 return (isProton) ? << 283 G4double res= 1.0; 241 G4NuclearRadii::CoulombFactor(Z, theA[Z], << 284 >> 285 // from G4ProtonInelasticCrossSection >> 286 if(isProton) { >> 287 >> 288 if (Z <= 1) { return kinEnergy*kinEnergy; } >> 289 >> 290 G4double elog = G4Log(kinEnergy/GeV)/llog10; >> 291 G4double aa = theA[Z]; >> 292 >> 293 G4double ff1 = 5.6 - 0.016*aa; // slope of the drop at medium energies. >> 294 G4double ff2 = 1.37 + 1.37/aa; // start of the slope. >> 295 G4double ff3 = 0.8 + 18./aa - 0.002*aa; // stephight >> 296 res = 1.0 + ff3*(1.0 - (1.0/(1+G4Exp(-ff1*(elog + ff2))))); >> 297 >> 298 ff1 = 8. - 8./aa - 0.008*aa; // slope of the rise >> 299 ff2 = 2.34 - 5.4/aa - 0.0028*aa; // start of the rise >> 300 res /= (1.0 + G4Exp(-ff1*(elog + ff2))); >> 301 >> 302 } >> 303 return res; 242 } 304 } 243 305 244 //....oooOO0OOooo........oooOO0OOooo........oo 306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 245 307 246 void G4BGGNucleonElasticXS::CrossSectionDescri 308 void G4BGGNucleonElasticXS::CrossSectionDescription(std::ostream& outFile) const 247 { 309 { 248 outFile << "The Barashenkov-Glauber-Gribov c 310 outFile << "The Barashenkov-Glauber-Gribov cross section handles elastic\n" 249 << "scattering of protons and neutro 311 << "scattering of protons and neutrons from nuclei using the\n" 250 << "Barashenkov parameterization bel 312 << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n" 251 << "parameterization above 91 GeV. n 313 << "parameterization above 91 GeV. n"; 252 } 314 } 253 315 254 //....oooOO0OOooo........oooOO0OOooo........oo 316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 255 317