source: trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc@ 953

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

import all except CVS
File size: 26.1 KB
Line 
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: G4BraggModel.cc,v 1.16 2007/07/28 13:30:53 vnivanch Exp $
27// GEANT4 tag $Name: $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4BraggModel
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 03.01.2002
39//
40// Modifications:
41//
42// 04-12-02 Fix problem of G4DynamicParticle constructor (V.Ivanchenko)
43// 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
44// 27-01-03 Make models region aware (V.Ivanchenko)
45// 13-02-03 Add name (V.Ivanchenko)
46// 04-06-03 Fix compilation warnings (V.Ivanchenko)
47// 12-09-04 Add lowestKinEnergy and change order of if in DEDX method (VI)
48// 11-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
49// 16-06-05 Fix problem of chemical formula (V.Ivantchenko)
50// 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
51// 25-04-06 Add stopping data from PSTAR (V.Ivanchenko)
52
53// Class Description:
54//
55// Implementation of energy loss and delta-electron production by
56// slow charged heavy particles
57
58// -------------------------------------------------------------------
59//
60
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64
65#include "G4BraggModel.hh"
66#include "Randomize.hh"
67#include "G4Electron.hh"
68#include "G4ParticleChangeForLoss.hh"
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72using namespace std;
73
74G4BraggModel::G4BraggModel(const G4ParticleDefinition* p, const G4String& nam)
75 : G4VEmModel(nam),
76 particle(0),
77 protonMassAMU(1.007276),
78 iMolecula(0),
79 isIon(false)
80{
81 if(p) SetParticle(p);
82 lowestKinEnergy = 1.0*keV;
83 theZieglerFactor = eV*cm2*1.0e-15;
84 theElectron = G4Electron::Electron();
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89G4BraggModel::~G4BraggModel()
90{}
91
92//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93
94G4double G4BraggModel::MinEnergyCut(const G4ParticleDefinition*,
95 const G4MaterialCutsCouple* couple)
96{
97 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
98}
99
100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102void G4BraggModel::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//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119G4double G4BraggModel::ComputeCrossSectionPerElectron(
120 const G4ParticleDefinition* p,
121 G4double kineticEnergy,
122 G4double cutEnergy,
123 G4double maxKinEnergy)
124{
125
126 G4double cross = 0.0;
127 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
128 G4double maxEnergy = min(tmax,maxKinEnergy);
129 if(cutEnergy < tmax) {
130
131 G4double energy = kineticEnergy + mass;
132 G4double energy2 = energy*energy;
133 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
134 cross = 1.0/cutEnergy - 1.0/maxEnergy - beta2*log(maxEnergy/cutEnergy)/tmax;
135
136 cross *= twopi_mc2_rcl2*chargeSquare/beta2;
137 }
138 // G4cout << "BR: e= " << kineticEnergy << " tmin= " << cutEnergy
139 // << " tmax= " << tmax << " cross= " << cross << G4endl;
140
141 return cross;
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145
146G4double G4BraggModel::ComputeCrossSectionPerAtom(
147 const G4ParticleDefinition* p,
148 G4double kineticEnergy,
149 G4double Z, G4double,
150 G4double cutEnergy,
151 G4double maxEnergy)
152{
153 G4double cross = Z*ComputeCrossSectionPerElectron
154 (p,kineticEnergy,cutEnergy,maxEnergy);
155 return cross;
156}
157
158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159
160G4double G4BraggModel::CrossSectionPerVolume(
161 const G4Material* material,
162 const G4ParticleDefinition* p,
163 G4double kineticEnergy,
164 G4double cutEnergy,
165 G4double maxEnergy)
166{
167 G4double eDensity = material->GetElectronDensity();
168 G4double cross = eDensity*ComputeCrossSectionPerElectron
169 (p,kineticEnergy,cutEnergy,maxEnergy);
170 return cross;
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
174
175G4double G4BraggModel::ComputeDEDXPerVolume(const G4Material* material,
176 const G4ParticleDefinition* p,
177 G4double kineticEnergy,
178 G4double cutEnergy)
179{
180 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
181 G4double tkin = kineticEnergy/massRate;
182 G4double dedx = 0.0;
183 if(tkin > lowestKinEnergy) dedx = DEDX(material, tkin);
184 else dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
185
186 if (cutEnergy < tmax) {
187
188 G4double tau = kineticEnergy/mass;
189 G4double gam = tau + 1.0;
190 G4double bg2 = tau * (tau+2.0);
191 G4double beta2 = bg2/(gam*gam);
192 G4double x = cutEnergy/tmax;
193
194 dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
195 * (material->GetElectronDensity())/beta2;
196 }
197
198 // now compute the total ionization loss
199
200 if (dedx < 0.0) dedx = 0.0 ;
201
202 dedx *= chargeSquare;
203
204 return dedx;
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208
209void G4BraggModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
210 const G4MaterialCutsCouple*,
211 const G4DynamicParticle* dp,
212 G4double xmin,
213 G4double maxEnergy)
214{
215 G4double tmax = MaxSecondaryKinEnergy(dp);
216 G4double xmax = std::min(tmax, maxEnergy);
217 if(xmin >= xmax) return;
218
219 G4double kineticEnergy = dp->GetKineticEnergy();
220 G4double energy = kineticEnergy + mass;
221 G4double energy2 = energy*energy;
222 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
223 G4double grej = 1.0;
224 G4double deltaKinEnergy, f;
225
226 G4ThreeVector direction = dp->GetMomentumDirection();
227
228 // sampling follows ...
229 do {
230 G4double q = G4UniformRand();
231 deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
232
233 f = 1.0 - beta2*deltaKinEnergy/tmax;
234
235 if(f > grej) {
236 G4cout << "G4BraggModel::SampleSecondary Warning! "
237 << "Majorant " << grej << " < "
238 << f << " for e= " << deltaKinEnergy
239 << G4endl;
240 }
241
242 } while( grej*G4UniformRand() >= f );
243
244 G4double deltaMomentum =
245 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
246 G4double totMomentum = energy*sqrt(beta2);
247 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
248 (deltaMomentum * totMomentum);
249 if(cost > 1.0) cost = 1.0;
250 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
251
252 G4double phi = twopi * G4UniformRand() ;
253
254 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
255 deltaDirection.rotateUz(direction);
256
257 // Change kinematics of primary particle
258 kineticEnergy -= deltaKinEnergy;
259 G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
260 finalP = finalP.unit();
261
262 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
263 fParticleChange->SetProposedMomentumDirection(finalP);
264
265 // create G4DynamicParticle object for delta ray
266 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
267 deltaKinEnergy);
268
269 vdp->push_back(delta);
270}
271
272//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
273
274G4bool G4BraggModel::HasMaterial(const G4Material* material)
275{
276 const size_t numberOfMolecula = 11 ;
277 SetMoleculaNumber(numberOfMolecula) ;
278 G4String chFormula = material->GetChemicalFormula() ;
279
280 // ICRU Report N49, 1993. Power's model for He.
281 static G4String molName[numberOfMolecula] = {
282 "Al_2O_3", "CO_2", "CH_4",
283 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
284 "C_3H_8", "SiO_2", "H_2O",
285 "H_2O-Gas", "Graphite" } ;
286
287 // Special treatment for water in gas state
288 const G4State theState = material->GetState() ;
289
290 if( theState == kStateGas && "H_2O" == chFormula) {
291 chFormula = G4String("H_2O-Gas");
292 }
293
294 // Search for the material in the table
295 for (size_t i=0; i<numberOfMolecula; i++) {
296 if (chFormula == molName[i]) {
297 SetMoleculaNumber(i) ;
298 return true ;
299 }
300 }
301 return false ;
302}
303
304//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
305
306G4double G4BraggModel::StoppingPower(const G4Material* material,
307 G4double kineticEnergy)
308{
309 G4double ionloss = 0.0 ;
310
311 if (iMolecula < 11) {
312
313 // The data and the fit from:
314 // ICRU Report N49, 1993. Ziegler's model for protons.
315 // Proton kinetic energy for parametrisation (keV/amu)
316
317 G4double T = kineticEnergy/(keV*protonMassAMU) ;
318
319 static G4double a[11][5] = {
320 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
321 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
322 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
323 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
324 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
325 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
326 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
327 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
328 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
329 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
330 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
331
332 static G4double atomicWeight[11] = {
333 101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
334 104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
335
336 if ( T < 10.0 ) {
337 ionloss = a[iMolecula][0] * sqrt(T) ;
338
339 } else if ( T < 10000.0 ) {
340 G4double slow = a[iMolecula][1] * pow(T, 0.45) ;
341 G4double shigh = log( 1.0 + a[iMolecula][3]/T
342 + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
343 ionloss = slow*shigh / (slow + shigh) ;
344 }
345
346 if ( ionloss < 0.0) ionloss = 0.0 ;
347 if ( 10 == iMolecula ) {
348 if (T < 100.0) {
349 ionloss *= (1.0+0.023+0.0066*log10(T));
350 }
351 else if (T < 700.0) {
352 ionloss *=(1.0+0.089-0.0248*log10(T-99.));
353 }
354 else if (T < 10000.0) {
355 ionloss *=(1.0+0.089-0.0248*log10(700.-99.));
356 }
357 }
358 ionloss /= atomicWeight[iMolecula];
359
360 // pure material (normally not the case for this function)
361 } else if(1 == (material->GetNumberOfElements())) {
362 G4double z = material->GetZ() ;
363 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
364 }
365
366 return ionloss;
367}
368
369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370
371G4double G4BraggModel::ElectronicStoppingPower(G4double z,
372 G4double kineticEnergy) const
373{
374 G4double ionloss ;
375 G4int i = G4int(z)-1 ; // index of atom
376 if(i < 0) i = 0 ;
377 if(i > 91) i = 91 ;
378
379 // The data and the fit from:
380 // ICRU Report 49, 1993. Ziegler's type of parametrisations.
381 // Proton kinetic energy for parametrisation (keV/amu)
382
383 G4double T = kineticEnergy/(keV*protonMassAMU) ;
384
385 static G4double a[92][5] = {
386 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
387 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
388 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
389 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
390 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
391 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
392 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
393 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
394 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
395 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
396 // Z= 11-20
397 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
398 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
399 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
400 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
401 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
402 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
403 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
404 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
405 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
406 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
407 // Z= 21-30
408 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
409 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
410 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
411 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
412 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
413 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
414 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
415 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
416 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
417 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
418 // Z= 31-40
419 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
420 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
421 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
422 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
423 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
424 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
425 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
426 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
427 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
428 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
429 // Z= 41-50
430 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
431 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
432 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
433 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
434 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
435 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
436 // {5.623, 6.354, 7160.0, 337.6, 0.013940}, // Ag Ziegler77
437 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2}, // Ag ICRU49
438 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
439 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
440 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
441 // Z= 51-60
442 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
443 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
444 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
445 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
446 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
447 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
448 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
449 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
450 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
451 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
452 // Z= 61-70
453 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
454 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
455 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
456 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
457 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
458 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
459 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
460 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
461 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
462 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
463 // Z= 71-80
464 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
465 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
466 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
467 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
468 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
469 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
470 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
471 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
472 // {4.856, 5.460, 18320.0, 438.5, 0.002542}, //Ziegler77
473 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2}, //ICRU49
474 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
475 // Z= 81-90
476 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
477 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
478 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
479 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
480 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
481 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
482 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
483 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
484 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
485 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
486 // Z= 91-92
487 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
488 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
489 };
490
491 G4double fac = 1.0 ;
492
493 // Carbon specific case for E < 40 keV
494 if ( T < 40.0 && 5 == i) {
495 fac = sqrt(T/40.0) ;
496 T = 40.0 ;
497
498 // Free electron gas model
499 } else if ( T < 10.0 ) {
500 fac = sqrt(T*0.1) ;
501 T =10.0 ;
502 }
503
504 // Main parametrisation
505 G4double slow = a[i][1] * pow(T, 0.45) ;
506 G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
507 ionloss = slow*shigh*fac / (slow + shigh) ;
508
509 if ( ionloss < 0.0) ionloss = 0.0 ;
510
511 return ionloss;
512}
513
514//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
515
516G4double G4BraggModel::DEDX(const G4Material* material,
517 G4double kineticEnergy)
518{
519 G4double eloss = 0.0;
520 const G4int numberOfElements = material->GetNumberOfElements();
521 const G4double* theAtomicNumDensityVector =
522 material->GetAtomicNumDensityVector();
523
524 // compaund material with parametrisation
525 G4int iNist = pstar.GetIndex(material);
526
527 if( iNist >= 0 ) {
528 return pstar.GetElectronicDEDX(iNist, kineticEnergy)*material->GetDensity();
529
530 } else if( HasMaterial(material) ) {
531
532 eloss = StoppingPower(material, kineticEnergy)*
533 material->GetDensity()/amu;
534
535 // pure material
536 } else if(1 == numberOfElements) {
537
538 G4double z = material->GetZ();
539 eloss = ElectronicStoppingPower(z, kineticEnergy)
540 * (material->GetTotNbOfAtomsPerVolume());
541
542
543 // Experimental data exist only for kinetic energy 125 keV
544 } else if( MolecIsInZiegler1988(material) ) {
545
546 // Cycle over elements - calculation based on Bragg's rule
547 G4double eloss125 = 0.0 ;
548 const G4ElementVector* theElementVector =
549 material->GetElementVector();
550
551 // loop for the elements in the material
552 for (G4int i=0; i<numberOfElements; i++) {
553 const G4Element* element = (*theElementVector)[i] ;
554 G4double z = element->GetZ() ;
555 eloss += ElectronicStoppingPower(z,kineticEnergy)
556 * theAtomicNumDensityVector[i] ;
557 eloss125 += ElectronicStoppingPower(z,125.0*keV)
558 * theAtomicNumDensityVector[i] ;
559 }
560
561 // Chemical factor is taken into account
562 eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
563
564 // Brugg's rule calculation
565 } else {
566 const G4ElementVector* theElementVector =
567 material->GetElementVector() ;
568
569 // loop for the elements in the material
570 for (G4int i=0; i<numberOfElements; i++)
571 {
572 const G4Element* element = (*theElementVector)[i] ;
573 eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
574 * theAtomicNumDensityVector[i];
575 }
576 }
577 return eloss*theZieglerFactor;
578}
579
580//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
581
582G4bool G4BraggModel::MolecIsInZiegler1988(const G4Material* material)
583{
584 // The list of molecules from
585 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
586 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
587
588 G4String myFormula = G4String(" ") ;
589 const G4String chFormula = material->GetChemicalFormula() ;
590 if (myFormula == chFormula ) return false ;
591
592 // There are no evidence for difference of stopping power depended on
593 // phase of the compound except for water. The stopping power of the
594 // water in gas phase can be predicted using Bragg's rule.
595 //
596 // No chemical factor for water-gas
597
598 myFormula = G4String("H_2O") ;
599 const G4State theState = material->GetState() ;
600 if( theState == kStateGas && myFormula == chFormula) return false ;
601
602 const size_t numberOfMolecula = 53 ;
603
604 // The coffecient from Table.4 of Ziegler & Manoyan
605 const G4double HeEff = 2.8735 ;
606
607 static G4String nameOfMol[numberOfMolecula] = {
608 "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
609 "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
610 "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
611 "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
612 "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
613 "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
614 "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
615 "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
616 "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
617 "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
618 "C_3H_6S", "C_4H_4S", "C_7H_8"
619 } ;
620
621 static G4double expStopping[numberOfMolecula] = {
622 66.1, 190.4, 258.7, 42.2, 141.5,
623 210.9, 279.6, 198.8, 31.0, 267.5,
624 122.8, 311.4, 260.3, 328.9, 391.3,
625 206.6, 374.0, 422.0, 432.0, 398.0,
626 554.0, 353.0, 326.0, 74.6, 220.5,
627 197.4, 362.0, 170.0, 330.5, 211.3,
628 262.3, 349.6, 51.3, 187.0, 236.9,
629 121.9, 35.8, 247.0, 292.6, 268.0,
630 262.3, 49.0, 398.9, 444.0, 22.91,
631 68.0, 155.0, 84.0, 74.2, 254.7,
632 306.8, 324.4, 420.0
633 } ;
634
635 static G4double expCharge[numberOfMolecula] = {
636 HeEff, HeEff, HeEff, 1.0, HeEff,
637 HeEff, HeEff, HeEff, 1.0, 1.0,
638 1.0, HeEff, HeEff, HeEff, HeEff,
639 HeEff, HeEff, HeEff, HeEff, HeEff,
640 HeEff, HeEff, HeEff, 1.0, HeEff,
641 HeEff, HeEff, HeEff, HeEff, HeEff,
642 HeEff, HeEff, 1.0, HeEff, HeEff,
643 HeEff, 1.0, HeEff, HeEff, HeEff,
644 HeEff, 1.0, HeEff, HeEff, 1.0,
645 1.0, 1.0, 1.0, 1.0, HeEff,
646 HeEff, HeEff, HeEff
647 } ;
648
649 static G4double numberOfAtomsPerMolecula[numberOfMolecula] = {
650 3.0, 7.0, 10.0, 4.0, 6.0,
651 9.0, 12.0, 7.0, 4.0, 24.0,
652 12.0, 14.0, 10.0, 13.0, 5.0,
653 5.0, 14.0, 18.0, 17.0, 17.0,
654 24.0, 15.0, 13.0, 9.0, 8.0,
655 6.0, 14.0, 8.0, 8.0, 9.0,
656 10.0, 15.0, 6.0, 7.0, 7.0,
657 3.0, 5.0, 5.0, 5.0, 5.0,
658 9.0, 3.0, 16.0, 14.0, 3.0,
659 9.0, 16.0, 11.0, 9.0, 10.0,
660 10.0, 9.0, 15.0
661 } ;
662
663 // Search for the compaund in the table
664 for (size_t i=0; i<numberOfMolecula; i++)
665 {
666 if(chFormula == nameOfMol[i]) {
667 G4double exp125 = expStopping[i] *
668 (material->GetTotNbOfAtomsPerVolume()) /
669 (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
670 SetExpStopPower125(exp125);
671 return true;
672 }
673 }
674
675 return false;
676}
677
678//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
679
680G4double G4BraggModel::ChemicalFactor(G4double kineticEnergy,
681 G4double eloss125) const
682{
683 // Approximation of Chemical Factor according to
684 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
685 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
686
687 G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
688 G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
689 G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
690 G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ;
691 G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
692 G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
693
694 G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
695 (1.0 + exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
696 (1.0 + exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
697
698 return factor ;
699}
700
701//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
702
Note: See TracBrowser for help on using the repository browser.