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