| Geant4 Cross Reference |
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 #include "G4BetheBlochModel.hh"
58 #include "Randomize.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4NistManager.hh"
62 #include "G4Electron.hh"
63 #include "G4LossTableManager.hh"
64 #include "G4EmCorrections.hh"
65 #include "G4EmParameters.hh"
66 #include "G4ParticleChangeForLoss.hh"
67 #include "G4ICRU90StoppingData.hh"
68 #include "G4Log.hh"
69 #include "G4DeltaAngle.hh"
70 #include <vector>
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 G4BetheBlochModel::G4BetheBlochModel(const G4ParticleDefinition*,
75 const G4String& nam)
76 : G4VEmModel(nam),
77 twoln10(2.0*G4Log(10.0)),
78 fAlphaTlimit(1*CLHEP::GeV),
79 fProtonTlimit(10*CLHEP::GeV)
80 {
81 theElectron = G4Electron::Electron();
82 corr = G4LossTableManager::Instance()->EmCorrections();
83 nist = G4NistManager::Instance();
84 SetLowEnergyLimit(2.0*CLHEP::MeV);
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89 G4BetheBlochModel::~G4BetheBlochModel() = default;
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92
93 void G4BetheBlochModel::Initialise(const G4ParticleDefinition* p,
94 const G4DataVector&)
95 {
96 if(p != particle) { SetupParameters(p); }
97
98 // always false before the run
99 SetDeexcitationFlag(false);
100
101 // initialisation once
102 if(nullptr == fParticleChange) {
103 const G4String& pname = particle->GetParticleName();
104 if(G4EmParameters::Instance()->UseICRU90Data() &&
105 (pname == "proton" || pname == "GenericIon" || pname == "alpha")) {
106 fICRU90 = nist->GetICRU90StoppingData();
107 }
108 if (pname == "GenericIon") {
109 isIon = true;
110 } else if (pname == "alpha") {
111 isAlpha = true;
112 } else if (particle->GetPDGCharge() > 1.1*CLHEP::eplus) {
113 isIon = true;
114 }
115
116 fParticleChange = GetParticleChangeForLoss();
117 if(UseAngularGeneratorFlag() && nullptr == GetAngularDistribution()) {
118 SetAngularDistribution(new G4DeltaAngle());
119 }
120 }
121 // initialisation for each new run
122 if(IsMaster() && nullptr != fICRU90) {
123 fICRU90->Initialise();
124 }
125 }
126
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128
129 G4double G4BetheBlochModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
130 const G4Material* mat,
131 G4double kinEnergy)
132 {
133 // this method is called only for ions, so no check if it is an ion
134 if(isAlpha) { return 1.0; }
135 chargeSquare = corr->EffectiveChargeSquareRatio(p, mat, kinEnergy);
136 return chargeSquare;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140
141 G4double G4BetheBlochModel::GetParticleCharge(const G4ParticleDefinition* p,
142 const G4Material* mat,
143 G4double kineticEnergy)
144 {
145 // this method is called only for ions, so no check if it is an ion
146 return corr->GetParticleCharge(p, mat, kineticEnergy);
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150
151 void G4BetheBlochModel::SetupParameters(const G4ParticleDefinition* p)
152 {
153 particle = p;
154 mass = particle->GetPDGMass();
155 spin = particle->GetPDGSpin();
156 G4double q = particle->GetPDGCharge()*inveplus;
157 chargeSquare = q*q;
158 ratio = electron_mass_c2/mass;
159 constexpr G4double aMag = 1./(0.5*eplus*CLHEP::hbar_Planck*CLHEP::c_squared);
160 G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
161 magMoment2 = magmom*magmom - 1.0;
162 formfact = 0.0;
163 tlimit = DBL_MAX;
164 if(particle->GetLeptonNumber() == 0) {
165 G4double x = 0.8426*CLHEP::GeV;
166 if(spin == 0.0 && mass < CLHEP::GeV) { x = 0.736*CLHEP::GeV; }
167 else if (mass > CLHEP::GeV) {
168 G4int iz = G4lrint(std::abs(q));
169 if(iz > 1) { x /= nist->GetA27(iz); }
170 }
171 formfact = 2.0*CLHEP::electron_mass_c2/(x*x);
172 tlimit = 2.0/formfact;
173 }
174 }
175
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177
178 G4double G4BetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
179 const G4MaterialCutsCouple* couple)
180 {
181 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185
186 G4double
187 G4BetheBlochModel::ComputeCrossSectionPerElectron(const G4ParticleDefinition* p,
188 G4double kineticEnergy,
189 G4double cut,
190 G4double maxKinEnergy)
191 {
192 G4double cross = 0.0;
193 const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
194 const G4double cutEnergy = std::min(std::min(cut,tmax), tlimit);
195 const G4double maxEnergy = std::min(tmax, maxKinEnergy);
196 if(cutEnergy < maxEnergy) {
197
198 G4double totEnergy = kineticEnergy + mass;
199 G4double energy2 = totEnergy*totEnergy;
200 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
201
202 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
204
205 // +term for spin=1/2 particle
206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
207
208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
209 }
210
211 // G4cout << "BB: e= " << kineticEnergy << " tmin= " << cutEnergy
212 // << " tmax= " << tmax << " cross= " << cross << G4endl;
213
214 return cross;
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218
219 G4double G4BetheBlochModel::ComputeCrossSectionPerAtom(
220 const G4ParticleDefinition* p,
221 G4double kinEnergy,
222 G4double Z, G4double,
223 G4double cutEnergy,
224 G4double maxEnergy)
225 {
226 return Z*ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,maxEnergy);
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230
231 G4double G4BetheBlochModel::CrossSectionPerVolume(
232 const G4Material* mat,
233 const G4ParticleDefinition* p,
234 G4double kinEnergy,
235 G4double cutEnergy,
236 G4double maxEnergy)
237 {
238 G4double sigma = mat->GetElectronDensity()
239 *ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,maxEnergy);
240 if(isAlpha) {
241 sigma *= corr->EffectiveChargeSquareRatio(p,mat,kinEnergy)/chargeSquare;
242 }
243 return sigma;
244 }
245
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247
248 G4double G4BetheBlochModel::ComputeDEDXPerVolume(const G4Material* material,
249 const G4ParticleDefinition* p,
250 G4double kineticEnergy,
251 G4double cut)
252 {
253 const G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
254 // projectile formfactor limit energy loss
255 const G4double cutEnergy = std::min(std::min(cut,tmax), tlimit);
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 /*
270 G4cout << "### DEDX ICRI90:" << (nullptr != fICRU90)
271 << " Ekin=" << kineticEnergy
272 << " " << p->GetParticleName()
273 << " q2=" << chargeSquare
274 << " inside " << material->GetName() << G4endl;
275 */
276 if(nullptr != fICRU90 && kineticEnergy < fProtonTlimit) {
277 if(material != currentMaterial) {
278 currentMaterial = material;
279 baseMaterial = material->GetBaseMaterial()
280 ? material->GetBaseMaterial() : material;
281 iICRU90 = fICRU90->GetIndex(baseMaterial);
282 }
283 if(iICRU90 >= 0) {
284 G4double dedx = 0.0;
285 // only for alpha
286 if(isAlpha) {
287 if(kineticEnergy <= fAlphaTlimit) {
288 dedx = fICRU90->GetElectronicDEDXforAlpha(iICRU90, kineticEnergy);
289 } else {
290 const G4double e = kineticEnergy*CLHEP::proton_mass_c2/mass;
291 dedx = fICRU90->GetElectronicDEDXforProton(iICRU90, e)*chargeSquare;
292 }
293 } else {
294 dedx = fICRU90->GetElectronicDEDXforProton(iICRU90, kineticEnergy)
295 *chargeSquare;
296 }
297 dedx *= material->GetDensity();
298 if(cutEnergy < tmax) {
299 dedx += (G4Log(xc) + (1.0 - xc)*beta2)*CLHEP::twopi_mc2_rcl2
300 *(eDensity*chargeSquare/beta2);
301 }
302 //G4cout << " iICRU90=" << iICRU90 << " dedx=" << dedx << G4endl;
303 if(dedx > 0.0) { return dedx; }
304 }
305 }
306 // general Bethe-Bloch formula
307 G4double dedx = G4Log(2.0*CLHEP::electron_mass_c2*bg2*cutEnergy/eexc2)
308 - (1.0 + xc)*beta2;
309
310 if(0.0 < spin) {
311 G4double del = 0.5*cutEnergy/(kineticEnergy + mass);
312 dedx += del*del;
313 }
314
315 // density correction
316 G4double x = G4Log(bg2)/twoln10;
317 dedx -= material->GetIonisation()->DensityCorrection(x);
318
319 // shell correction
320 dedx -= 2.0*corr->ShellCorrection(p,material,kineticEnergy);
321
322 // now compute the total ionization loss
323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*eDensity/beta2;
324
325 //High order correction different for hadrons and ions
326 if(isIon) {
327 dedx += corr->IonBarkasCorrection(p,material,kineticEnergy);
328 } else {
329 dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
330 }
331
332 dedx = std::max(dedx, 0.0);
333 /*
334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP::MeV << " dedx= " << dedx
335 << " " << material->GetName() << G4endl;
336 */
337 return dedx;
338 }
339
340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341
342 void G4BetheBlochModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
343 const G4DynamicParticle* dp,
344 const G4double& /*length*/,
345 G4double& eloss)
346 {
347 // no correction for alpha
348 if(isAlpha) { return; }
349
350 // no correction at the last step or at small step
351 const G4double preKinEnergy = dp->GetKineticEnergy();
352 if(eloss >= preKinEnergy || eloss < preKinEnergy*0.05) { return; }
353
354 // corrections for all charged particles with Q > 1
355 const G4ParticleDefinition* p = dp->GetDefinition();
356 if(p != particle) { SetupParameters(p); }
357 if(!isIon) { return; }
358
359 // effective energy and charge at a step
360 const G4double e = std::max(preKinEnergy - eloss*0.5, preKinEnergy*0.5);
361 const G4Material* mat = couple->GetMaterial();
362 const G4double q20 = corr->EffectiveChargeSquareRatio(p, mat, preKinEnergy);
363 const G4double q2 = corr->EffectiveChargeSquareRatio(p, mat, e);
364 const G4double qfactor = q2/q20;
365
366 /*
367 G4cout << "G4BetheBlochModel::CorrectionsAlongStep: Epre(MeV)="
368 << preKinEnergy << " Eeff(MeV)=" << e
369 << " eloss=" << eloss << " elossnew=" << eloss*qfactor
370 << " qfactor=" << qfactor << " Qpre=" << q20
371 << p->GetParticleName() <<G4endl;
372 */
373 eloss *= qfactor;
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377
378 void G4BetheBlochModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
379 const G4MaterialCutsCouple* couple,
380 const G4DynamicParticle* dp,
381 G4double cut,
382 G4double maxEnergy)
383 {
384 G4double kinEnergy = dp->GetKineticEnergy();
385 const G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(), kinEnergy);
386 const G4double minKinEnergy = std::min(cut, tmax);
387 const G4double maxKinEnergy = std::min(maxEnergy, tmax);
388 if(minKinEnergy >= maxKinEnergy) { return; }
389
390 //G4cout << "G4BetheBlochModel::SampleSecondaries Emin= " << minKinEnergy
391 // << " Emax= " << maxKinEnergy << G4endl;
392
393 const G4double totEnergy = kinEnergy + mass;
394 const G4double etot2 = totEnergy*totEnergy;
395 const G4double beta2 = kinEnergy*(kinEnergy + 2.0*mass)/etot2;
396
397 G4double deltaKinEnergy, f;
398 G4double f1 = 0.0;
399 G4double fmax = 1.0;
400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
401
402 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
403 G4double rndm[2];
404
405 // sampling without nuclear size effect
406 do {
407 rndmEngineMod->flatArray(2, rndm);
408 deltaKinEnergy = minKinEnergy*maxKinEnergy
409 /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
410
411 f = 1.0 - beta2*deltaKinEnergy/tmax;
412 if( 0.0 < spin ) {
413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
414 f += f1;
415 }
416
417 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
418 } while( fmax*rndm[1] > f);
419
420 // projectile formfactor - suppresion of high energy
421 // delta-electron production at high energy
422
423 G4double x = formfact*deltaKinEnergy;
424 if(x > 1.e-6) {
425
426 G4double x1 = 1.0 + x;
427 G4double grej = 1.0/(x1*x1);
428 if( 0.0 < spin ) {
429 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
431 }
432 if(grej > 1.1) {
433 G4cout << "### G4BetheBlochModel WARNING: grej= " << grej
434 << " " << dp->GetDefinition()->GetParticleName()
435 << " Ekin(MeV)= " << kinEnergy
436 << " delEkin(MeV)= " << deltaKinEnergy
437 << G4endl;
438 }
439 if(rndmEngineMod->flat() > grej) { return; }
440 }
441
442 G4ThreeVector deltaDirection;
443
444 if(UseAngularGeneratorFlag()) {
445 const G4Material* mat = couple->GetMaterial();
446 deltaDirection =
447 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy,
448 SelectRandomAtomNumber(mat),
449 mat);
450 } else {
451
452 G4double deltaMomentum =
453 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
454 G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
455 (deltaMomentum * dp->GetTotalMomentum());
456 cost = std::min(cost, 1.0);
457 const G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
458 const G4double phi = twopi*rndmEngineMod->flat();
459
460 deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
461 deltaDirection.rotateUz(dp->GetMomentumDirection());
462 }
463 /*
464 G4cout << "### G4BetheBlochModel "
465 << dp->GetDefinition()->GetParticleName()
466 << " Ekin(MeV)= " << kinEnergy
467 << " delEkin(MeV)= " << deltaKinEnergy
468 << " tmin(MeV)= " << minKinEnergy
469 << " tmax(MeV)= " << maxKinEnergy
470 << " dir= " << dp->GetMomentumDirection()
471 << " dirDelta= " << deltaDirection
472 << G4endl;
473 */
474 // create G4DynamicParticle object for delta ray
475 auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
476
477 vdp->push_back(delta);
478
479 // Change kinematics of primary particle
480 kinEnergy -= deltaKinEnergy;
481 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
482 finalP = finalP.unit();
483
484 fParticleChange->SetProposedKineticEnergy(kinEnergy);
485 fParticleChange->SetProposedMomentumDirection(finalP);
486 }
487
488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
489
490 G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
491 G4double kinEnergy)
492 {
493 // here particle type is checked for the case,
494 // when this model is shared between particles
495 if(pd != particle) { SetupParameters(pd); }
496 G4double tau = kinEnergy/mass;
497 return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) /
498 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
499 }
500
501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502