ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ICRU73QOModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ICRU73QOModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4ICRU73QOModel
33 //
34 // Author: Alexander Bagulya
35 //
36 // Creation date: 21.05.2010
37 //
38 // Modifications:
39 //
40 //
41 // -------------------------------------------------------------------
42 //
43 
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 #include "G4ICRU73QOModel.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "Randomize.hh"
52 #include "G4Electron.hh"
53 #include "G4ParticleChangeForLoss.hh"
54 #include "G4LossTableManager.hh"
55 #include "G4AntiProton.hh"
56 #include "G4DeltaAngle.hh"
57 #include "G4Log.hh"
58 #include "G4Exp.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 using namespace std;
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 G4ICRU73QOModel::G4ICRU73QOModel(const G4ParticleDefinition* p, const G4String& nam)
67  : G4VEmModel(nam),
68  particle(nullptr),
69  isInitialised(false)
70 {
71  mass = charge = chargeSquare = massRate = ratio = 0.0;
72  if(p) { SetParticle(p); }
73  SetHighEnergyLimit(10.0*MeV);
74 
75  lowestKinEnergy = 5.0*keV;
76 
77  sizeL0 = 67;
78  sizeL1 = 22;
79  sizeL2 = 14;
80 
81  theElectron = G4Electron::Electron();
82 
83  for (G4int i = 0; i < 100; ++i)
84  {
85  indexZ[i] = -1;
86  }
87  for(G4int i = 0; i < NQOELEM; ++i)
88  {
89  if(ZElementAvailable[i] > 0) {
90  indexZ[ZElementAvailable[i]] = i;
91  }
92  }
93  fParticleChange = nullptr;
94  denEffData = nullptr;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 void G4ICRU73QOModel::Initialise(const G4ParticleDefinition* p,
100  const G4DataVector&)
101 {
102  if(p != particle) SetParticle(p);
103 
104  // always false before the run
105  SetDeexcitationFlag(false);
106 
107  if(!isInitialised) {
108  isInitialised = true;
109 
110  if(UseAngularGeneratorFlag() && !GetAngularDistribution()) {
111  SetAngularDistribution(new G4DeltaAngle());
112  }
113 
114  G4String pname = particle->GetParticleName();
115  fParticleChange = GetParticleChangeForLoss();
116  const G4MaterialTable* mtab = G4Material::GetMaterialTable();
117  denEffData = (*mtab)[0]->GetIonisation()->GetDensityEffectData();
118  }
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
123 G4double G4ICRU73QOModel::ComputeCrossSectionPerElectron(
124  const G4ParticleDefinition* p,
125  G4double kineticEnergy,
126  G4double cutEnergy,
127  G4double maxKinEnergy)
128 {
129  G4double cross = 0.0;
130  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
131  G4double maxEnergy = std::min(tmax,maxKinEnergy);
132  if(cutEnergy < maxEnergy) {
133 
134  G4double energy = kineticEnergy + mass;
135  G4double energy2 = energy*energy;
136  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
137  cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
138  - beta2*G4Log(maxEnergy/cutEnergy)/tmax;
139 
140  cross *= CLHEP::twopi_mc2_rcl2*chargeSquare/beta2;
141  }
142 
143  return cross;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
148 G4double G4ICRU73QOModel::ComputeCrossSectionPerAtom(
149  const G4ParticleDefinition* p,
150  G4double kineticEnergy,
151  G4double Z, G4double,
152  G4double cutEnergy,
153  G4double maxEnergy)
154 {
155  G4double cross = Z*ComputeCrossSectionPerElectron
156  (p,kineticEnergy,cutEnergy,maxEnergy);
157  return cross;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 G4double G4ICRU73QOModel::CrossSectionPerVolume(
163  const G4Material* material,
164  const G4ParticleDefinition* p,
165  G4double kineticEnergy,
166  G4double cutEnergy,
167  G4double maxEnergy)
168 {
169  G4double eDensity = material->GetElectronDensity();
170  G4double cross = eDensity*ComputeCrossSectionPerElectron
171  (p,kineticEnergy,cutEnergy,maxEnergy);
172  return cross;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
177 G4double G4ICRU73QOModel::ComputeDEDXPerVolume(const G4Material* material,
178  const G4ParticleDefinition* p,
179  G4double kineticEnergy,
180  G4double cutEnergy)
181 {
182  SetParticle(p);
183  G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
184  G4double tkin = kineticEnergy/massRate;
185  G4double dedx = 0.0;
186  if(tkin > lowestKinEnergy) { dedx = DEDX(material, tkin); }
187  else { dedx = DEDX(material, lowestKinEnergy)*sqrt(tkin/lowestKinEnergy); }
188 
189  if (cutEnergy < tmax) {
190 
191  G4double tau = kineticEnergy/mass;
192  G4double gam = tau + 1.0;
193  G4double bg2 = tau * (tau+2.0);
194  G4double beta2 = bg2/(gam*gam);
195  G4double x = cutEnergy/tmax;
196 
197  dedx += chargeSquare*( G4Log(x) + (1.0 - x)*beta2 ) * twopi_mc2_rcl2
198  * material->GetElectronDensity()/beta2;
199  }
200  if(dedx < 0.0) { dedx = 0.0; }
201  return dedx;
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205 
206 G4double G4ICRU73QOModel::DEDX(const G4Material* material,
207  G4double kineticEnergy)
208 {
209  G4double eloss = 0.0;
210  const G4int numberOfElements = material->GetNumberOfElements();
211  const G4double* theAtomicNumDensityVector =
212  material->GetAtomicNumDensityVector();
213 
214  // Bragg's rule calculation
215  const G4ElementVector* theElementVector =
216  material->GetElementVector() ;
217 
218  // loop for the elements in the material
219  for (G4int i=0; i<numberOfElements; ++i)
220  {
221  const G4Element* element = (*theElementVector)[i] ;
222  eloss += DEDXPerElement(element->GetZasInt(), kineticEnergy)
223  * theAtomicNumDensityVector[i] * element->GetZ();
224  }
225  return eloss;
226 }
227 
228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229 
230 G4double G4ICRU73QOModel::DEDXPerElement(G4int AtomicNumber,
231  G4double kineticEnergy)
232 {
233  G4int Z = std::min(AtomicNumber, 97);
234  G4int nbOfShells = std::max(GetNumberOfShells(Z), 1);
235 
236  G4double v = CLHEP::c_light * std::sqrt( 2.0*kineticEnergy/proton_mass_c2 );
237 
238  G4double fBetheVelocity = CLHEP::fine_structure_const*CLHEP::c_light/v;
239 
240  G4double tau = kineticEnergy/proton_mass_c2;
241  G4double gam = tau + 1.0;
242  G4double bg2 = tau * (tau+2.0);
243  G4double beta2 = bg2/(gam*gam);
244 
245  G4double l0Term = 0, l1Term = 0, l2Term = 0;
246 
247  for (G4int nos = 0; nos < nbOfShells; ++nos){
248 
249  G4double NormalizedEnergy = (2.0*CLHEP::electron_mass_c2*beta2) /
250  GetShellEnergy(Z,nos);
251 
252  G4double shStrength = GetShellStrength(Z,nos);
253 
254  G4double l0 = GetL0(NormalizedEnergy);
255  l0Term += shStrength * l0;
256 
257  G4double l1 = GetL1(NormalizedEnergy);
258  l1Term += shStrength * l1;
259 
260  G4double l2 = GetL2(NormalizedEnergy);
261  l2Term += shStrength * l2;
262 
263  }
264  G4double dedx = 2*CLHEP::twopi_mc2_rcl2*chargeSquare*factorBethe[Z]*
265  (l0Term + charge*fBetheVelocity*l1Term
266  + chargeSquare*fBetheVelocity*fBetheVelocity*l2Term)/beta2;
267  return dedx;
268 }
269 
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272 
273 G4double G4ICRU73QOModel::GetOscillatorEnergy(G4int Z,
274  G4int nbOfTheShell) const
275 {
276  G4int idx = denEffData->GetElementIndex(Z, kStateUndefined);
277  if(idx == -1) { idx = denEffData->GetElementIndex(Z-1, kStateUndefined); }
278  G4double PlasmaEnergy = denEffData->GetPlasmaEnergy(idx);
279 
280  G4double PlasmaEnergy2 = PlasmaEnergy * PlasmaEnergy;
281 
282  G4double plasmonTerm = 0.66667
283  * G4AtomicShells::GetNumberOfElectrons(Z,nbOfTheShell)
284  * PlasmaEnergy2 / (Z*Z) ;
285 
286  static const G4double exphalf = G4Exp(0.5);
287  G4double ionTerm = exphalf *
288  (G4AtomicShells::GetBindingEnergy(Z,nbOfTheShell)) ;
289  G4double ionTerm2 = ionTerm*ionTerm ;
290 
291  G4double oscShellEnergy = std::sqrt( ionTerm2 + plasmonTerm );
292 
293  return oscShellEnergy;
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297 
298 G4int G4ICRU73QOModel::GetNumberOfShells(G4int Z) const
299 {
300  G4int nShell = 0;
301 
302  if(indexZ[Z] >= 0) {
303  nShell = nbofShellsForElement[indexZ[Z]];
304  } else {
305  nShell = G4AtomicShells::GetNumberOfShells(Z);
306  }
307  return nShell;
308 }
309 
310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
311 
312 G4double G4ICRU73QOModel::GetShellEnergy(G4int Z, G4int nbOfTheShell) const
313 {
314  G4double shellEnergy = 0.;
315 
316  G4int idx = indexZ[Z];
317 
318  if(idx >= 0) {
319  shellEnergy = ShellEnergy[startElemIndex[idx] + nbOfTheShell]*CLHEP::eV;
320  } else {
321  shellEnergy = GetOscillatorEnergy(Z, nbOfTheShell);
322  }
323 
324  return shellEnergy;
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328 
329 G4double G4ICRU73QOModel::GetShellStrength(G4int Z, G4int nbOfTheShell) const
330 {
331  G4double shellStrength = 0.;
332 
333  G4int idx = indexZ[Z];
334 
335  if(idx >= 0) {
336  shellStrength = SubShellOccupation[startElemIndex[idx] + nbOfTheShell] / Z;
337  } else {
338  shellStrength = G4double(G4AtomicShells::GetNumberOfElectrons(Z,nbOfTheShell))/Z;
339  }
340 
341  return shellStrength;
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345 
346 G4double G4ICRU73QOModel::GetL0(G4double normEnergy) const
347 {
348  G4int n;
349 
350  for(n = 0; n < sizeL0; n++) {
351  if( normEnergy < L0[n][0] ) break;
352  }
353  if(0 == n) { n = 1; }
354  if(n >= sizeL0) { n = sizeL0 - 1; }
355 
356  G4double l0 = L0[n][1];
357  G4double l0p = L0[n-1][1];
358  G4double bethe = l0p + (l0 - l0p) * ( normEnergy - L0[n-1][0]) /
359  (L0[n][0] - L0[n-1][0]);
360 
361  return bethe ;
362 }
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
365 
366 G4double G4ICRU73QOModel::GetL1(G4double normEnergy) const
367 {
368  G4int n;
369 
370  for(n = 0; n < sizeL1; n++) {
371  if( normEnergy < L1[n][0] ) break;
372  }
373  if(0 == n) n = 1 ;
374  if(n >= sizeL1) n = sizeL1 - 1 ;
375 
376  G4double l1 = L1[n][1];
377  G4double l1p = L1[n-1][1];
378  G4double barkas= l1p + (l1 - l1p) * ( normEnergy - L1[n-1][0]) /
379  (L1[n][0] - L1[n-1][0]);
380 
381  return barkas;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
386 G4double G4ICRU73QOModel::GetL2(G4double normEnergy) const
387 {
388  G4int n;
389  for(n = 0; n < sizeL2; n++) {
390  if( normEnergy < L2[n][0] ) break;
391  }
392  if(0 == n) n = 1 ;
393  if(n >= sizeL2) n = sizeL2 - 1 ;
394 
395  G4double l2 = L2[n][1];
396  G4double l2p = L2[n-1][1];
397  G4double bloch = l2p + (l2 - l2p) * ( normEnergy - L2[n-1][0]) /
398  (L2[n][0] - L2[n-1][0]);
399 
400  return bloch;
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 
405 void G4ICRU73QOModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
406  const G4DynamicParticle*,
407  G4double&,
408  G4double&,
409  G4double)
410 {}
411 
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
413 
414 void G4ICRU73QOModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
415  const G4MaterialCutsCouple* couple,
416  const G4DynamicParticle* dp,
417  G4double xmin,
418  G4double maxEnergy)
419 {
420  G4double tmax = MaxSecondaryKinEnergy(dp);
421  G4double xmax = std::min(tmax, maxEnergy);
422  if(xmin >= xmax) { return; }
423 
424  G4double kineticEnergy = dp->GetKineticEnergy();
425  G4double energy = kineticEnergy + mass;
426  G4double energy2 = energy*energy;
427  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/energy2;
428  G4double grej = 1.0;
429  G4double deltaKinEnergy, f;
430 
431  G4ThreeVector direction = dp->GetMomentumDirection();
432 
433  // sampling follows ...
434  do {
435  G4double x = G4UniformRand();
436  deltaKinEnergy = xmin*xmax/(xmin*(1.0 - x) + xmax*x);
437 
438  f = 1.0 - beta2*deltaKinEnergy/tmax;
439 
440  if(f > grej) {
441  G4cout << "G4ICRU73QOModel::SampleSecondary Warning! "
442  << "Majorant " << grej << " < "
443  << f << " for e= " << deltaKinEnergy
444  << G4endl;
445  }
446 
447  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
448  } while( grej*G4UniformRand() >= f );
449 
450  G4ThreeVector deltaDirection;
451 
452  if(UseAngularGeneratorFlag()) {
453  const G4Material* mat = couple->GetMaterial();
454  G4int Z = SelectRandomAtomNumber(mat);
455 
456  deltaDirection =
457  GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
458 
459  } else {
460 
461  G4double deltaMomentum =
462  sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
463  G4double totMomentum = energy*sqrt(beta2);
464  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
465  (deltaMomentum * totMomentum);
466  if(cost > 1.0) { cost = 1.0; }
467  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
468 
469  G4double phi = twopi * G4UniformRand() ;
470 
471  deltaDirection.set(sint*cos(phi),sint*sin(phi), cost) ;
472  deltaDirection.rotateUz(direction);
473  }
474  // create G4DynamicParticle object for delta ray
475  G4DynamicParticle* delta =
476  new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
477 
478  // Change kinematics of primary particle
479  kineticEnergy -= deltaKinEnergy;
480  G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
481  finalP = finalP.unit();
482 
483  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
484  fParticleChange->SetProposedMomentumDirection(finalP);
485 
486  vdp->push_back(delta);
487 }
488 
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
490 
491 G4double G4ICRU73QOModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
492  G4double kinEnergy)
493 {
494  if(pd != particle) { SetParticle(pd); }
495  G4double tau = kinEnergy/mass;
496  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
497  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
498  return tmax;
499 }
500 
501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502 
503 const G4int G4ICRU73QOModel::ZElementAvailable[NQOELEM] =
504  {1,2,4,6,7,8,10,13,14,-18,
505  22,26,28,29,32,36,42,47,
506  50,54,73,74,78,79,82,92};
507 
508 const G4int G4ICRU73QOModel::nbofShellsForElement[NQOELEM] =
509  {1,1,2,3,3,3,3,4,5,4,
510  5,5,5,5,6,4,6,6,
511  7,6,6,8,7,7,9,9};
512 
513 const G4int G4ICRU73QOModel::startElemIndex[NQOELEM] =
514  {0,1,2,4,7,10,13,16,20,25,
515  29,34,39,44,49,55,59,65,
516  71,78,84,90,98,105,112,121};
517 
518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
519 
520 // SubShellOccupation = Z * ShellStrength
521 const G4double G4ICRU73QOModel::SubShellOccupation[NQODATA] =
522  {
523  1.000, // H 0
524  2.000, // He 1
525  1.930, 2.070, // Be 2-3
526  1.992, 1.841, 2.167, // C 4-6
527  1.741, 1.680, 3.579, // N 7-9
528  1.802, 1.849, 4.349, // O 10-12
529  1.788, 2.028, 6.184, // Ne 13-15
530  1.623, 2.147, 6.259, 2.971, // Al 16-19
531  1.631, 2.094, 6.588, 2.041, 1.646, // Si 20-24
532  1.535, 8.655, 1.706, 6.104, // Ar 25-28
533  1.581, 8.358, 8.183, 2.000, 1.878, // Ti 29-33
534  1.516, 8.325, 8.461, 6.579, 1.119, // Fe 34-38
535  1.422, 7.81, 8.385, 8.216, 2.167, // Ni 39-43
536  1.458, 8.049, 8.79, 9.695, 1.008, // Cu 44-48
537  1.442, 7.791, 7.837, 10.122, 2.463, 2.345, // Ge 49-54
538  1.645, 7.765, 19.192, 7.398, // Kr 55-58
539  1.313, 6.409, 19.229, 8.633, 5.036, 1.380, // Mo 59-64
540  1.295, 6.219, 18.751, 8.748, 10.184, 1.803, // Ag 65-70
541  1.277, 6.099, 20.386, 8.011, 10.007, 2.272, 1.948, // Sn 71-77
542  1.563, 6.312, 21.868, 5.762, 11.245, 7.250, // Xe 78-83
543  0.9198, 6.5408, 18.9727, 24.9149, 15.0161, 6.6284, // Ta 84-89
544  1.202, 5.582, 19.527, 18.741, 8.411, 14.387, 4.042, 2.108, // W 90-97
545  1.159, 5.467, 18.802, 33.905, 8.300, 9.342, 1.025, // Pt 98-104
546  1.124, 5.331, 18.078, 34.604, 8.127, 10.414, 1.322, // Au 105-111
547  2.000, 8.000, 18.000, 18.000, 14.000, 8.000, 10.000, 2.000, 2.000, // Pb 112-120
548  2.000, 8.000, 18.000, 32.000, 18.000, 8.000, 2.000, 1.000, 3.000 // U 121-129
549 };
550 
551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
552 
553 // ShellEnergy in eV
554 const G4double G4ICRU73QOModel::ShellEnergy[NQODATA] =
555  {
556  19.2, // H
557  41.8, // He
558  209.11, 21.68, // Be
559  486.2, 60.95, 23.43, // C
560  732.61, 100.646, 23.550, // N
561  965.1, 129.85, 31.60, // O
562  1525.9, 234.9, 56.18, // Ne
563  2701, 476.5, 150.42, 16.89, // Al
564  3206.1, 586.4, 186.8, 23.52, 14.91, // Si
565  5551.6, 472.43, 124.85, 22.332, // Ar
566  8554.6, 850.58, 93.47, 39.19, 19.46, // Ti
567  12254.7, 1279.29, 200.35, 49.19, 17.66, // Fe
568  14346.9, 1532.28, 262.71, 74.37, 23.03, // Ni
569  15438.5, 1667.96, 294.1, 70.69, 16.447, // Cu
570  19022.1, 2150.79, 455.79, 179.87, 57.89, 20.95, // Ge
571  24643, 2906.4, 366.85, 22.24, // Kr
572  34394, 4365.3, 589.36, 129.42, 35.59, 18.42, // Mo
573  43664.3, 5824.91, 909.79, 175.47, 54.89, 19.63, // Ag
574  49948, 6818.2, 1036.1, 172.65, 70.89, 33.87, 14.54, // Sn
575  58987, 8159, 1296.6, 356.75, 101.03, 16.52, // Xe
576  88926, 18012, 3210, 575, 108.7, 30.8, // Ta
577  115025.9, 17827.44, 3214.36, 750.41, 305.21, 105.50, 38.09, 21.25, // W
578  128342, 20254, 3601.8, 608.1, 115.0, 42.75, 17.04, // Pt
579  131872, 20903, 3757.4, 682.1, 105.2, 44.89, 17.575, // Au
580  154449, 25067, 5105.0, 987.44, 247.59, 188.1, 40.61, 19.2, 15.17, // Pb
581  167282, 27868, 6022.7, 1020.4, 244.81, 51.33, 13, 11.06, 14.43 // U
582 };
583 
584 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
585 
586 // Data for L0 from: Sigmund P., Haagerup U. Phys. Rev. A34 (1986) 892-910
587 const G4double G4ICRU73QOModel::L0[67][2] =
588 {
589  {0.00, 0.000001},
590  {0.10, 0.000001},
591  {0.12, 0.00001},
592  {0.14, 0.00005},
593  {0.16, 0.00014},
594  {0.18, 0.00030},
595  {0.20, 0.00057},
596  {0.25, 0.00189},
597  {0.30, 0.00429},
598  {0.35, 0.00784},
599  {0.40, 0.01248},
600  {0.45, 0.01811},
601  {0.50, 0.02462},
602  {0.60, 0.03980},
603  {0.70, 0.05731},
604  {0.80, 0.07662},
605  {0.90, 0.09733},
606  {1.00, 0.11916},
607  {1.20, 0.16532},
608  {1.40, 0.21376},
609  {1.60, 0.26362},
610  {1.80, 0.31428},
611  {2.00, 0.36532},
612  {2.50, 0.49272},
613  {3.00, 0.61765},
614  {3.50, 0.73863},
615  {4.00, 0.85496},
616  {4.50, 0.96634},
617  {5.00, 1.07272},
618  {6.00, 1.27086},
619  {7.00, 1.45075},
620  {8.00, 1.61412},
621  {9.00, 1.76277},
622  {10.00, 1.89836},
623  {12.00, 2.13625},
624  {14.00, 2.33787},
625  {16.00, 2.51093},
626  {18.00, 2.66134},
627  {20.00, 2.79358},
628  {25.00, 3.06539},
629  {30.00, 3.27902},
630  {35.00, 3.45430},
631  {40.00, 3.60281},
632  {45.00, 3.73167},
633  {50.00, 3.84555},
634  {60.00, 4.04011},
635  {70.00, 4.20264},
636  {80.00, 4.34229},
637  {90.00, 4.46474},
638  {100.00, 4.57378},
639  {120.00, 4.76155},
640  {140.00, 4.91953},
641  {160.00, 5.05590},
642  {180.00, 5.17588},
643  {200.00, 5.28299},
644  {250.00, 5.50925},
645  {300.00, 5.69364},
646  {350.00, 5.84926},
647  {400.00, 5.98388},
648  {450.00, 6.10252},
649  {500.00, 6.20856},
650  {600.00, 6.39189},
651  {700.00, 6.54677},
652  {800.00, 6.68084},
653  {900.00, 6.79905},
654  {1000.00, 6.90474}
655 };
656 
657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658 
659 // Data for L1 from: Mikkelsen H.H., Sigmund P. Phys. Rev. A40 (1989) 101-116
660 const G4double G4ICRU73QOModel::L1[22][2] =
661 {
662  {0.00, -0.000001},
663  {0.10, -0.00001},
664  {0.20, -0.00049},
665  {0.30, -0.00084},
666  {0.40, 0.00085},
667  {0.50, 0.00519},
668  {0.60, 0.01198},
669  {0.70, 0.02074},
670  {0.80, 0.03133},
671  {0.90, 0.04369},
672  {1.00, 0.06035},
673  {2.00, 0.24023},
674  {3.00, 0.44284},
675  {4.00, 0.62012},
676  {5.00, 0.77031},
677  {6.00, 0.90390},
678  {7.00, 1.02705},
679  {8.00, 1.10867},
680  {9.00, 1.17546},
681  {10.00, 1.21599},
682  {15.00, 1.24349},
683  {20.00, 1.16752}
684 };
685 
686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
687 
688 // Data for L2 from: Mikkelsen H.H. Nucl. Instr. Meth. B58 (1991) 136-148
689 const G4double G4ICRU73QOModel::L2[14][2] =
690 {
691  {0.00, 0.000001},
692  {0.10, 0.00001},
693  {0.20, 0.00000},
694  {0.40, -0.00120},
695  {0.60, -0.00036},
696  {0.80, 0.00372},
697  {1.00, 0.01298},
698  {2.00, 0.08296},
699  {4.00, 0.21953},
700  {6.00, 0.23903},
701  {8.00, 0.20893},
702  {10.00, 0.10879},
703  {20.00, -0.88409},
704  {40.00, -1.13902}
705 };
706 
707 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
708 
709 // Correction obtained by V.Ivanchenko using G4BetheBlochModel
710 const G4double G4ICRU73QOModel::factorBethe[99] = { 1.0,
711 0.9637, 0.9872, 0.9469, 0.9875, 0.91, 0.989, 0.9507, 0.9773, 0.8621, 0.979, // 1 - 10
712 0.8357, 0.868, 0.9417, 0.9466, 0.8911, 0.905, 0.944, 0.9607, 0.928, 0.96, // 11 - 20
713 0.9098, 0.976, 0.8425, 0.8099, 0.7858, 0.947, 0.7248, 0.9106, 0.9246, 0.6821, // 21 - 30
714 0.7223, 0.9784, 0.774, 0.7953, 0.829, 0.9405, 0.8318, 0.8583, 0.8563, 0.8481, // 31 - 40
715 0.8207, 0.9033, 0.8063, 0.7837, 0.7818, 0.744, 0.875, 0.7693, 0.7871, 0.8459, // 41 - 50
716 0.8231, 0.8462, 0.853, 0.8736, 0.856, 0.8762, 0.8629, 0.8323, 0.8064, 0.7828, // 51 - 60
717 0.7533, 0.7273, 0.7093, 0.7157, 0.6823, 0.6612, 0.6418, 0.6395, 0.6323, 0.6221, // 61 - 70
718 0.6497, 0.6746, 0.8568, 0.8541, 0.6958, 0.6962, 0.7051, 0.863, 0.8588, 0.7226, // 71 - 80
719 0.7454, 0.78, 0.7783, 0.7996, 0.8216, 0.8632, 0.8558, 0.8792, 0.8745, 0.8676, // 81 - 90
720 0.8321, 0.8272, 0.7999, 0.7934, 0.7787, 0.7851, 0.7692, 0.7598};