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