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

Last change on this file since 1042 was 1007, checked in by garnier, 17 years ago

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