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

Last change on this file was 1340, checked in by garnier, 15 years ago

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