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