Geant4  10.02.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 93362 2015-10-19 13:45:19Z 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(nullptr),
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 = nullptr;
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.0 < 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.0 < 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.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
361 
362  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
363  G4double rndm[2];
364 
365  // sampling without nuclear size effect
366  do {
367  rndmEngineMod->flatArray(2, rndm);
368  deltaKinEnergy = minKinEnergy*maxKinEnergy
369  /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
370 
371  f = 1.0 - beta2*deltaKinEnergy/tmax;
372  if( 0.0 < spin ) {
373  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
374  f += f1;
375  }
376 
377  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
378  } while( fmax*rndm[1] > f);
379 
380  // projectile formfactor - suppresion of high energy
381  // delta-electron production at high energy
382 
383  G4double x = formfact*deltaKinEnergy;
384  if(x > 1.e-6) {
385 
386  G4double x1 = 1.0 + x;
387  G4double grej = 1.0/(x1*x1);
388  if( 0.0 < spin ) {
389  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
390  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
391  }
392  if(grej > 1.1) {
393  G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
394  << " " << dp->GetDefinition()->GetParticleName()
395  << " Ekin(MeV)= " << kineticEnergy
396  << " delEkin(MeV)= " << deltaKinEnergy
397  << G4endl;
398  }
399  if(rndmEngineMod->flat() > grej) { return; }
400  }
401 
402  G4ThreeVector deltaDirection;
403 
405 
406  const G4Material* mat = couple->GetMaterial();
407  G4int Z = SelectRandomAtomNumber(mat);
408 
409  deltaDirection =
410  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
411 
412  } else {
413 
414  G4double deltaMomentum =
415  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
416  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
417  (deltaMomentum * dp->GetTotalMomentum());
418  if(cost > 1.0) { cost = 1.0; }
419  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
420 
421  G4double phi = twopi*rndmEngineMod->flat();
422 
423  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
424  deltaDirection.rotateUz(dp->GetMomentumDirection());
425  }
426  /*
427  G4cout << "### G4BetheBlochModel "
428  << dp->GetDefinition()->GetParticleName()
429  << " Ekin(MeV)= " << kineticEnergy
430  << " delEkin(MeV)= " << deltaKinEnergy
431  << " tmin(MeV)= " << minKinEnergy
432  << " tmax(MeV)= " << maxKinEnergy
433  << " dir= " << dp->GetMomentumDirection()
434  << " dirDelta= " << deltaDirection
435  << G4endl;
436  */
437  // create G4DynamicParticle object for delta ray
438  G4DynamicParticle* delta =
439  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
440 
441  vdp->push_back(delta);
442 
443  // Change kinematics of primary particle
444  kineticEnergy -= deltaKinEnergy;
445  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
446  finalP = finalP.unit();
447 
450 }
451 
452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
453 
455  G4double kinEnergy)
456 {
457  // here particle type is checked for any method
458  SetParticle(pd);
459  G4double tau = kinEnergy/mass;
460  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
461  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
462  return std::min(tmax,tlimit);
463 }
464 
465 //....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:211
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:617
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4ParticleDefinition * GetDefinition() const
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:610
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
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:697
G4double IonBarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4EmCorrections * EmCorrections()
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
Definition: G4SIunits.hh:75
void SetProposedKineticEnergy(G4double proposedKinEnergy)
G4double IonHighOrderCorrections(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
static const double GeV
Definition: G4SIunits.hh:214
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:624
const G4double x[NPOINTSGL]
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
G4BetheBlochModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="BetheBloch")
G4double GetPDGMagneticMoment() const
G4ParticleDefinition * theElectron
double G4double
Definition: G4Types.hh:76
G4double GetChargeSquareRatio() const
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:781
static const double eplus
Definition: G4SIunits.hh:196
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:564