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 "G4mplIonisationModel.hh"
00057 #include "Randomize.hh"
00058 #include "G4PhysicalConstants.hh"
00059 #include "G4SystemOfUnits.hh"
00060 #include "G4LossTableManager.hh"
00061 #include "G4ParticleChangeForLoss.hh"
00062
00063
00064
00065 using namespace std;
00066
00067 G4mplIonisationModel::G4mplIonisationModel(G4double mCharge, const G4String& nam)
00068 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
00069 magCharge(mCharge),
00070 twoln10(log(100.0)),
00071 betalow(0.01),
00072 betalim(0.1),
00073 beta2lim(betalim*betalim),
00074 bg2lim(beta2lim*(1.0 + beta2lim))
00075 {
00076 nmpl = G4int(abs(magCharge) * 2 * fine_structure_const + 0.5);
00077 if(nmpl > 6) { nmpl = 6; }
00078 else if(nmpl < 1) { nmpl = 1; }
00079 pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
00080 chargeSquare = magCharge * magCharge;
00081 dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
00082 fParticleChange = 0;
00083 monopole = 0;
00084 mass = 0.0;
00085 }
00086
00087
00088
00089 G4mplIonisationModel::~G4mplIonisationModel()
00090 {}
00091
00092
00093
00094 void G4mplIonisationModel::SetParticle(const G4ParticleDefinition* p)
00095 {
00096 monopole = p;
00097 mass = monopole->GetPDGMass();
00098 G4double emin =
00099 std::min(LowEnergyLimit(),0.1*mass*(1/sqrt(1 - betalow*betalow) - 1));
00100 G4double emax =
00101 std::max(HighEnergyLimit(),10*mass*(1/sqrt(1 - beta2lim) - 1));
00102 SetLowEnergyLimit(emin);
00103 SetHighEnergyLimit(emax);
00104 }
00105
00106
00107
00108 void G4mplIonisationModel::Initialise(const G4ParticleDefinition* p,
00109 const G4DataVector&)
00110 {
00111 if(!monopole) { SetParticle(p); }
00112 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00113 }
00114
00115
00116
00117 G4double G4mplIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
00118 const G4ParticleDefinition* p,
00119 G4double kineticEnergy,
00120 G4double)
00121 {
00122 if(!monopole) { SetParticle(p); }
00123 G4double tau = kineticEnergy / mass;
00124 G4double gam = tau + 1.0;
00125 G4double bg2 = tau * (tau + 2.0);
00126 G4double beta2 = bg2 / (gam * gam);
00127 G4double beta = sqrt(beta2);
00128
00129
00130 G4double dedx = dedxlim*beta*material->GetDensity();
00131
00132
00133 if(beta > betalow) {
00134
00135
00136 if(beta >= betalim) {
00137 dedx = ComputeDEDXAhlen(material, bg2);
00138
00139 } else {
00140
00141 G4double dedx1 = dedxlim*betalow*material->GetDensity();
00142 G4double dedx2 = ComputeDEDXAhlen(material, bg2lim);
00143
00144
00145 G4double kapa2 = beta - betalow;
00146 G4double kapa1 = betalim - beta;
00147 dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
00148 }
00149 }
00150 return dedx;
00151 }
00152
00153
00154
00155 G4double G4mplIonisationModel::ComputeDEDXAhlen(const G4Material* material,
00156 G4double bg2)
00157 {
00158 G4double eDensity = material->GetElectronDensity();
00159 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
00160 G4double cden = material->GetIonisation()->GetCdensity();
00161 G4double mden = material->GetIonisation()->GetMdensity();
00162 G4double aden = material->GetIonisation()->GetAdensity();
00163 G4double x0den = material->GetIonisation()->GetX0density();
00164 G4double x1den = material->GetIonisation()->GetX1density();
00165
00166
00167 G4double dedx = log(2.0 * electron_mass_c2 * bg2 / eexc) - 0.5;
00168
00169
00170 G4double k = 0.406;
00171 if(nmpl > 1) k = 0.346;
00172
00173
00174 const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
00175
00176 dedx += 0.5 * k - B[nmpl];
00177
00178
00179 G4double deltam;
00180 G4double x = log(bg2) / twoln10;
00181 if ( x >= x0den ) {
00182 deltam = twoln10 * x - cden;
00183 if ( x < x1den ) deltam += aden * pow((x1den-x), mden);
00184 dedx -= 0.5 * deltam;
00185 }
00186
00187
00188 dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
00189
00190 if (dedx < 0.0) dedx = 0;
00191 return dedx;
00192 }
00193
00194
00195
00196 void G4mplIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
00197 const G4MaterialCutsCouple*,
00198 const G4DynamicParticle*,
00199 G4double,
00200 G4double)
00201 {}
00202
00203
00204
00205 G4double G4mplIonisationModel::SampleFluctuations(
00206 const G4Material* material,
00207 const G4DynamicParticle* dp,
00208 G4double& tmax,
00209 G4double& length,
00210 G4double& meanLoss)
00211 {
00212 G4double siga = Dispersion(material,dp,tmax,length);
00213 G4double loss = meanLoss;
00214 siga = sqrt(siga);
00215 G4double twomeanLoss = meanLoss + meanLoss;
00216
00217 if(twomeanLoss < siga) {
00218 G4double x;
00219 do {
00220 loss = twomeanLoss*G4UniformRand();
00221 x = (loss - meanLoss)/siga;
00222 } while (1.0 - 0.5*x*x < G4UniformRand());
00223 } else {
00224 do {
00225 loss = G4RandGauss::shoot(meanLoss,siga);
00226 } while (0.0 > loss || loss > twomeanLoss);
00227 }
00228 return loss;
00229 }
00230
00231
00232
00233 G4double G4mplIonisationModel::Dispersion(const G4Material* material,
00234 const G4DynamicParticle* dp,
00235 G4double& tmax,
00236 G4double& length)
00237 {
00238 G4double siga = 0.0;
00239 G4double tau = dp->GetKineticEnergy()/mass;
00240 if(tau > 0.0) {
00241 G4double electronDensity = material->GetElectronDensity();
00242 G4double gam = tau + 1.0;
00243 G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
00244 siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00245 * electronDensity * chargeSquare;
00246 }
00247 return siga;
00248 }
00249
00250