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