Geant4: G4BGGNucleonInelasticXS.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4BGGNucleonInelasticXS
00034 //
00035 // Author:        Vladimir Ivanchenko
00036 //
00037 // Creation date: 13.03.2007
00038 // Modifications:
00039 //
00040 //
00041 // -------------------------------------------------------------------
00042 //
00043 
00044 #include "G4BGGNucleonInelasticXS.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4GlauberGribovCrossSection.hh"
00047 #include "G4NucleonNuclearCrossSection.hh"
00048 #include "G4HadronNucleonXsc.hh"
00049 //#include "G4HadronInelasticDataSet.hh"
00050 #include "G4ComponentSAIDTotalXS.hh"
00051 #include "G4Proton.hh"
00052 #include "G4Neutron.hh"
00053 #include "G4NistManager.hh"
00054 #include "G4Material.hh"
00055 #include "G4Element.hh"
00056 #include "G4Isotope.hh"
00057 
00058 #include "G4CrossSectionDataSetRegistry.hh"
00059 
00060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00061 
00062 G4BGGNucleonInelasticXS::G4BGGNucleonInelasticXS(const G4ParticleDefinition* p)
00063  : G4VCrossSectionDataSet("Barashenkov-Glauber")
00064 {
00065   verboseLevel = 0;
00066   fGlauberEnergy = 91.*GeV;
00067   fLowEnergy = 14.*MeV;
00068   fHighEnergy = 5.*GeV;
00069   fSAIDHighEnergyLimit = 1.3*GeV;
00070   for (G4int i = 0; i < 93; ++i) {
00071     theGlauberFac[i] = 0.0;
00072     theCoulombFac[i] = 0.0;
00073     theA[i] = 1;
00074   }
00075   fNucleon = 0;
00076   fGlauber = 0;
00077   fHadron  = 0;
00078   //  fGHEISHA = 0;
00079   fSAID    = 0;
00080   particle = p;
00081   theProton= G4Proton::Proton();
00082   isProton = false;
00083   isInitialized = false;
00084 }
00085 
00086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00087 
00088 G4BGGNucleonInelasticXS::~G4BGGNucleonInelasticXS()
00089 {
00090   delete fHadron;
00091   delete fSAID;
00092 }
00093 
00094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00095 
00096 G4bool G4BGGNucleonInelasticXS::IsElementApplicable(const G4DynamicParticle*, 
00097                                                     G4int, const G4Material*)
00098 {
00099   return true;
00100 }
00101 
00102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00103 
00104 G4bool G4BGGNucleonInelasticXS::IsIsoApplicable(const G4DynamicParticle*, 
00105                                                 G4int Z, G4int A,  
00106                                                 const G4Element*,
00107                                                 const G4Material*)
00108 {
00109   return (1 == Z && 2 >= A);
00110 }
00111 
00112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00113 
00114 G4double
00115 G4BGGNucleonInelasticXS::GetElementCrossSection(const G4DynamicParticle* dp,
00116                                                 G4int ZZ, const G4Material*)
00117 {
00118   // this method should be called only for Z > 1
00119 
00120   G4double cross = 0.0;
00121   G4double ekin = dp->GetKineticEnergy();
00122   G4int Z = ZZ;
00123   if(1 == Z) {
00124     cross = 1.0115*GetIsoCrossSection(dp,1,1);
00125   } else if(2 == Z) {
00126     if(ekin > fGlauberEnergy) {
00127       cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
00128     } else {
00129       cross = fNucleon->GetElementCrossSection(dp, Z);
00130     }
00131 
00132   } else {
00133     if(Z > 92) { Z = 92; }
00134 
00135     if(ekin <= fLowEnergy) {
00136       cross = theCoulombFac[Z]*CoulombFactor(ekin, Z);
00137     } else if(ekin > fGlauberEnergy) {
00138       cross = theGlauberFac[Z]*fGlauber->GetInelasticGlauberGribov(dp, Z, theA[Z]);
00139     } else {
00140       cross = fNucleon->GetElementCrossSection(dp, Z);
00141     }
00142   }
00143 
00144   if(verboseLevel > 1) {
00145     G4cout << "G4BGGNucleonInelasticXS::GetCrossSection  for "
00146            << dp->GetDefinition()->GetParticleName()
00147            << "  Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
00148            << " in nucleus Z= " << Z << "  A= " << theA[Z]
00149            << " XS(b)= " << cross/barn 
00150            << G4endl;
00151   }
00152   return cross;
00153 }
00154 
00155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00156 
00157 G4double
00158 G4BGGNucleonInelasticXS::GetIsoCrossSection(const G4DynamicParticle* dp,
00159                                             G4int Z, G4int A,
00160                                             const G4Isotope*,
00161                                             const G4Element*,
00162                                             const G4Material*)
00163 {
00164   // this method should be called only for Z = 1
00165 
00166   G4double cross = 0.0;
00167   G4double ekin = dp->GetKineticEnergy();
00168 
00169   if(ekin <= fSAIDHighEnergyLimit) {
00170     cross = fSAID->GetInelasticIsotopeCrossSection(particle, ekin, 1, 1);
00171   } else if(ekin < fHighEnergy) {
00172     fHadron->GetHadronNucleonXscNS(dp, theProton);
00173     cross = (theCoulombFac[0]/ekin + 1)*fHadron->GetInelasticHadronNucleonXsc();
00174   } else {
00175     fHadron->GetHadronNucleonXscPDG(dp, theProton);
00176     cross = (theCoulombFac[1]/ekin + 1)*fHadron->GetInelasticHadronNucleonXsc();
00177   } 
00178   cross *= A; 
00179 
00180   if(verboseLevel > 1) {
00181     G4cout << "G4BGGNucleonInelasticXS::GetCrossSection  for "
00182            << dp->GetDefinition()->GetParticleName()
00183            << "  Ekin(GeV)= " << dp->GetKineticEnergy()/CLHEP::GeV
00184            << " in nucleus Z= " << Z << "  A= " << A
00185            << " XS(b)= " << cross/barn 
00186            << G4endl;
00187   }
00188   return cross;
00189 }
00190 
00191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00192 
00193 void G4BGGNucleonInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
00194 {
00195   if(&p == theProton || &p == G4Neutron::Neutron()) {
00196     particle = &p;
00197   } else {
00198     G4cout << "### G4BGGNucleonInelasticXS WARNING: is not applicable to " 
00199            << p.GetParticleName()
00200            << G4endl;
00201     throw G4HadronicException(__FILE__, __LINE__,
00202           "G4BGGNucleonElasticXS::BuildPhysicsTable is used for wrong particle");
00203     return;
00204   }
00205 
00206   if(isInitialized) { return; }
00207   isInitialized = true;
00208 
00209     fNucleon = (G4NucleonNuclearCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NucleonNuclearCrossSection::Default_Name());
00210     fGlauber = (G4GlauberGribovCrossSection*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4GlauberGribovCrossSection::Default_Name());
00211 
00212     fHadron  = new G4HadronNucleonXsc();
00213   //fGHEISHA = new G4HadronInelasticDataSet();
00214   fSAID    = new G4ComponentSAIDTotalXS();
00215 
00216   fNucleon->BuildPhysicsTable(*particle);
00217   fGlauber->BuildPhysicsTable(*particle);
00218 
00219   if(particle == theProton) { 
00220     isProton = true; 
00221     fSAIDHighEnergyLimit = 2*GeV;
00222     fHighEnergy = 2*GeV;
00223   }
00224 
00225   G4ParticleDefinition* part = const_cast<G4ParticleDefinition*>(particle);
00226   G4ThreeVector mom(0.0,0.0,1.0);
00227   G4DynamicParticle dp(part, mom, fGlauberEnergy);
00228 
00229   G4NistManager* nist = G4NistManager::Instance();
00230   G4int A;
00231 
00232   G4double csup, csdn;
00233 
00234   if(verboseLevel > 0) {
00235     G4cout << "### G4BGGNucleonInelasticXS::Initialise for "
00236            << particle->GetParticleName() << G4endl;
00237   }
00238   for(G4int iz=2; iz<93; iz++) {
00239 
00240     A = G4lrint(nist->GetAtomicMassAmu(iz));
00241     theA[iz] = A;
00242 
00243     csup = fGlauber->GetInelasticGlauberGribov(&dp, iz, A);
00244     csdn = fNucleon->GetElementCrossSection(&dp, iz);
00245 
00246     theGlauberFac[iz] = csdn/csup;
00247     if(verboseLevel > 0) {
00248       G4cout << "Z= " << iz <<  "  A= " << A 
00249              << " GlauberFactor= " << theGlauberFac[iz] << G4endl; 
00250     }
00251   }
00252   //const G4Material* mat = 0;
00253 
00254   dp.SetKineticEnergy(fSAIDHighEnergyLimit);
00255   fHadron->GetHadronNucleonXscNS(&dp, theProton);
00256   theCoulombFac[0] = fSAIDHighEnergyLimit*
00257     (fSAID->GetInelasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)
00258      /fHadron->GetInelasticHadronNucleonXsc() - 1);
00259   
00260   //G4cout << "Z=1 E(GeV)= " << fSAIDHighEnergyLimit/GeV
00261   //     << "  xsNS(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;  
00262   fHadron->GetHadronNucleonXscPDG(&dp, theProton);
00263   //G4cout << "  xsPDG(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;
00264   //G4cout << "  xsSAID(b)= " << fSAID->GetInelasticIsotopeCrossSection(particle,fSAIDHighEnergyLimit,1,1)/barn << G4endl;
00265 
00266   dp.SetKineticEnergy(fHighEnergy);
00267   fHadron->GetHadronNucleonXscPDG(&dp, theProton);
00268   G4double x = fHadron->GetInelasticHadronNucleonXsc();
00269 
00270   //G4cout << "Z=1 E(GeV)= " << fHighEnergy/GeV
00271   //     << "  xsPDG(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn;  
00272 
00273   fHadron->GetHadronNucleonXscNS(&dp, theProton);
00274   theCoulombFac[1] = fHighEnergy*((theCoulombFac[0]/fHighEnergy + 1)
00275                                   *fHadron->GetInelasticHadronNucleonXsc()/x - 1);
00276 
00277   fHadron->GetHadronNucleonXscNS(&dp, theProton);
00278   //G4cout << "  xsNS(b)= " << fHadron->GetInelasticHadronNucleonXsc()/barn << G4endl;
00279 
00280   if(verboseLevel > 0) {
00281     G4cout << "Z=1   A=1" << " CoulombFactor[0]= " << theCoulombFac[0]
00282            << " CoulombFactor[1]= " << theCoulombFac[1] << G4endl; 
00283   }
00284   theCoulombFac[2] = 1.0;
00285      
00286   dp.SetKineticEnergy(fLowEnergy);
00287   for(G4int iz=3; iz<93; iz++) {
00288     theCoulombFac[iz] = 
00289       fNucleon->GetElementCrossSection(&dp, iz)/CoulombFactor(fLowEnergy, iz);
00290 
00291     if(verboseLevel > 0) {
00292       G4cout << "Z= " << iz <<  "  A= " << theA[iz] 
00293              << " CoulombFactor= " << theCoulombFac[iz] << G4endl; 
00294     }
00295   }
00296 }
00297 
00298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00299 
00300 G4double G4BGGNucleonInelasticXS::CoulombFactor(G4double kinEnergy, G4int Z)
00301 {
00302   G4double res= 0.0;
00303   if(kinEnergy <= 0.0) { return res; }
00304   else if (Z <= 1) { return kinEnergy*kinEnergy; }
00305   
00306   G4double elog = std::log10(kinEnergy/GeV);
00307   G4double aa = theA[Z];
00308 
00309   // from G4ProtonInelasticCrossSection
00310   if(isProton) {
00311 
00312     G4double ff1 = 0.70 - 0.002*aa;           // slope of the drop at medium energies.
00313     G4double ff2 = 1.00 + 1/aa;               // start of the slope.
00314     G4double ff3 = 0.8  + 18/aa - 0.002*aa;   // stephight
00315     res = 1.0 + ff3*(1.0 - (1.0/(1+std::exp(-8*ff1*(elog + 1.37*ff2)))));
00316 
00317     ff1 = 1.   - 1./aa - 0.001*aa; // slope of the rise
00318     ff2 = 1.17 - 2.7/aa-0.0014*aa; // start of the rise
00319     res /= (1 + std::exp(-8.*ff1*(elog + 2*ff2)));
00320 
00321   } else {
00322 
00323     // from G4NeutronInelasticCrossSection
00324     G4double p3 = 0.6 + 13./aa - 0.0005*aa;
00325     G4double p4 = 7.2449 - 0.018242*aa;
00326     G4double p5 = 1.36 + 1.8/aa + 0.0005*aa;
00327     G4double p6 = 1. + 200./aa + 0.02*aa;
00328     G4double p7 = 3.0 - (aa-70.)*(aa-200.)/11000.;
00329 
00330     G4double firstexp  = std::exp(-p4*(elog + p5));
00331     G4double secondexp = std::exp(-p6*(elog + p7));
00332 
00333     res = (1.+p3*firstexp/(1. + firstexp))/(1. + secondexp);
00334 
00335   }
00336   return res;  
00337 }
00338 
00339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00340 
00341 void G4BGGNucleonInelasticXS::CrossSectionDescription(std::ostream& outFile) const
00342 {
00343   outFile << "The Barashenkov-Glauber-Gribov cross section calculates inelastic\n"
00344           << "scattering of protons and neutrons from nuclei using the\n"
00345           << "Barashenkov parameterization below 91 GeV and the Glauber-Gribov\n"
00346           << "parameterization above 91 GeV.  It uses the G4HadronNucleonXsc\n"
00347           << "cross section component for hydrogen targets, and the\n"
00348           << "G4GlauberGribovCrossSection component for other targets.\n";
00349 }
00350 
00351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......