Geant4  10.01.p03
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 83008 2014-07-24 14:49:52Z 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.5 == 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  // sampling follows ...
284  do {
285  G4double q = G4UniformRand();
286  deltaKinEnergy = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
287 
288  f = 1.0 - beta2*deltaKinEnergy/tmax;
289 
290  if(f > grej) {
291  G4cout << "G4BraggModel::SampleSecondary Warning! "
292  << "Majorant " << grej << " < "
293  << f << " for e= " << deltaKinEnergy
294  << G4endl;
295  }
296 
297  } while( grej*G4UniformRand() >= f );
298 
299  G4ThreeVector deltaDirection;
300 
302  const G4Material* mat = couple->GetMaterial();
303  G4int Z = SelectRandomAtomNumber(mat);
304 
305  deltaDirection =
306  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
307 
308  } else {
309 
310  G4double deltaMomentum =
311  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
312  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
313  (deltaMomentum * dp->GetTotalMomentum());
314  if(cost > 1.0) { cost = 1.0; }
315  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
316 
317  G4double phi = twopi * G4UniformRand() ;
318 
319  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
320  deltaDirection.rotateUz(dp->GetMomentumDirection());
321  }
322 
323  // create G4DynamicParticle object for delta ray
324  G4DynamicParticle* delta =
325  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
326 
327  // Change kinematics of primary particle
328  kineticEnergy -= deltaKinEnergy;
329  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
330  finalP = finalP.unit();
331 
334 
335  vdp->push_back(delta);
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339 
341  G4double kinEnergy)
342 {
343  if(pd != particle) { SetParticle(pd); }
344  G4double tau = kinEnergy/mass;
345  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
346  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
347  return tmax;
348 }
349 
350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
351 
353 {
354  return false;
355  /*
356  G4String chFormula = material->GetChemicalFormula();
357  if("" == chFormula) { return false; }
358 
359  // ICRU Report N49, 1993. Power's model for H
360  static const size_t numberOfMolecula = 11;
361  static const G4String molName[numberOfMolecula] = {
362  "Al_2O_3", "CO_2", "CH_4",
363  "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
364  "C_3H_8", "SiO_2", "H_2O",
365  "H_2O-Gas", "Graphite" } ;
366 
367  // Search for the material in the table
368  for (size_t i=0; i<numberOfMolecula; ++i) {
369  if (chFormula == molName[i]) {
370  iPSTAR = fPSTAR->GetIndex(matName[i]);
371  break;
372  }
373  }
374  return (iPSTAR >= 0);
375  */
376 }
377 
378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379 
381  G4double kineticEnergy)
382 {
383  G4double ionloss = 0.0 ;
384 
385  if (iMolecula >= 0) {
386 
387  // The data and the fit from:
388  // ICRU Report N49, 1993. Ziegler's model for protons.
389  // Proton kinetic energy for parametrisation (keV/amu)
390 
391  G4double T = kineticEnergy/(keV*protonMassAMU) ;
392 
393  static const G4double a[11][5] = {
394  {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
395  {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
396  {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
397  {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
398  {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
399  {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
400  {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
401  {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
402  {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
403  {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
404  {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
405 
406  static const G4double atomicWeight[11] = {
407  101.96128, 44.0098, 16.0426, 28.0536, 42.0804,
408  104.1512, 44.665, 60.0843, 18.0152, 18.0152, 12.0};
409 
410  if ( T < 10.0 ) {
411  ionloss = a[iMolecula][0] * sqrt(T) ;
412 
413  } else if ( T < 10000.0 ) {
414  G4double slow = a[iMolecula][1] * G4Exp(G4Log(T)* 0.45);
415  G4double shigh = G4Log( 1.0 + a[iMolecula][3]/T
416  + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
417  ionloss = slow*shigh / (slow + shigh) ;
418  }
419 
420  if ( ionloss < 0.0) ionloss = 0.0 ;
421  if ( 10 == iMolecula ) {
422  if (T < 100.0) {
423  ionloss *= (1.0+0.023+0.0066*G4Log(T)*invLog10);
424  }
425  else if (T < 700.0) {
426  ionloss *=(1.0+0.089-0.0248*G4Log(T-99.)*invLog10);
427  }
428  else if (T < 10000.0) {
429  ionloss *=(1.0+0.089-0.0248*G4Log(700.-99.)*invLog10);
430  }
431  }
432  ionloss /= atomicWeight[iMolecula];
433 
434  // pure material (normally not the case for this function)
435  } else if(1 == (material->GetNumberOfElements())) {
436  G4double z = material->GetZ() ;
437  ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
438  }
439 
440  return ionloss;
441 }
442 
443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
444 
446  G4double kineticEnergy) const
447 {
448  G4double ionloss ;
449  G4int i = G4int(z)-1 ; // index of atom
450  if(i < 0) i = 0 ;
451  if(i > 91) i = 91 ;
452 
453  // The data and the fit from:
454  // ICRU Report 49, 1993. Ziegler's type of parametrisations.
455  // Proton kinetic energy for parametrisation (keV/amu)
456 
457  G4double T = kineticEnergy/(keV*protonMassAMU) ;
458 
459  static const G4double a[92][5] = {
460  {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
461  {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
462  {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
463  {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
464  {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
465  {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
466  {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
467  {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
468  {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
469  {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
470  // Z= 11-20
471  {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
472  {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
473  {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
474  {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
475  {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
476  {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
477  {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
478  {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
479  {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
480  {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
481  // Z= 21-30
482  {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
483  {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
484  {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
485  {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
486  {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
487  {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
488  {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
489  {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
490  {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
491  {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
492  // Z= 31-40
493  {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
494  {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
495  {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
496  {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
497  {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
498  {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
499  {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
500  {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
501  {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
502  {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
503  // Z= 41-50
504  {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
505  {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
506  {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
507  {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
508  {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
509  {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
510  // {5.623, 6.354, 7160.0, 337.6, 0.013940}, // Ag Ziegler77
511  {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2}, // Ag ICRU49
512  {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
513  {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
514  {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
515  // Z= 51-60
516  {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
517  {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
518  {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
519  {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
520  {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
521  {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
522  {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
523  {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
524  {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
525  {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
526  // Z= 61-70
527  {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
528  {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
529  {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
530  {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
531  {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
532  {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
533  {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
534  {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
535  {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
536  {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
537  // Z= 71-80
538  {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
539  {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
540  {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
541  {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
542  {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
543  {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
544  {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
545  {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
546  // {4.856, 5.460, 18320.0, 438.5, 0.002542}, //Ziegler77
547  {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2}, //ICRU49
548  {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
549  // Z= 81-90
550  {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
551  {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
552  {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
553  {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
554  {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
555  {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
556  {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
557  {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
558  {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
559  {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
560  // Z= 91-92
561  {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
562  {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
563  };
564 
565  G4double fac = 1.0 ;
566 
567  // Carbon specific case for E < 40 keV
568  if ( T < 40.0 && 5 == i) {
569  fac = sqrt(T/40.0) ;
570  T = 40.0 ;
571 
572  // Free electron gas model
573  } else if ( T < 10.0 ) {
574  fac = sqrt(T*0.1) ;
575  T =10.0 ;
576  }
577 
578  // Main parametrisation
579  G4double slow = a[i][1] * G4Exp(G4Log(T) * 0.45) ;
580  G4double shigh = G4Log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
581  ionloss = slow*shigh*fac / (slow + shigh) ;
582 
583  if ( ionloss < 0.0) { ionloss = 0.0; }
584 
585  return ionloss;
586 }
587 
588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
589 
590 G4double G4BraggModel::DEDX(const G4Material* material, G4double kineticEnergy)
591 {
592  G4double eloss = 0.0;
593 
594  // check DB
595  if(material != currentMaterial) {
596  currentMaterial = material;
597  iPSTAR = -1;
598  iMolecula = -1;
599  if( !HasMaterial(material) ) { iPSTAR = fPSTAR->GetIndex(material); }
600 
601  //G4cout << "%%% " <<material->GetName() << " iMolecula= "
602  // << iMolecula << " iPSTAR= " << iPSTAR << G4endl;
603 
604  }
605 
606  const G4int numberOfElements = material->GetNumberOfElements();
607  const G4double* theAtomicNumDensityVector =
608  material->GetAtomicNumDensityVector();
609 
610  if( iPSTAR >= 0 ) {
611  return
612  fPSTAR->GetElectronicDEDX(iPSTAR, kineticEnergy)*material->GetDensity();
613 
614  } else if(iMolecula >= 0) {
615 
616  eloss = StoppingPower(material, kineticEnergy)*
617  material->GetDensity()/amu;
618 
619  // Pure material ICRU49 paralmeterisation
620  } else if(1 == numberOfElements) {
621 
622  G4double z = material->GetZ();
623  eloss = ElectronicStoppingPower(z, kineticEnergy)
624  * (material->GetTotNbOfAtomsPerVolume());
625 
626 
627  // Experimental data exist only for kinetic energy 125 keV
628  } else if( MolecIsInZiegler1988(material) ) {
629 
630  // Loop over elements - calculation based on Bragg's rule
631  G4double eloss125 = 0.0 ;
632  const G4ElementVector* theElementVector =
633  material->GetElementVector();
634 
635  // Loop for the elements in the material
636  for (G4int i=0; i<numberOfElements; i++) {
637  const G4Element* element = (*theElementVector)[i] ;
638  G4double z = element->GetZ() ;
639  eloss += ElectronicStoppingPower(z,kineticEnergy)
640  * theAtomicNumDensityVector[i] ;
641  eloss125 += ElectronicStoppingPower(z,125.0*keV)
642  * theAtomicNumDensityVector[i] ;
643  }
644 
645  // Chemical factor is taken into account
646  eloss *= ChemicalFactor(kineticEnergy, eloss125) ;
647 
648  // Brugg's rule calculation
649  } else {
650  const G4ElementVector* theElementVector =
651  material->GetElementVector() ;
652 
653  // loop for the elements in the material
654  for (G4int i=0; i<numberOfElements; i++)
655  {
656  const G4Element* element = (*theElementVector)[i] ;
657  eloss += ElectronicStoppingPower(element->GetZ(), kineticEnergy)
658  * theAtomicNumDensityVector[i];
659  }
660  }
661  return eloss*theZieglerFactor;
662 }
663 
664 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665 
667 {
668  // The list of molecules from
669  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
670  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
671 
672  G4String myFormula = G4String(" ") ;
673  const G4String chFormula = material->GetChemicalFormula() ;
674  if (myFormula == chFormula ) { return false; }
675 
676  // There are no evidence for difference of stopping power depended on
677  // phase of the compound except for water. The stopping power of the
678  // water in gas phase can be predicted using Bragg's rule.
679  //
680  // No chemical factor for water-gas
681 
682  myFormula = G4String("H_2O") ;
683  const G4State theState = material->GetState() ;
684  if( theState == kStateGas && myFormula == chFormula) return false ;
685 
686 
687  // The coffecient from Table.4 of Ziegler & Manoyan
688  static const G4double HeEff = 2.8735 ;
689 
690  static const size_t numberOfMolecula = 53;
691  static const G4String nameOfMol[53] = {
692  "H_2O", "C_2H_4O", "C_3H_6O", "C_2H_2", "C_H_3OH",
693  "C_2H_5OH", "C_3H_7OH", "C_3H_4", "NH_3", "C_14H_10",
694  "C_6H_6", "C_4H_10", "C_4H_6", "C_4H_8O", "CCl_4",
695  "CF_4", "C_6H_8", "C_6H_12", "C_6H_10O", "C_6H_10",
696  "C_8H_16", "C_5H_10", "C_5H_8", "C_3H_6-Cyclopropane","C_2H_4F_2",
697  "C_2H_2F_2", "C_4H_8O_2", "C_2H_6", "C_2F_6", "C_2H_6O",
698  "C_3H_6O", "C_4H_10O", "C_2H_4", "C_2H_4O", "C_2H_4S",
699  "SH_2", "CH_4", "CCLF_3", "CCl_2F_2", "CHCl_2F",
700  "(CH_3)_2S", "N_2O", "C_5H_10O", "C_8H_6", "(CH_2)_N",
701  "(C_3H_6)_N","(C_8H_8)_N", "C_3H_8", "C_3H_6-Propylene", "C_3H_6O",
702  "C_3H_6S", "C_4H_4S", "C_7H_8"
703  };
704 
705  static const G4double expStopping[numberOfMolecula] = {
706  66.1, 190.4, 258.7, 42.2, 141.5,
707  210.9, 279.6, 198.8, 31.0, 267.5,
708  122.8, 311.4, 260.3, 328.9, 391.3,
709  206.6, 374.0, 422.0, 432.0, 398.0,
710  554.0, 353.0, 326.0, 74.6, 220.5,
711  197.4, 362.0, 170.0, 330.5, 211.3,
712  262.3, 349.6, 51.3, 187.0, 236.9,
713  121.9, 35.8, 247.0, 292.6, 268.0,
714  262.3, 49.0, 398.9, 444.0, 22.91,
715  68.0, 155.0, 84.0, 74.2, 254.7,
716  306.8, 324.4, 420.0
717  } ;
718 
719  static const G4double expCharge[53] = {
720  HeEff, HeEff, HeEff, 1.0, HeEff,
721  HeEff, HeEff, HeEff, 1.0, 1.0,
722  1.0, HeEff, HeEff, HeEff, HeEff,
723  HeEff, HeEff, HeEff, HeEff, HeEff,
724  HeEff, HeEff, HeEff, 1.0, HeEff,
725  HeEff, HeEff, HeEff, HeEff, HeEff,
726  HeEff, HeEff, 1.0, HeEff, HeEff,
727  HeEff, 1.0, HeEff, HeEff, HeEff,
728  HeEff, 1.0, HeEff, HeEff, 1.0,
729  1.0, 1.0, 1.0, 1.0, HeEff,
730  HeEff, HeEff, HeEff
731  } ;
732 
733  static const G4double numberOfAtomsPerMolecula[53] = {
734  3.0, 7.0, 10.0, 4.0, 6.0,
735  9.0, 12.0, 7.0, 4.0, 24.0,
736  12.0, 14.0, 10.0, 13.0, 5.0,
737  5.0, 14.0, 18.0, 17.0, 17.0,
738  24.0, 15.0, 13.0, 9.0, 8.0,
739  6.0, 14.0, 8.0, 8.0, 9.0,
740  10.0, 15.0, 6.0, 7.0, 7.0,
741  3.0, 5.0, 5.0, 5.0, 5.0,
742  9.0, 3.0, 16.0, 14.0, 3.0,
743  9.0, 16.0, 11.0, 9.0, 10.0,
744  10.0, 9.0, 15.0
745  } ;
746 
747  // Search for the compaund in the table
748  for (size_t i=0; i<numberOfMolecula; i++)
749  {
750  if(chFormula == nameOfMol[i]) {
751  G4double exp125 = expStopping[i] *
752  (material->GetTotNbOfAtomsPerVolume()) /
753  (expCharge[i] * numberOfAtomsPerMolecula[i]) ;
754  SetExpStopPower125(exp125);
755  return true;
756  }
757  }
758 
759  return false;
760 }
761 
762 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
763 
765  G4double eloss125) const
766 {
767  // Approximation of Chemical Factor according to
768  // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds,
769  // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228.
770 
771  G4double gamma = 1.0 + kineticEnergy/proton_mass_c2 ;
772  G4double gamma25 = 1.0 + 25.0*keV /proton_mass_c2 ;
773  G4double gamma125 = 1.0 + 125.0*keV/proton_mass_c2 ;
774  G4double beta = sqrt(1.0 - 1.0/(gamma*gamma)) ;
775  G4double beta25 = sqrt(1.0 - 1.0/(gamma25*gamma25)) ;
776  G4double beta125 = sqrt(1.0 - 1.0/(gamma125*gamma125)) ;
777 
778  G4double factor = 1.0 + (expStopPower125/eloss125 - 1.0) *
779  (1.0 + G4Exp( 1.48 * ( beta125/beta25 - 7.0 ) ) ) /
780  (1.0 + G4Exp( 1.48 * ( beta/beta25 - 7.0 ) ) ) ;
781 
782  return factor ;
783 }
784 
785 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
786 
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:464
G4double ratio
static const double MeV
Definition: G4SIunits.hh:193
G4double GetZ() const
Definition: G4Material.cc:632
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:120
static const double cm2
Definition: G4SIunits.hh:107
std::vector< G4Element * > G4ElementVector
G4bool HasMaterial(const G4Material *material)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4BraggModel(const G4ParticleDefinition *p=0, const G4String &nam="Bragg")
Definition: G4BraggModel.cc:85
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:599
G4double GetDensity() const
Definition: G4Material.hh:180
G4bool IsMaster() const
Definition: G4VEmModel.hh:700
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:592
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:707
G4double GetTotalMomentum() const
G4double expStopPower125
G4double chargeSquare
#define G4UniformRand()
Definition: Randomize.hh:93
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:679
bool G4bool
Definition: G4Types.hh:79
G4EmCorrections * EmCorrections()
const G4ThreeVector & GetMomentumDirection() const
G4bool isInitialised
virtual G4double GetParticleCharge(const G4ParticleDefinition *p, const G4Material *mat, G4double kineticEnergy)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static const double GeV
Definition: G4SIunits.hh:196
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:194
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)
G4double energy(const ThreeVector &p, const G4double m)
static const G4double invLog10
Definition: G4BraggModel.cc:82
G4EmCorrections * corr
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:606
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:195
G4double theZieglerFactor
G4State GetState() const
Definition: G4Material.hh:181
G4double protonMassAMU
double G4double
Definition: G4Types.hh:76
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:763
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:546
G4double ElectronicStoppingPower(G4double z, G4double kineticEnergy) const