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

Last change on this file since 969 was 961, checked in by garnier, 17 years ago

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