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