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

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
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 $
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.