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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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