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