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 #include "G4NeutronInelasticCrossSection.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4DynamicParticle.hh"
00036 #include "G4HadTmpUtil.hh"
00037 #include "G4NistManager.hh"
00038 #include "G4Pow.hh"
00039
00040 G4NeutronInelasticCrossSection::G4NeutronInelasticCrossSection()
00041 : G4VCrossSectionDataSet("Wellisch-Laidlaw"),
00042 minEnergy(19.9*MeV), maxEnergy(19.9*GeV)
00043 {}
00044
00045 G4NeutronInelasticCrossSection::~G4NeutronInelasticCrossSection()
00046 {}
00047
00048 void
00049 G4NeutronInelasticCrossSection::CrossSectionDescription(std::ostream& outFile) const
00050 {
00051 outFile << "G4NeutronInelasticCrossSection calculates the inelastic neutron\n"
00052 << "scattering cross section for nuclei using the Wellisch-Laidlaw\n"
00053 << "parameterization between 19.9 MeV and 19.9 GeV. Above 19.9 GeV\n"
00054 << "the cross section is assumed to be constant.\n";
00055 }
00056
00057 G4bool
00058 G4NeutronInelasticCrossSection::IsElementApplicable(
00059 const G4DynamicParticle* part, G4int Z, const G4Material*)
00060 {
00061 G4double e = part->GetKineticEnergy();
00062 return (1 < Z && e > minEnergy);
00063 }
00064
00065 G4double G4NeutronInelasticCrossSection::
00066 GetElementCrossSection(const G4DynamicParticle* aPart,
00067 G4int Z, const G4Material*)
00068 {
00069 G4int A = G4int(G4NistManager::Instance()->GetAtomicMassAmu(Z));
00070 return GetCrossSection(aPart->GetKineticEnergy(), Z, A);
00071 }
00072
00073 G4double
00074 G4NeutronInelasticCrossSection::GetCrossSection(G4double anEnergy,
00075 G4int Z, G4int A)
00076 {
00077 if(anEnergy > maxEnergy) { anEnergy = maxEnergy; }
00078 G4double cross_section = 0.0;
00079 if(anEnergy < keV) { return cross_section; }
00080
00081 G4Pow* g4pow = G4Pow::GetInstance();
00082 G4double A13 = g4pow->Z13(A);
00083
00084 G4double elog = std::log10(anEnergy/MeV);
00085 G4int nOfNeutrons = A - Z;
00086 G4double atomicNumber = G4double(A);
00087 const G4double p1=1.3773;
00088 G4double p2 = 1. + 10./atomicNumber - 0.0006*atomicNumber;
00089 G4double p3 = 0.6+ 13./atomicNumber - 0.0005*atomicNumber;
00090 G4double p4 = 7.2449 - 0.018242*atomicNumber;
00091 G4double p5 = 1.64 - 1.8/atomicNumber - 0.0005*atomicNumber;
00092 G4double p6 = 1. + 200./atomicNumber + 0.02*atomicNumber;
00093 G4double p7 = (atomicNumber-70.)*(atomicNumber-200.)/11000.;
00094
00095 G4double logN = g4pow->logZ(nOfNeutrons);
00096 G4double part1 = pi*p1*p1*logN;
00097 G4double part2 = 1.+ A13 - p2*(1.-1./A13);
00098
00099 G4double firstexp = -p4*(elog-p5);
00100 G4double first = 1. + std::exp(firstexp);
00101 G4double corr = 1. + p3*(1.-1./first);
00102
00103 G4double secondexp= -p6*(elog-p7);
00104 G4double secondv = 1.+std::exp(secondexp);
00105 G4double corr2 = 1./secondv;
00106
00107 G4double xsec = corr*corr2*part1*part2*10.*millibarn;
00108 if(xsec < 0.0) { xsec = 0.0; }
00109 return xsec;
00110 }