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

Last change on this file since 1109 was 1055, checked in by garnier, 17 years ago

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