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

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

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