source: trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc @ 1197

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

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