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

Last change on this file since 825 was 819, checked in by garnier, 16 years ago

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