Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4BraggIonModel.cc
Go to the documentation of this file.
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 file
30//
31//
32// File name: G4BraggIonModel
33//
34// Author: Vladimir Ivanchenko
35//
36// Creation date: 13.10.2004
37//
38// Modifications:
39// 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko)
40// 29-11-05 Do not use G4Alpha class (V.Ivantchenko)
41// 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
42// 25-04-06 Add stopping data from ASTAR (V.Ivanchenko)
43// 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko)
44// 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
45// CorrectionsAlongStep needed for ions(V.Ivanchenko)
46//
47
48// Class Description:
49//
50// Implementation of energy loss and delta-electron production by
51// slow charged heavy particles
52
53// -------------------------------------------------------------------
54//
55
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59#include "G4BraggIonModel.hh"
61#include "G4SystemOfUnits.hh"
62#include "Randomize.hh"
63#include "G4Electron.hh"
65#include "G4LossTableManager.hh"
66#include "G4EmCorrections.hh"
67#include "G4DeltaAngle.hh"
69#include "G4NistManager.hh"
70#include "G4Log.hh"
71#include "G4Exp.hh"
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75using namespace std;
76
77G4ASTARStopping* G4BraggIonModel::fASTAR = nullptr;
78
80 const G4String& nam)
81 : G4VEmModel(nam),
82 corr(nullptr),
83 particle(nullptr),
84 fParticleChange(nullptr),
85 fICRU90(nullptr),
86 currentMaterial(nullptr),
87 baseMaterial(nullptr),
88 iMolecula(-1),
89 iASTAR(-1),
90 iICRU90(-1),
91 isIon(false)
92{
93 SetHighEnergyLimit(2.0*MeV);
94
95 HeMass = 3.727417*GeV;
96 rateMassHe2p = HeMass/proton_mass_c2;
97 lowestKinEnergy = 1.0*keV/rateMassHe2p;
98 massFactor = 1000.*amu_c2/HeMass;
99 theZieglerFactor = eV*cm2*1.0e-15;
100 theElectron = G4Electron::Electron();
101 corrFactor = 1.0;
102 if(p) { SetParticle(p); }
103 else { SetParticle(theElectron); }
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
109{
110 if(IsMaster()) { delete fASTAR; fASTAR = nullptr; }
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114
116 const G4DataVector&)
117{
118 if(p != particle) { SetParticle(p); }
119
120 corrFactor = chargeSquare;
121
122 // always false before the run
123 SetDeexcitationFlag(false);
124
125 if(IsMaster()) {
126 if(nullptr == fASTAR) { fASTAR = new G4ASTARStopping(); }
127 if(particle->GetPDGMass() < GeV) { fASTAR->Initialise(); }
128 if(G4EmParameters::Instance()->UseICRU90Data()) {
129 if(!fICRU90) {
131 } else if(particle->GetPDGMass() < GeV) { fICRU90->Initialise(); }
132 }
133 }
134
135 if(nullptr == fParticleChange) {
136
139 }
140 G4String pname = particle->GetParticleName();
141 if(particle->GetParticleType() == "nucleus" &&
142 pname != "deuteron" && pname != "triton" &&
143 pname != "alpha+" && pname != "helium" &&
144 pname != "hydrogen") { isIon = true; }
145
147
148 fParticleChange = GetParticleChangeForLoss();
149 }
150}
151
152//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
153
155 const G4MaterialCutsCouple* couple)
156{
158}
159
160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161
163 const G4Material* mat,
164 G4double kineticEnergy)
165{
166 //G4cout<<"G4BraggIonModel::GetChargeSquareRatio e= "<<kineticEnergy<<G4endl;
167 // this method is called only for ions
168 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
169 corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
170 return corrFactor;
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
176 const G4Material* mat,
177 G4double kineticEnergy)
178{
179 //G4cout<<"G4BraggIonModel::GetParticleCharge e= "<<kineticEnergy <<
180 // " q= " << corr->GetParticleCharge(p,mat,kineticEnergy) <<G4endl;
181 // this method is called only for ions
182 return corr->GetParticleCharge(p,mat,kineticEnergy);
183}
184
185//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
188 const G4ParticleDefinition* p,
189 G4double kineticEnergy,
190 G4double cutEnergy,
191 G4double maxKinEnergy)
192{
193 G4double cross = 0.0;
194 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
195 G4double maxEnergy = std::min(tmax,maxKinEnergy);
196 if(cutEnergy < tmax) {
197
198 G4double energy = kineticEnergy + mass;
199 G4double energy2 = energy*energy;
200 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
201 cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
202 - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
203
204 if( 0.0 < spin ) { cross += 0.5*(maxEnergy - cutEnergy)/energy2; }
205
206 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
207 }
208 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
209 // << " tmax= " << tmax << " cross= " << cross << G4endl;
210
211 return cross;
212}
213
214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
215
217 const G4ParticleDefinition* p,
218 G4double kineticEnergy,
220 G4double cutEnergy,
221 G4double maxEnergy)
222{
223 return
224 Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
225}
226
227//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228
230 const G4Material* material,
231 const G4ParticleDefinition* p,
232 G4double kineticEnergy,
233 G4double cutEnergy,
234 G4double maxEnergy)
235{
236 return material->GetElectronDensity()
237 *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
238}
239
240//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241
243 const G4ParticleDefinition* p,
244 G4double kineticEnergy,
245 G4double cutEnergy)
246{
247 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
248 G4double tmin = min(cutEnergy, tmax);
249 G4double tkin = kineticEnergy/massRate;
250 G4double dedx = 0.0;
251
252 if(tkin < lowestKinEnergy) {
253 dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
254 } else {
255 dedx = DEDX(material, tkin);
256 }
257
258 if (cutEnergy < tmax) {
259
260 G4double tau = kineticEnergy/mass;
261 G4double gam = tau + 1.0;
262 G4double bg2 = tau * (tau+2.0);
263 G4double beta2 = bg2/(gam*gam);
264 G4double x = tmin/tmax;
265
266 dedx += (G4Log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
267 * (material->GetElectronDensity())/beta2;
268 }
269
270 // now compute the total ionization loss
271
272 dedx = std::max(dedx, 0.0);
273 dedx *= chargeSquare;
274 /*
275 G4cout << "Bragg: tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = "
276 << dedx*gram/(MeV*cm2*material->GetDensity())
277 << " q2 = " << chargeSquare << G4endl;
278 */
279 return dedx;
280}
281
282//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
283
285 const G4DynamicParticle* dp,
286 G4double& eloss,
287 G4double&,
288 G4double /*length*/)
289{
290 // this method is called only for ions
291 const G4ParticleDefinition* p = dp->GetDefinition();
292 const G4Material* mat = couple->GetMaterial();
293 G4double preKinEnergy = dp->GetKineticEnergy();
294 G4double e = preKinEnergy - eloss*0.5;
295 if(e < 0.0) { e = preKinEnergy*0.5; }
296
297 G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
299 G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
300 eloss *= qfactor;
301
302 //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " << e
303 // << " qfactor= " << qfactor << " " << p->GetParticleName() <<G4endl;
304}
305
306//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307
308void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
309 const G4MaterialCutsCouple* couple,
310 const G4DynamicParticle* dp,
311 G4double xmin,
312 G4double maxEnergy)
313{
315 G4double xmax = std::min(tmax, maxEnergy);
316 if(xmin >= xmax) { return; }
317
318 G4double kineticEnergy = dp->GetKineticEnergy();
319 G4double energy = kineticEnergy + mass;
320 G4double energy2 = energy*energy;
321 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
322 G4double grej = 1.0;
323 G4double deltaKinEnergy, f;
324
325 CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
326 G4double rndm[2];
327
328 // sampling follows ...
329 do {
330 rndmEngineMod->flatArray(2, rndm);
331 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
332
333 f = 1.0 - beta2*deltaKinEnergy/tmax;
334
335 if(f > grej) {
336 G4cout << "G4BraggIonModel::SampleSecondary Warning! "
337 << "Majorant " << grej << " < "
338 << f << " for e= " << deltaKinEnergy
339 << G4endl;
340 }
341
342 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
343 } while( grej*rndm[1] >= f );
344
345 G4ThreeVector deltaDirection;
346
348 const G4Material* mat = couple->GetMaterial();
350
351 deltaDirection =
352 GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
353
354 } else {
355
356 G4double deltaMomentum =
357 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
358 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
359 (deltaMomentum * dp->GetTotalMomentum());
360 if(cost > 1.0) { cost = 1.0; }
361 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
362
363 G4double phi = twopi*rndmEngineMod->flat();
364
365 deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
366 deltaDirection.rotateUz(dp->GetMomentumDirection());
367 }
368
369 // create G4DynamicParticle object for delta ray
370 G4DynamicParticle* delta =
371 new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
372
373 vdp->push_back(delta);
374
375 // Change kinematics of primary particle
376 kineticEnergy -= deltaKinEnergy;
377 G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
378 finalP = finalP.unit();
379
380 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
381 fParticleChange->SetProposedMomentumDirection(finalP);
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385
387 G4double kinEnergy)
388{
389 if(pd != particle) { SetParticle(pd); }
390 G4double tau = kinEnergy/mass;
391 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
392 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
393 return tmax;
394}
395
396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
397
398void G4BraggIonModel::HasMaterial(const G4Material* mat)
399{
400 const G4String& chFormula = mat->GetChemicalFormula();
401 if(chFormula.empty()) { return; }
402
403 // ICRU Report N49, 1993. Ziegler model for He.
404
405 static const size_t numberOfMolecula = 11;
406 static const G4String molName[numberOfMolecula] = {
407 "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate",
408 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate",
409 "Polysterene", "SiO_2", "NaI", "H_2O",
410 "Graphite" };
411
412 // Search for the material in the table
413 for (size_t i=0; i<numberOfMolecula; ++i) {
414 if (chFormula == molName[i]) {
415 iMolecula = i;
416 return;
417 }
418 }
419 return;
420}
421
422//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
423
424G4double G4BraggIonModel::StoppingPower(const G4Material* material,
425 G4double kineticEnergy)
426{
427 G4double ionloss = 0.0 ;
428
429 if (iMolecula >= 0) {
430
431 // The data and the fit from:
432 // ICRU Report N49, 1993. Ziegler's model for alpha
433 // He energy in internal units of parametrisation formula (MeV)
434
435 G4double T = kineticEnergy*rateMassHe2p/MeV ;
436
437 static const G4float a[11][5] = {
438 {9.43672f, 0.54398f, 84.341f, 1.3705f, 57.422f},
439 {67.1503f, 0.41409f, 404.512f, 148.97f, 20.99f},
440 {5.11203f, 0.453f, 36.718f, 50.6f, 28.058f},
441 {61.793f, 0.48445f, 361.537f, 57.889f, 50.674f},
442 {7.83464f, 0.49804f, 160.452f, 3.192f, 0.71922f},
443 {19.729f, 0.52153f, 162.341f, 58.35f, 25.668f},
444 {26.4648f, 0.50112f, 188.913f, 30.079f, 16.509f},
445 {7.8655f, 0.5205f, 63.96f, 51.32f, 67.775f},
446 {8.8965f, 0.5148f, 339.36f, 1.7205f, 0.70423f},
447 {2.959f, 0.53255f, 34.247f, 60.655f, 15.153f},
448 {3.80133f, 0.41590f, 12.9966f, 117.83f, 242.28f} };
449
450 static const G4double atomicWeight[11] = {
451 101.96128f, 44.0098f, 16.0426f, 28.0536f, 42.0804f,
452 104.1512f, 44.665f, 60.0843f, 18.0152f, 18.0152f, 12.0f};
453
454 G4int i = iMolecula;
455
456 G4double slow = (G4double)(a[i][0]);
457
458 G4double x1 = (G4double)(a[i][1]);
459 G4double x2 = (G4double)(a[i][2]);
460 G4double x3 = (G4double)(a[i][3]);
461 G4double x4 = (G4double)(a[i][4]);
462
463 // Free electron gas model
464 if ( T < 0.001 ) {
465 G4double shigh = G4Log( 1.0 + x3*1000.0 + x4*0.001 ) *x2*1000.0;
466 ionloss = slow*shigh / (slow + shigh) ;
467 ionloss *= sqrt(T*1000.0) ;
468
469 // Main parametrisation
470 } else {
471 slow *= G4Exp(G4Log(T*1000.0)*x1) ;
472 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T ;
473 ionloss = slow*shigh / (slow + shigh) ;
474 /*
475 G4cout << "## " << i << ". T= " << T << " slow= " << slow
476 << " a0= " << a[i][0] << " a1= " << a[i][1]
477 << " shigh= " << shigh
478 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
479 << G4endl;
480 */
481 }
482 ionloss = std::max(ionloss, 0.0);
483
484 // He effective charge
485 G4double aa = (G4double)atomicWeight[iMolecula];
486 ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
487
488 // pure material (normally not the case for this function)
489 } else if(1 == (material->GetNumberOfElements())) {
490 G4double z = material->GetZ() ;
491 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
492 }
493
494 return ionloss;
495}
496
497//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
498
499G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
500 G4double kineticEnergy) const
501{
502 G4double ionloss ;
503 G4int i = std::min(std::max(G4lrint(z)-1,0),91); // index of atom
504
505 // The data and the fit from:
506 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
507 // Proton kinetic energy for parametrisation (keV/amu)
508
509 // He energy in internal units of parametrisation formula (MeV)
510 G4double T = kineticEnergy*rateMassHe2p/MeV ;
511
512 static const G4float a[92][5] = {
513 { 0.35485f, 0.6456f, 6.01525f, 20.8933f, 4.3515f
514 },{ 0.58f, 0.59f, 6.3f, 130.0f, 44.07f
515 },{ 1.42f, 0.49f, 12.25f, 32.0f, 9.161f
516 },{ 2.206f, 0.51f, 15.32f, 0.25f, 8.995f //Be Ziegler77
517 // },{ 2.1895f, 0.47183,7.2362f, 134.30f, 197.96f //Be from ICRU
518 },{ 3.691f, 0.4128f, 18.48f, 50.72f, 9.0f
519 },{ 3.83523f, 0.42993f,12.6125f, 227.41f, 188.97f
520 // },{ 1.9259f, 0.5550f, 27.15125f, 26.0665f, 6.2768f //too many digits
521 },{ 1.9259f, 0.5550f, 27.1513f, 26.0665f, 6.2768f
522 },{ 2.81015f, 0.4759f, 50.0253f, 10.556f, 1.0382f
523 },{ 1.533f, 0.531f, 40.44f, 18.41f, 2.718f
524 },{ 2.303f, 0.4861f, 37.01f, 37.96f, 5.092f
525 // Z= 11-20
526 },{ 9.894f, 0.3081f, 23.65f, 0.384f, 92.93f
527 },{ 4.3f, 0.47f, 34.3f, 3.3f, 12.74f
528 },{ 2.5f, 0.625f, 45.7f, 0.1f, 4.359f
529 },{ 2.1f, 0.65f, 49.34f, 1.788f, 4.133f
530 },{ 1.729f, 0.6562f, 53.41f, 2.405f, 3.845f
531 },{ 1.402f, 0.6791f, 58.98f, 3.528f, 3.211f
532 },{ 1.117f, 0.7044f, 69.69f, 3.705f, 2.156f
533 },{ 2.291f, 0.6284f, 73.88f, 4.478f, 2.066f
534 },{ 8.554f, 0.3817f, 83.61f, 11.84f, 1.875f
535 },{ 6.297f, 0.4622f, 65.39f, 10.14f, 5.036f
536 // Z= 21-30
537 },{ 5.307f, 0.4918f, 61.74f, 12.4f, 6.665f
538 },{ 4.71f, 0.5087f, 65.28f, 8.806f, 5.948f
539 },{ 6.151f, 0.4524f, 83.0f, 18.31f, 2.71f
540 },{ 6.57f, 0.4322f, 84.76f, 15.53f, 2.779f
541 },{ 5.738f, 0.4492f, 84.6f, 14.18f, 3.101f
542 },{ 5.013f, 0.4707f, 85.8f, 16.55f, 3.211f
543 },{ 4.32f, 0.4947f, 76.14f, 10.85f, 5.441f
544 },{ 4.652f, 0.4571f, 80.73f, 22.0f, 4.952f
545 },{ 3.114f, 0.5236f, 76.67f, 7.62f, 6.385f
546 },{ 3.114f, 0.5236f, 76.67f, 7.62f, 7.502f
547 // Z= 31-40
548 },{ 3.114f, 0.5236f, 76.67f, 7.62f, 8.514f
549 },{ 5.746f, 0.4662f, 79.24f, 1.185f, 7.993f
550 },{ 2.792f, 0.6346f, 106.1f, 0.2986f, 2.331f
551 },{ 4.667f, 0.5095f, 124.3f, 2.102f, 1.667f
552 },{ 2.44f, 0.6346f, 105.0f, 0.83f, 2.851f
553 },{ 1.413f, 0.7377f, 147.9f, 1.466f, 1.016f
554 },{ 11.72f, 0.3826f, 102.8f, 9.231f, 4.371f
555 },{ 7.126f, 0.4804f, 119.3f, 5.784f, 2.454f
556 },{ 11.61f, 0.3955f, 146.7f, 7.031f, 1.423f
557 },{ 10.99f, 0.41f, 163.9f, 7.1f, 1.052f
558 // Z= 41-50
559 },{ 9.241f, 0.4275f, 163.1f, 7.954f, 1.102f
560 },{ 9.276f, 0.418f, 157.1f, 8.038f, 1.29f
561 },{ 3.999f, 0.6152f, 97.6f, 1.297f, 5.792f
562 },{ 4.306f, 0.5658f, 97.99f, 5.514f, 5.754f
563 },{ 3.615f, 0.6197f, 86.26f, 0.333f, 8.689f
564 },{ 5.8f, 0.49f, 147.2f, 6.903f, 1.289f
565 },{ 5.6f, 0.49f, 130.0f, 10.0f, 2.844f
566 },{ 3.55f, 0.6068f, 124.7f, 1.112f, 3.119f
567 },{ 3.6f, 0.62f, 105.8f, 0.1692f, 6.026f
568 },{ 5.4f, 0.53f, 103.1f, 3.931f, 7.767f
569 // Z= 51-60
570 },{ 3.97f, 0.6459f, 131.8f, 0.2233f, 2.723f
571 },{ 3.65f, 0.64f, 126.8f, 0.6834f, 3.411f
572 },{ 3.118f, 0.6519f, 164.9f, 1.208f, 1.51f
573 },{ 3.949f, 0.6209f, 200.5f, 1.878f, 0.9126f
574 },{ 14.4f, 0.3923f, 152.5f, 8.354f, 2.597f
575 },{ 10.99f, 0.4599f, 138.4f, 4.811f, 3.726f
576 },{ 16.6f, 0.3773f, 224.1f, 6.28f, 0.9121f
577 },{ 10.54f, 0.4533f, 159.3f, 4.832f, 2.529f
578 },{ 10.33f, 0.4502f, 162.0f, 5.132f, 2.444f
579 },{ 10.15f, 0.4471f, 165.6f, 5.378f, 2.328f
580 // Z= 61-70
581 },{ 9.976f, 0.4439f, 168.0f, 5.721f, 2.258f
582 },{ 9.804f, 0.4408f, 176.2f, 5.675f, 1.997f
583 },{ 14.22f, 0.363f, 228.4f, 7.024f, 1.016f
584 },{ 9.952f, 0.4318f, 233.5f, 5.065f, 0.9244f
585 },{ 9.272f, 0.4345f, 210.0f, 4.911f, 1.258f
586 },{ 10.13f, 0.4146f, 225.7f, 5.525f, 1.055f
587 },{ 8.949f, 0.4304f, 213.3f, 5.071f, 1.221f
588 },{ 11.94f, 0.3783f, 247.2f, 6.655f, 0.849f
589 },{ 8.472f, 0.4405f, 195.5f, 4.051f, 1.604f
590 },{ 8.301f, 0.4399f, 203.7f, 3.667f, 1.459f
591 // Z= 71-80
592 },{ 6.567f, 0.4858f, 193.0f, 2.65f, 1.66f
593 },{ 5.951f, 0.5016f, 196.1f, 2.662f, 1.589f
594 },{ 7.495f, 0.4523f, 251.4f, 3.433f, 0.8619f
595 },{ 6.335f, 0.4825f, 255.1f, 2.834f, 0.8228f
596 },{ 4.314f, 0.5558f, 214.8f, 2.354f, 1.263f
597 },{ 4.02f, 0.5681f, 219.9f, 2.402f, 1.191f
598 },{ 3.836f, 0.5765f, 210.2f, 2.742f, 1.305f
599 },{ 4.68f, 0.5247f, 244.7f, 2.749f, 0.8962f
600 },{ 2.892f, 0.6204f, 208.6f, 2.415f, 1.416f //Au Z77
601 // },{ 3.223f, 0.5883f, 232.7f, 2.954f, 1.05 //Au ICRU
602 },{ 2.892f, 0.6204f, 208.6f, 2.415f, 1.416f
603 // Z= 81-90
604 },{ 4.728f, 0.5522f, 217.0f, 3.091f, 1.386f
605 },{ 6.18f, 0.52f, 170.0f, 4.0f, 3.224f
606 },{ 9.0f, 0.47f, 198.0f, 3.8f, 2.032f
607 },{ 2.324f, 0.6997f, 216.0f, 1.599f, 1.399f
608 },{ 1.961f, 0.7286f, 223.0f, 1.621f, 1.296f
609 },{ 1.75f, 0.7427f, 350.1f, 0.9789f, 0.5507f
610 },{ 10.31f, 0.4613f, 261.2f, 4.738f, 0.9899f
611 },{ 7.962f, 0.519f, 235.7f, 4.347f, 1.313f
612 },{ 6.227f, 0.5645f, 231.9f, 3.961f, 1.379f
613 },{ 5.246f, 0.5947f, 228.6f, 4.027f, 1.432f
614 // Z= 91-92
615 },{ 5.408f, 0.5811f, 235.7f, 3.961f, 1.358f
616 },{ 5.218f, 0.5828f, 245.0f, 3.838f, 1.25f}
617 };
618
619 G4double slow = (G4double)(a[i][0]);
620
621 G4double x1 = (G4double)(a[i][1]);
622 G4double x2 = (G4double)(a[i][2]);
623 G4double x3 = (G4double)(a[i][3]);
624 G4double x4 = (G4double)(a[i][4]);
625
626 // Free electron gas model
627 if ( T < 0.001 ) {
628 G4double shigh = G4Log( 1.0 + x3*1000.0 + x4*0.001 )* x2*1000.0;
629 ionloss = slow*shigh*sqrt(T*1000.0) / (slow + shigh) ;
630
631 // Main parametrisation
632 } else {
633 slow *= G4Exp(G4Log(T*1000.0)*x1);
634 G4double shigh = G4Log( 1.0 + x3/T + x4*T ) * x2/T;
635 ionloss = slow*shigh / (slow + shigh) ;
636 /*
637 G4cout << "## " << i << ". T= " << T << " slow= " << slow
638 << " a0= " << a[i][0] << " a1= " << a[i][1]
639 << " shigh= " << shigh
640 << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
641 << G4endl;
642 */
643 }
644 ionloss = std::max(ionloss, 0.0);
645
646 // He effective charge
647 ionloss /= HeEffChargeSquare(z, T);
648
649 return ionloss;
650}
651
652//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
653
654G4double G4BraggIonModel::DEDX(const G4Material* material,
655 G4double kineticEnergy)
656{
657 G4double eloss = 0.0;
658 // check DB
659 if(material != currentMaterial) {
660 currentMaterial = material;
661 baseMaterial = material->GetBaseMaterial()
662 ? material->GetBaseMaterial() : material;
663 iASTAR = -1;
664 iMolecula = -1;
665 iICRU90 = fICRU90 ? fICRU90->GetIndex(baseMaterial) : -1;
666
667 if(iICRU90 < 0) {
668 iASTAR = fASTAR->GetIndex(baseMaterial);
669 if(iASTAR < 0) { HasMaterial(baseMaterial); }
670 }
671 //G4cout << "%%% " <<material->GetName() << " iMolecula= "
672 // << iMolecula << " iPSTAR= " << iPSTAR
673 // << " iICRU90= " << iICRU90<< G4endl;
674 }
675 if(iICRU90 >= 0) {
676 return fICRU90->GetElectronicDEDXforAlpha(iICRU90, kineticEnergy)
677 *material->GetDensity()/chargeSquare;
678 }
679 if( iASTAR >= 0 ) {
680 G4double T = kineticEnergy*rateMassHe2p;
681 G4int zeff = G4lrint(material->GetTotNbOfElectPerVolume()/
682 material->GetTotNbOfAtomsPerVolume());
683 return fASTAR->GetElectronicDEDX(iASTAR, T)*material->GetDensity()/
684 HeEffChargeSquare(zeff, T/MeV);
685 }
686
687 const G4int numberOfElements = material->GetNumberOfElements();
688 const G4double* theAtomicNumDensityVector =
689 material->GetAtomicNumDensityVector();
690
691 if(iMolecula >= 0) {
692
693 eloss = StoppingPower(baseMaterial, kineticEnergy)*
694 material->GetDensity()/amu;
695
696 // pure material
697 } else if(1 == numberOfElements) {
698
699 G4double z = material->GetZ();
700 eloss = ElectronicStoppingPower(z, kineticEnergy)
701 * (material->GetTotNbOfAtomsPerVolume());
702
703 // Brugg's rule calculation
704 } else {
705 const G4ElementVector* theElementVector =
706 material->GetElementVector() ;
707
708 // loop for the elements in the material
709 for (G4int i=0; i<numberOfElements; i++)
710 {
711 const G4Element* element = (*theElementVector)[i] ;
712 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
713 * theAtomicNumDensityVector[i];
714 }
715 }
716 return eloss*theZieglerFactor;
717}
718
719//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
720
721G4double G4BraggIonModel::HeEffChargeSquare(G4double z,
722 G4double kinEnergyHeInMeV) const
723{
724 // The aproximation of He effective charge from:
725 // J.F.Ziegler, J.P. Biersack, U. Littmark
726 // The Stopping and Range of Ions in Matter,
727 // Vol.1, Pergamon Press, 1985
728
729 static const G4double c[6] = {0.2865, 0.1266, -0.001429,
730 0.02402,-0.01135, 0.001475};
731
732 G4double e = std::max(0.0, G4Log(kinEnergyHeInMeV*massFactor));
733 G4double x = c[0] ;
734 G4double y = 1.0 ;
735 for (G4int i=1; i<6; ++i) {
736 y *= e;
737 x += y * c[i];
738 }
739
740 G4double w = 7.6 - e ;
741 w = 1.0 + (0.007 + 0.00005*z) * G4Exp( -w*w ) ;
742 w = 4.0 * (1.0 - G4Exp(-x)) * w * w ;
743
744 return w;
745}
746
747//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
748
std::vector< G4Element * > G4ElementVector
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4double G4Log(G4double x)
Definition: G4Log.hh:226
float G4float
Definition: G4Types.hh:84
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
Hep3Vector unit() const
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:33
virtual double flat()=0
virtual void flatArray(const int size, double *vect)=0
G4int GetIndex(const G4Material *) const
G4double GetElectronicDEDX(G4int idx, G4double energy) const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *couple) override
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy) override
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length) override
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy) override
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy) override
G4BraggIonModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="BraggIon")
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy, G4double maxEnergy) override
virtual ~G4BraggIonModel()
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4double GetZ() const
Definition: G4Element.hh:130
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetParticleCharge(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
static G4EmParameters * Instance()
G4int GetIndex(const G4Material *) const
G4double GetElectronicDEDXforAlpha(const G4Material *, G4double scaledKinEnergy) const
G4double GetMeanExcitationEnergy() const
static G4LossTableManager * Instance()
G4EmCorrections * EmCorrections()
const G4Material * GetMaterial() const
G4double GetDensity() const
Definition: G4Material.hh:178
const G4String & GetChemicalFormula() const
Definition: G4Material.hh:176
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:231
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
G4double GetZ() const
Definition: G4Material.cc:701
G4double GetTotNbOfElectPerVolume() const
Definition: G4Material.hh:210
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:224
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:214
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4ICRU90StoppingData * GetICRU90StoppingData()
static G4NistManager * Instance()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
const G4String & GetParticleType() const
const G4String & GetParticleName() const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
virtual void SetParticleAndCharge(const G4ParticleDefinition *, G4double q2)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757
G4VEmFluctuationModel * GetModelOfFluctuations()
Definition: G4VEmModel.hh:604
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:611
G4bool IsMaster() const
Definition: G4VEmModel.hh:736
G4int SelectRandomAtomNumber(const G4Material *)
Definition: G4VEmModel.cc:315
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:813
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618
G4bool UseAngularGeneratorFlag() const
Definition: G4VEmModel.hh:708
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:510
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
int G4lrint(double ad)
Definition: templates.hh:134