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