Geant4  10.01.p02
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: G4BetheBlochModel.cc 88979 2015-03-17 10:10:21Z gcosmo $
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"
68 #include "G4Log.hh"
69 #include "G4DeltaAngle.hh"
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
73 using namespace std;
74 
76  const G4String& nam)
77  : G4VEmModel(nam),
78  particle(0),
79  tlimit(DBL_MAX),
80  twoln10(2.0*G4Log(10.0)),
81  bg2lim(0.0169),
82  taulim(8.4146e-3),
83  isIon(false),
84  isInitialised(false)
85 {
86  fParticleChange = 0;
88  if(p) {
89  SetGenericIon(p);
90  SetParticle(p);
91  } else {
93  }
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
102 {}
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107  const G4DataVector&)
108 {
109  SetGenericIon(p);
110  SetParticle(p);
111 
112  //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName()
113  // << " isIon= " << isIon
114  // << G4endl;
115 
116  // always false before the run
117  SetDeexcitationFlag(false);
118 
119  if(!isInitialised) {
120  isInitialised = true;
124  }
125  }
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 
131  const G4Material* mat,
132  G4double kineticEnergy)
133 {
134  // this method is called only for ions
135  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
136  corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
137  return corrFactor;
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
143  const G4Material* mat,
144  G4double kineticEnergy)
145 {
146  //G4cout<<"G4BetheBlochModel::GetParticleCharge e= "<<kineticEnergy <<
147  // " q= " << corr->GetParticleCharge(p,mat,kineticEnergy) <<G4endl;
148  // this method is called only for ions, so no check if it is an ion
149  return corr->GetParticleCharge(p,mat,kineticEnergy);
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 
155 {
156  mass = particle->GetPDGMass();
157  spin = particle->GetPDGSpin();
159  chargeSquare = q*q;
161  ratio = electron_mass_c2/mass;
162  G4double magmom =
163  particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared);
164  magMoment2 = magmom*magmom - 1.0;
165  formfact = 0.0;
166  if(particle->GetLeptonNumber() == 0) {
167  G4double x = 0.8426*GeV;
168  if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
169  else if(mass > GeV) {
170  x /= nist->GetZ13(mass/proton_mass_c2);
171  // tlimit = 51.2*GeV*A13[iz]*A13[iz];
172  }
173  formfact = 2.0*electron_mass_c2/(x*x);
174  tlimit = 2.0/formfact;
175  }
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 
181  const G4MaterialCutsCouple* couple)
182 {
183  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
187 
188 G4double
190  G4double kineticEnergy,
191  G4double cutEnergy,
192  G4double maxKinEnergy)
193 {
194  G4double cross = 0.0;
195  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
196  G4double maxEnergy = min(tmax,maxKinEnergy);
197  if(cutEnergy < maxEnergy) {
198 
199  G4double totEnergy = kineticEnergy + mass;
200  G4double energy2 = totEnergy*totEnergy;
201  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
202 
203  cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
204  - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
205 
206  // +term for spin=1/2 particle
207  if( 0.5 == spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
208 
209  // High order correction different for hadrons and ions
210  // nevetheless they are applied to reduce high energy transfers
211  // if(!isIon)
212  //cross += corr->FiniteSizeCorrectionXS(p,currentMaterial,
213  // kineticEnergy,cutEnergy);
214 
215  cross *= twopi_mc2_rcl2*chargeSquare/beta2;
216  }
217 
218  // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
219  // << " tmax= " << tmax << " cross= " << cross << G4endl;
220 
221  return cross;
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
225 
227  const G4ParticleDefinition* p,
228  G4double kineticEnergy,
229  G4double Z, G4double,
230  G4double cutEnergy,
231  G4double maxEnergy)
232 {
234  (p,kineticEnergy,cutEnergy,maxEnergy);
235  return cross;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
241  const G4Material* material,
242  const G4ParticleDefinition* p,
243  G4double kineticEnergy,
244  G4double cutEnergy,
245  G4double maxEnergy)
246 {
247  G4double eDensity = material->GetElectronDensity();
249  (p,kineticEnergy,cutEnergy,maxEnergy);
250  return cross;
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
256  const G4ParticleDefinition* p,
257  G4double kineticEnergy,
258  G4double cut)
259 {
260  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
261  G4double cutEnergy = std::min(cut,tmax);
262 
263  G4double tau = kineticEnergy/mass;
264  G4double gam = tau + 1.0;
265  G4double bg2 = tau * (tau+2.0);
266  G4double beta2 = bg2/(gam*gam);
267 
268  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
269  G4double eexc2 = eexc*eexc;
270 
271  G4double eDensity = material->GetElectronDensity();
272 
273  G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
274  - (1.0 + cutEnergy/tmax)*beta2;
275 
276  if(0.5 == spin) {
277  G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
278  dedx += del*del;
279  }
280 
281  // density correction
282  G4double x = G4Log(bg2)/twoln10;
283  dedx -= material->GetIonisation()->DensityCorrection(x);
284 
285  // shell correction
286  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
287 
288  // now compute the total ionization loss
289  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
290 
291  //High order correction different for hadrons and ions
292  if(isIon) {
293  dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
294  } else {
295  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
296  }
297 
298  if (dedx < 0.0) { dedx = 0.0; }
299 
300  //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
301  // << " " << material->GetName() << G4endl;
302 
303  return dedx;
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307 
309  const G4DynamicParticle* dp,
310  G4double& eloss,
311  G4double&,
312  G4double length)
313 {
314  if(isIon) {
315  const G4ParticleDefinition* p = dp->GetDefinition();
316  const G4Material* mat = couple->GetMaterial();
317  G4double preKinEnergy = dp->GetKineticEnergy();
318  G4double e = preKinEnergy - eloss*0.5;
319  if(e < preKinEnergy*0.75) { e = preKinEnergy*0.75; }
320 
323  G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
324  G4double highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
325  G4double elossnew = eloss*qfactor + highOrder;
326  if(elossnew > preKinEnergy) { elossnew = preKinEnergy; }
327  else if(elossnew < eloss*0.5) { elossnew = eloss*0.5; }
328  eloss = elossnew;
329  //G4cout << "G4BetheBlochModel::CorrectionsAlongStep: e= " << preKinEnergy
330  // << " qfactor= " << qfactor
331  // << " highOrder= " << highOrder << " ("
332  // << highOrder/eloss << ")" << G4endl;
333  }
334 }
335 
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337 
338 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
339  const G4MaterialCutsCouple* couple,
340  const G4DynamicParticle* dp,
341  G4double minKinEnergy,
342  G4double maxEnergy)
343 {
344  G4double kineticEnergy = dp->GetKineticEnergy();
345  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
346 
347  G4double maxKinEnergy = std::min(maxEnergy,tmax);
348  if(minKinEnergy >= maxKinEnergy) { return; }
349 
350  //G4cout << "G4BetheBlochModel::SampleSecondaries Emin= " << minKinEnergy
351  // << " Emax= " << maxKinEnergy << G4endl;
352 
353  G4double totEnergy = kineticEnergy + mass;
354  G4double etot2 = totEnergy*totEnergy;
355  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
356 
357  G4double deltaKinEnergy, f;
358  G4double f1 = 0.0;
359  G4double fmax = 1.0;
360  if( 0.5 == spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
361 
362  // sampling without nuclear size effect
363  do {
364  G4double q = G4UniformRand();
365  deltaKinEnergy = minKinEnergy*maxKinEnergy
366  /(minKinEnergy*(1.0 - q) + maxKinEnergy*q);
367 
368  f = 1.0 - beta2*deltaKinEnergy/tmax;
369  if( 0.5 == spin ) {
370  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
371  f += f1;
372  }
373 
374  } while( fmax*G4UniformRand() > f);
375 
376  // projectile formfactor - suppresion of high energy
377  // delta-electron production at high energy
378 
379  G4double x = formfact*deltaKinEnergy;
380  if(x > 1.e-6) {
381 
382  G4double x1 = 1.0 + x;
383  G4double grej = 1.0/(x1*x1);
384  if( 0.5 == spin ) {
385  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
386  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
387  }
388  if(grej > 1.1) {
389  G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
390  << " " << dp->GetDefinition()->GetParticleName()
391  << " Ekin(MeV)= " << kineticEnergy
392  << " delEkin(MeV)= " << deltaKinEnergy
393  << G4endl;
394  }
395  if(G4UniformRand() > grej) return;
396  }
397 
398  G4ThreeVector deltaDirection;
399 
401 
402  const G4Material* mat = couple->GetMaterial();
403  G4int Z = SelectRandomAtomNumber(mat);
404 
405  deltaDirection =
406  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
407 
408  } else {
409 
410  G4double deltaMomentum =
411  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
412  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
413  (deltaMomentum * dp->GetTotalMomentum());
414  if(cost > 1.0) { cost = 1.0; }
415  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
416 
417  G4double phi = twopi * G4UniformRand() ;
418 
419  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
420  deltaDirection.rotateUz(dp->GetMomentumDirection());
421  }
422  /*
423  G4cout << "### G4BetheBlochModel "
424  << dp->GetDefinition()->GetParticleName()
425  << " Ekin(MeV)= " << kineticEnergy
426  << " delEkin(MeV)= " << deltaKinEnergy
427  << " tmin(MeV)= " << minKinEnergy
428  << " tmax(MeV)= " << maxKinEnergy
429  << " dir= " << dp->GetMomentumDirection()
430  << " dirDelta= " << deltaDirection
431  << G4endl;
432  */
433  // create G4DynamicParticle object for delta ray
434  G4DynamicParticle* delta =
435  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
436 
437  vdp->push_back(delta);
438 
439  // Change kinematics of primary particle
440  kineticEnergy -= deltaKinEnergy;
441  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
442  finalP = finalP.unit();
443 
446 }
447 
448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
449 
451  G4double kinEnergy)
452 {
453  // here particle type is checked for any method
454  SetParticle(pd);
455  G4double tau = kinEnergy/mass;
456  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
457  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
458  return std::min(tmax,tlimit);
459 }
460 
461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4double HighOrderCorrections(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy, G4double cutEnergy)
static const double MeV
Definition: G4SIunits.hh:193
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:120
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetZ13(G4double Z)
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:599
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4ParticleDefinition * GetDefinition() const
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:592
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4ParticleDefinition * particle
const G4String & GetParticleName() const
static const G4double bg2lim
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
void SetGenericIon(const G4ParticleDefinition *p)
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple)
G4NistManager * nist
G4EmCorrections * corr
G4double GetTotalMomentum() const
G4BetheBlochModel(const G4ParticleDefinition *p=0, const G4String &nam="BetheBloch")
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy)
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4bool UseAngularGeneratorFlag() const
Definition: G4VEmModel.hh:679
G4double IonBarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4EmCorrections * EmCorrections()
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double IonHighOrderCorrections(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
static const double GeV
Definition: G4SIunits.hh:196
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static const G4double f1
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, G4double &eloss, G4double &, G4double length)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4double ShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double DensityCorrection(G4double x)
G4double GetPDGMass() const
void SetParticle(const G4ParticleDefinition *p)
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:606
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ParticleChangeForLoss * fParticleChange
G4double GetPDGSpin() const
G4double GetMeanExcitationEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
static const G4double twoln10
#define G4endl
Definition: G4ios.hh:61
G4double GetPDGMagneticMoment() const
G4ParticleDefinition * theElectron
double G4double
Definition: G4Types.hh:76
G4double GetChargeSquareRatio() const
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:763
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
#define DBL_MAX
Definition: templates.hh:83
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.hh:546