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

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

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