Geant4 11.1.1
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// 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