00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4BGGNucleonElasticXS.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4GlauberGribovCrossSection.hh"
00047 #include "G4NucleonNuclearCrossSection.hh"
00048 #include "G4HadronNucleonXsc.hh"
00049 #include "G4ComponentSAIDTotalXS.hh"
00050 #include "G4Proton.hh"
00051 #include "G4Neutron.hh"
00052 #include "G4NistManager.hh"
00053
00054 #include "G4CrossSectionDataSetRegistry.hh"
00055
00056 G4BGGNucleonElasticXS::G4BGGNucleonElasticXS(const G4ParticleDefinition* p)
00057 : G4VCrossSectionDataSet("Barashenkov-Glauber")
00058 {
00059 verboseLevel = 0;
00060 fGlauberEnergy = 91.*GeV;
00061 fLowEnergy = 14.*MeV;
00062 fSAIDHighEnergyLimit = 1.3*GeV;
00063 for (G4int i = 0; i < 93; ++i) {
00064 theGlauberFac[i] = 0.0;
00065 theCoulombFac[i] = 0.0;
00066 theA[i] = 1;
00067 }
00068 fNucleon = 0;
00069 fGlauber = 0;
00070 fHadron = 0;
00071 fSAID = 0;
00072 particle = p;
00073 theProton= G4Proton::Proton();
00074 isProton = false;
00075 isInitialized = false;
00076 }
00077
00078
00079
00080 G4BGGNucleonElasticXS::~G4BGGNucleonElasticXS()
00081 {
00082 delete fHadron;
00083 delete fSAID;
00084 }
00085
00086
00087
00088 G4bool
00089 G4BGGNucleonElasticXS::IsElementApplicable(const G4DynamicParticle*, G4int,
00090 const G4Material*)
00091 {
00092 return true;
00093 }
00094
00095
00096
00097 G4bool G4BGGNucleonElasticXS::IsIsoApplicable(const G4DynamicParticle*,
00098 G4int Z, G4int A,
00099 const G4Element*,
00100 const G4Material*)
00101 {
00102 return (1 == Z && 2 >= A);
00103 }
00104
00105
00106
00107 G4double
00108 G4BGGNucleonElasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
00109 G4int ZZ, const G4Material*)
00110 {
00111
00112
00113 G4double cross = 0.0;
00114 G4double ekin = dp->GetKineticEnergy();
00115 G4int Z = ZZ;
00116 if(1 == Z) {
00117 cross = 1.0115*GetIsoCrossSection(dp,1,1);
00118 } else {
00119 if(Z > 92) { Z = 92; }
00120
00121 if(ekin <= fLowEnergy) {
00122 cross = theCoulombFac[Z];
00123
00124 } else if(ekin > fGlauberEnergy) {
00125 cross = theGlauberFac[Z]*fGlauber->GetElasticGlauberGribov(dp, Z, theA[Z]);
00126 } else {
00127 cross = fNucleon->GetElasticCrossSection(dp, Z);
00128 }
00129 }
00130
00131 if(verboseLevel > 1) {
00132 G4cout << "G4BGGNucleonElasticXS::GetElementCrossSection for "
00133 << dp->GetDefinition()->GetParticleName()
00134 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
00135 << " in nucleus Z= " << Z << " A= " << theA[Z]
00136 << " XS(b)= " << cross/barn
00137 << G4endl;
00138 }
00139 return cross;
00140 }
00141
00142
00143
00144 G4double
00145 G4BGGNucleonElasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
00146 G4int Z, G4int A,
00147 const G4Isotope*,
00148 const G4Element*,
00149 const G4Material*)
00150 {
00151
00152
00153 G4double cross = 0.0;
00154 G4double ekin = dp->GetKineticEnergy();
00155
00156 if(ekin <= fSAIDHighEnergyLimit) {
00157 cross = fSAID->GetElasticIsotopeCrossSection(particle, ekin, 1, 1);
00158 } else {
00159 fHadron->GetHadronNucleonXscPDG(dp, theProton);
00160 cross = theCoulombFac[1]*fHadron->GetElasticHadronNucleonXsc();
00161 }
00162
00163
00164
00165 cross *= A;
00166
00167 if(verboseLevel > 1) {
00168 G4cout << "G4BGGNucleonElasticXS::GetIsoCrossSection for "
00169 << dp->GetDefinition()->GetParticleName()
00170 << " Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
00171 << " in nucleus Z= " << Z << " A= " << A
00172 << " XS(b)= " << cross/barn
00173 << G4endl;
00174 }
00175 return cross;
00176 }
00177
00178
00179
00180 void G4BGGNucleonElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
00181 {
00182 if(&p == theProton || &p == G4Neutron::Neutron()) {
00183 particle = &p;
00184
00185 } else {
00186 G4cout << "### G4BGGNucleonElasticXS WARNING: is not applicable to "
00187 << p.GetParticleName()
00188 << G4endl;
00189 throw G4HadronicException(__FILE__, __LINE__,
00190 "G4BGGNucleonElasticXS::BuildPhysicsTable is used for wrong particle");
00191 return;
00192 }
00193
00194 if(isInitialized) { return; }
00195 isInitialized = true;
00196
00197 fNucleon = (G4NucleonNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NucleonNuclearCrossSection::Default_Name());
00198 fGlauber = (G4GlauberGribovCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4GlauberGribovCrossSection::Default_Name());
00199
00200 fHadron = new G4HadronNucleonXsc();
00201 fSAID = new G4ComponentSAIDTotalXS();
00202 fNucleon->BuildPhysicsTable(*particle);
00203 fGlauber->BuildPhysicsTable(*particle);
00204 if(particle == theProton) {
00205 isProton = true;
00206 fSAIDHighEnergyLimit = 3*GeV;
00207 }
00208
00209 G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
00210 G4ThreeVector mom(0.0,0.0,1.0);
00211 G4DynamicParticle dp(part, mom, fGlauberEnergy);
00212
00213 G4NistManager* nist = G4NistManager::Instance();
00214
00215 G4double csup, csdn;
00216 G4int A;
00217
00218 if(verboseLevel > 0) {
00219 G4cout << "### G4BGGNucleonElasticXS::Initialise for "
00220 << particle->GetParticleName() << G4endl;
00221 }
00222
00223 for(G4int iz=2; iz<93; iz++) {
00224
00225 A = G4lrint(nist->GetAtomicMassAmu(iz));
00226 theA[iz] = A;
00227
00228 csup = fGlauber->GetElasticGlauberGribov(&dp, iz, A);
00229 csdn = fNucleon->GetElasticCrossSection(&dp, iz);
00230
00231 theGlauberFac[iz] = csdn/csup;
00232 if(verboseLevel > 0) {
00233 G4cout << "Z= " << iz << " A= " << A
00234 << " factor= " << theGlauberFac[iz] << G4endl;
00235 }
00236 }
00237 dp.SetKineticEnergy(fSAIDHighEnergyLimit);
00238 fHadron->GetHadronNucleonXscPDG(&dp, theProton);
00239 theCoulombFac[1] =
00240 fSAID->GetElasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)
00241 /fHadron->GetElasticHadronNucleonXsc();
00242 if(verboseLevel > 0) {
00243 G4cout << "Z=1 A=1" << " CoulombFactor= " << theCoulombFac[1] << G4endl;
00244 }
00245
00246 dp.SetKineticEnergy(fLowEnergy);
00247 for(G4int iz=2; iz<93; iz++) {
00248 theCoulombFac[iz] = fNucleon->GetElasticCrossSection(&dp, iz);
00249 if(verboseLevel > 0) {
00250 G4cout << "Z= " << iz << " A= " << theA[iz]
00251 << " factor= " << theCoulombFac[iz] << G4endl;
00252 }
00253 }
00254 }
00255
00256
00257
00258 void G4BGGNucleonElasticXS::CrossSectionDescription(std::ostream& outFile) const
00259 {
00260 outFile << "The Barashenkov-Glauber-Gribov cross section handles elastic\n"
00261 << "scattering of protons and neutrons from nuclei using the\n"
00262 << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
00263 << "parameterization above 91 GeV. n";
00264 }
00265
00266