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