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 // By JPW, working, but to be cleaned up. @@@ 00027 // G.Folger, 29-sept-2006: extend to 1TeV, using a constant above 20GeV 00028 // 22 Dec 2006 - DHW added isotope dependence 00029 // G.Folger, 25-Nov-2009: extend to 100TeV, using a constant above 20GeV 00030 // V.Ivanchenko, 18-Aug-2011: migration to new design and cleanup; 00031 // make it applicable for Z>1 00032 // 00033 00034 #include "G4ProtonInelasticCrossSection.hh" 00035 #include "G4SystemOfUnits.hh" 00036 #include "G4DynamicParticle.hh" 00037 #include "G4Proton.hh" 00038 #include "G4HadTmpUtil.hh" 00039 #include "G4NistManager.hh" 00040 00041 G4ProtonInelasticCrossSection::G4ProtonInelasticCrossSection() 00042 : G4VCrossSectionDataSet("Axen-Wellisch"),thEnergy(19.8*CLHEP::GeV) 00043 { 00044 nist = G4NistManager::Instance(); 00045 theProton = G4Proton::Proton(); 00046 } 00047 00048 G4ProtonInelasticCrossSection::~G4ProtonInelasticCrossSection() 00049 {} 00050 00051 G4bool 00052 G4ProtonInelasticCrossSection::IsElementApplicable( 00053 const G4DynamicParticle* aPart, 00054 G4int Z, const G4Material*) 00055 { 00056 return ((1 < Z) && (aPart->GetDefinition() == theProton)); 00057 } 00058 00059 G4double G4ProtonInelasticCrossSection::GetElementCrossSection( 00060 const G4DynamicParticle* aPart, 00061 G4int Z, const G4Material*) 00062 { 00063 return GetProtonCrossSection(aPart->GetKineticEnergy(), Z); 00064 } 00065 00066 G4double G4ProtonInelasticCrossSection::GetProtonCrossSection( 00067 G4double kineticEnergy, G4int Z) 00068 { 00069 if(kineticEnergy <= 0.0) { return 0.0; } 00070 00071 // constant cross section above ~20GeV 00072 if (kineticEnergy > thEnergy) { kineticEnergy = thEnergy; } 00073 00074 G4double a = nist->GetAtomicMassAmu(Z); 00075 G4double a13 = std::pow(a,-0.3333333333); 00076 G4int nOfNeutrons = G4lrint(a) - Z; 00077 kineticEnergy /=GeV; 00078 G4double alog10E = std::log10(kineticEnergy); 00079 00080 const G4double nuleonRadius=1.36E-15; 00081 const G4double fac=CLHEP::pi*nuleonRadius*nuleonRadius; 00082 00083 G4double b0 = 2.247-0.915*(1 - a13); 00084 G4double fac1 = b0*(1 - a13); 00085 G4double fac2 = 1.; 00086 if(nOfNeutrons > 1) { fac2=std::log((G4double(nOfNeutrons))); } 00087 G4double crossSection = 1.0E31*fac*fac2*(1. + 1./a13 - fac1); 00088 00089 // high energy correction 00090 crossSection *= (1 - 0.15*std::exp(-kineticEnergy))/(1.0 - 0.0007*a); 00091 00092 // first try on low energies: rise 00093 G4double ff1= 0.70-0.002*a; // slope of the drop at medium energies. 00094 G4double ff2= 1.00+1/a; // start of the slope. 00095 G4double ff3= 0.8+18/a-0.002*a; // stephight 00096 00097 G4double ff4= 1.0 - (1.0/(1+std::exp(-8*ff1*(alog10E + 1.37*ff2)))); 00098 00099 crossSection *= (1 + ff3*ff4); 00100 00101 // low energy return to zero 00102 00103 ff1=1. - 1./a - 0.001*a; // slope of the rise 00104 ff2=1.17 - 2.7/a - 0.0014*a; // start of the rise 00105 00106 ff4=-8.*ff1*(alog10E + 2.0*ff2); 00107 00108 crossSection *= millibarn/(1. + std::exp(ff4)); 00109 return crossSection; 00110 }