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 #include "G4IonsKoxCrossSection.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4NucleiProperties.hh"
00038 #include "G4HadTmpUtil.hh"
00039 #include "G4NistManager.hh"
00040
00041 G4IonsKoxCrossSection::G4IonsKoxCrossSection()
00042 : G4VCrossSectionDataSet("IonsKox"), lowerLimit ( 10*MeV ),
00043 r0 ( 1.1*fermi ), rc ( 1.3*fermi )
00044 {}
00045
00046 G4IonsKoxCrossSection::~G4IonsKoxCrossSection()
00047 {}
00048
00049 void
00050 G4IonsKoxCrossSection::CrossSectionDescription(std::ostream& outFile) const
00051 {
00052 outFile << "G4IonsKoxCrossSection calculates the total reaction cross\n"
00053 << "section for nucleus-nucleus scattering using the Kox\n"
00054 << "parameterization. It is valid for projectiles and targets\n"
00055 << "of all Z, at projectile energies up to 10 GeV/n. If the\n"
00056 << "projectile energy is less than 10 MeV/n, a zero cross section\n"
00057 << "is returned.\n";
00058 }
00059
00060 G4bool G4IonsKoxCrossSection::IsElementApplicable(const G4DynamicParticle* aDP,
00061 G4int, const G4Material*)
00062 {
00063 return (1 <= aDP->GetDefinition()->GetBaryonNumber());
00064 }
00065
00066 G4double
00067 G4IonsKoxCrossSection::GetElementCrossSection(
00068 const G4DynamicParticle* aParticle, G4int ZZ, const G4Material*)
00069 {
00070 G4double xsection = 0.0;
00071
00072 G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
00073 G4int Zp = G4int(aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5);
00074 G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
00075
00076
00077
00078
00079 G4int At = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(ZZ));
00080 G4int Zt = ZZ;
00081
00082 G4double one_third = 1.0 / 3.0;
00083
00084 G4double cubicrAt = std::pow ( G4double(At) , G4double(one_third) );
00085 G4double cubicrAp = std::pow ( G4double(Ap) , G4double(one_third) );
00086
00087
00088 G4double Bc = Zt * Zp / ( (rc/fermi) * (cubicrAp+cubicrAt) );
00089
00090 G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
00091 G4double proj_mass = aParticle->GetMass();
00092 G4double proj_momentum = aParticle->GetMomentum().mag();
00093
00094 G4double Ecm = calEcm ( proj_mass , targ_mass , proj_momentum );
00095 if( Ecm <= Bc) return xsection;
00096
00097 G4double Rvol = r0 * ( cubicrAp + cubicrAt );
00098
00099
00100 G4double c = calCeValue ( ke_per_N / MeV );
00101
00102 G4double a = 1.85;
00103 G4double Rsurf = r0 * (a*cubicrAp * cubicrAt/(cubicrAp + cubicrAt) - c);
00104 G4double D = 5.0 * ( At - 2 * Zt ) * Zp / ( Ap * At );
00105 Rsurf = Rsurf + D * fermi;
00106
00107 G4double Rint = Rvol + Rsurf;
00108 xsection = pi * Rint * Rint * ( 1 - Bc / ( Ecm / MeV ) );
00109
00110 return xsection;
00111 }
00112
00113 G4double
00114 G4IonsKoxCrossSection::calEcm(G4double mp, G4double mt, G4double Plab)
00115 {
00116 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
00117 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
00118 G4double Pcm = Plab * mt / Ecm;
00119 G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
00120 return KEcm;
00121 }
00122
00123
00124 G4double G4IonsKoxCrossSection::calCeValue(const G4double ke)
00125 {
00126
00127
00128
00129
00130
00131
00132 G4double Ce;
00133 G4double log10_ke = std::log10 ( ke );
00134 if (log10_ke > 1.5)
00135 {
00136 Ce = - 10.0 / std::pow ( G4double(log10_ke) , G4double(5) ) + 2.0;
00137 }
00138 else
00139 {
00140 Ce = (-10.0/std::pow(G4double(1.5), G4double(5) ) + 2.0) /
00141 std::pow(G4double(1.5), G4double(3)) * std::pow(G4double(log10_ke), G4double(3) );
00142
00143 }
00144 return Ce;
00145 }
00146