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 #include "G4BGGPionInelasticXS.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4GlauberGribovCrossSection.hh"
00046 #include "G4UPiNuclearCrossSection.hh"
00047 #include "G4HadronNucleonXsc.hh"
00048
00049 #include "G4ComponentSAIDTotalXS.hh"
00050
00051 #include "G4Proton.hh"
00052 #include "G4PionPlus.hh"
00053 #include "G4PionMinus.hh"
00054 #include "G4NistManager.hh"
00055
00056
00057 G4BGGPionInelasticXS::G4BGGPionInelasticXS(const G4ParticleDefinition* p)
00058 : G4VCrossSectionDataSet("Barashenkov-Glauber-Gribov")
00059 {
00060 verboseLevel = 0;
00061 fGlauberEnergy = 91.*GeV;
00062 fLowEnergy = 20.*MeV;
00063 fSAIDHighEnergyLimit = 2.6*GeV;
00064 SetMinKinEnergy(0.0);
00065 SetMaxKinEnergy(100*TeV);
00066
00067 for (G4int i = 0; i < 93; i++) {
00068 theGlauberFac[i] = 0.0;
00069 theCoulombFac[i] = 0.0;
00070 theA[i] = 1;
00071 }
00072 fPion = 0;
00073 fGlauber = 0;
00074 fHadron = 0;
00075
00076 fSAID = 0;
00077 particle = p;
00078 theProton= G4Proton::Proton();
00079 isPiplus = false;
00080 isInitialized = false;
00081 }
00082
00083
00084 G4BGGPionInelasticXS::~G4BGGPionInelasticXS()
00085 {
00086 delete fGlauber;
00087 delete fPion;
00088 delete fHadron;
00089 delete fSAID;
00090 }
00091
00092
00093
00094 G4bool
00095 G4BGGPionInelasticXS::IsElementApplicable(const G4DynamicParticle*, G4int,
00096 const G4Material*)
00097 {
00098 return true;
00099 }
00100
00101
00102
00103 G4bool G4BGGPionInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
00104 G4int Z, G4int A,
00105 const G4Element*,
00106 const G4Material*)
00107 {
00108 return (1 == Z && 2 >= A);
00109 }
00110
00111
00112
00113 G4double
00114 G4BGGPionInelasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
00115 G4int ZZ, const G4Material*)
00116 {
00117
00118
00119 G4double cross = 0.0;
00120 G4double ekin = dp->GetKineticEnergy();
00121 G4int Z = ZZ;
00122 if(1 == Z) {
00123 cross = 1.0115*GetIsoCrossSection(dp,1,1);
00124 } else {
00125 if(Z > 92) { Z = 92; }
00126
00127 if(ekin <= fLowEnergy && !isPiplus) {
00128 cross = theCoulombFac[Z];
00129 } else if(ekin <= 2*MeV && isPiplus) {
00130 cross = theCoulombFac[Z]*CoulombFactor(ekin, Z);
00131 } else if(ekin > fGlauberEnergy) {
00132 cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
00133 } else {
00134 cross = fPion->GetInelasticCrossSection(dp, Z, theA[Z]);
00135 }
00136 }
00137 if(verboseLevel > 1) {
00138 G4cout << "G4BGGPionInelasticXS::GetCrossSection for "
00139 << dp->GetDefinition()->GetParticleName()
00140 << " Ekin(GeV)= " << dp->GetKineticEnergy()
00141 << " in nucleus Z= " << Z << " A= " << theA[Z]
00142 << " XS(b)= " << cross/barn
00143 << G4endl;
00144 }
00145 return cross;
00146 }
00147
00148
00149
00150 G4double
00151 G4BGGPionInelasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
00152 G4int Z, G4int A,
00153 const G4Isotope*,
00154 const G4Element*,
00155 const G4Material*)
00156 {
00157
00158
00159 G4double cross = 0.0;
00160 G4double ekin = dp->GetKineticEnergy();
00161
00162 if(ekin <= fSAIDHighEnergyLimit) {
00163 cross = fSAID->GetInelasticIsotopeCrossSection(particle, ekin, 1, 1);
00164 } else {
00165 fHadron->GetHadronNucleonXscPDG(dp, theProton);
00166 cross = (theCoulombFac[1]/ekin + 1)*fHadron->GetInelasticHadronNucleonXsc();
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 cross *= A;
00189
00190 if(verboseLevel > 1) {
00191 G4cout << "G4BGGPionInelasticXS::GetCrossSection for "
00192 << dp->GetDefinition()->GetParticleName()
00193 << " Ekin(GeV)= " << dp->GetKineticEnergy()
00194 << " in nucleus Z= " << Z << " A= " << A
00195 << " XS(b)= " << cross/barn
00196 << G4endl;
00197 }
00198 return cross;
00199 }
00200
00201
00202
00203 void G4BGGPionInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
00204 {
00205 if(&p == G4PionPlus::PionPlus() || &p == G4PionMinus::PionMinus()) {
00206 particle = &p;
00207 } else {
00208 G4cout << "### G4BGGPionInelasticXS WARNING: is not applicable to "
00209 << p.GetParticleName()
00210 << G4endl;
00211 throw G4HadronicException(__FILE__, __LINE__,
00212 "G4BGGPionInelasticXS::BuildPhysicsTable is used for wrong particle");
00213 return;
00214 }
00215
00216 if(isInitialized) { return; }
00217 isInitialized = true;
00218
00219 fPion = new G4UPiNuclearCrossSection();
00220 fGlauber = new G4GlauberGribovCrossSection();
00221 fHadron = new G4HadronNucleonXsc();
00222
00223 fSAID = new G4ComponentSAIDTotalXS();
00224
00225 fPion->BuildPhysicsTable(*particle);
00226 fGlauber->BuildPhysicsTable(*particle);
00227 if(particle == G4PionPlus::PionPlus()) { isPiplus = true; }
00228
00229 G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
00230 G4ThreeVector mom(0.0,0.0,1.0);
00231 G4DynamicParticle dp(part, mom, fGlauberEnergy);
00232
00233 G4NistManager* nist = G4NistManager::Instance();
00234
00235 G4double csup, csdn;
00236 G4int A;
00237
00238 if(verboseLevel > 0) {
00239 G4cout << "### G4BGGPionInelasticXS::Initialise for "
00240 << particle->GetParticleName()
00241 << " isPiplus: " << isPiplus
00242 << G4endl;
00243 }
00244
00245 for(G4int iz=2; iz<93; iz++) {
00246
00247 A = G4lrint(nist->GetAtomicMassAmu(iz));
00248 theA[iz] = A;
00249
00250 csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
00251 csdn = fPion->GetInelasticCrossSection(&dp, iz, theA[iz]);
00252
00253 theGlauberFac[iz] = csdn/csup;
00254 if(verboseLevel > 0) {
00255 G4cout << "Z= " << iz << " A= " << A
00256 << " factor= " << theGlauberFac[iz] << G4endl;
00257 }
00258 }
00259 dp.SetKineticEnergy(fSAIDHighEnergyLimit);
00260 fHadron->GetHadronNucleonXscPDG(&dp, theProton);
00261 theCoulombFac[1] = fSAIDHighEnergyLimit*
00262 (fSAID->GetInelasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)
00263 /fHadron->GetInelasticHadronNucleonXsc() - 1);
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279 if(isPiplus) {
00280 dp.SetKineticEnergy(2*MeV);
00281 for(G4int iz=2; iz<93; iz++) {
00282 theCoulombFac[iz] = fPion->GetInelasticCrossSection(&dp, iz, theA[iz])
00283 /CoulombFactor(2*MeV,iz);
00284 }
00285
00286 } else {
00287 dp.SetKineticEnergy(fLowEnergy);
00288
00289
00290 for(G4int iz=2; iz<93; iz++) {
00291 theCoulombFac[iz] = fPion->GetInelasticCrossSection(&dp, iz, theA[iz]);
00292 }
00293 }
00294 }
00295
00296
00297
00298 G4double G4BGGPionInelasticXS::CoulombFactor(G4double kinEnergy, G4int Z)
00299 {
00300 G4int A = theA[Z];
00301 G4double res= 0.0;
00302 if(kinEnergy <= DBL_MIN) { return res; }
00303 else if(A < 2) { return kinEnergy*kinEnergy; }
00304
00305 G4double elog = std::log10(6.7*kinEnergy/GeV);
00306 G4double aa = A;
00307
00308
00309 G4double ff1 = 0.70 - 0.002*aa;
00310 G4double ff2 = 1.00 + 1/aa;
00311 G4double ff3 = 0.8 + 18/aa - 0.002*aa;
00312 res = 1.0 + ff3*(1.0 - (1.0/(1+std::exp(-8*ff1*(elog + 1.37*ff2)))));
00313
00314 ff1 = 1. - 1./aa - 0.001*aa;
00315 ff2 = 1.17 - 2.7/aa-0.0014*aa;
00316 res /= (1 + std::exp(-8.*ff1*(elog + 2*ff2)));
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327 return res;
00328 }
00329
00330
00331
00332 void
00333 G4BGGPionInelasticXS::CrossSectionDescription(std::ostream& outFile) const
00334 {
00335 outFile << "The Barashenkov-Glauber-Gribov cross section handles inelastic\n"
00336 << "pion scattering from nuclei at all energies. The Barashenkov\n"
00337 << "parameterization is used below 91 GeV and the Glauber-Gribov\n"
00338 << "parameterization is used above 91 GeV.\n";
00339 }
00340
00341