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 "G4mplIonisationWithDeltaModel.hh"
00057 #include "Randomize.hh"
00058 #include "G4PhysicalConstants.hh"
00059 #include "G4SystemOfUnits.hh"
00060 #include "G4LossTableManager.hh"
00061 #include "G4ParticleChangeForLoss.hh"
00062 #include "G4Electron.hh"
00063 #include "G4DynamicParticle.hh"
00064
00065
00066
00067 using namespace std;
00068
00069 G4mplIonisationWithDeltaModel::G4mplIonisationWithDeltaModel(G4double mCharge,
00070 const G4String& nam)
00071 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
00072 magCharge(mCharge),
00073 twoln10(log(100.0)),
00074 betalow(0.01),
00075 betalim(0.1),
00076 beta2lim(betalim*betalim),
00077 bg2lim(beta2lim*(1.0 + beta2lim))
00078 {
00079 nmpl = G4lrint(std::fabs(magCharge) * 2 * fine_structure_const);
00080 if(nmpl > 6) { nmpl = 6; }
00081 else if(nmpl < 1) { nmpl = 1; }
00082 pi_hbarc2_over_mc2 = pi * hbarc * hbarc / electron_mass_c2;
00083 chargeSquare = magCharge * magCharge;
00084 dedxlim = 45.*nmpl*nmpl*GeV*cm2/g;
00085 fParticleChange = 0;
00086 theElectron = G4Electron::Electron();
00087 G4cout << "### Monopole ionisation model with d-electron production, Gmag= "
00088 << magCharge/eplus << G4endl;
00089 monopole = 0;
00090 mass = 0.0;
00091 }
00092
00093
00094
00095 G4mplIonisationWithDeltaModel::~G4mplIonisationWithDeltaModel()
00096 {}
00097
00098
00099
00100 void G4mplIonisationWithDeltaModel::SetParticle(const G4ParticleDefinition* p)
00101 {
00102 monopole = p;
00103 mass = monopole->GetPDGMass();
00104 G4double emin =
00105 std::min(LowEnergyLimit(),0.1*mass*(1/sqrt(1 - betalow*betalow) - 1));
00106 G4double emax =
00107 std::max(HighEnergyLimit(),10*mass*(1/sqrt(1 - beta2lim) - 1));
00108 SetLowEnergyLimit(emin);
00109 SetHighEnergyLimit(emax);
00110 }
00111
00112
00113
00114 void
00115 G4mplIonisationWithDeltaModel::Initialise(const G4ParticleDefinition* p,
00116 const G4DataVector&)
00117 {
00118 if(!monopole) { SetParticle(p); }
00119 if(!fParticleChange) { fParticleChange = GetParticleChangeForLoss(); }
00120 }
00121
00122
00123
00124 G4double
00125 G4mplIonisationWithDeltaModel::ComputeDEDXPerVolume(const G4Material* material,
00126 const G4ParticleDefinition* p,
00127 G4double kineticEnergy,
00128 G4double maxEnergy)
00129 {
00130 if(!monopole) { SetParticle(p); }
00131 G4double tmax = MaxSecondaryEnergy(p,kineticEnergy);
00132 G4double cutEnergy = std::min(tmax, maxEnergy);
00133 cutEnergy = std::max(LowEnergyLimit(), cutEnergy);
00134 G4double tau = kineticEnergy / mass;
00135 G4double gam = tau + 1.0;
00136 G4double bg2 = tau * (tau + 2.0);
00137 G4double beta2 = bg2 / (gam * gam);
00138 G4double beta = sqrt(beta2);
00139
00140
00141 G4double dedx = dedxlim*beta*material->GetDensity();
00142
00143
00144 if(beta > betalow) {
00145
00146
00147 if(beta >= betalim) {
00148 dedx = ComputeDEDXAhlen(material, bg2, cutEnergy);
00149
00150 } else {
00151
00152 G4double dedx1 = dedxlim*betalow*material->GetDensity();
00153 G4double dedx2 = ComputeDEDXAhlen(material, bg2lim, cutEnergy);
00154
00155
00156 G4double kapa2 = beta - betalow;
00157 G4double kapa1 = betalim - beta;
00158 dedx = (kapa1*dedx1 + kapa2*dedx2)/(kapa1 + kapa2);
00159 }
00160 }
00161 return dedx;
00162 }
00163
00164
00165
00166 G4double
00167 G4mplIonisationWithDeltaModel::ComputeDEDXAhlen(const G4Material* material,
00168 G4double bg2,
00169 G4double cutEnergy)
00170 {
00171 G4double eDensity = material->GetElectronDensity();
00172 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
00173
00174
00175 G4double dedx =
00176 0.5*(log(2.0 * electron_mass_c2 * bg2*cutEnergy / (eexc*eexc)) - 1.0);
00177
00178
00179 G4double k = 0.406;
00180 if(nmpl > 1) { k = 0.346; }
00181
00182
00183 const G4double B[7] = { 0.0, 0.248, 0.672, 1.022, 1.243, 1.464, 1.685};
00184
00185 dedx += 0.5 * k - B[nmpl];
00186
00187
00188 G4double x = log(bg2)/twoln10;
00189 dedx -= material->GetIonisation()->DensityCorrection(x);
00190
00191
00192 dedx *= pi_hbarc2_over_mc2 * eDensity * nmpl * nmpl;
00193
00194 if (dedx < 0.0) { dedx = 0; }
00195 return dedx;
00196 }
00197
00198
00199
00200 G4double
00201 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerElectron(
00202 const G4ParticleDefinition* p,
00203 G4double kineticEnergy,
00204 G4double cut,
00205 G4double maxKinEnergy)
00206 {
00207 if(!monopole) { SetParticle(p); }
00208 G4double cross = 0.0;
00209 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
00210 G4double maxEnergy = std::min(tmax,maxKinEnergy);
00211 G4double cutEnergy = std::max(LowEnergyLimit(), cut);
00212 if(cutEnergy < maxEnergy) {
00213 cross = (0.5/cutEnergy - 0.5/maxEnergy)*pi_hbarc2_over_mc2 * nmpl * nmpl;
00214 }
00215 return cross;
00216 }
00217
00218
00219
00220 G4double
00221 G4mplIonisationWithDeltaModel::ComputeCrossSectionPerAtom(
00222 const G4ParticleDefinition* p,
00223 G4double kineticEnergy,
00224 G4double Z, G4double,
00225 G4double cutEnergy,
00226 G4double maxEnergy)
00227 {
00228 G4double cross =
00229 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
00230 return cross;
00231 }
00232
00233
00234
00235 void
00236 G4mplIonisationWithDeltaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00237 const G4MaterialCutsCouple*,
00238 const G4DynamicParticle* dp,
00239 G4double minKinEnergy,
00240 G4double maxEnergy)
00241 {
00242 G4double kineticEnergy = dp->GetKineticEnergy();
00243 G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
00244
00245 G4double maxKinEnergy = std::min(maxEnergy,tmax);
00246 if(minKinEnergy >= maxKinEnergy) { return; }
00247
00248
00249
00250
00251
00252 G4double totEnergy = kineticEnergy + mass;
00253 G4double etot2 = totEnergy*totEnergy;
00254 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
00255
00256
00257 G4double q = G4UniformRand();
00258 G4double deltaKinEnergy = minKinEnergy*maxKinEnergy
00259 /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
00260
00261
00262 G4double totMomentum = totEnergy*sqrt(beta2);
00263 G4double deltaMomentum =
00264 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
00265 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
00266 (deltaMomentum * totMomentum);
00267 if(cost > 1.0) { cost = 1.0; }
00268
00269 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
00270
00271 G4double phi = twopi * G4UniformRand() ;
00272
00273 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
00274 G4ThreeVector direction = dp->GetMomentumDirection();
00275 deltaDirection.rotateUz(direction);
00276
00277
00278 G4DynamicParticle* delta =
00279 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
00280
00281 vdp->push_back(delta);
00282
00283
00284 kineticEnergy -= deltaKinEnergy;
00285 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
00286 finalP = finalP.unit();
00287
00288 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
00289 fParticleChange->SetProposedMomentumDirection(finalP);
00290 }
00291
00292
00293
00294 G4double G4mplIonisationWithDeltaModel::SampleFluctuations(
00295 const G4Material* material,
00296 const G4DynamicParticle* dp,
00297 G4double& tmax,
00298 G4double& length,
00299 G4double& meanLoss)
00300 {
00301 G4double siga = Dispersion(material,dp,tmax,length);
00302 G4double loss = meanLoss;
00303 siga = sqrt(siga);
00304 G4double twomeanLoss = meanLoss + meanLoss;
00305
00306 if(twomeanLoss < siga) {
00307 G4double x;
00308 do {
00309 loss = twomeanLoss*G4UniformRand();
00310 x = (loss - meanLoss)/siga;
00311 } while (1.0 - 0.5*x*x < G4UniformRand());
00312 } else {
00313 do {
00314 loss = G4RandGauss::shoot(meanLoss,siga);
00315 } while (0.0 > loss || loss > twomeanLoss);
00316 }
00317 return loss;
00318 }
00319
00320
00321
00322 G4double
00323 G4mplIonisationWithDeltaModel::Dispersion(const G4Material* material,
00324 const G4DynamicParticle* dp,
00325 G4double& tmax,
00326 G4double& length)
00327 {
00328 G4double siga = 0.0;
00329 G4double tau = dp->GetKineticEnergy()/mass;
00330 if(tau > 0.0) {
00331 G4double electronDensity = material->GetElectronDensity();
00332 G4double gam = tau + 1.0;
00333 G4double invbeta2 = (gam*gam)/(tau * (tau+2.0));
00334 siga = (invbeta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
00335 * electronDensity * chargeSquare;
00336 }
00337 return siga;
00338 }
00339
00340
00341
00342 G4double
00343 G4mplIonisationWithDeltaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
00344 G4double kinEnergy)
00345 {
00346 G4double tau = kinEnergy/mass;
00347 return 2.0*electron_mass_c2*tau*(tau + 2.);
00348 }
00349
00350