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 "G4TripathiCrossSection.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4DynamicParticle.hh"
00037 #include "G4ParticleTable.hh"
00038 #include "G4IonTable.hh"
00039 #include "G4HadTmpUtil.hh"
00040 #include "G4Proton.hh"
00041 #include "G4NistManager.hh"
00042
00043 G4TripathiCrossSection::G4TripathiCrossSection()
00044 : G4VCrossSectionDataSet("Tripathi")
00045 {}
00046
00047 G4TripathiCrossSection::~G4TripathiCrossSection()
00048 {}
00049
00050 G4bool
00051 G4TripathiCrossSection::IsElementApplicable(const G4DynamicParticle* aPart,
00052 G4int, const G4Material*)
00053 {
00054 G4bool result = false;
00055 if ( (aPart->GetDefinition()->GetBaryonNumber()>2.5) &&
00056 ( aPart->GetKineticEnergy()/aPart->GetDefinition()->GetBaryonNumber()<1*GeV) ) {
00057 result = true;
00058 }
00059 return result;
00060 }
00061
00062 G4double G4TripathiCrossSection::
00063 GetElementCrossSection(const G4DynamicParticle* aPart, G4int ZZ,
00064 const G4Material*)
00065 {
00066 G4double result = 0.;
00067 G4double targetAtomicNumber = G4NistManager::Instance()->GetAtomicMassAmu(ZZ);
00068 G4double nTargetProtons = ZZ;
00069
00070 G4double kineticEnergy = aPart->GetKineticEnergy()/MeV;
00071 G4double nProjProtons = aPart->GetDefinition()->GetPDGCharge();
00072 G4double projectileAtomicNumber =
00073 aPart->GetDefinition()->GetBaryonNumber();
00074
00075 const G4double nuleonRadius=1.1E-15;
00076 const G4double myNuleonRadius=1.36E-15;
00077
00078
00079 G4double targetMass =
00080 G4ParticleTable::GetParticleTable()->GetIonTable()
00081 ->GetIonMass(G4lrint(nTargetProtons), G4lrint(targetAtomicNumber));
00082 G4LorentzVector pTarget(0,0,0,targetMass);
00083 G4LorentzVector pProjectile(aPart->Get4Momentum());
00084 pTarget = pTarget+pProjectile;
00085 G4double E_cm = (pTarget.mag()-targetMass-pProjectile.m())/MeV;
00086 if(E_cm <= DBL_MIN) { return result; }
00087
00088 G4double r_rms_p = 0.6 * myNuleonRadius *
00089 std::pow(projectileAtomicNumber, 1./3.);
00090 G4double r_rms_t = 0.6 * myNuleonRadius *
00091 std::pow(targetAtomicNumber, 1./3.);
00092
00093
00094 G4double r_p = 1.29*r_rms_p/nuleonRadius ;
00095 G4double r_t = 1.29*r_rms_t/nuleonRadius;
00096
00097
00098 G4double Radius = r_p + r_t +
00099 1.2*(std::pow(targetAtomicNumber, 1./3.) +
00100 std::pow(projectileAtomicNumber, 1./3.))/std::pow(E_cm, 1./3.);
00101
00102
00103 G4double B = 1.44*nProjProtons*nTargetProtons/Radius;
00104 if(E_cm <= B) return result;
00105
00106 G4double Energy = kineticEnergy/projectileAtomicNumber;
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117 G4double D;
00118 if (nProjProtons==1 && projectileAtomicNumber==1)
00119 {
00120 D = 2.05;
00121 }
00122 else if (nProjProtons==2 && projectileAtomicNumber==4)
00123 {
00124 D = 2.77-(8.0E-3*targetAtomicNumber)+
00125 (1.8E-5*targetAtomicNumber*targetAtomicNumber)
00126 - 0.8/(1+std::exp((250.-Energy)/75.));
00127 }
00128 else
00129 {
00130
00131
00132
00133
00134
00135
00136 D = 1.75;
00137 }
00138
00139 G4double C_E = D * (1-std::exp(-Energy/40.)) -
00140 0.292*std::exp(-Energy/792.)*std::cos(0.229*std::pow(Energy, 0.453));
00141
00142
00143 G4double S = std::pow(projectileAtomicNumber, 1./3.)*
00144 std::pow(targetAtomicNumber, 1./3.)/
00145 (std::pow(projectileAtomicNumber, 1./3.) +
00146 std::pow(targetAtomicNumber, 1./3.));
00147
00148
00149 G4double deltaE = 1.85*S + 0.16*S/std::pow(E_cm,1./3.) - C_E +
00150 0.91*(targetAtomicNumber-2.*nTargetProtons)*nProjProtons/
00151 (targetAtomicNumber*projectileAtomicNumber);
00152
00153
00154 result = pi * nuleonRadius*nuleonRadius *
00155 std::pow(( std::pow(targetAtomicNumber, 1./3.) +
00156 std::pow(projectileAtomicNumber, 1./3.) + deltaE),2.) *
00157 (1-B/E_cm);
00158
00159 if(result < 0.) { result = 0.; }
00160 return result*m2;
00161
00162 }
00163