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 #include "G4IonCoulombCrossSection.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "Randomize.hh"
00059 #include "G4Proton.hh"
00060 #include "G4LossTableManager.hh"
00061
00062
00063
00064
00065 using namespace std;
00066
00067 G4IonCoulombCrossSection::G4IonCoulombCrossSection():
00068 cosThetaMin(1.0),
00069 cosThetaMax(-1.0),
00070 alpha2(fine_structure_const*fine_structure_const)
00071 {
00072 fNistManager = G4NistManager::Instance();
00073 theProton = G4Proton::Proton();
00074 particle=0;
00075
00076 G4double p0 = electron_mass_c2*classic_electr_radius;
00077 coeff = twopi*p0*p0;
00078
00079 cosTetMinNuc=0;
00080 cosTetMaxNuc=0;
00081 nucXSection =0;
00082
00083
00084 chargeSquare = spin = mass =0;
00085 tkinLab = momLab2 = invbetaLab2=0;
00086 tkin = mom2 = invbeta2=0;
00087
00088 targetZ = targetMass = screenZ =0;
00089 ScreenRSquare=0.;
00090
00091 etag = ecut = 0.0;
00092
00093 }
00094
00095
00096 G4IonCoulombCrossSection::~G4IonCoulombCrossSection()
00097 {}
00098
00099
00100 void G4IonCoulombCrossSection::Initialise(const G4ParticleDefinition* p,
00101 G4double CosThetaLim)
00102 {
00103
00104 SetupParticle(p);
00105 nucXSection = 0.0;
00106 tkin = targetZ = mom2 = DBL_MIN;
00107 ecut = etag = DBL_MAX;
00108 particle = p;
00109
00110
00111 cosThetaMin = CosThetaLim;
00112
00113 }
00114
00115
00116 void G4IonCoulombCrossSection::SetupKinematic(G4double ekin,
00117 G4double cut,G4int iz )
00118 {
00119 if(ekin != tkinLab || ecut != cut) {
00120
00121
00122 tkinLab = ekin;
00123 momLab2 = tkinLab*(tkinLab + 2.0*mass);
00124 invbetaLab2 = 1.0 + mass*mass/momLab2;
00125
00126 G4double etot = tkinLab + mass;
00127 G4double ptot = sqrt(momLab2);
00128 G4double m12 = mass*mass;
00129
00130 targetMass=fNistManager->GetAtomicMassAmu(iz)*amu_c2;
00131
00132
00133
00134
00135
00136 G4double Ecm=sqrt(m12 + targetMass*targetMass + 2.0*etot*targetMass);
00137 G4double mu_rel=mass*targetMass/Ecm;
00138 G4double momCM= ptot*targetMass/Ecm;
00139
00140 mom2 = momCM*momCM;
00141 invbeta2 = 1.0 + mu_rel*mu_rel/mom2;
00142 tkin = momCM*sqrt(invbeta2) - mu_rel;
00143
00144
00145 cosTetMinNuc = cosThetaMin;
00146 cosTetMaxNuc = cosThetaMax;
00147
00148 }
00149
00150 }
00151
00152
00153
00154 void G4IonCoulombCrossSection::SetupTarget(G4double Z, G4double e, G4int heavycorr)
00155 {
00156 if(Z != targetZ || e != etag) {
00157 etag = e;
00158 targetZ = Z;
00159 G4int iz= G4int(Z);
00160
00161 SetScreenRSquare(iz);
00162 screenZ =0;
00163 screenZ = ScreenRSquare/mom2;
00164
00165
00166
00167 if(heavycorr!=0 && particle != theProton){
00168 G4double corr=5.*twopi*Z*std::sqrt(chargeSquare*alpha2);
00169 corr=std::pow(corr,0.12);
00170 screenZ *=(1.13 + corr*3.76*Z*Z*chargeSquare*invbeta2*alpha2)/2.;
00171
00172
00173
00174 }else{ screenZ *=(1.13 + 3.76*Z*Z*chargeSquare*invbeta2*alpha2)/2.;
00175
00176
00177 }
00178
00179 if(1 == iz && particle == theProton && cosTetMaxNuc < 0.0) {
00180 cosTetMaxNuc = 0.0;
00181 }
00182
00183 }
00184 }
00185
00186
00187 void G4IonCoulombCrossSection::SetScreenRSquare(G4int iz){
00188
00189 G4double a0 = electron_mass_c2/0.88534;
00190
00191
00192 G4double Z1=std::sqrt(chargeSquare);
00193 G4double Z2=targetZ;
00194
00195 G4double Z1023=std::pow(Z1,0.23);
00196 G4double Z2023=std::pow(Z2,0.23);
00197
00198
00199 G4double x=a0*(Z1023+Z2023);
00200
00201
00202 if(particle == theProton){
00203 x = a0*fNistManager->GetZ13(iz);
00204 }
00205 ScreenRSquare = alpha2*x*x;
00206
00207
00208 }
00209
00210
00211 G4double G4IonCoulombCrossSection::NuclearCrossSection()
00212 {
00213
00214
00215
00216 G4double fac = coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2;
00217
00218 nucXSection = 0.0;
00219
00220 G4double x = 1.0 - cosTetMinNuc;
00221 G4double x1 = x + screenZ;
00222
00223
00224 if(cosTetMaxNuc < cosTetMinNuc) {
00225
00226 nucXSection =fac*(cosTetMinNuc - cosTetMaxNuc)/
00227 (x1*(1.0 - cosTetMaxNuc + screenZ));
00228
00229 }
00230
00231 return nucXSection;
00232 }
00233
00234
00235
00236 G4double G4IonCoulombCrossSection::SampleCosineTheta()
00237 {
00238
00239 if(cosTetMaxNuc >= cosTetMinNuc) return 0.0;
00240
00241 G4double x1 = 1. - cosTetMinNuc + screenZ;
00242 G4double x2 = 1. - cosTetMaxNuc + screenZ;
00243 G4double dx = cosTetMinNuc - cosTetMaxNuc;
00244 G4double z1;
00245
00246 z1 = x1*x2/(x1 + G4UniformRand()*dx) - screenZ;
00247
00248 return z1;
00249 }
00250
00251
00252
00253
00254
00255