| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // -------------------------------------------
27 //
28 // GEANT4 Class header file
29 //
30 //
31 // File name: G4BetheBlochModel
32 //
33 // Author: Vladimir Ivanchenko on base
34 //
35 // Creation date: 03.01.2002
36 //
37 // Modifications:
38 //
39 // 04-12-02 Fix problem of G4DynamicParticle c
40 // 23-12-02 Change interface in order to move
41 // 27-01-03 Make models region aware (V.Ivanch
42 // 13-02-03 Add name (V.Ivanchenko)
43 // 24-03-05 Add G4EmCorrections (V.Ivanchenko)
44 // 11-04-05 Major optimisation of internal int
45 // 11-02-06 ComputeCrossSectionPerElectron, Co
46 // 12-02-06 move G4LossTableManager::Instance(
47 // in constructor (mma)
48 // 12-08-08 Added methods GetParticleCharge, G
49 // CorrectionsAlongStep needed for io
50 //
51 // -------------------------------------------
52 //
53
54 //....oooOO0OOooo........oooOO0OOooo........oo
55 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
73
74 G4BetheBlochModel::G4BetheBlochModel(const G4P
75 const G4S
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()->EmCor
83 nist = G4NistManager::Instance();
84 SetLowEnergyLimit(2.0*CLHEP::MeV);
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4BetheBlochModel::~G4BetheBlochModel() = defa
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93 void G4BetheBlochModel::Initialise(const G4Par
94 const G4Dat
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->GetParti
104 if(G4EmParameters::Instance()->UseICRU90Da
105 (pname == "proton" || pname == "Generic
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*
113 isIon = true;
114 }
115
116 fParticleChange = GetParticleChangeForLoss
117 if(UseAngularGeneratorFlag() && nullptr ==
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........oo
128
129 G4double G4BetheBlochModel::GetChargeSquareRat
130
131
132 {
133 // this method is called only for ions, so n
134 if(isAlpha) { return 1.0; }
135 chargeSquare = corr->EffectiveChargeSquareRa
136 return chargeSquare;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oo
140
141 G4double G4BetheBlochModel::GetParticleCharge(
142
143
144 {
145 // this method is called only for ions, so n
146 return corr->GetParticleCharge(p, mat, kinet
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oo
150
151 void G4BetheBlochModel::SetupParameters(const
152 {
153 particle = p;
154 mass = particle->GetPDGMass();
155 spin = particle->GetPDGSpin();
156 G4double q = particle->GetPDGCharge()*invepl
157 chargeSquare = q*q;
158 ratio = electron_mass_c2/mass;
159 constexpr G4double aMag = 1./(0.5*eplus*CLHE
160 G4double magmom = particle->GetPDGMagneticMo
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 =
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*
172 tlimit = 2.0/formfact;
173 }
174 }
175
176 //....oooOO0OOooo........oooOO0OOooo........oo
177
178 G4double G4BetheBlochModel::MinEnergyCut(const
179 const
180 {
181 return couple->GetMaterial()->GetIonisation(
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186 G4double
187 G4BetheBlochModel::ComputeCrossSectionPerElect
188
189
190
191 {
192 G4double cross = 0.0;
193 const G4double tmax = MaxSecondaryEnergy(p,
194 const G4double cutEnergy = std::min(std::min
195 const G4double maxEnergy = std::min(tmax, ma
196 if(cutEnergy < maxEnergy) {
197
198 G4double totEnergy = kineticEnergy + mass;
199 G4double energy2 = totEnergy*totEnergy;
200 G4double beta2 = kineticEnergy*(kineti
201
202 cross = (maxEnergy - cutEnergy)/(cutEnergy
203 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
204
205 // +term for spin=1/2 particle
206 if( 0.0 < spin ) { cross += 0.5*(maxEnergy
207
208 cross *= CLHEP::twopi_mc2_rcl2*chargeSquar
209 }
210
211 // G4cout << "BB: e= " << kineticEnergy <<
212 // << " tmax= " << tmax << " cross=
213
214 return cross;
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oo
218
219 G4double G4BetheBlochModel::ComputeCrossSectio
220 con
221
222
223
224
225 {
226 return Z*ComputeCrossSectionPerElectron(p,ki
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oo
230
231 G4double G4BetheBlochModel::CrossSectionPerVol
232 con
233 con
234
235
236
237 {
238 G4double sigma = mat->GetElectronDensity()
239 *ComputeCrossSectionPerElectron(p,kinEnerg
240 if(isAlpha) {
241 sigma *= corr->EffectiveChargeSquareRatio(
242 }
243 return sigma;
244 }
245
246 //....oooOO0OOooo........oooOO0OOooo........oo
247
248 G4double G4BetheBlochModel::ComputeDEDXPerVolu
249
250
251
252 {
253 const G4double tmax = MaxSecondaryEnergy(p,
254 // projectile formfactor limit energy loss
255 const G4double cutEnergy = std::min(std::min
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()->
264 G4double eexc2 = eexc*eexc;
265
266 G4double eDensity = material->GetElectronDen
267
268 // added ICRU90 stopping data for limited li
269 /*
270 G4cout << "### DEDX ICRI90:" << (nullptr !=
271 << " Ekin=" << kineticEnergy
272 << " " << p->GetParticleName()
273 << " q2=" << chargeSquare
274 << " inside " << material->GetName() << G4
275 */
276 if(nullptr != fICRU90 && kineticEnergy < fPr
277 if(material != currentMaterial) {
278 currentMaterial = material;
279 baseMaterial = material->GetBaseMaterial
280 ? material->GetBaseMaterial() : materi
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(
289 } else {
290 const G4double e = kineticEnergy*CLH
291 dedx = fICRU90->GetElectronicDEDXforProton
292 }
293 } else {
294 dedx = fICRU90->GetElectronicDEDXforPr
295 *chargeSquare;
296 }
297 dedx *= material->GetDensity();
298 if(cutEnergy < tmax) {
299 dedx += (G4Log(xc) + (1.0 - xc)*beta2)
300 *(eDensity*chargeSquare/beta2);
301 }
302 //G4cout << " iICRU90=" << iICRU90 <<
303 if(dedx > 0.0) { return dedx; }
304 }
305 }
306 // general Bethe-Bloch formula
307 G4double dedx = G4Log(2.0*CLHEP::electron_ma
308 - (1.0 + xc)*beta2;
309
310 if(0.0 < spin) {
311 G4double del = 0.5*cutEnergy/(kineticEnerg
312 dedx += del*del;
313 }
314
315 // density correction
316 G4double x = G4Log(bg2)/twoln10;
317 dedx -= material->GetIonisation()->DensityCo
318
319 // shell correction
320 dedx -= 2.0*corr->ShellCorrection(p,material
321
322 // now compute the total ionization loss
323 dedx *= CLHEP::twopi_mc2_rcl2*chargeSquare*e
324
325 //High order correction different for hadron
326 if(isIon) {
327 dedx += corr->IonBarkasCorrection(p,materi
328 } else {
329 dedx += corr->HighOrderCorrections(p,mater
330 }
331
332 dedx = std::max(dedx, 0.0);
333 /*
334 G4cout << "E(MeV)= " << kineticEnergy/CLHEP:
335 << " " << material->GetName() << G
336 */
337 return dedx;
338 }
339
340 //....oooOO0OOooo........oooOO0OOooo........oo
341
342 void G4BetheBlochModel::CorrectionsAlongStep(c
343 c
344 c
345 G
346 {
347 // no correction for alpha
348 if(isAlpha) { return; }
349
350 // no correction at the last step or at smal
351 const G4double preKinEnergy = dp->GetKinetic
352 if(eloss >= preKinEnergy || eloss < preKinEn
353
354 // corrections for all charged particles wit
355 const G4ParticleDefinition* p = dp->GetDefin
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 - e
361 const G4Material* mat = couple->GetMaterial(
362 const G4double q20 = corr->EffectiveChargeSq
363 const G4double q2 = corr->EffectiveChargeSqu
364 const G4double qfactor = q2/q20;
365
366 /*
367 G4cout << "G4BetheBlochModel::CorrectionsA
368 << preKinEnergy << " Eeff(MeV)=" << e
369 << " eloss=" << eloss << " elossnew=" << e
370 << " qfactor=" << qfactor << " Qpre=" << q
371 << p->GetParticleName() <<G4endl;
372 */
373 eloss *= qfactor;
374 }
375
376 //....oooOO0OOooo........oooOO0OOooo........oo
377
378 void G4BetheBlochModel::SampleSecondaries(std:
379 cons
380 cons
381 G4do
382 G4do
383 {
384 G4double kinEnergy = dp->GetKineticEnergy();
385 const G4double tmax = MaxSecondaryEnergy(dp-
386 const G4double minKinEnergy = std::min(cut,
387 const G4double maxKinEnergy = std::min(maxEn
388 if(minKinEnergy >= maxKinEnergy) { return; }
389
390 //G4cout << "G4BetheBlochModel::SampleSecond
391 // << " Emax= " << maxKinEnergy << G
392
393 const G4double totEnergy = kinEnergy + mass;
394 const G4double etot2 = totEnergy*totEnergy;
395 const G4double beta2 = kinEnergy*(kinEnergy
396
397 G4double deltaKinEnergy, f;
398 G4double f1 = 0.0;
399 G4double fmax = 1.0;
400 if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*
401
402 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
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
410
411 f = 1.0 - beta2*deltaKinEnergy/tmax;
412 if( 0.0 < spin ) {
413 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/e
414 f += f1;
415 }
416
417 // Loop checking, 03-Aug-2015, Vladimir Iv
418 } while( fmax*rndm[1] > f);
419
420 // projectile formfactor - suppresion of hig
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*delta
430 grej *= (1.0 + magMoment2*(x2 - f1/f)/(1
431 }
432 if(grej > 1.1) {
433 G4cout << "### G4BetheBlochModel WARNING
434 << " " << dp->GetDefinition()->G
435 << " Ekin(MeV)= " << kinEnergy
436 << " delEkin(MeV)= " << deltaKinE
437 << G4endl;
438 }
439 if(rndmEngineMod->flat() > grej) { return;
440 }
441
442 G4ThreeVector deltaDirection;
443
444 if(UseAngularGeneratorFlag()) {
445 const G4Material* mat = couple->GetMateria
446 deltaDirection =
447 GetAngularDistribution()->SampleDirectio
448 SelectRandomAtomNumber(mat),
449 mat);
450 } else {
451
452 G4double deltaMomentum =
453 std::sqrt(deltaKinEnergy * (deltaKinEner
454 G4double cost = deltaKinEnergy * (totEnerg
455 (deltaMomentum * dp->GetTotalMomentum())
456 cost = std::min(cost, 1.0);
457 const G4double sint = std::sqrt((1.0 - cos
458 const G4double phi = twopi*rndmEngineMod->
459
460 deltaDirection.set(sint*std::cos(phi),sint
461 deltaDirection.rotateUz(dp->GetMomentumDir
462 }
463 /*
464 G4cout << "### G4BetheBlochModel "
465 << dp->GetDefinition()->GetParticle
466 << " Ekin(MeV)= " << kinEnergy
467 << " delEkin(MeV)= " << deltaKinEne
468 << " tmin(MeV)= " << minKinEnergy
469 << " tmax(MeV)= " << maxKinEnergy
470 << " dir= " << dp->GetMomentumDirec
471 << " dirDelta= " << deltaDirection
472 << G4endl;
473 */
474 // create G4DynamicParticle object for delta
475 auto delta = new G4DynamicParticle(theElectr
476
477 vdp->push_back(delta);
478
479 // Change kinematics of primary particle
480 kinEnergy -= deltaKinEnergy;
481 G4ThreeVector finalP = dp->GetMomentum() - d
482 finalP = finalP.unit();
483
484 fParticleChange->SetProposedKineticEnergy(ki
485 fParticleChange->SetProposedMomentumDirectio
486 }
487
488 //....oooOO0OOooo........oooOO0OOooo........oo
489
490 G4double G4BetheBlochModel::MaxSecondaryEnergy
491
492 {
493 // here particle type is checked for the cas
494 // when this model is shared between particl
495 if(pd != particle) { SetupParameters(pd); }
496 G4double tau = kinEnergy/mass;
497 return 2.0*CLHEP::electron_mass_c2*tau*(tau
498 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
499 }
500
501 //....oooOO0OOooo........oooOO0OOooo........oo
502