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