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

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

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