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