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

Last change on this file since 1317 was 1315, checked in by garnier, 15 years ago

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