Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BraggIonModel.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: G4BraggIonModel
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 13.10.2004
38 //
39 // Modifications:
40 // 11-05-05 Major optimisation of internal interfaces (V.Ivantchenko)
41 // 29-11-05 Do not use G4Alpha class (V.Ivantchenko)
42 // 15-02-06 ComputeCrossSectionPerElectron, ComputeCrossSectionPerAtom (mma)
43 // 25-04-06 Add stopping data from ASTAR (V.Ivanchenko)
44 // 23-10-06 Reduce lowestKinEnergy to 0.25 keV (V.Ivanchenko)
45 // 12-08-08 Added methods GetParticleCharge, GetChargeSquareRatio,
46 // CorrectionsAlongStep needed for ions(V.Ivanchenko)
47 //
48 
49 // Class Description:
50 //
51 // Implementation of energy loss and delta-electron production by
52 // slow charged heavy particles
53 
54 // -------------------------------------------------------------------
55 //
56 
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 #include "G4BraggIonModel.hh"
62 #include "G4PhysicalConstants.hh"
63 #include "G4SystemOfUnits.hh"
64 #include "Randomize.hh"
65 #include "G4Electron.hh"
66 #include "G4ParticleChangeForLoss.hh"
67 #include "G4LossTableManager.hh"
68 #include "G4EmCorrections.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 using namespace std;
73 
74 G4BraggIonModel::G4BraggIonModel(const G4ParticleDefinition* p,
75  const G4String& nam)
76  : G4VEmModel(nam),
77  corr(0),
78  particle(0),
79  fParticleChange(0),
80  currentMaterial(0),
81  iMolecula(-1),
82  iASTAR(-1),
83  isIon(false),
84  isInitialised(false)
85 {
86  SetHighEnergyLimit(2.0*MeV);
87 
88  HeMass = 3.727417*GeV;
89  rateMassHe2p = HeMass/proton_mass_c2;
90  lowestKinEnergy = 1.0*keV/rateMassHe2p;
91  massFactor = 1000.*amu_c2/HeMass;
92  theZieglerFactor = eV*cm2*1.0e-15;
93  theElectron = G4Electron::Electron();
94  corrFactor = 1.0;
95  if(p) { SetParticle(p); }
96  else { SetParticle(theElectron); }
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 G4BraggIonModel::~G4BraggIonModel()
102 {}
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
106 void G4BraggIonModel::Initialise(const G4ParticleDefinition* p,
107  const G4DataVector&)
108 {
109  if(p != particle) { SetParticle(p); }
110 
111  corrFactor = chargeSquare;
112 
113  // always false before the run
114  SetDeexcitationFlag(false);
115 
116  if(!isInitialised) {
117  isInitialised = true;
118 
119  G4String pname = particle->GetParticleName();
120  if(particle->GetParticleType() == "nucleus" &&
121  pname != "deuteron" && pname != "triton" &&
122  pname != "alpha+" && pname != "helium" &&
123  pname != "hydrogen") { isIon = true; }
124 
125  corr = G4LossTableManager::Instance()->EmCorrections();
126 
127  fParticleChange = GetParticleChangeForLoss();
128  }
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 G4double G4BraggIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
134  const G4Material* mat,
135  G4double kineticEnergy)
136 {
137  //G4cout << "G4BraggIonModel::GetChargeSquareRatio e= " << kineticEnergy << G4endl;
138  // this method is called only for ions
139  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
140  corrFactor = q2*corr->EffectiveChargeCorrection(p,mat,kineticEnergy);
141  return corrFactor;
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145 
146 G4double G4BraggIonModel::GetParticleCharge(const G4ParticleDefinition* p,
147  const G4Material* mat,
148  G4double kineticEnergy)
149 {
150  //G4cout << "G4BraggIonModel::GetParticleCharge e= " << kineticEnergy << G4endl;
151  // this method is called only for ions
152  return corr->GetParticleCharge(p,mat,kineticEnergy);
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
156 
157 G4double G4BraggIonModel::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 < tmax) {
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 
183 G4double G4BraggIonModel::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 
197 G4double G4BraggIonModel::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 
212 G4double G4BraggIonModel::ComputeDEDXPerVolume(const G4Material* material,
213  const G4ParticleDefinition* p,
214  G4double kineticEnergy,
215  G4double cutEnergy)
216 {
217  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
218  G4double tmin = min(cutEnergy, tmax);
219  G4double tkin = kineticEnergy/massRate;
220  G4double dedx = 0.0;
221 
222  if(tkin < lowestKinEnergy) {
223  dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy);
224  } else {
225  dedx = DEDX(material, tkin);
226  }
227 
228  if (cutEnergy < tmax) {
229 
230  G4double tau = kineticEnergy/mass;
231  G4double gam = tau + 1.0;
232  G4double bg2 = tau * (tau+2.0);
233  G4double beta2 = bg2/(gam*gam);
234  G4double x = tmin/tmax;
235 
236  dedx += (log(x) + (1.0 - x)*beta2) * twopi_mc2_rcl2
237  * (material->GetElectronDensity())/beta2;
238  }
239 
240  // now compute the total ionization loss
241 
242  if (dedx < 0.0) dedx = 0.0 ;
243 
244  dedx *= chargeSquare;
245 
246  //G4cout << " tkin(MeV) = " << tkin/MeV << " dedx(MeVxcm^2/g) = "
247  // << dedx*gram/(MeV*cm2*material->GetDensity())
248  // << " q2 = " << chargeSquare << G4endl;
249 
250  return dedx;
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
255 void G4BraggIonModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
256  const G4DynamicParticle* dp,
257  G4double& eloss,
258  G4double&,
259  G4double /*length*/)
260 {
261  // this method is called only for ions
262  const G4ParticleDefinition* p = dp->GetDefinition();
263  const G4Material* mat = couple->GetMaterial();
264  G4double preKinEnergy = dp->GetKineticEnergy();
265  G4double e = preKinEnergy - eloss*0.5;
266  if(e < 0.0) { e = preKinEnergy*0.5; }
267 
268  G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
269  GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
270  G4double qfactor = q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor;
271  eloss *= qfactor;
272 
273  //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " << e
274  // << " qfactor= " << qfactor << " " << p->GetParticleName() <<G4endl;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 void G4BraggIonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
280  const G4MaterialCutsCouple*,
281  const G4DynamicParticle* dp,
282  G4double xmin,
283  G4double maxEnergy)
284 {
285  G4double tmax = MaxSecondaryKinEnergy(dp);
286  G4double xmax = std::min(tmax, maxEnergy);
287  if(xmin >= xmax) { return; }
288 
289  G4double kineticEnergy = dp->GetKineticEnergy();
290  G4double energy = kineticEnergy + mass;
291  G4double energy2 = energy*energy;
292  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
293  G4double grej = 1.0;
294  G4double deltaKinEnergy, f;
295 
296  G4ThreeVector direction = dp->GetMomentumDirection();
297 
298  // sampling follows ...
299  do {
300  G4double q = G4UniformRand();
301  deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
302 
303  f = 1.0 - beta2*deltaKinEnergy/tmax;
304 
305  if(f > grej) {
306  G4cout << "G4BraggIonModel::SampleSecondary Warning! "
307  << "Majorant " << grej << " < "
308  << f << " for e= " << deltaKinEnergy
309  << G4endl;
310  }
311 
312  } while( grej*G4UniformRand() >= f );
313 
314  G4double deltaMomentum =
315  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
316  G4double totMomentum = energy*sqrt(beta2);
317  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
318  (deltaMomentum * totMomentum);
319  if(cost > 1.0) { cost = 1.0; }
320  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
321 
322  G4double phi = twopi * G4UniformRand() ;
323 
324  G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
325  deltaDirection.rotateUz(direction);
326 
327  // create G4DynamicParticle object for delta ray
328  G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,
329  deltaKinEnergy);
330 
331  vdp->push_back(delta);
332 
333  // Change kinematics of primary particle
334  kineticEnergy -= deltaKinEnergy;
335  G4ThreeVector finalP = direction*totMomentum - deltaDirection*deltaMomentum;
336  finalP = finalP.unit();
337 
338  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
339  fParticleChange->SetProposedMomentumDirection(finalP);
340 }
341 
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
343 
344 G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
345  G4double kinEnergy)
346 {
347  if(pd != particle) { SetParticle(pd); }
348  G4double tau = kinEnergy/mass;
349  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
350  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
351  return tmax;
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355 
356 G4bool G4BraggIonModel::HasMaterial(const G4Material* material)
357 {
358  G4String chFormula = material->GetChemicalFormula();
359  if("" == chFormula) { return false; }
360 
361  // ICRU Report N49, 1993. Ziegler model for He.
362  const size_t numberOfMolecula = 11;
363  const G4String molName[numberOfMolecula] = {
364  "CaF_2", "Cellulose_Nitrate", "LiF", "Policarbonate",
365  "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polymethly_Methacralate",
366  "Polysterene", "SiO_2", "NaI", "H_2O",
367  "Graphite" };
368  const G4int idxASTAR[numberOfMolecula] = {
369  17, 19, 33, 51,
370  52, 54,
371  56, 62, 43, 71,
372  13};
373 
374  // Search for the material in the table
375  for (size_t i=0; i<numberOfMolecula; ++i) {
376  if (chFormula == molName[i]) {
377  iMolecula = -1;
378  iASTAR = idxASTAR[i];
379  return true;
380  }
381  }
382  return false ;
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
386 
387 G4double G4BraggIonModel::StoppingPower(const G4Material* material,
388  G4double kineticEnergy)
389 {
390  G4double ionloss = 0.0 ;
391 
392  if (iMolecula >= 0) {
393 
394  // The data and the fit from:
395  // ICRU Report N49, 1993. Ziegler's model for alpha
396  // He energy in internal units of parametrisation formula (MeV)
397 
398  G4double T = kineticEnergy*rateMassHe2p/MeV ;
399 
400  static G4double a[11][5] = {
401  {9.43672, 0.54398, 84.341, 1.3705, 57.422},
402  {67.1503, 0.41409, 404.512, 148.97, 20.99},
403  {5.11203, 0.453, 36.718, 50.6, 28.058},
404  {61.793, 0.48445, 361.537, 57.889, 50.674},
405  {7.83464, 0.49804, 160.452, 3.192, 0.71922},
406  {19.729, 0.52153, 162.341, 58.35, 25.668},
407  {26.4648, 0.50112, 188.913, 30.079, 16.509},
408  {7.8655, 0.5205, 63.96, 51.32, 67.775},
409  {8.8965, 0.5148, 339.36, 1.7205, 0.70423},
410  {2.959, 0.53255, 34.247, 60.655, 15.153},
411  {3.80133, 0.41590, 12.9966, 117.83, 242.28} };
412 
413  static G4double atomicWeight[11] = {
414  101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
415  104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
416 
417  G4int i = iMolecula;
418 
419  // Free electron gas model
420  if ( T < 0.001 ) {
421  G4double slow = a[i][0] ;
422  G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
423  * a[i][2]*1000.0 ;
424  ionloss = slow*shigh / (slow + shigh) ;
425  ionloss *= sqrt(T*1000.0) ;
426 
427  // Main parametrisation
428  } else {
429  G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
430  G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
431  ionloss = slow*shigh / (slow + shigh) ;
432  /*
433  G4cout << "## " << i << ". T= " << T << " slow= " << slow
434  << " a0= " << a[i][0] << " a1= " << a[i][1]
435  << " shigh= " << shigh
436  << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
437  << G4endl;
438  */
439  }
440  if ( ionloss < 0.0) ionloss = 0.0 ;
441 
442  // He effective charge
443  G4double aa = atomicWeight[iMolecula];
444  ionloss /= (HeEffChargeSquare(0.5*aa, T)*aa);
445 
446  // pure material (normally not the case for this function)
447  } else if(1 == (material->GetNumberOfElements())) {
448  G4double z = material->GetZ() ;
449  ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
450  }
451 
452  return ionloss;
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456 
457 G4double G4BraggIonModel::ElectronicStoppingPower(G4double z,
458  G4double kineticEnergy) const
459 {
460  G4double ionloss ;
461  G4int i = G4int(z)-1 ; // index of atom
462  if(i < 0) i = 0 ;
463  if(i > 91) i = 91 ;
464 
465  // The data and the fit from:
466  // ICRU Report 49, 1993. Ziegler's type of parametrisations.
467  // Proton kinetic energy for parametrisation (keV/amu)
468 
469  // He energy in internal units of parametrisation formula (MeV)
470  G4double T = kineticEnergy*rateMassHe2p/MeV ;
471 
472  static G4double a[92][5] = {
473  {0.35485, 0.6456, 6.01525, 20.8933, 4.3515
474  },{ 0.58, 0.59, 6.3, 130.0, 44.07
475  },{ 1.42, 0.49, 12.25, 32.0, 9.161
476  },{ 2.206, 0.51, 15.32, 0.25, 8.995 //Be Ziegler77
477  // },{ 2.1895, 0.47183,7.2362, 134.30, 197.96 //Be from ICRU
478  },{ 3.691, 0.4128, 18.48, 50.72, 9.0
479  },{ 3.83523, 0.42993,12.6125, 227.41, 188.97
480  },{ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768
481  },{ 2.81015, 0.4759, 50.0253, 10.556, 1.0382
482  },{ 1.533, 0.531, 40.44, 18.41, 2.718
483  },{ 2.303, 0.4861, 37.01, 37.96, 5.092
484  // Z= 11-20
485  },{ 9.894, 0.3081, 23.65, 0.384, 92.93
486  },{ 4.3, 0.47, 34.3, 3.3, 12.74
487  },{ 2.5, 0.625, 45.7, 0.1, 4.359
488  },{ 2.1, 0.65, 49.34, 1.788, 4.133
489  },{ 1.729, 0.6562, 53.41, 2.405, 3.845
490  },{ 1.402, 0.6791, 58.98, 3.528, 3.211
491  },{ 1.117, 0.7044, 69.69, 3.705, 2.156
492  },{ 2.291, 0.6284, 73.88, 4.478, 2.066
493  },{ 8.554, 0.3817, 83.61, 11.84, 1.875
494  },{ 6.297, 0.4622, 65.39, 10.14, 5.036
495  // Z= 21-30
496  },{ 5.307, 0.4918, 61.74, 12.4, 6.665
497  },{ 4.71, 0.5087, 65.28, 8.806, 5.948
498  },{ 6.151, 0.4524, 83.0, 18.31, 2.71
499  },{ 6.57, 0.4322, 84.76, 15.53, 2.779
500  },{ 5.738, 0.4492, 84.6, 14.18, 3.101
501  },{ 5.013, 0.4707, 85.8, 16.55, 3.211
502  },{ 4.32, 0.4947, 76.14, 10.85, 5.441
503  },{ 4.652, 0.4571, 80.73, 22.0, 4.952
504  },{ 3.114, 0.5236, 76.67, 7.62, 6.385
505  },{ 3.114, 0.5236, 76.67, 7.62, 7.502
506  // Z= 31-40
507  },{ 3.114, 0.5236, 76.67, 7.62, 8.514
508  },{ 5.746, 0.4662, 79.24, 1.185, 7.993
509  },{ 2.792, 0.6346, 106.1, 0.2986, 2.331
510  },{ 4.667, 0.5095, 124.3, 2.102, 1.667
511  },{ 2.44, 0.6346, 105.0, 0.83, 2.851
512  },{ 1.413, 0.7377, 147.9, 1.466, 1.016
513  },{ 11.72, 0.3826, 102.8, 9.231, 4.371
514  },{ 7.126, 0.4804, 119.3, 5.784, 2.454
515  },{ 11.61, 0.3955, 146.7, 7.031, 1.423
516  },{ 10.99, 0.41, 163.9, 7.1, 1.052
517  // Z= 41-50
518  },{ 9.241, 0.4275, 163.1, 7.954, 1.102
519  },{ 9.276, 0.418, 157.1, 8.038, 1.29
520  },{ 3.999, 0.6152, 97.6, 1.297, 5.792
521  },{ 4.306, 0.5658, 97.99, 5.514, 5.754
522  },{ 3.615, 0.6197, 86.26, 0.333, 8.689
523  },{ 5.8, 0.49, 147.2, 6.903, 1.289
524  },{ 5.6, 0.49, 130.0, 10.0, 2.844
525  },{ 3.55, 0.6068, 124.7, 1.112, 3.119
526  },{ 3.6, 0.62, 105.8, 0.1692, 6.026
527  },{ 5.4, 0.53, 103.1, 3.931, 7.767
528  // Z= 51-60
529  },{ 3.97, 0.6459, 131.8, 0.2233, 2.723
530  },{ 3.65, 0.64, 126.8, 0.6834, 3.411
531  },{ 3.118, 0.6519, 164.9, 1.208, 1.51
532  },{ 3.949, 0.6209, 200.5, 1.878, 0.9126
533  },{ 14.4, 0.3923, 152.5, 8.354, 2.597
534  },{ 10.99, 0.4599, 138.4, 4.811, 3.726
535  },{ 16.6, 0.3773, 224.1, 6.28, 0.9121
536  },{ 10.54, 0.4533, 159.3, 4.832, 2.529
537  },{ 10.33, 0.4502, 162.0, 5.132, 2.444
538  },{ 10.15, 0.4471, 165.6, 5.378, 2.328
539  // Z= 61-70
540  },{ 9.976, 0.4439, 168.0, 5.721, 2.258
541  },{ 9.804, 0.4408, 176.2, 5.675, 1.997
542  },{ 14.22, 0.363, 228.4, 7.024, 1.016
543  },{ 9.952, 0.4318, 233.5, 5.065, 0.9244
544  },{ 9.272, 0.4345, 210.0, 4.911, 1.258
545  },{ 10.13, 0.4146, 225.7, 5.525, 1.055
546  },{ 8.949, 0.4304, 213.3, 5.071, 1.221
547  },{ 11.94, 0.3783, 247.2, 6.655, 0.849
548  },{ 8.472, 0.4405, 195.5, 4.051, 1.604
549  },{ 8.301, 0.4399, 203.7, 3.667, 1.459
550  // Z= 71-80
551  },{ 6.567, 0.4858, 193.0, 2.65, 1.66
552  },{ 5.951, 0.5016, 196.1, 2.662, 1.589
553  },{ 7.495, 0.4523, 251.4, 3.433, 0.8619
554  },{ 6.335, 0.4825, 255.1, 2.834, 0.8228
555  },{ 4.314, 0.5558, 214.8, 2.354, 1.263
556  },{ 4.02, 0.5681, 219.9, 2.402, 1.191
557  },{ 3.836, 0.5765, 210.2, 2.742, 1.305
558  },{ 4.68, 0.5247, 244.7, 2.749, 0.8962
559  },{ 2.892, 0.6204, 208.6, 2.415, 1.416 //Au Z77
560  // },{ 3.223, 0.5883, 232.7, 2.954, 1.05 //Au ICRU
561  },{ 2.892, 0.6204, 208.6, 2.415, 1.416
562  // Z= 81-90
563  },{ 4.728, 0.5522, 217.0, 3.091, 1.386
564  },{ 6.18, 0.52, 170.0, 4.0, 3.224
565  },{ 9.0, 0.47, 198.0, 3.8, 2.032
566  },{ 2.324, 0.6997, 216.0, 1.599, 1.399
567  },{ 1.961, 0.7286, 223.0, 1.621, 1.296
568  },{ 1.75, 0.7427, 350.1, 0.9789, 0.5507
569  },{ 10.31, 0.4613, 261.2, 4.738, 0.9899
570  },{ 7.962, 0.519, 235.7, 4.347, 1.313
571  },{ 6.227, 0.5645, 231.9, 3.961, 1.379
572  },{ 5.246, 0.5947, 228.6, 4.027, 1.432
573  // Z= 91-92
574  },{ 5.408, 0.5811, 235.7, 3.961, 1.358
575  },{ 5.218, 0.5828, 245.0, 3.838, 1.25}
576  };
577 
578  // Free electron gas model
579  if ( T < 0.001 ) {
580  G4double slow = a[i][0] ;
581  G4double shigh = log( 1.0 + a[i][3]*1000.0 + a[i][4]*0.001 )
582  * a[i][2]*1000.0 ;
583  ionloss = slow*shigh / (slow + shigh) ;
584  ionloss *= sqrt(T*1000.0) ;
585 
586  // Main parametrisation
587  } else {
588  G4double slow = a[i][0] * pow((T*1000.0), a[i][1]) ;
589  G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
590  ionloss = slow*shigh / (slow + shigh) ;
591  /*
592  G4cout << "## " << i << ". T= " << T << " slow= " << slow
593  << " a0= " << a[i][0] << " a1= " << a[i][1]
594  << " shigh= " << shigh
595  << " dedx= " << ionloss << " q^2= " << HeEffChargeSquare(z, T*MeV)
596  << G4endl;
597  */
598  }
599  if ( ionloss < 0.0) { ionloss = 0.0; }
600 
601  // He effective charge
602  ionloss /= HeEffChargeSquare(z, T);
603 
604  return ionloss;
605 }
606 
607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
608 
609 G4double G4BraggIonModel::DEDX(const G4Material* material,
610  G4double kineticEnergy)
611 {
612  G4double eloss = 0.0;
613  // check DB
614  if(material != currentMaterial) {
615  currentMaterial = material;
616  iASTAR = -1;
617  iMolecula = -1;
618  if( !HasMaterial(material) ) { iASTAR = astar.GetIndex(material); }
619  }
620 
621  const G4int numberOfElements = material->GetNumberOfElements();
622  const G4double* theAtomicNumDensityVector =
623  material->GetAtomicNumDensityVector();
624 
625  if( iASTAR >= 0 ) {
626  G4double T = kineticEnergy*rateMassHe2p;
627  return astar.GetElectronicDEDX(iASTAR, T)*material->GetDensity()/
628  HeEffChargeSquare(astar.GetEffectiveZ(iASTAR), T/MeV);
629 
630  } else if(iMolecula >= 0) {
631 
632  eloss = StoppingPower(material, kineticEnergy)*
633  material->GetDensity()/amu;
634 
635  // pure material
636  } else if(1 == numberOfElements) {
637 
638  G4double z = material->GetZ();
639  eloss = ElectronicStoppingPower(z, kineticEnergy)
640  * (material->GetTotNbOfAtomsPerVolume());
641 
642  // Brugg's rule calculation
643  } else {
644  const G4ElementVector* theElementVector =
645  material->GetElementVector() ;
646 
647  // loop for the elements in the material
648  for (G4int i=0; i<numberOfElements; i++)
649  {
650  const G4Element* element = (*theElementVector)[i] ;
651  eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
652  * theAtomicNumDensityVector[i];
653  }
654  }
655  return eloss*theZieglerFactor;
656 }
657 
658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
659 
660 G4double G4BraggIonModel::HeEffChargeSquare(G4double z,
661  G4double kinEnergyHeInMeV) const
662 {
663  // The aproximation of He effective charge from:
664  // J.F.Ziegler, J.P. Biersack, U. Littmark
665  // The Stopping and Range of Ions in Matter,
666  // Vol.1, Pergamon Press, 1985
667 
668  static G4double c[6] = {0.2865, 0.1266, -0.001429,
669  0.02402,-0.01135, 0.001475};
670 
671  G4double e = std::max(0.0,std::log(kinEnergyHeInMeV*massFactor));
672  G4double x = c[0] ;
673  G4double y = 1.0 ;
674  for (G4int i=1; i<6; i++) {
675  y *= e ;
676  x += y * c[i] ;
677  }
678 
679  G4double w = 7.6 - e ;
680  w = 1.0 + (0.007 + 0.00005*z) * exp( -w*w ) ;
681  w = 4.0 * (1.0 - exp(-x)) * w * w ;
682 
683  return w;
684 }
685 
686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
687 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4