source: trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc@ 899

Last change on this file since 899 was 819, checked in by garnier, 17 years ago

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