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