Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4LindhardSorensenIonModel.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: G4LindhardSorensenIonModel.cc 108805 2018-03-07 18:58:45Z vnivanch $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class header file
31 //
32 //
33 // File name: G4LindhardSorensenIonModel
34 //
35 // Author: Alexander Bagulya & Vladimir Ivanchenko
36 //
37 // Creation date: 16.04.2018
38 //
39 //
40 // -------------------------------------------------------------------
41 //
42 
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
48 #include "Randomize.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4Electron.hh"
52 #include "G4LossTableManager.hh"
53 #include "G4EmCorrections.hh"
55 #include "G4Log.hh"
56 #include "G4DeltaAngle.hh"
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 using namespace std;
62 
64 
66  const G4String& nam)
67  : G4VEmModel(nam),
68  particle(nullptr),
69  tlimit(DBL_MAX),
70  twoln10(2.0*G4Log(10.0))
71 {
72  fParticleChange = nullptr;
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 
83 {}
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 
88  const G4DataVector&)
89 {
90  SetParticle(p);
91 
92  //G4cout << "G4LindhardSorensenIonModel::Initialise for " << p->GetParticleName()
93  // << G4endl;
94 
95  // always false before the run
96  SetDeexcitationFlag(false);
97 
98  if(nullptr == fParticleChange) {
102  }
103  }
104  if(IsMaster() && !lsdata) {
106  }
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 G4double
113  const G4Material*,
114  G4double)
115 {
116  return chargeSquare;
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 G4double
123  const G4Material*,
124  G4double)
125 {
126  return charge;
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130 
132 {
133  mass = particle->GetPDGMass();
134  spin = particle->GetPDGSpin();
136  Zin = G4lrint(charge);
139  static const G4double aMag = 1./(0.5*eplus*hbar_Planck*c_squared);
140  G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
141  magMoment2 = magmom*magmom - 1.0;
142  if(Zin <= 1) {
143  formfact = (spin == 0.0 && mass < GeV) ? 1.181e-6 : 1.548e-6;
144  } else {
145  G4double x = nist->GetA27(Zin);
146  formfact = 3.969e-6*x*x;
147  }
148  tlimit = std::sqrt(0.414/formfact +
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153 
155  const G4MaterialCutsCouple* couple)
156 {
157  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 G4double
164  const G4ParticleDefinition* p,
166  G4double cutEnergy,
167  G4double maxKinEnergy)
168 {
169  G4double cross = 0.0;
170  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
171  G4double maxEnergy = std::min(tmax,maxKinEnergy);
172  if(cutEnergy < maxEnergy) {
173 
174  G4double totEnergy = kineticEnergy + mass;
175  G4double energy2 = totEnergy*totEnergy;
176  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
177 
178  cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
179  - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
180 
181  // +term for spin=1/2 particle
182  if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
183 
184  cross *= twopi_mc2_rcl2*chargeSquare/beta2;
185  }
186 
187  // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
188  // << " tmax= " << tmax << " cross= " << cross << G4endl;
189 
190  return cross;
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194 
196  const G4ParticleDefinition* p,
199  G4double cutEnergy,
200  G4double maxEnergy)
201 {
203  (p,kineticEnergy,cutEnergy,maxEnergy);
204  return cross;
205 }
206 
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208 
210  const G4Material* material,
211  const G4ParticleDefinition* p,
213  G4double cutEnergy,
214  G4double maxEnergy)
215 {
216  G4double eDensity = material->GetElectronDensity();
218  (p,kineticEnergy,cutEnergy,maxEnergy);
219  return cross;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
223 
224 G4double
226  const G4ParticleDefinition* p,
228  G4double cut)
229 {
230  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
231  G4double cutEnergy = std::min(cut,tmax);
232 
233  G4double tau = kineticEnergy/mass;
234  G4double gam = tau + 1.0;
235  G4double bg2 = tau * (tau+2.0);
236  G4double beta2 = bg2/(gam*gam);
237 
238  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
239  G4double eexc2 = eexc*eexc;
240 
241  G4double eDensity = material->GetElectronDensity();
242 
243  G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
244  - (1.0 + cutEnergy/tmax)*beta2;
245 
246  if(0.0 < spin) {
247  G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
248  dedx += del*del;
249  }
250 
251  // density correction
252  G4double x = G4Log(bg2)/twoln10;
253  dedx -= material->GetIonisation()->DensityCorrection(x);
254 
255  // shell correction
256  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
257 
258  //High order correction different for hadrons and ions
259  dedx += 2.0*corr->BarkasCorrection(p,material,kineticEnergy);
260  dedx = std::max(dedx, 0.0);
261 
262  // now compute the total ionization loss
263  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
264 
265  //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
266  // << " " << material->GetName() << G4endl;
267 
268  return dedx;
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272 
273 void
275  const G4DynamicParticle* dp,
276  G4double& eloss,
277  G4double&,
278  G4double length)
279 {
280  const G4ParticleDefinition* p = dp->GetDefinition();
281  SetParticle(p);
282  G4double eDensity = couple->GetMaterial()->GetElectronDensity();
283  G4double preKinEnergy = dp->GetKineticEnergy();
284  G4double e = preKinEnergy - eloss*0.5;
286 
287  G4double tau = e/mass;
288  G4double gam = tau + 1.0;
289  G4double beta2 = tau * (tau+2.0)/(gam*gam);
290  G4double deltaL0 =
291  2.0*corr->BarkasCorrection (p, couple->GetMaterial(), e)*(charge-1.)/charge;
292  G4double deltaL = lsdata->GetDeltaL(Zin, tau);
293 
294  G4double elossnew =
295  eloss + twopi_mc2_rcl2*chargeSquare*eDensity*(deltaL+deltaL0)*length/beta2;
296  /*
297  G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)= "
298  << preKinEnergy/GeV << " eloss(MeV)= " << eloss
299  << " L= " << eloss*beta2/(twopi_mc2_rcl2*chargeSquare*eDensity*length)
300  << " dL0= " << deltaL0
301  << " dL= " << deltaL << G4endl;
302  */
303  if(elossnew > preKinEnergy) { elossnew = preKinEnergy; }
304  else if(elossnew < 0.0) { elossnew = eloss*0.5; }
305 
306  eloss = elossnew;
307 }
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310 
312  vector<G4DynamicParticle*>* vdp,
313  const G4MaterialCutsCouple* couple,
314  const G4DynamicParticle* dp,
315  G4double minKinEnergy,
316  G4double maxEnergy)
317 {
320 
321  G4double maxKinEnergy = std::min(maxEnergy,tmax);
322  if(minKinEnergy >= maxKinEnergy) { return; }
323 
324  //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= " << minKinEnergy
325  // << " Emax= " << maxKinEnergy << G4endl;
326 
327  G4double totEnergy = kineticEnergy + mass;
328  G4double etot2 = totEnergy*totEnergy;
329  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
330 
331  G4double deltaKinEnergy, f;
332  G4double f1 = 0.0;
333  G4double fmax = 1.0;
334  if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
335 
336  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
337  G4double rndm[2];
338 
339  // sampling without nuclear size effect
340  do {
341  rndmEngineMod->flatArray(2, rndm);
342  deltaKinEnergy = minKinEnergy*maxKinEnergy
343  /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
344 
345  f = 1.0 - beta2*deltaKinEnergy/tmax;
346  if( 0.0 < spin ) {
347  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
348  f += f1;
349  }
350 
351  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
352  } while( fmax*rndm[1] > f);
353 
354  // projectile formfactor - suppresion of high energy
355  // delta-electron production at high energy
356 
357  G4double x = formfact*deltaKinEnergy*(deltaKinEnergy + 2*electron_mass_c2);
358  if(x > 1.e-6) {
359 
360  G4double x1 = 1.0 + x;
361  G4double grej = 1.0/(x1*x1);
362  if( 0.0 < spin ) {
363  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
364  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
365  }
366  if(grej > 1.1) {
367  G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej
368  << " " << dp->GetDefinition()->GetParticleName()
369  << " Ekin(MeV)= " << kineticEnergy
370  << " delEkin(MeV)= " << deltaKinEnergy
371  << G4endl;
372  }
373  if(rndmEngineMod->flat() > grej) { return; }
374  }
375 
376  G4ThreeVector deltaDirection;
377 
379 
380  const G4Material* mat = couple->GetMaterial();
382 
383  deltaDirection =
384  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
385 
386  } else {
387 
388  G4double deltaMomentum =
389  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
390  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
391  (deltaMomentum * dp->GetTotalMomentum());
392  if(cost > 1.0) { cost = 1.0; }
393  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
394 
395  G4double phi = twopi*rndmEngineMod->flat();
396 
397  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
398  deltaDirection.rotateUz(dp->GetMomentumDirection());
399  }
400  /*
401  G4cout << "### G4LindhardSorensenIonModel "
402  << dp->GetDefinition()->GetParticleName()
403  << " Ekin(MeV)= " << kineticEnergy
404  << " delEkin(MeV)= " << deltaKinEnergy
405  << " tmin(MeV)= " << minKinEnergy
406  << " tmax(MeV)= " << maxKinEnergy
407  << " dir= " << dp->GetMomentumDirection()
408  << " dirDelta= " << deltaDirection
409  << G4endl;
410  */
411  // create G4DynamicParticle object for delta ray
412  G4DynamicParticle* delta =
413  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
414 
415  vdp->push_back(delta);
416 
417  // Change kinematics of primary particle
418  kineticEnergy -= deltaKinEnergy;
419  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
420  finalP = finalP.unit();
421 
424 }
425 
426 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
427 
428 G4double
430  G4double kinEnergy)
431 {
432  // here particle type is checked for any method
433  SetParticle(pd);
434  G4double tau = kinEnergy/mass;
435  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
436  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
437  return std::min(tmax,tlimit);
438 }
439 
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Float_t x
Definition: compare.C:6
void set(double x, double y, double z)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetDeltaL(G4int Z, G4double gamma) const
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
G4LindhardSorensenIonModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="LindhardSorensen")
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
G4double BarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
const char * p
Definition: xmltok.h:285
static G4LindhardSorensenData * lsdata
const G4String & GetParticleName() const
G4double DensityCorrection(G4double x)
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)
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
G4double GetPDGMass() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
G4ParticleChangeForLoss * fParticleChange
Float_t Z
void SetParticle(const G4ParticleDefinition *p)
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:116
static constexpr double c_squared
const G4ParticleDefinition * particle
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
static constexpr double twopi
Definition: G4SIunits.hh:76
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Float_t f1
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:582
Float_t mat
G4double GetMeanExcitationEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4double GetPDGMagneticMoment() const
Hep3Vector unit() const
G4ThreeVector GetMomentum() const
static constexpr double eplus
Definition: G4SIunits.hh:199
int G4lrint(double ad)
Definition: templates.hh:151
int G4int
Definition: G4Types.hh:78
G4double GetA27(G4int Z) const
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) override
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()
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *couple, const G4DynamicParticle *dp, G4double &eloss, G4double &, G4double length) override
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:575
const G4Material * GetMaterial() const
G4double GetTotalMomentum() const
Float_t x2[n_points_geant4]
Definition: compare.C:26
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy) override
#define DBL_MAX
Definition: templates.hh:83
static constexpr double GeV
Definition: G4SIunits.hh:217
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) override
G4double GetElectronDensity() const
Definition: G4Material.hh:218
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4NistManager * Instance()