ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BetheBlochModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class header file
29 //
30 //
31 // File name: G4BetheBlochModel
32 //
33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
34 //
35 // Creation date: 03.01.2002
36 //
37 // Modifications:
38 //
39 // 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
40 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
41 // 27-01-03 Make models region aware (V.Ivanchenko)
42 // 13-02-03 Add name (V.Ivanchenko)
43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)
44 // 11-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
45 // 11-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
46 // 12-02-06 move G4LossTableManager::Instance()->EmCorrections()
47 // in constructor (mma)
48 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
49 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
50 //
51 // -------------------------------------------------------------------
52 //
53 
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
58 #include "G4BetheBlochModel.hh"
59 #include "Randomize.hh"
60 #include "G4PhysicalConstants.hh"
61 #include "G4SystemOfUnits.hh"
62 #include "G4Electron.hh"
63 #include "G4LossTableManager.hh"
64 #include "G4EmCorrections.hh"
65 #include "G4ParticleChangeForLoss.hh"
66 #include "G4ICRU90StoppingData.hh"
67 #include "G4Log.hh"
68 #include "G4DeltaAngle.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 using namespace std;
73 
74 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition*,
75  const G4String& nam)
76  : G4VEmModel(nam),
77  particle(nullptr),
78  fICRU90(nullptr),
79  currentMaterial(nullptr),
80  baseMaterial(nullptr),
81  tlimit(DBL_MAX),
82  twoln10(2.0*G4Log(10.0)),
83  fAlphaTlimit(CLHEP::GeV),
84  fProtonTlimit(10*CLHEP::GeV),
85  iICRU90(-1),
86  isIon(false)
87 {
88  fParticleChange = nullptr;
89  theElectron = G4Electron::Electron();
90  SetParticle(theElectron);
91  corr = G4LossTableManager::Instance()->EmCorrections();
92  nist = G4NistManager::Instance();
93  SetLowEnergyLimit(2.0*MeV);
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 
98 G4BetheBlochModel::~G4BetheBlochModel()
99 {}
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
103 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
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(IsMaster() && G4EmParameters::Instance()->UseICRU90Data()) {
117  if(!fICRU90) { fICRU90 = nist->GetICRU90StoppingData(); }
118  else if(particle->GetPDGMass() < GeV) { fICRU90->Initialise(); }
119  }
120 
121  if(nullptr == fParticleChange) {
122  fParticleChange = GetParticleChangeForLoss();
123  if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
124  SetAngularDistribution(new G4DeltaAngle());
125  }
126  }
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
131 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
132  const G4Material* mat,
133  G4double kineticEnergy)
134 {
135  // this method is called only for ions
136  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
137  corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
138  return corrFactor;
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 
143 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
144  const G4Material* mat,
145  G4double kineticEnergy)
146 {
147  //G4cout<<"G4BetheBlochModel::GetParticleCharge e= "<<kineticEnergy <<
148  // " q= " << corr->GetParticleCharge(p,mat,kineticEnergy) <<G4endl;
149  // this method is called only for ions, so no check if it is an ion
150  return corr->GetParticleCharge(p,mat,kineticEnergy);
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
154 
155 void G4BetheBlochModel::SetupParameters()
156 {
157  mass = particle->GetPDGMass();
158  spin = particle->GetPDGSpin();
159  G4double q = particle->GetPDGCharge()*inveplus;
160  chargeSquare = q*q;
161  corrFactor = chargeSquare;
162  ratio = electron_mass_c2/mass;
163  static const G4double aMag = 1./(0.5*eplus*hbar_Planck*c_squared);
164  G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
165  magMoment2 = magmom*magmom - 1.0;
166  formfact = 0.0;
167  tlimit = DBL_MAX;
168  if(particle->GetLeptonNumber() == 0) {
169  G4int iz = G4lrint(q);
170  if(iz <= 1) {
171  formfact = (spin == 0.0 && mass < GeV) ? 1.181e-6 : 1.548e-6;
172  } else {
173  G4double x = nist->GetA27(iz);
174  formfact = 3.969e-6*x*x;
175  }
176  tlimit = std::sqrt(0.414/formfact +
177  electron_mass_c2*electron_mass_c2) - electron_mass_c2;
178  }
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
184  const G4MaterialCutsCouple* couple)
185 {
186  return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
191 G4double
192 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
193  G4double kineticEnergy,
194  G4double cutEnergy,
195  G4double maxKinEnergy)
196 {
197  G4double cross = 0.0;
198  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
199  G4double maxEnergy = std::min(tmax,maxKinEnergy);
200  if(cutEnergy < maxEnergy) {
201 
202  G4double totEnergy = kineticEnergy + mass;
203  G4double energy2 = totEnergy*totEnergy;
204  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
205 
206  cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
207  - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
208 
209  // +term for spin=1/2 particle
210  if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
211 
212  cross *= twopi_mc2_rcl2*chargeSquare/beta2;
213  }
214 
215  // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
216  // << " tmax= " << tmax << " cross= " << cross << G4endl;
217 
218  return cross;
219 }
220 
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
222 
223 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
224  const G4ParticleDefinition* p,
225  G4double kineticEnergy,
226  G4double Z, G4double,
227  G4double cutEnergy,
228  G4double maxEnergy)
229 {
230  return
231  Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 
236 G4double G4BetheBlochModel::CrossSectionPerVolume(
237  const G4Material* material,
238  const G4ParticleDefinition* p,
239  G4double kineticEnergy,
240  G4double cutEnergy,
241  G4double maxEnergy)
242 {
243  return material->GetElectronDensity()
244  *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
245 }
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248 
249 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
250  const G4ParticleDefinition* p,
251  G4double kineticEnergy,
252  G4double cut)
253 {
254  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
255  G4double cutEnergy = std::min(cut,tmax);
256 
257  G4double tau = kineticEnergy/mass;
258  G4double gam = tau + 1.0;
259  G4double bg2 = tau * (tau+2.0);
260  G4double beta2 = bg2/(gam*gam);
261  G4double xc = cutEnergy/tmax;
262 
263  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
264  G4double eexc2 = eexc*eexc;
265 
266  G4double eDensity = material->GetElectronDensity();
267 
268  // added ICRU90 stopping data for limited list of materials
269  if(fICRU90) {
270  if(material != currentMaterial) {
271  currentMaterial = material;
272  baseMaterial = material->GetBaseMaterial()
273  ? material->GetBaseMaterial() : material;
274  iICRU90 = fICRU90->GetIndex(baseMaterial);
275  }
276  if(iICRU90 >= 0) {
277  G4double e = kineticEnergy*proton_mass_c2/mass;
278  G4double dedx = 0.0;
279  if(chargeSquare > 1.1 && e < fAlphaTlimit) {
280  dedx = fICRU90->GetElectronicDEDXforAlpha(iICRU90, e)
281  *material->GetDensity()*0.25;
282  } else if(chargeSquare < 1.1 && e < fProtonTlimit) {
283  dedx = fICRU90->GetElectronicDEDXforProton(iICRU90, e)
284  *material->GetDensity();
285  }
286  if(cutEnergy < tmax) {
287  dedx += (G4Log(xc) + (1.0 - xc)*beta2)*twopi_mc2_rcl2
288  *eDensity/beta2;
289  return std::max(chargeSquare*dedx, 0.0);
290  }
291  }
292  }
293  // general Bethe-Bloch formula
294  G4double dedx = G4Log(2.0*electron_mass_c2*bg2*cutEnergy/eexc2)
295  - (1.0 + xc)*beta2;
296 
297  if(0.0 < spin) {
298  G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
299  dedx += del*del;
300  }
301 
302  // density correction
303  G4double x = G4Log(bg2)/twoln10;
304  dedx -= material->GetIonisation()->DensityCorrection(x);
305 
306  // shell correction
307  dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
308 
309  // now compute the total ionization loss
310  dedx *= twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
311 
312  //High order correction different for hadrons and ions
313  if(isIon) {
314  dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
315  } else {
316  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
317  }
318 
319  dedx = std::max(dedx, 0.0);
320 
321  //G4cout << "E(MeV)= " << kineticEnergy/MeV << " dedx= " << dedx
322  // << " " << material->GetName() << G4endl;
323 
324  return dedx;
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328 
329 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
330  const G4DynamicParticle* dp,
331  G4double& eloss,
332  G4double&,
333  G4double length)
334 {
335  if(isIon) {
336  const G4Material* mat = couple->GetMaterial();
337  const G4ParticleDefinition* p = dp->GetDefinition();
338  G4double preKinEnergy = dp->GetKineticEnergy();
339  G4double e = preKinEnergy - eloss*0.5;
340  if(e < preKinEnergy*0.75) { e = preKinEnergy*0.75; }
341 
342  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
343  GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
344  G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
345 
346  // no high order correction for ICRU90 data
347  baseMaterial = mat->GetBaseMaterial() ? mat->GetBaseMaterial() : mat;
348  G4double highOrder = 0.0;
349  if(!fICRU90 || fICRU90->GetIndex(baseMaterial) < 0) {
350  highOrder = length*corr->IonHighOrderCorrections(p,couple,e);
351  }
352  G4double elossnew = eloss*qfactor + highOrder;
353  eloss = std::max(std::min(elossnew,preKinEnergy),eloss*0.5);
354  //G4cout << "G4BetheBlochModel::CorrectionsAlongStep: e= " << preKinEnergy
355  // << " qfactor= " << qfactor
356  // << " highOrder= " << highOrder << " ("
357  // << highOrder/eloss << ")" << G4endl;
358  }
359 }
360 
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
362 
363 void G4BetheBlochModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
364  const G4MaterialCutsCouple* couple,
365  const G4DynamicParticle* dp,
366  G4double minKinEnergy,
367  G4double maxEnergy)
368 {
369  G4double kineticEnergy = dp->GetKineticEnergy();
370  G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(),kineticEnergy);
371 
372  G4double maxKinEnergy = std::min(maxEnergy,tmax);
373  if(minKinEnergy >= maxKinEnergy) { return; }
374 
375  //G4cout << "G4BetheBlochModel::SampleSecondaries Emin= " << minKinEnergy
376  // << " Emax= " << maxKinEnergy << G4endl;
377 
378  G4double totEnergy = kineticEnergy + mass;
379  G4double etot2 = totEnergy*totEnergy;
380  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
381 
382  G4double deltaKinEnergy, f;
383  G4double f1 = 0.0;
384  G4double fmax = 1.0;
385  if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
386 
387  CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
388  G4double rndm[2];
389 
390  // sampling without nuclear size effect
391  do {
392  rndmEngineMod->flatArray(2, rndm);
393  deltaKinEnergy = minKinEnergy*maxKinEnergy
394  /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
395 
396  f = 1.0 - beta2*deltaKinEnergy/tmax;
397  if( 0.0 < spin ) {
398  f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
399  f += f1;
400  }
401 
402  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
403  } while( fmax*rndm[1] > f);
404 
405  // projectile formfactor - suppresion of high energy
406  // delta-electron production at high energy
407 
408  G4double x = formfact*deltaKinEnergy*(deltaKinEnergy + 2*electron_mass_c2);
409  if(x > 1.e-6) {
410 
411  G4double x1 = 1.0 + x;
412  G4double grej = 1.0/(x1*x1);
413  if( 0.0 < spin ) {
414  G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
415  grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
416  }
417  if(grej > 1.1) {
418  G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
419  << " " << dp->GetDefinition()->GetParticleName()
420  << " Ekin(MeV)= " << kineticEnergy
421  << " delEkin(MeV)= " << deltaKinEnergy
422  << G4endl;
423  }
424  if(rndmEngineMod->flat() > grej) { return; }
425  }
426 
427  G4ThreeVector deltaDirection;
428 
429  if(UseAngularGeneratorFlag()) {
430 
431  const G4Material* mat = couple->GetMaterial();
432  G4int Z = SelectRandomAtomNumber(mat);
433 
434  deltaDirection =
435  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
436 
437  } else {
438 
439  G4double deltaMomentum =
440  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
441  G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
442  (deltaMomentum * dp->GetTotalMomentum());
443  if(cost > 1.0) { cost = 1.0; }
444  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
445 
446  G4double phi = twopi*rndmEngineMod->flat();
447 
448  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
449  deltaDirection.rotateUz(dp->GetMomentumDirection());
450  }
451  /*
452  G4cout << "### G4BetheBlochModel "
453  << dp->GetDefinition()->GetParticleName()
454  << " Ekin(MeV)= " << kineticEnergy
455  << " delEkin(MeV)= " << deltaKinEnergy
456  << " tmin(MeV)= " << minKinEnergy
457  << " tmax(MeV)= " << maxKinEnergy
458  << " dir= " << dp->GetMomentumDirection()
459  << " dirDelta= " << deltaDirection
460  << G4endl;
461  */
462  // create G4DynamicParticle object for delta ray
463  G4DynamicParticle* delta =
464  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
465 
466  vdp->push_back(delta);
467 
468  // Change kinematics of primary particle
469  kineticEnergy -= deltaKinEnergy;
470  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
471  finalP = finalP.unit();
472 
473  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
474  fParticleChange->SetProposedMomentumDirection(finalP);
475 }
476 
477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
478 
479 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
480  G4double kinEnergy)
481 {
482  // here particle type is checked for any method
483  SetParticle(pd);
484  G4double tau = kinEnergy/mass;
485  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
486  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
487  return std::min(tmax,tlimit);
488 }
489 
490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......