ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LindhardSorensenIonModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class header file
30 //
31 //
32 // File name: G4LindhardSorensenIonModel
33 //
34 // Author: Alexander Bagulya & Vladimir Ivanchenko
35 //
36 // Creation date: 16.04.2018
37 //
38 //
39 // -------------------------------------------------------------------
40 //
41 
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
46 #include "G4LindhardSorensenIonModel.hh"
47 #include "Randomize.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Electron.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4EmCorrections.hh"
53 #include "G4ParticleChangeForLoss.hh"
54 #include "G4Log.hh"
55 #include "G4DeltaAngle.hh"
56 #include "G4LindhardSorensenData.hh"
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 
60 using namespace std;
61 
62 G4LindhardSorensenData* G4LindhardSorensenIonModel::lsdata = nullptr;
63 
64 G4LindhardSorensenIonModel::G4LindhardSorensenIonModel(const G4ParticleDefinition*,
65  const G4String& nam)
66  : G4VEmModel(nam),
67  particle(nullptr),
68  tlimit(DBL_MAX),
69  twoln10(2.0*G4Log(10.0))
70 {
71  fParticleChange = nullptr;
72  theElectron = G4Electron::Electron();
73  SetParticle(theElectron);
74  corr = G4LossTableManager::Instance()->EmCorrections();
75  nist = G4NistManager::Instance();
76  SetLowEnergyLimit(2.0*MeV);
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
81 G4LindhardSorensenIonModel::~G4LindhardSorensenIonModel()
82 {}
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
86 void G4LindhardSorensenIonModel::Initialise(const G4ParticleDefinition* p,
87  const G4DataVector&)
88 {
89  SetParticle(p);
90 
91  //G4cout << "G4LindhardSorensenIonModel::Initialise for " << p->GetParticleName()
92  // << G4endl;
93 
94  // always false before the run
95  SetDeexcitationFlag(false);
96 
97  if(nullptr == fParticleChange) {
98  fParticleChange = GetParticleChangeForLoss();
99  if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
100  SetAngularDistribution(new G4DeltaAngle());
101  }
102  }
103  if(IsMaster() && !lsdata) {
104  lsdata = new G4LindhardSorensenData();
105  }
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 
110 G4double
111 G4LindhardSorensenIonModel::GetChargeSquareRatio(const G4ParticleDefinition*,
112  const G4Material*,
113  G4double)
114 {
115  return chargeSquare;
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
120 G4double
121 G4LindhardSorensenIonModel::GetParticleCharge(const G4ParticleDefinition*,
122  const G4Material*,
123  G4double)
124 {
125  return charge;
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
130 void G4LindhardSorensenIonModel::SetupParameters()
131 {
132  mass = particle->GetPDGMass();
133  spin = particle->GetPDGSpin();
134  charge = particle->GetPDGCharge()*inveplus;
135  Zin = G4lrint(charge);
136  chargeSquare = charge*charge;
137  ratio = electron_mass_c2/mass;
138  static const G4double aMag = 1./(0.5*eplus*hbar_Planck*c_squared);
139  G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
140  magMoment2 = magmom*magmom - 1.0;
141  if(Zin <= 1) {
142  formfact = (spin == 0.0 && mass < GeV) ? 1.181e-6 : 1.548e-6;
143  } else {
144  G4double x = nist->GetA27(Zin);
145  formfact = 3.969e-6*x*x;
146  }
147  tlimit = std::sqrt(0.414/formfact +
148  electron_mass_c2*electron_mass_c2) - electron_mass_c2;
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152 
153 G4double G4LindhardSorensenIonModel::MinEnergyCut(const G4ParticleDefinition*,
154  const G4MaterialCutsCouple* couple)
155 {
156  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
157 }
158 
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
160 
161 G4double
162 G4LindhardSorensenIonModel::ComputeCrossSectionPerElectron(
163  const G4ParticleDefinition* p,
164  G4double kineticEnergy,
165  G4double cutEnergy,
166  G4double maxKinEnergy)
167 {
168  G4double cross = 0.0;
169  // take into account formfactor
170  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy),tlimit);
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 
195 G4double G4LindhardSorensenIonModel::ComputeCrossSectionPerAtom(
196  const G4ParticleDefinition* p,
197  G4double kineticEnergy,
198  G4double Z, G4double,
199  G4double cutEnergy,
200  G4double maxEnergy)
201 {
202  return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
206 
207 G4double G4LindhardSorensenIonModel::CrossSectionPerVolume(
208  const G4Material* material,
209  const G4ParticleDefinition* p,
210  G4double kineticEnergy,
211  G4double cutEnergy,
212  G4double maxEnergy)
213 {
214  return material->GetElectronDensity()
215  *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 G4double
221 G4LindhardSorensenIonModel::ComputeDEDXPerVolume(const G4Material* material,
222  const G4ParticleDefinition* p,
223  G4double kineticEnergy,
224  G4double cut)
225 {
226  // formfactor is taken into account in CorrectionsAlongStep(..)
227  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
228  G4double cutEnergy = std::min(cut,tmax);
229 
230  G4double tau = kineticEnergy/mass;
231  G4double gam = tau + 1.0;
232  G4double bg2 = tau * (tau+2.0);
233  G4double beta2 = bg2/(gam*gam);
234 
235  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
236  G4double eexc2 = eexc*eexc;
237 
238  G4double eDensity = material->GetElectronDensity();
239 
240  G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
241  - (1.0 + cutEnergy/tmax)*beta2;
242 
243  if(0.0 < spin) {
244  G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
245  dedx += del*del;
246  }
247 
248  // density correction
249  G4double x = G4Log(bg2)/twoln10;
250  dedx -= material->GetIonisation()->DensityCorrection(x);
251 
252  // shell correction
253  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
254 
255  //High order correction different for hadrons and ions
256  dedx += 2.0*corr->BarkasCorrection(p,material,kineticEnergy);
257  dedx = std::max(dedx, 0.0);
258 
259  // now compute the total ionization loss
260  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
261 
262  //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
263  // << " " << material->GetName() << G4endl;
264 
265  return dedx;
266 }
267 
268 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
269 
270 void
271 G4LindhardSorensenIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
272  const G4DynamicParticle* dp,
273  G4double& eloss,
274  G4double&,
275  G4double length)
276 {
277  const G4ParticleDefinition* p = dp->GetDefinition();
278  SetParticle(p);
279  G4double eDensity = couple->GetMaterial()->GetElectronDensity();
280  G4double preKinEnergy = dp->GetKineticEnergy();
281  G4double e = preKinEnergy - eloss*0.5;
282  GetModelOfFluctuations()->SetParticleAndCharge(p, chargeSquare);
283 
284  G4double tau = e/mass;
285  G4double gam = tau + 1.0;
286  G4double beta2 = tau * (tau+2.0)/(gam*gam);
287  G4double deltaL0 =
288  2.0*corr->BarkasCorrection (p, couple->GetMaterial(), e)*(charge-1.)/charge;
289  G4double deltaL = lsdata->GetDeltaL(Zin, gam);
290 
291  G4double elossnew =
292  eloss + twopi_mc2_rcl2*chargeSquare*eDensity*(deltaL+deltaL0)*length/beta2;
293  /*
294  G4cout << "G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)= "
295  << preKinEnergy/GeV << " eloss(MeV)= " << eloss
296  << " L= " << eloss*beta2/(twopi_mc2_rcl2*chargeSquare*eDensity*length)
297  << " dL0= " << deltaL0
298  << " dL= " << deltaL << G4endl;
299  */
300  if(elossnew > preKinEnergy) { elossnew = preKinEnergy; }
301  else if(elossnew < 0.0) { elossnew = eloss*0.5; }
302 
303  eloss = elossnew;
304 }
305 
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307 
308 void G4LindhardSorensenIonModel::SampleSecondaries(
309  vector<G4DynamicParticle*>* vdp,
310  const G4MaterialCutsCouple* couple,
311  const G4DynamicParticle* dp,
312  G4double minKinEnergy,
313  G4double maxEnergy)
314 {
315  G4double kineticEnergy = dp->GetKineticEnergy();
316  // take into account formfactor
317  G4double tmax =
318  std::min(MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy),tlimit);
319 
320  G4double maxKinEnergy = std::min(maxEnergy,tmax);
321  if(minKinEnergy >= maxKinEnergy) { return; }
322 
323  //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= "
324  // << minKinEnergy << " Emax= " << maxKinEnergy << G4endl;
325 
326  G4double totEnergy = kineticEnergy + mass;
327  G4double etot2 = totEnergy*totEnergy;
328  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
329 
330  G4double deltaKinEnergy, f;
331  G4double f1 = 0.0;
332  G4double fmax = 1.0;
333  if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
334 
335  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
336  G4double rndm[2];
337 
338  // sampling without nuclear size effect
339  do {
340  rndmEngineMod->flatArray(2, rndm);
341  deltaKinEnergy = minKinEnergy*maxKinEnergy
342  /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
343 
344  f = 1.0 - beta2*deltaKinEnergy/tmax;
345  if( 0.0 < spin ) {
346  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
347  f += f1;
348  }
349 
350  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
351  } while( fmax*rndm[1] > f);
352 
353  // projectile formfactor - suppresion of high energy
354  // delta-electron production at high energy
355 
356  G4double x = formfact*deltaKinEnergy*(deltaKinEnergy + 2*electron_mass_c2);
357  if(x > 1.e-6) {
358 
359  G4double x1 = 1.0 + x;
360  G4double grej = 1.0/(x1*x1);
361  if( 0.0 < spin ) {
362  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
363  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
364  }
365  if(grej > 1.1) {
366  G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej
367  << " " << dp->GetDefinition()->GetParticleName()
368  << " Ekin(MeV)= " << kineticEnergy
369  << " delEkin(MeV)= " << deltaKinEnergy
370  << G4endl;
371  }
372  if(rndmEngineMod->flat() > grej) { return; }
373  }
374 
375  G4ThreeVector deltaDirection;
376 
377  if(UseAngularGeneratorFlag()) {
378 
379  const G4Material* mat = couple->GetMaterial();
380  G4int Z = SelectRandomAtomNumber(mat);
381 
382  deltaDirection =
383  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
384 
385  } else {
386 
387  G4double deltaMomentum =
388  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
389  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
390  (deltaMomentum * dp->GetTotalMomentum());
391  if(cost > 1.0) { cost = 1.0; }
392  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
393 
394  G4double phi = twopi*rndmEngineMod->flat();
395 
396  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
397  deltaDirection.rotateUz(dp->GetMomentumDirection());
398  }
399  /*
400  G4cout << "### G4LindhardSorensenIonModel "
401  << dp->GetDefinition()->GetParticleName()
402  << " Ekin(MeV)= " << kineticEnergy
403  << " delEkin(MeV)= " << deltaKinEnergy
404  << " tmin(MeV)= " << minKinEnergy
405  << " tmax(MeV)= " << maxKinEnergy
406  << " dir= " << dp->GetMomentumDirection()
407  << " dirDelta= " << deltaDirection
408  << G4endl;
409  */
410  // create G4DynamicParticle object for delta ray
411  G4DynamicParticle* delta =
412  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
413 
414  vdp->push_back(delta);
415 
416  // Change kinematics of primary particle
417  kineticEnergy -= deltaKinEnergy;
418  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
419  finalP = finalP.unit();
420 
421  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
422  fParticleChange->SetProposedMomentumDirection(finalP);
423 }
424 
425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426 
427 G4double
428 G4LindhardSorensenIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
429  G4double kinEnergy)
430 {
431  // here particle type is checked for any method
432  SetParticle(pd);
433  G4double tau = kinEnergy/mass;
434  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
435  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
436  // formfactor is not taken into account
437  return tmax;
438 }
439 
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......