Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4EmCorrections.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 // $Id: G4EmCorrections.cc 110572 2018-05-30 13:08:12Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 //
32 // File name: G4EmCorrections
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 13.01.2005
37 //
38 // Modifications:
39 // 05.05.2005 V.Ivanchenko Fix misprint in Mott term
40 // 26.11.2005 V.Ivanchenko Fix effective charge for heavy ions using original paper
41 // 28.04.2006 V.Ivanchenko General cleanup, add finite size corrections
42 // 13.05.2006 V.Ivanchenko Add corrections for ion stopping
43 // 08.05.2007 V.Ivanchenko Use G4IonTable for ion mass instead of NistTable to avoid
44 // division by zero
45 // 29.02.2008 V.Ivanchenko use expantions for log and power function
46 // 21.04.2008 Updated computations for ions (V.Ivanchenko)
47 // 20.05.2008 Removed Finite Size correction (V.Ivanchenko)
48 // 19.04.2012 Fix reproducibility problem (A.Ribon)
49 //
50 //
51 // Class Description:
52 //
53 // This class provides calculation of EM corrections to ionisation
54 //
55 
56 // -------------------------------------------------------------------
57 //
58 
59 #include "G4EmCorrections.hh"
60 #include "Randomize.hh"
61 #include "G4PhysicalConstants.hh"
62 #include "G4SystemOfUnits.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4IonTable.hh"
65 #include "G4VEmModel.hh"
66 #include "G4Proton.hh"
67 #include "G4GenericIon.hh"
68 #include "G4LPhysicsFreeVector.hh"
69 #include "G4PhysicsLogVector.hh"
70 #include "G4ProductionCutsTable.hh"
71 #include "G4MaterialCutsCouple.hh"
72 #include "G4AtomicShells.hh"
73 #include "G4LPhysicsFreeVector.hh"
74 #include "G4Log.hh"
75 #include "G4Exp.hh"
76 #include "G4Pow.hh"
77 #include "G4Threading.hh"
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
82 
83 const G4double G4EmCorrections::ZD[11] =
84  {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15};
85 const G4double G4EmCorrections::UK[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662,
86  2.0817, 2.0945, 2.0999, 2.1049, 2.1132,
87  2.1197, 2.1246, 2.1280, 2.1292, 2.1301,
88  2.1310, 2.1310, 2.1300, 2.1283, 2.1271};
89 const G4double G4EmCorrections::VK[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247,
90  8.3219, 8.3201, 8.3194, 8.3191, 8.3188,
91  8.3191, 8.3199, 8.3211, 8.3218, 8.3226,
92  8.3244, 8.3264, 8.3285, 8.3308, 8.3320};
94 const G4double G4EmCorrections::Eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02,
95  0.03, 0.04, 0.05, 0.06, 0.08,
96  0.1, 0.15, 0.2, 0.3, 0.4,
97  0.5, 0.6, 0.7, 0.8, 1.0,
98  1.2, 1.4, 1.5, 1.7, 2.0, 3.0, 5.0, 7.0, 10.};
101 const G4double G4EmCorrections::UL[] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828,
102  1.4379, 1.5032, 1.5617, 1.6608, 1.7401,
103  1.8036, 1.8543, 1.8756, 1.8945, 1.9262,
104  1.9508, 1.9696, 1.9836, 1.9890, 1.9935,
105  2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028};
106 G4double G4EmCorrections::VL[] = {0.0};
107 
111 
113 {
114  particle = nullptr;
115  curParticle= nullptr;
116  material = nullptr;
117  curMaterial= nullptr;
118  theElementVector = nullptr;
119  atomDensity= nullptr;
120  curVector = nullptr;
121  ionLEModel = nullptr;
122  ionHEModel = nullptr;
123 
124  kinEnergy = 0.0;
125  verbose = verb;
126  massFactor = 1.0;
127  eth = 2.0*CLHEP::MeV;
128  nbinCorr = 20;
129  eCorrMin = 25.*CLHEP::keV;
130  eCorrMax = 250.*CLHEP::MeV;
131 
134 
136  mass = tau = gamma = bg2 = beta2 = beta = ba2 = tmax = charge = q2 = 0.0;
137 
138  // Constants
140 
141  // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111.
142  // "Shell corrections for K- and L- electrons
143 
144  nK = 20;
145  nL = 26;
146  nEtaK = 29;
147  nEtaL = 28;
148 
149  isMaster = false;
150 
151  // fill vectors
152  if(BarkasCorr == nullptr) { Initialise(); }
153 }
154 
155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156 
158 {
159  for(G4int i=0; i<nIons; ++i) {delete stopData[i];}
160  if(isMaster) {
161  delete BarkasCorr;
162  delete ThetaK;
163  delete ThetaL;
164  BarkasCorr = ThetaK = ThetaL = nullptr;
165  }
166 }
167 
169  const G4Material* mat,
171 {
172  if(kineticEnergy != kinEnergy || p != particle) {
173  particle = p;
175  mass = p->GetPDGMass();
176  tau = kineticEnergy / mass;
177  gamma = 1.0 + tau;
178  bg2 = tau * (tau+2.0);
179  beta2 = bg2/(gamma*gamma);
180  beta = std::sqrt(beta2);
181  ba2 = beta2/alpha2;
184  /(1. + 2.0*gamma*ratio + ratio*ratio);
186  if(charge > 1.5) { charge = effCharge.EffectiveCharge(p,mat,kinEnergy); }
187  q2 = charge*charge;
188  }
189  if(mat != material) {
190  material = mat;
194  }
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 
200  const G4Material* mat,
202 {
203  // . Z^3 Barkas effect in the stopping power of matter for charged particles
204  // J.C Ashley and R.H.Ritchie
205  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
206  // and ICRU49 report
207  // valid for kineticEnergy < 0.5 MeV
208  // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
209 
210  SetupKinematics(p, mat, e);
211  if(tau <= 0.0) { return 0.0; }
212 
213  G4double Barkas = BarkasCorrection (p, mat, e);
214  G4double Bloch = BlochCorrection (p, mat, e);
215  G4double Mott = MottCorrection (p, mat, e);
216 
217  G4double sum = (2.0*(Barkas + Bloch) + Mott);
218 
219  if(verbose > 1) {
220  G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
221  << " Bloch= " << Bloch << " Mott= " << Mott
222  << " Sum= " << sum << " q2= " << q2 << G4endl;
223  G4cout << " ShellCorrection: " << ShellCorrection(p, mat, e)
224  << " Kshell= " << KShellCorrection(p, mat, e)
225  << " Lshell= " << LShellCorrection(p, mat, e)
226  << " " << mat->GetName() << G4endl;
227  }
229  return sum;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
235  const G4Material* mat,
236  G4double e)
237 {
238  return 2.0*BarkasCorrection(p, mat, e)*
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243 
245  const G4Material* mat,
246  G4double e)
247 {
248  // . Z^3 Barkas effect in the stopping power of matter for charged particles
249  // J.C Ashley and R.H.Ritchie
250  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
251  // and ICRU49 report
252  // valid for kineticEnergy < 0.5 MeV
253  // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
254  SetupKinematics(p, mat, e);
255  if(tau <= 0.0) { return 0.0; }
256 
257  G4double Barkas = BarkasCorrection (p, mat, e);
258  G4double Bloch = BlochCorrection (p, mat, e);
259  G4double Mott = MottCorrection (p, mat, e);
260 
261  G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott;
262 
263  if(verbose > 1) {
264  G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
265  << " Bloch= " << Bloch << " Mott= " << Mott
266  << " Sum= " << sum << G4endl;
267  }
269 
270  if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; }
271  return sum;
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275 
277  const G4MaterialCutsCouple* couple,
278  G4double e)
279 {
280  // . Z^3 Barkas effect in the stopping power of matter for charged particles
281  // J.C Ashley and R.H.Ritchie
282  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
283  // and ICRU49 report
284  // valid for kineticEnergy < 0.5 MeV
285  // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
286 
287  G4double sum = 0.0;
288 
289  if(ionHEModel) {
291  if(Z >= 100) Z = 99;
292  else if(Z < 1) Z = 1;
293 
294  G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2;
295  G4int ionPDG = p->GetPDGEncoding();
296  if(thcorr.find(ionPDG)==thcorr.end()) { // Not found: fill the map
297  std::vector<G4double> v;
298  for(size_t i=0; i<ncouples; ++i){
299  v.push_back(ethscaled*ComputeIonCorrections(p,currmat[i],ethscaled));
300  }
301  thcorr.insert(std::pair< G4int, std::vector<G4double> >(ionPDG,v));
302  }
303 
304  //G4cout << " map size=" << thcorr.size() << G4endl;
305  //for(std::map< G4int, std::vector<G4double> >::iterator
306  // it = thcorr.begin(); it != thcorr.end(); ++it){
307  // G4cout << "\t map element: first (key)=" << it->first
308  // << "\t second (vector): vec size=" << (it->second).size() << G4endl;
309  // for(size_t i=0; i<(it->second).size(); ++i){
310  // G4cout << "\t \t vec element: [" << i << "]=" << (it->second)[i]
311  //<< G4endl; } }
312 
313  G4double rest = (thcorr.find(ionPDG)->second)[couple->GetIndex()];
314 
315  sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
316 
317  if(verbose > 1) {
318  G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl;
319  }
320  }
321  return sum;
322 }
323 
324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
325 
327  const G4Material* mat,
328  G4double e)
329 {
330  SetupKinematics(p, mat, e);
332  G4double eexc2 = eexc*eexc;
333  G4double dedx = 0.5*G4Log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2;
334  return dedx;
335 }
336 
337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
338 
340  const G4Material* mat,
341  G4double e)
342 {
343  SetupKinematics(p, mat, e);
344  G4double dedx = 0.5*tmax/(kinEnergy + mass);
345  return 0.5*dedx*dedx;
346 }
347 
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349 
351  const G4Material* mat,
352  G4double e)
353 {
354  SetupKinematics(p, mat, e);
355  G4double term = 0.0;
356  for (G4int i = 0; i<numberOfElements; ++i) {
357 
358  G4double Z = (*theElementVector)[i]->GetZ();
359  G4int iz = (*theElementVector)[i]->GetZasInt();
360  G4double f = 1.0;
361  G4double Z2= (Z-0.3)*(Z-0.3);
362  if(1 == iz) {
363  f = 0.5;
364  Z2 = 1.0;
365  }
366  G4double eta = ba2/Z2;
367  G4double tet = Z2*(1. + Z2*0.25*alpha2);
368  if(11 < iz) { tet = ThetaK->Value(Z); }
369  term += f*atomDensity[i]*KShell(tet,eta)/Z;
370  }
371 
373 
374  return term;
375 }
376 
377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378 
380  const G4Material* mat,
381  G4double e)
382 {
383  SetupKinematics(p, mat, e);
384  G4double term = 0.0;
385  for (G4int i = 0; i<numberOfElements; ++i) {
386 
387  G4double Z = (*theElementVector)[i]->GetZ();
388  G4int iz = (*theElementVector)[i]->GetZasInt();
389  if(2 < iz) {
390  G4double Zeff = Z - ZD[10];
391  if(iz < 10) { Zeff = Z - ZD[iz]; }
392  G4double Z2= Zeff*Zeff;
393  G4double f = 0.125;
394  G4double eta = ba2/Z2;
395  G4double tet = ThetaL->Value(Z);
397  for(G4int j=1; j<nmax; ++j) {
399  if(15 >= iz) {
400  if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
401  else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
402  }
403  //G4cout << " LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
404  // << " ThetaL= " << tet << G4endl;
405  term += f*ne*atomDensity[i]*LShell(tet,eta)/Z;
406  }
407  }
408  }
409 
411 
412  return term;
413 }
414 
415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
416 
418 {
419  G4double corr = 0.0;
420 
421  static const G4double TheK[20] =
422  {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78,
423  0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95};
424 
425 
426  G4double x = tet;
427  G4int itet = 0;
428  G4int ieta = 0;
429  if(tet < TheK[0]) {
430  x = TheK[0];
431  } else if(tet > TheK[nK-1]) {
432  x = TheK[nK-1];
433  itet = nK-2;
434  } else {
435  itet = Index(x, TheK, nK);
436  }
437  // assimptotic case
438  if(eta >= Eta[nEtaK-1]) {
439  corr =
440  (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) +
441  Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta +
442  Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta;
443  } else {
444  G4double y = eta;
445  if(eta < Eta[0]) {
446  y = Eta[0];
447  } else {
448  ieta = Index(y, Eta, nEtaK);
449  }
450  corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
451  CK[itet][ieta], CK[itet+1][ieta],
452  CK[itet][ieta+1], CK[itet+1][ieta+1]);
453  //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]
454  // <<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
455  // <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta]
456  // <<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl;
457  }
458  //G4cout << "Kshell: tet= " << tet << " eta= " << eta << " C= " << corr
459  // << " itet= " << itet << " ieta= " << ieta <<G4endl;
460  return corr;
461 }
462 
463 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
464 
466 {
467  G4double corr = 0.0;
468 
469  static const G4double TheL[26] =
470  {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40,
471  0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.52, 0.54, 0.55, 0.56,
472  0.58, 0.60, 0.62, 0.64, 0.65, 0.66};
473 
474  G4double x = tet;
475  G4int itet = 0;
476  G4int ieta = 0;
477  if(tet < TheL[0]) {
478  x = TheL[0];
479  } else if(tet > TheL[nL-1]) {
480  x = TheL[nL-1];
481  itet = nL-2;
482  } else {
483  itet = Index(x, TheL, nL);
484  }
485 
486  // assimptotic case
487  if(eta >= Eta[nEtaL-1]) {
488  corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1])
489  + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta
490  )/eta;
491  } else {
492  G4double y = eta;
493  if(eta < Eta[0]) {
494  y = Eta[0];
495  } else {
496  ieta = Index(y, Eta, nEtaL);
497  }
498  corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1],
499  CL[itet][ieta], CL[itet+1][ieta],
500  CL[itet][ieta+1], CL[itet+1][ieta+1]);
501  //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet]
502  // <<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
503  // <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta]
504  // <<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl;
505  }
506  //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<eta<<" itet= "<<itet
507  // <<" ieta= "<<ieta<<" Corr= "<<corr<<G4endl;
508  return corr;
509 }
510 
511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
512 
514  const G4Material* mat,
515  G4double e)
516 {
517  SetupKinematics(p, mat, e);
518  G4double taulim= 8.0*MeV/mass;
519  G4double bg2lim= taulim * (taulim+2.0);
520 
521  G4double* shellCorrectionVector =
523  G4double sh = 0.0;
524  G4double x = 1.0;
525  G4double taul = material->GetIonisation()->GetTaul();
526 
527  if ( bg2 >= bg2lim ) {
528  for (G4int k=0; k<3; ++k) {
529  x *= bg2 ;
530  sh += shellCorrectionVector[k]/x;
531  }
532 
533  } else {
534  for (G4int k=0; k<3; ++k) {
535  x *= bg2lim ;
536  sh += shellCorrectionVector[k]/x;
537  }
538  sh *= G4Log(tau/taul)/G4Log(taulim/taul);
539  }
540  sh *= 0.5;
541  return sh;
542 }
543 
544 
545 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
546 
548  const G4Material* mat,
549  G4double ekin)
550 {
551  SetupKinematics(p, mat, ekin);
552 
553  G4double term = 0.0;
554  //G4cout << "### G4EmCorrections::ShellCorrection " << mat->GetName()
555  // << " " << ekin/MeV << " MeV " << G4endl;
556  for (G4int i = 0; i<numberOfElements; ++i) {
557 
558  G4double res = 0.0;
559  G4double res0 = 0.0;
560  G4double Z = (*theElementVector)[i]->GetZ();
561  G4int iz = (*theElementVector)[i]->GetZasInt();
562  G4double Z2= (Z-0.3)*(Z-0.3);
563  G4double f = 1.0;
564  if(1 == iz) {
565  f = 0.5;
566  Z2 = 1.0;
567  }
568  G4double eta = ba2/Z2;
569  G4double tet = Z2*(1. + Z2*0.25*alpha2);
570  if(11 < iz) { tet = ThetaK->Value(Z); }
571  res0 = f*KShell(tet,eta);
572  res += res0;
573  //G4cout << " Z= " << iz << " Shell 0" << " tet= " << tet
574  // << " eta= " << eta << " resK= " << res0 << G4endl;
575  if(2 < iz) {
576  G4double Zeff = Z - ZD[10];
577  if(iz < 10) { Zeff = Z - ZD[iz]; }
578  Z2= Zeff*Zeff;
579  eta = ba2/Z2;
580  f = 0.125;
581  tet = ThetaL->Value(Z);
583  G4int nmax = std::min(4, ntot);
584  G4double norm = 0.0;
585  G4double eshell = 0.0;
586  for(G4int j=1; j<nmax; ++j) {
588  if(15 >= iz) {
589  if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
590  else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
591  }
592  norm += ne;
593  eshell += tet*ne;
594  res0 = f*ne*LShell(tet,eta);
595  res += res0;
596  //G4cout << " Z= " << iz << " Shell " << j << " Ne= " << ne
597  // << " tet= " << tet << " eta= " << eta
598  // << " resL= " << res0 << G4endl;
599  }
600  if(ntot > nmax) {
601  eshell /= norm;
602 
603  static const G4double HM[53] = {
604  12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4,
605  10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87,
606  5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38,
607  4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 3.95, 3.93, 3.91,
608  3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89,
609  3.90, 3.92, 3.93 };
610  static const G4double HN[31] = {
611  75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8,
612  24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7,
613  19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4,
614  18.2};
615 
616  // Add M-shell
617  if(28 > iz) {
618  res += f*(iz - 10)*LShell(eshell,HM[iz-11]*eta);
619  } else if(63 > iz) {
620  res += f*18*LShell(eshell,HM[iz-11]*eta);
621  } else {
622  res += f*18*LShell(eshell,HM[52]*eta);
623  }
624  // Add N-shell
625  if(32 < iz) {
626  if(60 > iz) {
627  res += f*(iz - 28)*LShell(eshell,HN[iz-33]*eta);
628  } else if(63 > iz) {
629  res += 4*LShell(eshell,HN[iz-33]*eta);
630  } else {
631  res += 4*LShell(eshell,HN[30]*eta);
632  }
633  // Add O-P-shells
634  if(60 < iz) {
635  res += f*(iz - 60)*LShell(eshell,150*eta);
636  }
637  }
638  }
639  }
640  term += res*atomDensity[i]/Z;
641  }
642 
644  //G4cout << "# Shell Correction= " << term << G4endl;
645  return term;
646 }
647 
648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
649 
651  const G4Material* mat,
652  G4double e)
653 {
654  SetupKinematics(p, mat, e);
655 
661 
662  G4double dedx = 0.0;
663 
664  // density correction
665  static const G4double twoln10 = 2.0*G4Log(10.0);
666  G4double x = G4Log(bg2)/twoln10;
667  if ( x >= x0den ) {
668  dedx = twoln10*x - cden ;
669  if ( x < x1den ) dedx += aden*G4Exp(G4Log(x1den-x)*mden) ;
670  }
671 
672  return dedx;
673 }
674 
675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
676 
678  const G4Material* mat,
679  G4double e)
680 {
681  // . Z^3 Barkas effect in the stopping power of matter for charged particles
682  // J.C Ashley and R.H.Ritchie
683  // Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397
684  // valid for kineticEnergy > 0.5 MeV
685 
686  SetupKinematics(p, mat, e);
687  G4double BarkasTerm = 0.0;
688 
689  for (G4int i = 0; i<numberOfElements; ++i) {
690 
691  G4double Z = (*theElementVector)[i]->GetZ();
692  G4int iz = (*theElementVector)[i]->GetZasInt();
693  if(iz == 47) {
694  BarkasTerm += atomDensity[i]*0.006812*G4Exp(-G4Log(beta)*0.9);
695  } else if(iz >= 64) {
696  BarkasTerm += atomDensity[i]*0.002833*G4Exp(-G4Log(beta)*1.2);
697  } else {
698 
699  G4double X = ba2 / Z;
700  G4double b = 1.3;
701  if(1 == iz) {
702  if(material->GetName() == "G4_lH2") { b = 0.6; }
703  else { b = 1.8; }
704  }
705  else if(2 == iz) { b = 0.6; }
706  else if(10 >= iz) { b = 1.8; }
707  else if(17 >= iz) { b = 1.4; }
708  else if(18 == iz) { b = 1.8; }
709  else if(25 >= iz) { b = 1.4; }
710  else if(50 >= iz) { b = 1.35;}
711 
712  G4double W = b/std::sqrt(X);
713 
714  G4double val = BarkasCorr->Value(W);
715  if(W > BarkasCorr->Energy(46)) {
716  val *= BarkasCorr->Energy(46)/W;
717  }
718  // G4cout << "i= " << i << " b= " << b << " W= " << W
719  // << " Z= " << Z << " X= " << X << " val= " << val<< G4endl;
720  BarkasTerm += val*atomDensity[i] / (std::sqrt(Z*X)*X);
721  }
722  }
723 
724  BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
725 
726  return BarkasTerm;
727 }
728 
729 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
730 
732  const G4Material* mat,
733  G4double e)
734 {
735  SetupKinematics(p, mat, e);
736 
737  G4double y2 = q2/ba2;
738 
739  G4double term = 1.0/(1.0 + y2);
740  G4double del;
741  G4double j = 1.0;
742  do {
743  j += 1.0;
744  del = 1.0/(j* (j*j + y2));
745  term += del;
746  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
747  } while (del > 0.01*term);
748 
749  G4double res = -y2*term;
750  return res;
751 }
752 
753 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
754 
756  const G4Material* mat,
757  G4double e)
758 {
759  SetupKinematics(p, mat, e);
761  return mterm;
762 }
763 
764 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
765 
766 G4double
768  const G4Material* mat,
769  G4double ekin)
770 {
771  G4double factor = 1.0;
772  if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || nIons <= 0) { return factor; }
773 
774  if(verbose > 1) {
775  G4cout << "EffectiveChargeCorrection: " << p->GetParticleName()
776  << " in " << mat->GetName()
777  << " ekin(MeV)= " << ekin/MeV << G4endl;
778  }
779 
780  if(p != curParticle || mat != curMaterial) {
781  curParticle = p;
782  curMaterial = mat;
783  curVector = nullptr;
784  currentZ = p->GetAtomicNumber();
785  if(verbose > 1) {
786  G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= "
787  << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl;
788  }
790  idx = -1;
791 
792  for(G4int i=0; i<nIons; ++i) {
793  if(materialList[i] == mat && currentZ == Zion[i]) {
794  idx = i;
795  break;
796  }
797  }
798  //G4cout << " idx= " << idx << " dz= " << G4endl;
799  if(idx >= 0) {
800  if(!ionList[idx]) { BuildCorrectionVector(); }
801  if(ionList[idx]) { curVector = stopData[idx]; }
802  } else { return factor; }
803  }
804  if(curVector) {
805  factor = curVector->Value(ekin*massFactor);
806  if(verbose > 1) {
807  G4cout << "E= " << ekin << " factor= " << factor << " massfactor= "
808  << massFactor << G4endl;
809  }
810  }
811  return factor;
812 }
813 
814 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
815 
817  const G4String& mname,
818  G4PhysicsVector* dVector)
819 {
820  G4int i = 0;
821  for(; i<nIons; ++i) {
822  if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break;
823  }
824  if(i == nIons) {
825  Zion.push_back(Z);
826  Aion.push_back(A);
827  materialName.push_back(mname);
828  materialList.push_back(nullptr);
829  ionList.push_back(nullptr);
830  stopData.push_back(dVector);
831  nIons++;
832  if(verbose > 1) {
833  G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname
834  << " idx= " << i << G4endl;
835  }
836  }
837 }
838 
839 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840 
842 {
843  if(!ionLEModel || !ionHEModel) {
844  return;
845  }
846 
848  G4int Z = Zion[idx];
849  if(currentZ != Z) {
850  ion = ionTable->GetIon(Z, Aion[idx], 0);
851  }
852  //G4cout << "BuildCorrectionVector: idx= " << idx << " Z= " << Z
853  // << " curZ= " << currentZ << G4endl;
854 
857 
859  G4double massRatio = proton_mass_c2/ion->GetPDGMass();
860 
861  if(verbose > 1) {
862  G4cout << "BuildCorrectionVector: Stopping for "
863  << curParticle->GetParticleName() << " in "
864  << materialName[idx] << " Ion Z= " << Z << " A= " << A
865  << " massRatio= " << massRatio << G4endl;
866  }
867 
868  G4PhysicsLogVector* vv =
870  vv->SetSpline(true);
871  G4double e, eion, dedx, dedx1;
872  G4double eth0 = v->Energy(0);
873  G4double escal = eth/massRatio;
874  G4double qe =
876  G4double dedxt =
878  G4double dedx1t =
881  G4double rest = escal*(dedxt - dedx1t);
882  //G4cout << "Escal(MeV)= "<<escal<<" dedxt0= " <<dedxt
883  // << " dedxt1= " << dedx1t << G4endl;
884 
885  for(G4int i=0; i<=nbinCorr; ++i) {
886  e = vv->Energy(i);
887  escal = e/massRatio;
888  eion = escal/A;
889  if(eion <= eth0) {
890  dedx = v->Value(eth0)*std::sqrt(eion/eth0);
891  } else {
892  dedx = v->Value(eion);
893  }
895  if(e <= eth) {
896  dedx1 = ionLEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe;
897  } else {
898  dedx1 = ionHEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe +
899  ComputeIonCorrections(curParticle, curMaterial, escal) + rest/escal;
900  }
901  vv->PutValue(i, dedx/dedx1);
902  if(verbose > 1) {
903  G4cout << " E(meV)= " << e/MeV << " Correction= " << dedx/dedx1
904  << " " << dedx << " " << dedx1
905  << " massF= " << massFactor << G4endl;
906  }
907  }
908  delete v;
909  ionList[idx] = ion;
910  stopData[idx] = vv;
911  if(verbose > 1) { G4cout << "End data set " << G4endl; }
912 }
913 
914 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
915 
917 {
919  ncouples = tb->GetTableSize();
920  if(currmat.size() != ncouples) {
921  currmat.resize(ncouples);
922  for(std::map< G4int, std::vector<G4double> >::iterator it =
923  thcorr.begin(); it != thcorr.end(); ++it){
924  (it->second).clear();
925  }
926  thcorr.clear();
927  for(size_t i=0; i<ncouples; ++i) {
929  G4String nam = currmat[i]->GetName();
930  for(G4int j=0; j<nIons; ++j) {
931  if(nam == materialName[j]) { materialList[j] = currmat[i]; }
932  }
933  }
934  }
935 }
936 
937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
938 
940 {
941  if(G4Threading::IsMasterThread()) { isMaster = true; }
942 
943  // Z^3 Barkas effect in the stopping power of matter for charged particles
944  // J.C Ashley and R.H.Ritchie
945  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
946  G4int i, j;
947  static const G4double fTable[47][2] = {
948  { 0.02, 21.5},
949  { 0.03, 20.0},
950  { 0.04, 18.0},
951  { 0.05, 15.6},
952  { 0.06, 15.0},
953  { 0.07, 14.0},
954  { 0.08, 13.5},
955  { 0.09, 13.},
956  { 0.1, 12.2},
957  { 0.2, 9.25},
958  { 0.3, 7.0},
959  { 0.4, 6.0},
960  { 0.5, 4.5},
961  { 0.6, 3.5},
962  { 0.7, 3.0},
963  { 0.8, 2.5},
964  { 0.9, 2.0},
965  { 1.0, 1.7},
966  { 1.2, 1.2},
967  { 1.3, 1.0},
968  { 1.4, 0.86},
969  { 1.5, 0.7},
970  { 1.6, 0.61},
971  { 1.7, 0.52},
972  { 1.8, 0.5},
973  { 1.9, 0.43},
974  { 2.0, 0.42},
975  { 2.1, 0.3},
976  { 2.4, 0.2},
977  { 3.0, 0.13},
978  { 3.08, 0.1},
979  { 3.1, 0.09},
980  { 3.3, 0.08},
981  { 3.5, 0.07},
982  { 3.8, 0.06},
983  { 4.0, 0.051},
984  { 4.1, 0.04},
985  { 4.8, 0.03},
986  { 5.0, 0.024},
987  { 5.1, 0.02},
988  { 6.0, 0.013},
989  { 6.5, 0.01},
990  { 7.0, 0.009},
991  { 7.1, 0.008},
992  { 8.0, 0.006},
993  { 9.0, 0.0032},
994  { 10.0, 0.0025} };
995 
996  BarkasCorr = new G4LPhysicsFreeVector(47, 0.02, 10.);
997  for(i=0; i<47; ++i) { BarkasCorr->PutValues(i, fTable[i][0], fTable[i][1]); }
998  BarkasCorr->SetSpline(true);
999 
1000  static const G4double SK[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137,
1001  1.7754, 1.7396, 1.7223, 1.7063, 1.6752,
1002  1.6461, 1.6189, 1.5933, 1.5811, 1.5693,
1003  1.5467, 1.5254, 1.5053, 1.4863, 1.4772};
1004  static const G4double TK[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608,
1005  2.4388, 2.4163, 2.4044, 2.3933, 2.3701,
1006  2.3466, 2.3229, 2.2992, 2.2872, 2.2753,
1007  2.2515, 2.2277, 2.2040, 2.1804, 2.1686};
1008 
1009  static const G4double SL[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283,
1010  10.3424, 10.0371, 9.7537, 9.2443, 8.8005,
1011  8.4114, 8.0683, 7.9117, 7.7641, 7.4931,
1012  7.2506, 7.0327, 6.8362, 6.7452, 6.6584,
1013  6.4969, 6.3498, 6.2154, 6.0923, 6.0345, 5.9792};
1014  static const G4double TL[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226,
1015  28.6128, 28.1449, 27.6991, 26.8674, 26.1061,
1016  25.4058, 24.7587, 24.4531, 24.1583, 23.5992,
1017  23.0771, 22.5880, 22.1285, 21.9090, 21.6958,
1018  21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546};
1019 
1020  static const G4double bk1[29][11] = {
1021  {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, 2.03521E-8, 2.41370E-8, 2.87247E-8, 3.13778E-8, 3.43072E-8, 4.11274E-8, 4.94946E-8},
1022  {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8, 1.03022E-7, 1.21760E-7, 1.44370E-7, 1.57398E-7, 1.71747E-7, 2.05023E-7, 2.45620E-7},
1023  {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5.60760E-7, 6.59478E-7, 7.77847E-7, 8.45709E-7, 9.20187E-7, 1.09192E-6, 1.29981E-6},
1024  {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6, 3.68791E-6, 4.30423E-6, 5.03578E-6, 5.45200E-6, 5.90633E-6, 6.94501E-6, 8.18757E-6},
1025  {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1.35313E-5, 1.56829E-5, 1.82138E-5, 1.96439E-5, 2.11973E-5, 2.47216E-5, 2.88935E-5},
1026  {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7.90324E-5, 9.04832E-5, 1.03744E-4, 1.11147E-4, 1.19122E-4, 1.36980E-4, 1.57744E-4},
1027  {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2.60584E-4, 2.95248E-4, 3.34870E-4, 3.56771E-4, 3.80200E-4, 4.32104E-4, 4.91584E-4},
1028  {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6.31597E-4, 7.09244E-4, 7.97023E-4, 8.45134E-4, 8.96410E-4, 1.00867E-3, 1.13590E-3},
1029  {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1.26476E-3, 1.46928E-3, 1.57113E-3, 1.65921E-3, 1.75244E-3, 1.95562E-3, 2.18336E-3},
1030  {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3.57027E-3, 3.92793E-3, 4.32246E-3, 4.53473E-3, 4.75768E-3, 5.23785E-3, 5.76765E-3},
1031  {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7.59781E-3, 8.27424E-3, 9.01187E-3, 9.40534E-3, 9.81623E-3, 1.06934E-2, 1.16498E-2},
1032  {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2.64300E-2, 2.82832E-2, 3.02661E-2, 3.13090E-2, 3.23878E-2, 3.46580E-2, 3.70873E-2},
1033  {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5.80617E-2, 6.14403E-2, 6.50152E-2, 6.68798E-2, 6.87981E-2, 7.28020E-2, 7.70407E-2},
1034  {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1.53743E-1, 1.60536E-1, 1.67641E-1, 1.71315E-1, 1.75074E-1, 1.82852E-1, 1.90997E-1},
1035  {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2.79776E-1, 2.89946E-1, 3.00525E-1, 3.05974E-1, 3.11533E-1, 3.22994E-1, 3.34935E-1},
1036  {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4.23427E-1, 4.36726E-1, 4.50519E-1, 4.57610E-1, 4.64834E-1, 4.79700E-1, 4.95148E-1},
1037  {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5.75948E-1, 5.92092E-1, 6.08811E-1, 6.17396E-1, 6.26136E-1, 6.44104E-1, 6.62752E-1},
1038  {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7.31650E-1, 7.50373E-1, 7.69748E-1, 7.79591E-1, 7.89811E-1, 8.10602E-1, 8.32167E-1},
1039  {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8.86910E-1, 9.07979E-1, 9.29772E-1, 9.40953E-1, 9.52331E-1, 9.75701E-1, 9.99934E-1},
1040  {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.21317, 1.23921, 1.25257, 1.26616, 1.29408, 1.32303},
1041  {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.49875, 1.52842, 1.54364, 1.55913, 1.59095, 1.62396},
1042  {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.76244, 1.79526, 1.81210, 1.82925, 1.86448, 1.90104},
1043  {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.88617, 1.92042, 1.93800, 1.95590, 1.99269, 2.03087},
1044  {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.11827, 2.15555, 2.17409, 2.19337, 2.23302, 2.27419},
1045  {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.42395, 2.46429, 2.48401, 2.50612, 2.54955, 2.59468},
1046  {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.24748, 3.29664, 3.32192, 3.34770, 3.40081, 3.45611},
1047  {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.27010, 4.33056, 4.36172, 4.39353, 4.45918, 4.52772},
1048  {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.91402, 4.98200, 5.01707, 5.05290, 5.12695, 5.20436},
1049  {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359}
1050  };
1051 
1052  static const G4double bk2[29][11] = {
1053  {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, 8.84294E-8, 1.08253E-7, 1.33148E-7, 1.64573E-7, 2.04459E-7, 2.28346E-7, 2.55370E-7},
1054  {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7, 4.32017E-7, 5.25688E-7, 6.42391E-7, 7.88464E-7, 9.72171E-7, 1.08140E-6, 1.20435E-6},
1055  {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6, 2.23662E-6, 2.69889E-6, 3.26860E-6, 3.26860E-6, 4.84882E-6, 5.36428E-6, 5.94048E-6},
1056  {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5, 1.36329E-5, 1.62480E-5, 1.94200E-5, 2.32783E-5, 2.79850E-5, 3.07181E-5, 3.37432E-5},
1057  {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5, 4.67351E-5, 5.51033E-5, 6.51154E-5, 7.71154E-5, 9.15431E-5, 9.98212E-5, 1.08909E-4},
1058  {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4, 2.43007E-4, 2.81460E-4, 3.26458E-4, 3.79175E-4, 4.41006E-4, 4.75845E-4, 5.13606E-4},
1059  {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4, 7.28042E-4, 8.31425E-4, 9.50341E-4, 1.08721E-3, 1.24485E-3, 1.33245E-3, 1.42650E-3},
1060  {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3, 1.62855E-3, 1.83861E-3, 2.07396E-3, 2.34750E-3, 2.65469E-3, 2.82358E-3, 3.00358E-3},
1061  {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3, 3.04636E-3, 3.40681E-3, 3.81132E-3, 4.26536E-3, 4.77507E-3, 5.05294E-3, 5.34739E-3},
1062  {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3, 7.70916E-3, 8.49478E-3, 9.36187E-3, 1.03189E-2, 1.13754E-2, 1.19441E-2, 1.25417E-2},
1063  {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1.50707E-2, 1.64235E-2, 1.78989E-2, 1.95083E-2, 2.12640E-2, 2.22009E-2, 2.31795E-2},
1064  {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2, 4.54488E-2, 4.86383E-2, 5.20542E-2, 5.57135E-2, 5.96350E-2, 6.17003E-2, 6.38389E-2},
1065  {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9.13200E-2, 9.66589E-2, 1.02320E-1, 1.08326E-1, 1.14701E-1, 1.18035E-1, 1.21472E-1},
1066  {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2.17848E-1, 2.27689E-1, 2.38022E-1, 2.48882E-1, 2.60304E-1, 2.66239E-1, 2.72329E-1},
1067  {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3.73925E-1, 3.88089E-1, 4.02900E-1, 4.18404E-1, 4.34647E-1, 4.43063E-1, 4.51685E-1},
1068  {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5.45354E-1, 5.63515E-1, 5.82470E-1, 6.02273E-1, 6.22986E-1, 6.33705E-1, 6.44677E-1},
1069  {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7.23214E-1, 7.45041E-1, 7.67800E-1, 7.91559E-1, 8.16391E-1, 8.29235E-1, 8.42380E-1},
1070  {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9.02008E-1, 9.27198E-1, 9.53454E-1, 9.80856E-1, 1.00949, 1.02430, 1.03945},
1071  {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1.10667, 1.13615, 1.16692, 1.19907, 1.21570, 1.23272},
1072  {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1.45057, 1.48582, 1.52263, 1.56111, 1.58102, 1.60140},
1073  {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1.76954, 1.80983, 1.85192, 1.89596, 1.91876, 1.94211},
1074  {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2.06251, 2.10727, 2.15406, 2.20304, 2.22842, 2.25442},
1075  {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2.19963, 2.24644, 2.29539, 2.34666, 2.37323, 2.40045},
1076  {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2.45641, 2.50703, 2.56000, 2.61552, 2.64430, 2.67381},
1077  {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2.79481, 2.85052, 2.90887, 2.97009, 3.00185, 3.03442},
1078  {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3.70268, 3.77170, 3.84418, 3.92040, 3.96003, 4.00073},
1079  {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4.83543, 4.92217, 5.01353, 5.10992, 5.16014, 5.21181},
1080  {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5.55344, 5.65226, 5.79496, 5.90517, 5.96269, 6.02191},
1081  {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535}
1082  };
1083 
1084  static const G4double bls1[28][10] = {
1085  {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3.0658E-4, 3.4511E-4, 3.8795E-4, 4.3542E-4, 4.6100E-4, 4.8786E-4},
1086  {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7.7167E-4, 8.4592E-4, 9.2605E-4, 1.0125E-3, 1.0583E-3, 1.1058E-3},
1087  {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7856E-3, 1.9181E-3, 2.1615E-3, 2.3178E-3, 2.4019E-3, 2.4904E-3},
1088  {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4.2354E-3, 4.5396E-3, 4.8853E-3, 5.2820E-3, 5.5031E-3, 5.7414E-3},
1089  {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0336E-3, 8.6984E-3, 9.4638E-3, 1.0348E-2, 1.0841E-2, 1.1372E-2},
1090  {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0784E-2, 2.2797E-2, 2.5066E-2, 2.7622E-2, 2.9020E-2, 3.0503E-2},
1091  {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5300E-2, 4.9254E-2, 5.3619E-2, 5.8436E-2, 6.1028E-2, 6.3752E-2},
1092  {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8282E-2, 8.4470E-2, 9.1206E-2, 9.8538E-2, 1.0244E-1, 1.0652E-1},
1093  {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1943E-1, 1.2796E-1, 1.3715E-1, 1.4706E-1, 1.5231E-1, 1.5776E-1},
1094  {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2038E-1, 2.3351E-1, 2.4751E-1, 2.6244E-1, 2.7027E-1, 2.7837E-1},
1095  {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.3796E-1, 3.5537E-1, 3.7381E-1, 3.9336E-1, 4.0357E-1, 4.1410E-1},
1096  {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6330E-1, 6.8974E-1, 7.1753E-1, 7.4678E-1, 7.6199E-1, 7.7761E-1},
1097  {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.0037, 1.0372, 1.0723, 1.1092, 1.1284, 1.1480},
1098  {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.6438, 1.6899, 1.7382, 1.7632, 1.7889},
1099  {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.1659, 2.2199, 2.2765, 2.3059, 2.3360},
1100  {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.5587, 2.6190, 2.6821, 2.7148, 2.7484},
1101  {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.8884, 2.9538, 3.0222, 3.0577, 3.0941},
1102  {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.1652, 3.2349, 3.3079, 3.3457, 3.3845},
1103  {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.4232, 3.4965, 3.5734, 3.6133, 3.6542},
1104  {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.8866, 3.9661, 4.0495, 4.0928, 4.1371},
1105  {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.2130, 4.2974, 4.3861, 4.4321, 4.4794},
1106  {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.5170, 4.6056, 4.6988, 4.7471, 4.7968},
1107  {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.6521, 4.7426, 4.8378, 4.8872, 4.9379},
1108  {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.8949, 4.9888, 5.0876, 5.1388, 5.1915},
1109  {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.2053, 5.3036, 5.4070, 5.4607, 5.5159},
1110  {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.9489, 6.0581, 6.1730, 6.2328, 6.2943},
1111  {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.8003, 6.9231, 7.0525, 7.1199, 7.1892},
1112  {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384}
1113  };
1114 
1115  static const G4double bls2[28][10] = {
1116  {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7.5494E-4, 7.9587E-4, 8.3883E-4, 9.3160E-4, 1.0352E-3, 1.1529E-3},
1117  {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1.5719E-3, 1.6451E-3, 1.7231E-3, 1.8969E-3, 2.1009E-3, 2.3459E-3},
1118  {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4479E-3, 3.6149E-3, 3.7976E-3, 4.2187E-3, 4.7320E-3, 5.3636E-3},
1119  {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8.5370E-2, 9.0407E-3, 9.5910E-3, 1.0850E-2, 1.2358E-2, 1.4165E-2},
1120  {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7327E-2, 1.8478E-2, 1.9612E-2, 2.2160E-2, 2.5130E-2, 2.8594E-2},
1121  {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6170E-2, 4.8708E-2, 5.1401E-2, 5.7297E-2, 6.3943E-2, 7.1441E-2},
1122  {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1150E-2, 9.5406E-2, 9.9881E-2, 1.0954E-1, 1.2023E-1, 1.3208E-1},
1123  {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4632E-1, 1.5234E-1, 1.5864E-1, 1.7211E-1, 1.8686E-1, 2.0304E-1},
1124  {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0991E-1, 2.1767E-1, 2.2576E-1, 2.4295E-1, 2.6165E-1, 2.8201E-1},
1125  {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5403E-1, 3.6506E-1, 3.7650E-1, 4.0067E-1, 4.2673E-1, 4.5488E-1},
1126  {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1115E-1, 5.2514E-1, 5.3961E-1, 5.7008E-1, 6.0277E-1, 6.3793E-1},
1127  {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1954E-1, 9.3973E-1, 9.6056E-1, 1.0043, 1.0509, 1.1008},
1128  {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.3506, 1.3765, 1.4308, 1.4886, 1.5504},
1129  {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.0528, 2.0864, 2.1569, 2.2319, 2.3120},
1130  {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.6452, 2.6847, 2.7674, 2.8554, 2.9494},
1131  {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.0937, 3.1378, 3.2301, 3.3285, 3.4337},
1132  {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.4691, 3.5171, 3.6175, 3.7246, 3.8391},
1133  {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.7848, 3.8360, 3.9433, 4.0578, 4.1804},
1134  {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.0764, 4.1304, 4.2438, 4.3648, 4.4944},
1135  {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.5962, 4.6551, 4.7786, 4.9106, 5.0520},
1136  {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.9685, 5.0314, 5.1633, 5.3043, 5.4555},
1137  {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.3114, 5.3776, 5.5166, 5.6653, 5.8249},
1138  {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.4640, 5.5317, 5.6739, 5.8260, 5.9893},
1139  {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.7373, 5.8088, 5.9568, 6.1152, 6.2853},
1140  {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.0896, 6.1636, 6.3192, 6.4857, 6.6647},
1141  {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.9351, 7.0180, 7.1925, 7.3795, 7.5808},
1142  {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.9146, 8.0087, 8.2071, 8.4200, 8.6496},
1143  {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133}
1144  };
1145 
1146  static const G4double bls3[28][9] = {
1147  {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1.6524E-3, 1.9078E-3, 2.2414E-3, 2.6889E-3, 3.3006E-3},
1148  {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3.5045E-3, 4.1260E-3, 4.9376E-3, 6.0050E-3, 7.4152E-3},
1149  {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3491E-3, 9.8871E-3, 1.1822E-2, 1.4261E-2, 1.7335E-2},
1150  {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2.2053E-2, 2.5803E-2, 3.0308E-2, 3.5728E-2, 4.2258E-2},
1151  {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2850E-2, 4.9278E-2, 5.6798E-2, 6.5610E-2, 7.5963E-2},
1152  {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0032E-1, 1.1260E-1, 1.2656E-1, 1.4248E-1, 1.6071E-1},
1153  {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7614E-1, 1.9434E-1, 2.1473E-1, 2.3766E-1, 2.6357E-1},
1154  {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6199E-1, 2.8590E-1, 3.1248E-1, 3.4212E-1, 3.7536E-1},
1155  {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5522E-1, 3.8459E-1, 4.1704E-1, 4.5306E-1, 4.9326E-1},
1156  {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5453E-1, 5.9397E-1, 6.3728E-1, 6.8507E-1, 7.3810E-1},
1157  {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.6141E-1, 8.0992E-1, 8.6301E-1, 9.2142E-1, 9.8604E-1},
1158  {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.3424, 1.4163, 1.4974, 1.5868},
1159  {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.8484, 1.9394, 2.0390, 2.1489},
1160  {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.6980, 2.8159, 2.9451, 3.0876},
1161  {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.4030, 3.5416, 3.6938, 3.8620},
1162  {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.9421, 4.0978, 4.2688, 4.4580},
1163  {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.3935, 4.5636, 4.7506, 4.9576},
1164  {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.7644, 4.9470, 5.1478, 5.3703},
1165  {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.1233, 5.3169, 5.5300, 5.7661},
1166  {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.7400, 5.9523, 6.1863, 6.4458},
1167  {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.1927, 6.4206, 6.6719, 6.9510},
1168  {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.6038, 6.8451, 7.1113, 7.4071},
1169  {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.7871, 7.0344, 7.3073, 7.6107},
1170  {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.1175, 7.3757, 7.6609, 7.9782},
1171  {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.5418, 7.8143, 8.1156, 8.4510},
1172  {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.5708, 8.8796, 9.2214, 9.6027},
1173  {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.7834, 10.1384, 10.5323, 10.9722},
1174  {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356}
1175  };
1176 
1177  static const G4double bll1[28][10] = {
1178  {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5.6969E-5, 7.1625E-5, 9.0279E-5, 1.1407E-4, 1.2834E-4, 1.4447E-4},
1179  {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2.7006E-4, 3.3049E-4, 4.0498E-4, 4.9688E-4, 5.5061E-4, 6.1032E-4},
1180  {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2178E-3, 1.4459E-3, 1.7174E-3, 2.0405E-3, 2.2245E-3, 2.4252E-3},
1181  {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5.5731E-3, 6.3968E-3, 7.3414E-3, 8.4242E-3, 9.0236E-3, 9.6652E-3},
1182  {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4486E-2, 1.6264E-2, 1.8256E-2, 2.0487E-2, 2.1702E-2, 2.2989E-2},
1183  {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7163E-2, 5.1543E-2, 5.6423E-2, 6.1540E-2, 6.4326E-2, 6.7237E-2},
1184  {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7747E-2, 1.0516E-1, 1.1313E-1, 1.2171E-1, 1.2625E-1, 1.3096E-1},
1185  {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6253E-1, 1.7306E-1, 1.8430E-1, 1.9630E-1, 2.0261E-1, 2.0924E-1},
1186  {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3794E-1, 2.5153E-1, 2.6596E-1, 2.8130E-1, 2.8934E-1, 2.9763E-1},
1187  {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0955E-1, 4.2888E-1, 4.4928E-1, 4.7086E-1, 4.8212E-1, 4.9371E-1},
1188  {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.9601E-1, 6.2049E-1, 6.4627E-1, 6.7346E-1, 6.8762E-1, 7.0218E-1},
1189  {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.1073, 1.1443, 1.1832, 1.2035, 1.2243},
1190  {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.5624, 1.6083, 1.6566, 1.6817, 1.7076},
1191  {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.3322, 2.3915, 2.4542, 2.4869, 2.5205},
1192  {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.9377, 3.0072, 3.0807, 3.1192, 3.1587},
1193  {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.4122, 3.4898, 3.5721, 3.6151, 3.6595},
1194  {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.8303, 3.9146, 4.0042, 4.0511, 4.0995},
1195  {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.1520, 4.2421, 4.3380, 4.3882, 4.4401},
1196  {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.4333, 4.5285, 4.6298, 4.6830, 4.7380},
1197  {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.9106, 5.0144, 5.1250, 5.1831, 5.2432},
1198  {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.2729, 5.3837, 5.5050, 5.5642, 5.6287},
1199  {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.5782, 5.6950, 5.8198, 5.8855, 5.9554},
1200  {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.7128, 5.8323, 5.9601, 6.0274, 6.0972},
1201  {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.9541, 6.0785, 6.2116, 6.2818, 6.3546},
1202  {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.2618, 6.3925, 6.5327, 6.6066, 6.6833},
1203  {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.0062, 7.1529, 7.3107, 7.3941, 7.4807},
1204  {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.8982, 8.0653, 8.2454, 8.3409, 8.4402},
1205  {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927}
1206  };
1207 
1208  static const G4double bll2[28][10] = {
1209  {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3.7977E-4, 4.2945E-4, 4.8582E-4, 6.2244E-4, 7.9858E-4, 1.0258E-3},
1210  {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1.4021E-3, 1.5570E-3, 1.7292E-3, 2.1335E-3, 2.6335E-3, 3.2515E-3},
1211  {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8457E-3, 5.2839E-3, 5.7617E-3, 6.8504E-3, 8.1442E-3, 9.6816E-3},
1212  {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1.6717E-2, 1.7898E-2, 1.9163E-2, 2.1964E-2, 2.5173E-2, 2.8851E-2},
1213  {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6371E-2, 3.8514E-2, 4.0784E-2, 4.5733E-2, 5.1288E-2, 5.7531E-2},
1214  {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5852E-2, 1.0021E-1, 1.0478E-1, 1.1458E-1, 1.2535E-1, 1.3721E-1},
1215  {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7596E-1, 1.8265E-1, 1.8962E-1, 2.0445E-1, 2.2058E-1, 2.3818E-1},
1216  {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7045E-1, 2.7944E-1, 2.8877E-1, 3.0855E-1, 3.2995E-1, 3.5318E-1},
1217  {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7473E-1, 3.8594E-1, 3.9756E-1, 4.2212E-1, 4.4861E-1, 4.7727E-1},
1218  {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0032E-1, 6.1569E-1, 6.3159E-1, 6.6512E-1, 7.0119E-1, 7.4012E-1},
1219  {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.3556E-1, 8.5472E-1, 8.7454E-1, 9.1630E-1, 9.6119E-1, 1.0096},
1220  {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.4421, 1.4703, 1.5300, 1.5942, 1.6636},
1221  {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.9788, 2.0142, 2.0889, 2.1694, 2.2567},
1222  {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.8762, 2.9231, 3.0222, 3.1295, 3.2463},
1223  {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.5803, 3.6361, 3.7546, 3.8835, 4.0242},
1224  {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.1350, 4.1983, 4.3330, 4.4799, 4.6408},
1225  {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.6206, 4.6904, 4.8390, 5.0013, 5.1796},
1226  {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.0009, 5.0762, 5.2370, 5.4129, 5.6066},
1227  {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.3338, 5.4141, 5.5857, 5.7738, 5.9811},
1228  {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.8989, 5.9878, 6.1780, 6.3870, 6.6179},
1229  {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.3343, 6.4303, 6.6361, 6.8626, 7.1137},
1230  {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.7021, 6.8043, 7.0237, 7.2655, 7.5338},
1231  {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.8650, 6.9700, 7.1954, 7.4442, 7.7203},
1232  {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.1578, 7.2679, 7.5045, 7.7660, 8.0565},
1233  {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.5328, 7.6496, 7.9010, 8.1791, 8.4886},
1234  {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.4475, 8.5814, 8.8702, 9.1908, 9.5488},
1235  {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.5572, 9.7132, 10.0504, 10.4259, 10.8466},
1236  {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250}
1237  };
1238 
1239  static const G4double bll3[28][9] = {
1240  {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2.1858E-3, 2.8163E-3, 3.6302E-3, 4.6814E-3, 6.0395E-3},
1241  {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6.1257E-3, 7.5675E-3, 9.3502E-3, 1.1556E-2, 1.4290E-2},
1242  {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6263E-2, 1.9336E-2, 2.2999E-2, 2.7370E-2, 3.2603E-2},
1243  {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4.3483E-2, 4.9898E-2, 5.7304E-2, 6.5884E-2, 7.5861E-2},
1244  {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1413E-2, 9.1539E-2, 1.0304E-1, 1.1617E-1, 1.3121E-1},
1245  {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7451E-1, 1.9244E-1, 2.1244E-1, 2.3496E-1, 2.6044E-1},
1246  {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0180E-1, 3.2751E-1, 3.5608E-1, 3.8803E-1, 4.2401E-1},
1247  {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3635E-1, 4.6973E-1, 5.0672E-1, 5.4798E-1, 5.9436E-1},
1248  {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7943E-1, 6.2028E-1, 6.6549E-1, 7.1589E-1, 7.7252E-1},
1249  {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7836E-1, 9.3355E-1, 9.9462E-1, 1.0627, 1.1394},
1250  {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.2503, 1.3265, 1.4116, 1.5076},
1251  {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.0108, 2.1217, 2.2462, 2.3876},
1252  {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.6971, 2.8391, 2.9994, 3.1822},
1253  {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.8445, 4.0404, 4.2631, 4.5193},
1254  {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.7530, 4.9944, 5.2703, 5.5895},
1255  {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.4811, 5.7618, 6.0840, 6.4583},
1256  {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.1173, 6.4326, 6.7958, 7.2191},
1257  {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.6306, 6.9769, 7.3767, 7.8440},
1258  {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.0827, 7.4570, 7.8900, 8.3972},
1259  {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.8546, 8.2778, 8.7690, 9.3464},
1260  {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.4652, 8.9302, 9.4713, 10.1090},
1261  {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.9851, 9.4866, 10.0713, 10.7619},
1262  {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.2181, 9.7352, 10.3399, 11.0546},
1263  {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.6367, 10.1856, 10.8270, 11.5863},
1264  {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.1793, 10.7688, 11.4590, 12.2775},
1265  {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11.5229, 12.2172, 13.0332, 14.0048},
1266  {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 13.1898, 14.0213, 15.0024, 16.1752},
1267  {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420}
1268  };
1269 
1270  G4double b, bs;
1271  for(i=0; i<nEtaK; ++i) {
1272 
1273  G4double et = Eta[i];
1274  G4double loget = G4Log(et);
1275 
1276  for(j=0; j<nK; ++j) {
1277 
1278  if(j < 10) { b = bk2[i][10-j]; }
1279  else { b = bk1[i][20-j]; }
1280 
1281  CK[j][i] = SK[j]*loget + TK[j] - b;
1282 
1283  if(i == nEtaK-1) {
1284  ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]);
1285  //G4cout << "i= " << i << " j= " << j
1286  // << " CK[j][i]= " << CK[j][i]
1287  // << " ZK[j]= " << ZK[j] << " b= " << b << G4endl;
1288  }
1289  }
1290  //G4cout << G4endl;
1291  if(i < nEtaL) {
1292  //G4cout << " LShell:" <<G4endl;
1293  for(j=0; j<nL; ++j) {
1294 
1295  if(j < 8) {
1296  bs = bls3[i][8-j];
1297  b = bll3[i][8-j];
1298  } else if(j < 17) {
1299  bs = bls2[i][17-j];
1300  b = bll2[i][17-j];
1301  } else {
1302  bs = bls1[i][26-j];
1303  b = bll1[i][26-j];
1304  }
1305  G4double c = SL[j]*loget + TL[j];
1306  CL[j][i] = c - bs - 3.0*b;
1307  if(i == nEtaL-1) {
1308  VL[j] = et*(et*CL[j][i] - UL[j]);
1309  //G4cout << "i= " << i << " j= " << j
1310  // << " CL[j][i]= " << CL[j][i]
1311  // << " VL[j]= " << VL[j] << " b= " << b << " bs= " << bs
1312  // << " et= " << et << G4endl;
1313  //" UL= " << UL[j] << " TL= " << TL[j] << " SL= " << SL[j] <<G4endl;
1314  }
1315  }
1316  //G4cout << G4endl;
1317  }
1318  }
1319 
1320  static const G4double xzk[34] = { 11.7711,
1321  13.3669, 15.5762, 17.1715, 18.7667, 20.8523, 23.0606, 24.901, 26.9861, 29.4394, 31.77,
1322  34.3457, 37.4119, 40.3555, 42.3177, 44.7705, 47.2234, 50.78, 53.8458, 56.4214, 58.3834,
1323  60.9586, 63.6567, 66.5998, 68.807, 71.8728, 74.5706, 77.3911, 81.8056, 85.7297, 89.8988,
1324  93.4549, 96.2753, 99.709};
1325  static const G4double yzk[34] = { 0.70663,
1326  0.72033, 0.73651, 0.74647, 0.75518, 0.76388, 0.77258, 0.78129, 0.78625, 0.7937, 0.79991,
1327  0.80611, 0.8123, 0.8185, 0.82097, 0.82467, 0.82838, 0.83457, 0.83702, 0.84198, 0.8432,
1328  0.84565, 0.84936, 0.85181, 0.85303, 0.85548, 0.85794, 0.8604, 0.86283, 0.86527, 0.86646,
1329  0.86891, 0.87011, 0.87381};
1330 
1331  static const G4double xzl[36] = { 15.5102,
1332  16.7347, 17.9592, 19.551, 21.0204, 22.6122, 24.9388, 27.3878, 29.5918, 31.3061, 32.898,
1333  34.4898, 36.2041, 38.4082, 40.3674, 42.5714, 44.898, 47.4694, 49.9184, 52.7347, 55.9184,
1334  59.3469, 61.9184, 64.6122, 67.4286, 71.4694, 75.2653, 78.3265, 81.2653, 85.551, 88.7347,
1335  91.551, 94.2449, 96.449, 98.4082, 99.7551};
1336  static const G4double yzl[36] = { 0.29875,
1337  0.31746, 0.33368, 0.35239, 0.36985, 0.38732, 0.41102, 0.43472, 0.45343, 0.4659, 0.47713,
1338  0.4896, 0.50083, 0.51331, 0.52328, 0.53077, 0.54075, 0.54823, 0.55572, 0.56445, 0.57193,
1339  0.58191, 0.5869, 0.59189, 0.60062, 0.60686, 0.61435, 0.61809, 0.62183, 0.62931, 0.6343,
1340  0.6368, 0.64054, 0.64304, 0.64428, 0.64678};
1341 
1342  ThetaK = new G4LPhysicsFreeVector(34, xzk[0], xzk[33]);
1343  ThetaL = new G4LPhysicsFreeVector(36, xzl[0], xzl[35]);
1344  for(i=0; i<34; ++i) { ThetaK->PutValues(i, xzk[i], yzk[i]); }
1345  for(i=0; i<36; ++i) { ThetaL->PutValues(i, xzl[i], yzl[i]); }
1346  ThetaK->SetSpline(true);
1347  ThetaL->SetSpline(true);
1348 }
1349 
1350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1351 
static const G4double UL[26]
Float_t x
Definition: compare.C:6
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
static G4double CL[26][28]
G4double Energy(size_t index) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
const G4ParticleDefinition * particle
G4double HighOrderCorrections(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy, G4double cutEnergy)
static G4double CK[20][29]
std::map< G4int, std::vector< G4double > > thcorr
const G4ElementVector * theElementVector
G4int GetAtomicNumber() const
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double MottCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
std::vector< const G4Material * > materialList
static const G4double ZD[11]
static constexpr double twopi_mc2_rcl2
#define G4endl
Definition: G4ios.hh:61
G4double ShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Float_t y
Definition: compare.C:6
G4double BlochCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double BarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
const char * p
Definition: xmltok.h:285
static G4LPhysicsFreeVector * BarkasCorr
G4double SpinCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
const G4String & GetParticleName() const
G4VEmModel * ionHEModel
std::vector< G4String > materialName
const G4Material * material
G4IonTable * GetIonTable() const
G4double GetX1density() const
G4double GetPDGCharge() const
G4double Bethe(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4bool IsMasterThread()
Definition: G4Threading.cc:130
G4IonTable * ionTable
const G4double inveplus
std::vector< G4PhysicsVector * > stopData
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double * GetShellCorrectionVector() const
G4double LShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetPDGMass() const
G4double ComputeIonCorrections(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
static const G4double UK[20]
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double Value(G4double xv, G4double x1, G4double x2, G4double y1, G4double y2) const
void PutValues(size_t index, G4double e, G4double dataValue)
Float_t y2[n_points_geant4]
Definition: compare.C:26
static constexpr double proton_mass_c2
const G4int nmax
const G4String & GetName() const
Definition: G4Material.hh:179
Float_t Z
std::vector< G4int > Zion
G4double KShell(G4double theta, G4double eta)
double G4double
Definition: G4Types.hh:76
static constexpr double amu_c2
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
void BuildCorrectionVector()
static constexpr double electron_mass_c2
static constexpr double MeV
G4double GetTaul() const
void SetSpline(G4bool)
const G4double * atomDensity
static G4int GetNumberOfShells(G4int Z)
static const G4double VK[20]
double A(double temperature)
G4double IonHighOrderCorrections(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
static const G4double * SL[nLA]
G4double GetAdensity() const
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
static G4LPhysicsFreeVector * ThetaK
static G4LPhysicsFreeVector * ThetaL
Float_t mat
G4double GetMeanExcitationEnergy() const
G4double KShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetX0density() const
void AddStoppingData(G4int Z, G4int A, const G4String &materialName, G4PhysicsVector *dVector)
G4EmCorrections(G4int verb)
const G4Material * curMaterial
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4VEmModel * ionLEModel
Double_t Z2
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ParticleDefinition * curParticle
G4double IonBarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
int G4lrint(double ad)
Definition: templates.hh:151
G4double LShell(G4double theta, G4double eta)
G4double GetCdensity() const
static G4double VL[26]
int G4int
Definition: G4Types.hh:78
G4ionEffectiveCharge effCharge
G4int Index(G4double x, const G4double *y, G4int n) const
Float_t norm
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:227
G4double GetMdensity() const
std::vector< G4int > Aion
G4double ShellCorrectionSTD(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
static G4double ZK[20]
G4double EffectiveCharge(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4double Value2(G4double xv, G4double yv, G4double x1, G4double x2, G4double y1, G4double y2, G4double z11, G4double z21, G4double z12, G4double z22) const
G4double DensityCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4PhysicsVector * curVector
static const G4double Eta[29]
const G4Material * GetMaterial() const
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:217
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:504
std::vector< const G4Material * > currmat
void SetupKinematics(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:210
static constexpr double fine_structure_const
Float_t X
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.cc:243
static constexpr double keV
std::vector< const G4ParticleDefinition * > ionList
G4double GetElectronDensity() const
Definition: G4Material.hh:218
void PutValue(size_t index, G4double theValue)
static constexpr double pi
Definition: SystemOfUnits.h:54
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
vec_iX clear()
static constexpr double eplus
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual ~G4EmCorrections()
static G4double tet[DIM]