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
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00066
00067 #include "G4TripathiLightCrossSection.hh"
00068 #include "G4PhysicalConstants.hh"
00069 #include "G4SystemOfUnits.hh"
00070 #include "G4DynamicParticle.hh"
00071 #include "G4WilsonRadius.hh"
00072 #include "G4NucleiProperties.hh"
00073 #include "G4HadTmpUtil.hh"
00074 #include "G4NistManager.hh"
00075 #include "G4Pow.hh"
00076
00077 G4TripathiLightCrossSection::G4TripathiLightCrossSection ()
00078 : G4VCrossSectionDataSet("TripathiLightIons")
00079 {
00080
00081
00082
00083
00084 theWilsonRadius = new G4WilsonRadius();
00085 r_0 = 1.1 * fermi;
00086
00087
00088
00089
00090
00091
00092 lowEnergyCheck = false;
00093 }
00095
00096 G4TripathiLightCrossSection::~G4TripathiLightCrossSection ()
00097 {
00098
00099
00100
00101 delete theWilsonRadius;
00102 }
00104
00105 G4bool
00106 G4TripathiLightCrossSection::IsElementApplicable(const G4DynamicParticle* theProjectile,
00107 G4int ZT, const G4Material*)
00108 {
00109 G4bool result = false;
00110 G4int AT = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(ZT));
00111 G4int ZP = G4lrint(theProjectile->GetDefinition()->GetPDGCharge()/eplus);
00112 G4int AP = theProjectile->GetDefinition()->GetBaryonNumber();
00113 if (theProjectile->GetKineticEnergy()/AP < 10.0*GeV &&
00114 ((AT==1 && ZT==1) || (AP==1 && ZP==1) ||
00115 (AT==1 && ZT==0) || (AP==1 && ZP==0) ||
00116 (AT==2 && ZT==1) || (AP==2 && ZP==1) ||
00117 (AT==3 && ZT==2) || (AP==3 && ZP==2) ||
00118 (AT==4 && ZT==2) || (AP==4 && ZP==2))) { result = true; }
00119 return result;
00120 }
00121
00123
00124 G4double
00125 G4TripathiLightCrossSection::GetElementCrossSection(const G4DynamicParticle* theProjectile,
00126 G4int ZT, const G4Material*)
00127 {
00128
00129 G4double result = 0.0;
00130
00131
00132
00133
00134 G4double xAT= G4NistManager::Instance()->GetAtomicMassAmu(ZT);
00135 G4int AT = G4lrint(xAT);
00136 G4double EA = theProjectile->GetKineticEnergy()/MeV;
00137 G4int AP = theProjectile->GetDefinition()->GetBaryonNumber();
00138 G4double xAP= G4double(AP);
00139 G4double ZP = G4lrint(theProjectile->GetDefinition()->GetPDGCharge()/eplus);
00140 G4double E = EA / xAP;
00141
00142 G4Pow* g4pow = G4Pow::GetInstance();
00143
00144 G4double AT13 = g4pow->Z13(AT);
00145 G4double AP13 = g4pow->Z13(AP);
00146
00147
00148
00149 G4double mT = G4NucleiProperties::GetNuclearMass(AT, ZT);
00150 G4LorentzVector pT(0.0, 0.0, 0.0, mT);
00151 G4LorentzVector pP(theProjectile->Get4Momentum());
00152 pT += pP;
00153 G4double E_cm = (pT.mag()-mT-pP.m())/MeV;
00154
00155
00156
00157
00158
00159
00160 if (E_cm <= 0.0) { return 0.; }
00161 if (E_cm <= (0.8 + 0.04*ZT)*xAP && !lowEnergyCheck) { return 0.; }
00162
00163 G4double E_cm13 = g4pow->A13(E_cm);
00164
00165
00166
00167
00168 G4double r_rms_p = theWilsonRadius->GetWilsonRMSRadius(xAP);
00169 G4double r_rms_t = theWilsonRadius->GetWilsonRMSRadius(xAT);
00170
00171 G4double r_p = 1.29*r_rms_p;
00172 G4double r_t = 1.29*r_rms_t;
00173
00174 G4double Radius = (r_p + r_t)/fermi + 1.2*(AT13 + AP13)/E_cm13;
00175
00176 G4double B = 1.44 * ZP * ZT / Radius;
00177
00178
00179
00180
00181 G4double T1 = 0.0;
00182 G4double D = 0.0;
00183 G4double G = 0.0;
00184
00185 if ((AT==1 && ZT==1) || (AP==1 && ZP==1)) {
00186 T1 = 23.0;
00187 D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
00188
00189 } else if ((AT==1 && ZT==0) || (AP==1 && ZP==0)) {
00190 T1 = 18.0;
00191 D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
00192
00193 } else if ((AT==2 && ZT==1) || (AP==2 && ZP==1)) {
00194 T1 = 23.0;
00195 D = 1.65 + 0.1/(1+std::exp((500.0-E)/200.0));
00196
00197 } else if ((AT==3 && ZT==2) || (AP==3 && ZP==2)) {
00198 T1 = 40.0;
00199 D = 1.55;
00200
00201 } else if (AP==4 && ZP==2) {
00202 if (AT==4 && ZT==2) {T1 = 40.0; G = 300.0;}
00203 else if (ZT==4) {T1 = 25.0; G = 300.0;}
00204 else if (ZT==7) {T1 = 40.0; G = 500.0;}
00205 else if (ZT==13) {T1 = 25.0; G = 300.0;}
00206 else if (ZT==26) {T1 = 40.0; G = 300.0;}
00207 else {T1 = 40.0; G = 75.0;}
00208 D = 2.77 - 8.0E-3*AT + 1.8E-5*AT*AT-0.8/(1.0+std::exp((250.0-E)/G));
00209 }
00210 else if (AT==4 && ZT==2) {
00211 if (AP==4 && ZP==2) {T1 = 40.0; G = 300.0;}
00212 else if (ZP==4) {T1 = 25.0; G = 300.0;}
00213 else if (ZP==7) {T1 = 40.0; G = 500.0;}
00214 else if (ZP==13) {T1 = 25.0; G = 300.0;}
00215 else if (ZP==26) {T1 = 40.0; G = 300.0;}
00216 else {T1 = 40.0; G = 75.0;}
00217 D = 2.77 - 8.0E-3*AP + 1.8E-5*AP*AP-0.8/(1.0+std::exp((250.0-E)/G));
00218 }
00219
00220
00221
00222
00223
00224 G4double C_E = D*(1.0-std::exp(-E/T1)) -
00225 0.292*std::exp(-E/792.0)*std::cos(0.229*std::pow(E,0.453));
00226
00227 G4double S = AP13*AT13/(AP13 + AT13);
00228
00229 G4double deltaE = 0.0;
00230 G4double X1 = 0.0;
00231 if (AT >= AP)
00232 {
00233 deltaE = 1.85*S + 0.16*S/E_cm13 - C_E + 0.91*(AT-2*ZT)*ZP/(xAT*xAP);
00234 X1 = 2.83 - 3.1E-2*AT + 1.7E-4*AT*AT;
00235 }
00236 else
00237 {
00238 deltaE = 1.85*S + 0.16*S/E_cm13 - C_E + 0.91*(AP-2*ZP)*ZT/(xAT*xAP);
00239 X1 = 2.83 - 3.1E-2*AP + 1.7E-4*AP*AP;
00240 }
00241 G4double S_L = 1.2 + 1.6*(1.0-std::exp(-E/15.0));
00242
00243 G4double X_m = 1.0 - X1*std::exp(-E/(X1*S_L));
00244
00245
00246
00247
00248
00249
00250 G4double R_c = 1.0;
00251 if (AP==1 && ZP==1)
00252 {
00253 if (AT==2 && ZT==1) R_c = 13.5;
00254 else if (AT==3 && ZT==2) R_c = 21.0;
00255 else if (AT==4 && ZT==2) R_c = 27.0;
00256 else if (ZT==3) R_c = 2.2;
00257 }
00258 else if (AT==1 && ZT==1)
00259 {
00260 if (AP==2 && ZP==1) R_c = 13.5;
00261 else if (AP==3 && ZP==2) R_c = 21.0;
00262 else if (AP==4 && ZP==2) R_c = 27.0;
00263 else if (ZP==3) R_c = 2.2;
00264 }
00265 else if (AP==2 && ZP==1)
00266 {
00267 if (AT==2 && ZT==1) R_c = 13.5;
00268 else if (AT==4 && ZT==2) R_c = 13.5;
00269 else if (AT==12 && ZT==6) R_c = 6.0;
00270 }
00271 else if (AT==2 && ZT==1)
00272 {
00273 if (AP==2 && ZP==1) R_c = 13.5;
00274 else if (AP==4 && ZP==2) R_c = 13.5;
00275 else if (AP==12 && ZP==6) R_c = 6.0;
00276 }
00277 else if ((AP==4 && ZP==2 && (ZT==73 || ZT==79)) ||
00278 (AT==4 && ZT==2 && (ZP==73 || ZP==79))) R_c = 0.6;
00279
00280
00281
00282
00283
00284
00285
00286 G4double xr = r_0*(AT13 + AP13 + deltaE);
00287 result = pi * xr * xr * (1.0 - R_c*B/E_cm) * X_m;
00288
00289 if (result < 0.0) {
00290 result = 0.0;
00291
00292 } else if (!lowEnergyCheck && E < 6.0) {
00293 G4double f = 0.95;
00294 G4DynamicParticle slowerProjectile = *theProjectile;
00295 slowerProjectile.SetKineticEnergy(f * EA * MeV);
00296
00297 G4bool savelowenergy = lowEnergyCheck;
00298 SetLowEnergyCheck(true);
00299 G4double resultp = GetElementCrossSection(&slowerProjectile, ZT);
00300 SetLowEnergyCheck(savelowenergy);
00301
00302 if (resultp > result) { result = 0.0; }
00303 }
00304
00305 return result;
00306 }
00307