Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetheBlochModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class header file
31 //
32 //
33 // File name: G4BetheBlochModel
34 //
35 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
42 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanchenko)
44 // 13-02-03 Add name (V.Ivanchenko)
45 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)
46 // 11-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
47 // 11-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
48 // 12-02-06 move G4LossTableManager::Instance()->EmCorrections()
49 // in constructor (mma)
50 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
51 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
52 //
53 // -------------------------------------------------------------------
54 //
55 
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 #include "G4BetheBlochModel.hh"
61 #include "Randomize.hh"
62 #include "G4PhysicalConstants.hh"
63 #include "G4SystemOfUnits.hh"
64 #include "G4Electron.hh"
65 #include "G4LossTableManager.hh"
66 #include "G4EmCorrections.hh"
67 #include "G4ParticleChangeForLoss.hh"
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 using namespace std;
72 
73 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition* p,
74  const G4String& nam)
75  : G4VEmModel(nam),
76  particle(0),
77  tlimit(DBL_MAX),
78  twoln10(2.0*log(10.0)),
79  bg2lim(0.0169),
80  taulim(8.4146e-3),
81  isIon(false),
82  isInitialised(false)
83 {
84  fParticleChange = 0;
85  theElectron = G4Electron::Electron();
86  if(p) {
87  SetGenericIon(p);
88  SetParticle(p);
89  } else {
90  SetParticle(theElectron);
91  }
92  corr = G4LossTableManager::Instance()->EmCorrections();
93  nist = G4NistManager::Instance();
94  SetLowEnergyLimit(2.0*MeV);
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 G4BetheBlochModel::~G4BetheBlochModel()
100 {}
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
105  const G4DataVector&)
106 {
107  SetGenericIon(p);
108  SetParticle(p);
109 
110  //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName()
111  // << " isIon= " << isIon
112  // << G4endl;
113 
114  // always false before the run
115  SetDeexcitationFlag(false);
116 
117  if(!isInitialised) {
118  isInitialised = true;
119  fParticleChange = GetParticleChangeForLoss();
120  }
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
125 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
126  const G4Material* mat,
127  G4double kineticEnergy)
128 {
129  // this method is called only for ions
130  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
131  corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
132  return corrFactor;
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 
137 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
138  const G4Material* mat,
139  G4double kineticEnergy)
140 {
141  // this method is called only for ions, so no check if it is an ion
142  return corr->GetParticleCharge(p,mat,kineticEnergy);
143 }
144 
145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146 
147 void G4BetheBlochModel::SetupParameters()
148 {
149  mass = particle->GetPDGMass();
150  spin = particle->GetPDGSpin();
151  G4double q = particle->GetPDGCharge()/eplus;
152  chargeSquare = q*q;
153  corrFactor = chargeSquare;
154  ratio = electron_mass_c2/mass;
155  G4double magmom =
156  particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared);
157  magMoment2 = magmom*magmom - 1.0;
158  formfact = 0.0;
159  if(particle->GetLeptonNumber() == 0) {
160  G4double x = 0.8426*GeV;
161  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
162  else if(mass > GeV) {
163  x /= nist->GetZ13(mass/proton_mass_c2);
164  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
165  }
166  formfact = 2.0*electron_mass_c2/(x*x);
167  tlimit = 2.0/formfact;
168  }
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 
173 G4double
174 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
175  G4double kineticEnergy,
176  G4double cutEnergy,
177  G4double maxKinEnergy)
178 {
179  G4double cross = 0.0;
180  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
181  G4double maxEnergy = min(tmax,maxKinEnergy);
182  if(cutEnergy < maxEnergy) {
183 
184  G4double totEnergy = kineticEnergy + mass;
185  G4double energy2 = totEnergy*totEnergy;
186  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
187 
188  cross = 1.0/cutEnergy - 1.0/maxEnergy
189  - beta2*log(maxEnergy/cutEnergy)/tmax;
190 
191  // +term for spin=1/2 particle
192  if( 0.5 == spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
193 
194  // High order correction different for hadrons and ions
195  // nevetheless they are applied to reduce high energy transfers
196  // if(!isIon)
197  //cross += corr->FiniteSizeCorrectionXS(p,currentMaterial,
198  // kineticEnergy,cutEnergy);
199 
200  cross *= twopi_mc2_rcl2*chargeSquare/beta2;
201  }
202 
203  // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
204  // << " tmax= " << tmax << " cross= " << cross << G4endl;
205 
206  return cross;
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 
211 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
212  const G4ParticleDefinition* p,
213  G4double kineticEnergy,
214  G4double Z, G4double,
215  G4double cutEnergy,
216  G4double maxEnergy)
217 {
218  G4double cross = Z*ComputeCrossSectionPerElectron
219  (p,kineticEnergy,cutEnergy,maxEnergy);
220  return cross;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
224 
225 G4double G4BetheBlochModel::CrossSectionPerVolume(
226  const G4Material* material,
227  const G4ParticleDefinition* p,
228  G4double kineticEnergy,
229  G4double cutEnergy,
230  G4double maxEnergy)
231 {
232  G4double eDensity = material->GetElectronDensity();
233  G4double cross = eDensity*ComputeCrossSectionPerElectron
234  (p,kineticEnergy,cutEnergy,maxEnergy);
235  return cross;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
240 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
241  const G4ParticleDefinition* p,
242  G4double kineticEnergy,
243  G4double cut)
244 {
245  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
246  G4double cutEnergy = std::min(cut,tmax);
247 
248  G4double tau = kineticEnergy/mass;
249  G4double gam = tau + 1.0;
250  G4double bg2 = tau * (tau+2.0);
251  G4double beta2 = bg2/(gam*gam);
252 
253  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
254  G4double eexc2 = eexc*eexc;
255 
256  G4double eDensity = material->GetElectronDensity();
257 
258  G4double dedx = log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
259  - (1.0 + cutEnergy/tmax)*beta2;
260 
261  if(0.5 == spin) {
262  G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
263  dedx += del*del;
264  }
265 
266  // density correction
267  G4double x = log(bg2)/twoln10;
268  dedx -= material->GetIonisation()->DensityCorrection(x);
269 
270  // shell correction
271  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
272 
273  // now compute the total ionization loss
274  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
275 
276  //High order correction different for hadrons and ions
277  if(isIon) {
278  dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
279  } else {
280  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
281  }
282 
283  if (dedx < 0.0) { dedx = 0.0; }
284 
285  //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
286  // << " " << material->GetName() << G4endl;
287 
288  return dedx;
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
293 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
294  const G4DynamicParticle* dp,
295  G4double& eloss,
296  G4double&,
297  G4double length)
298 {
299  if(isIon) {
300  const G4ParticleDefinition* p = dp->GetDefinition();
301  const G4Material* mat = couple->GetMaterial();
302  G4double preKinEnergy = dp->GetKineticEnergy();
303  G4double e = preKinEnergy - eloss*0.5;
304  if(e < preKinEnergy*0.75) { e = preKinEnergy*0.75; }
305 
306  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
307  GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
308  G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
309  G4double highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
310  G4double elossnew = eloss*qfactor + highOrder;
311  if(elossnew > preKinEnergy) { elossnew = preKinEnergy; }
312  else if(elossnew < eloss*0.5) { elossnew = eloss*0.5; }
313  eloss = elossnew;
314  //G4cout << "G4BetheBlochModel::CorrectionsAlongStep: e= " << preKinEnergy
315  // << " qfactor= " << qfactor
316  // << " highOrder= " << highOrder << " (" << highOrder/eloss << ")" << G4endl;
317  }
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321 
322 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
323  const G4MaterialCutsCouple*,
324  const G4DynamicParticle* dp,
325  G4double minKinEnergy,
326  G4double maxEnergy)
327 {
328  G4double kineticEnergy = dp->GetKineticEnergy();
329  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
330 
331  G4double maxKinEnergy = std::min(maxEnergy,tmax);
332  if(minKinEnergy >= maxKinEnergy) { return; }
333 
334  G4double totEnergy = kineticEnergy + mass;
335  G4double etot2 = totEnergy*totEnergy;
336  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
337 
338  G4double deltaKinEnergy, f;
339  G4double f1 = 0.0;
340  G4double fmax = 1.0;
341  if( 0.5 == spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
342 
343  // sampling without nuclear size effect
344  do {
345  G4double q = G4UniformRand();
346  deltaKinEnergy = minKinEnergy*maxKinEnergy
347  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
348 
349  f = 1.0 - beta2*deltaKinEnergy/tmax;
350  if( 0.5 == spin ) {
351  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
352  f += f1;
353  }
354 
355  } while( fmax*G4UniformRand() > f);
356 
357  // projectile formfactor - suppresion of high energy
358  // delta-electron production at high energy
359 
360  G4double x = formfact*deltaKinEnergy;
361  if(x > 1.e-6) {
362 
363  G4double x1 = 1.0 + x;
364  G4double grej = 1.0/(x1*x1);
365  if( 0.5 == spin ) {
366  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
367  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
368  }
369  if(grej > 1.1) {
370  G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
371  << " " << dp->GetDefinition()->GetParticleName()
372  << " Ekin(MeV)= " << kineticEnergy
373  << " delEkin(MeV)= " << deltaKinEnergy
374  << G4endl;
375  }
376  if(G4UniformRand() > grej) return;
377  }
378 
379  // delta-electron is produced
380  G4double totMomentum = totEnergy*sqrt(beta2);
381  G4double deltaMomentum =
382  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
383  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
384  (deltaMomentum * totMomentum);
385  /*
386  if(cost > 1.0) {
387  G4cout << "### G4BetheBlochModel WARNING: cost= "
388  << cost << " > 1 for "
389  << dp->GetDefinition()->GetParticleName()
390  << " Ekin(MeV)= " << kineticEnergy
391  << " p(MeV/c)= " << totMomentum
392  << " delEkin(MeV)= " << deltaKinEnergy
393  << " delMom(MeV/c)= " << deltaMomentum
394  << " tmin(MeV)= " << minKinEnergy
395  << " tmax(MeV)= " << maxKinEnergy
396  << " dir= " << dp->GetMomentumDirection()
397  << G4endl;
398  cost = 1.0;
399  }
400  */
401  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
402 
403  G4double phi = twopi * G4UniformRand() ;
404 
405 
406  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost);
407  G4ThreeVector direction = dp->GetMomentumDirection();
408  deltaDirection.rotateUz(direction);
409 
410  // create G4DynamicParticle object for delta ray
411  G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
412  deltaDirection,deltaKinEnergy);
413 
414  vdp->push_back(delta);
415 
416  // Change kinematics of primary particle
417  kineticEnergy -= deltaKinEnergy;
418  G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
419  finalP = finalP.unit();
420 
421  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
422  fParticleChange->SetProposedMomentumDirection(finalP);
423 }
424 
425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426 
427 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
428  G4double kinEnergy)
429 {
430  // here particle type is checked for any method
431  SetParticle(pd);
432  G4double tau = kinEnergy/mass;
433  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
434  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
435  return std::min(tmax,tlimit);
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4