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

Last change on this file since 1183 was 1055, checked in by garnier, 15 years ago

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