Geant4  10.02.p01
G4EmCorrections.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: G4EmCorrections.cc 93566 2015-10-26 14:50:09Z 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 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
80 const G4double G4EmCorrections::inveplus = 1.0/CLHEP::eplus;
81 static const G4double twoln10 = 2.0*G4Log(10.0);
82 
83 G4EmCorrections::G4EmCorrections(G4int verb)
84 {
85  particle = nullptr;
86  curParticle= nullptr;
87  material = nullptr;
88  curMaterial= nullptr;
89  curVector = nullptr;
90  kinEnergy = 0.0;
91  ionLEModel = nullptr;
92  ionHEModel = nullptr;
93  nIons = 0;
94  verbose = verb;
95  ncouples = 0;
96  massFactor = 1.0;
97  eth = 2.0*MeV;
98  nbinCorr = 20;
99  eCorrMin = 25.*keV;
100  eCorrMax = 250.*MeV;
101  nist = G4NistManager::Instance();
102  ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
103  BarkasCorr = ThetaK = ThetaL = 0;
104  Initialise();
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
109 G4EmCorrections::~G4EmCorrections()
110 {
111  for(G4int i=0; i<nIons; ++i) {delete stopData[i];}
112  delete BarkasCorr;
113  delete ThetaK;
114  delete ThetaL;
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
119 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p,
120  const G4Material* mat,
121  G4double e, G4double)
122 {
123  // . Z^3 Barkas effect in the stopping power of matter for charged particles
124  // J.C Ashley and R.H.Ritchie
125  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
126  // and ICRU49 report
127  // valid for kineticEnergy < 0.5 MeV
128  // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
129 
130  SetupKinematics(p, mat, e);
131  if(tau <= 0.0) { return 0.0; }
132 
133  G4double Barkas = BarkasCorrection (p, mat, e);
134  G4double Bloch = BlochCorrection (p, mat, e);
135  G4double Mott = MottCorrection (p, mat, e);
136 
137  G4double sum = (2.0*(Barkas + Bloch) + Mott);
138 
139  if(verbose > 1) {
140  G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
141  << " Bloch= " << Bloch << " Mott= " << Mott
142  << " Sum= " << sum << " q2= " << q2 << G4endl;
143  G4cout << " ShellCorrection: " << ShellCorrection(p, mat, e)
144  << " Kshell= " << KShellCorrection(p, mat, e)
145  << " Lshell= " << LShellCorrection(p, mat, e)
146  << " " << mat->GetName() << G4endl;
147  }
148  sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
149  return sum;
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 
154 G4double G4EmCorrections::IonBarkasCorrection(const G4ParticleDefinition* p,
155  const G4Material* mat,
156  G4double e)
157 {
158  // . Z^3 Barkas effect in the stopping power of matter for charged particles
159  // J.C Ashley and R.H.Ritchie
160  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
161  // and ICRU49 report
162  // valid for kineticEnergy < 0.5 MeV
163 
164  SetupKinematics(p, mat, e);
165  G4double res = 0.0;
166  if(tau > 0.0)
167  res = 2.0*BarkasCorrection(p, mat, e)*
168  material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
169  return res;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 
174 G4double G4EmCorrections::ComputeIonCorrections(const G4ParticleDefinition* p,
175  const G4Material* mat,
176  G4double e)
177 {
178  // . Z^3 Barkas effect in the stopping power of matter for charged particles
179  // J.C Ashley and R.H.Ritchie
180  // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
181  // and ICRU49 report
182  // valid for kineticEnergy < 0.5 MeV
183  // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
184  SetupKinematics(p, mat, e);
185  if(tau <= 0.0) { return 0.0; }
186 
187  G4double Barkas = BarkasCorrection (p, mat, e);
188  G4double Bloch = BlochCorrection (p, mat, e);
189  G4double Mott = MottCorrection (p, mat, e);
190 
191  G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott;
192 
193  if(verbose > 1) {
194  G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
195  << " Bloch= " << Bloch << " Mott= " << Mott
196  << " Sum= " << sum << G4endl;
197  }
198  sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
199 
200  if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; }
201  return sum;
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205 
206 G4double G4EmCorrections::IonHighOrderCorrections(const G4ParticleDefinition* p,
207  const G4MaterialCutsCouple* couple,
208  G4double e)
209 {
210 // . Z^3 Barkas effect in the stopping power of matter for charged particles
211 // J.C Ashley and R.H.Ritchie
212 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
213 // and ICRU49 report
214 // valid for kineticEnergy < 0.5 MeV
215 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
216 
217  G4double sum = 0.0;
218 
219  if(ionHEModel) {
220  G4int Z = G4lrint(p->GetPDGCharge()*inveplus);
221  if(Z >= 100) Z = 99;
222  else if(Z < 1) Z = 1;
223 
224  G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2;
225  G4int ionPDG = p->GetPDGEncoding();
226  if(thcorr.find(ionPDG)==thcorr.end()) { // Not found: fill the map
227  std::vector<G4double> v;
228  for(size_t i=0; i<ncouples; ++i){
229  v.push_back(ethscaled*ComputeIonCorrections(p,currmat[i],ethscaled));
230  }
231  thcorr.insert(std::pair< G4int, std::vector<G4double> >(ionPDG,v));
232  }
233 
234  //G4cout << " map size=" << thcorr.size() << G4endl;
235  //for(std::map< G4int, std::vector<G4double> >::iterator
236  // it = thcorr.begin(); it != thcorr.end(); ++it){
237  // G4cout << "\t map element: first (key)=" << it->first
238  // << "\t second (vector): vec size=" << (it->second).size() << G4endl;
239  // for(size_t i=0; i<(it->second).size(); ++i){
240  // G4cout << "\t \t vec element: [" << i << "]=" << (it->second)[i]
241  //<< G4endl; } }
242 
243  G4double rest = (thcorr.find(ionPDG)->second)[couple->GetIndex()];
244 
245  sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
246 
247  if(verbose > 1) {
248  G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl;
249  }
250  }
251  return sum;
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
256 G4double G4EmCorrections::Bethe(const G4ParticleDefinition* p,
257  const G4Material* mat,
258  G4double e)
259 {
260  SetupKinematics(p, mat, e);
261  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
262  G4double eexc2 = eexc*eexc;
263  G4double dedx = 0.5*G4Log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2;
264  return dedx;
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
268 
269 G4double G4EmCorrections::SpinCorrection(const G4ParticleDefinition* p,
270  const G4Material* mat,
271  G4double e)
272 {
273  SetupKinematics(p, mat, e);
274  G4double dedx = 0.5*tmax/(kinEnergy + mass);
275  return 0.5*dedx*dedx;
276 }
277 
278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 
280 G4double G4EmCorrections:: KShellCorrection(const G4ParticleDefinition* p,
281  const G4Material* mat,
282  G4double e)
283 {
284  SetupKinematics(p, mat, e);
285  G4double term = 0.0;
286  for (G4int i = 0; i<numberOfElements; ++i) {
287 
288  G4double Z = (*theElementVector)[i]->GetZ();
289  G4int iz = G4lrint(Z);
290  G4double f = 1.0;
291  G4double Z2= (Z-0.3)*(Z-0.3);
292  if(1 == iz) {
293  f = 0.5;
294  Z2 = 1.0;
295  }
296  G4double eta = ba2/Z2;
297  G4double tet = Z2*(1. + Z2*0.25*alpha2);
298  if(11 < iz) { tet = ThetaK->Value(Z); }
299  term += f*atomDensity[i]*KShell(tet,eta)/Z;
300  }
301 
302  term /= material->GetTotNbOfAtomsPerVolume();
303 
304  return term;
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 
309 G4double G4EmCorrections:: LShellCorrection(const G4ParticleDefinition* p,
310  const G4Material* mat,
311  G4double e)
312 {
313  SetupKinematics(p, mat, e);
314  G4double term = 0.0;
315  for (G4int i = 0; i<numberOfElements; ++i) {
316 
317  G4double Z = (*theElementVector)[i]->GetZ();
318  G4int iz = G4lrint(Z);
319  if(2 < iz) {
320  G4double Zeff = Z - ZD[10];
321  if(iz < 10) { Zeff = Z - ZD[iz]; }
322  G4double Z2= Zeff*Zeff;
323  G4double f = 0.125;
324  G4double eta = ba2/Z2;
325  G4double tet = ThetaL->Value(Z);
326  G4int nmax = std::min(4,G4AtomicShells::GetNumberOfShells(iz));
327  for(G4int j=1; j<nmax; ++j) {
328  G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
329  if(15 >= iz) {
330  if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
331  else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
332  }
333  //G4cout << " LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
334  // << " ThetaL= " << tet << G4endl;
335  term += f*ne*atomDensity[i]*LShell(tet,eta)/Z;
336  }
337  }
338  }
339 
340  term /= material->GetTotNbOfAtomsPerVolume();
341 
342  return term;
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346 
347 G4double G4EmCorrections::KShell(G4double tet, G4double eta)
348 {
349  G4double corr = 0.0;
350 
351  G4double x = tet;
352  G4int itet = 0;
353  G4int ieta = 0;
354  if(tet < TheK[0]) {
355  x = TheK[0];
356  } else if(tet > TheK[nK-1]) {
357  x = TheK[nK-1];
358  itet = nK-2;
359  } else {
360  itet = Index(x, TheK, nK);
361  }
362  // assimptotic case
363  if(eta >= Eta[nEtaK-1]) {
364  corr =
365  (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) +
366  Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta +
367  Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta;
368  } else {
369  G4double y = eta;
370  if(eta < Eta[0]) {
371  y = Eta[0];
372  } else {
373  ieta = Index(y, Eta, nEtaK);
374  }
375  corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
376  CK[itet][ieta], CK[itet+1][ieta],
377  CK[itet][ieta+1], CK[itet+1][ieta+1]);
378  //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]
379  // <<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
380  // <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta]
381  // <<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl;
382  }
383  //G4cout << "Kshell: tet= " << tet << " eta= " << eta << " C= " << corr
384  // << " itet= " << itet << " ieta= " << ieta <<G4endl;
385  return corr;
386 }
387 
388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 
390 G4double G4EmCorrections::LShell(G4double tet, G4double eta)
391 {
392  G4double corr = 0.0;
393 
394  G4double x = tet;
395  G4int itet = 0;
396  G4int ieta = 0;
397  if(tet < TheL[0]) {
398  x = TheL[0];
399  } else if(tet > TheL[nL-1]) {
400  x = TheL[nL-1];
401  itet = nL-2;
402  } else {
403  itet = Index(x, TheL, nL);
404  }
405 
406  // assimptotic case
407  if(eta >= Eta[nEtaL-1]) {
408  corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1])
409  + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta
410  )/eta;
411  } else {
412  G4double y = eta;
413  if(eta < Eta[0]) {
414  y = Eta[0];
415  } else {
416  ieta = Index(y, Eta, nEtaL);
417  }
418  corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1],
419  CL[itet][ieta], CL[itet+1][ieta],
420  CL[itet][ieta+1], CL[itet+1][ieta+1]);
421  //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet]
422  // <<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
423  // <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta]
424  // <<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl;
425  }
426  //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<eta<<" itet= "<<itet
427  // <<" ieta= "<<ieta<<" Corr= "<<corr<<G4endl;
428  return corr;
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432 
433 G4double G4EmCorrections::ShellCorrectionSTD(const G4ParticleDefinition* p,
434  const G4Material* mat,
435  G4double e)
436 {
437  SetupKinematics(p, mat, e);
438  G4double taulim= 8.0*MeV/mass;
439  G4double bg2lim= taulim * (taulim+2.0);
440 
441  G4double* shellCorrectionVector =
442  material->GetIonisation()->GetShellCorrectionVector();
443  G4double sh = 0.0;
444  G4double x = 1.0;
445  G4double taul = material->GetIonisation()->GetTaul();
446 
447  if ( bg2 >= bg2lim ) {
448  for (G4int k=0; k<3; ++k) {
449  x *= bg2 ;
450  sh += shellCorrectionVector[k]/x;
451  }
452 
453  } else {
454  for (G4int k=0; k<3; ++k) {
455  x *= bg2lim ;
456  sh += shellCorrectionVector[k]/x;
457  }
458  sh *= G4Log(tau/taul)/G4Log(taulim/taul);
459  }
460  sh *= 0.5;
461  return sh;
462 }
463 
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 
467 G4double G4EmCorrections::ShellCorrection(const G4ParticleDefinition* p,
468  const G4Material* mat,
469  G4double ekin)
470 {
471  SetupKinematics(p, mat, ekin);
472 
473  G4double term = 0.0;
474  //G4cout << "### G4EmCorrections::ShellCorrection " << mat->GetName()
475  // << " " << ekin/MeV << " MeV " << G4endl;
476  for (G4int i = 0; i<numberOfElements; ++i) {
477 
478  G4double res = 0.0;
479  G4double res0 = 0.0;
480  G4double Z = (*theElementVector)[i]->GetZ();
481  G4int iz = G4lrint(Z);
482  G4double Z2= (Z-0.3)*(Z-0.3);
483  G4double f = 1.0;
484  if(1 == iz) {
485  f = 0.5;
486  Z2 = 1.0;
487  }
488  G4double eta = ba2/Z2;
489  G4double tet = Z2*(1. + Z2*0.25*alpha2);
490  if(11 < iz) { tet = ThetaK->Value(Z); }
491  res0 = f*KShell(tet,eta);
492  res += res0;
493  //G4cout << " Z= " << iz << " Shell 0" << " tet= " << tet
494  // << " eta= " << eta << " resK= " << res0 << G4endl;
495  if(2 < iz) {
496  G4double Zeff = Z - ZD[10];
497  if(iz < 10) { Zeff = Z - ZD[iz]; }
498  Z2= Zeff*Zeff;
499  eta = ba2/Z2;
500  f = 0.125;
501  tet = ThetaL->Value(Z);
502  G4int ntot = G4AtomicShells::GetNumberOfShells(iz);
503  G4int nmax = std::min(4, ntot);
504  G4double norm = 0.0;
505  G4double eshell = 0.0;
506  for(G4int j=1; j<nmax; ++j) {
507  G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
508  if(15 >= iz) {
509  if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
510  else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
511  }
512  norm += ne;
513  eshell += tet*ne;
514  res0 = f*ne*LShell(tet,eta);
515  res += res0;
516  //G4cout << " Z= " << iz << " Shell " << j << " Ne= " << ne
517  // << " tet= " << tet << " eta= " << eta
518  // << " resL= " << res0 << G4endl;
519  }
520  if(ntot > nmax) {
521  eshell /= norm;
522  // Add M-shell
523  if(28 > iz) {
524  res += f*(iz - 10)*LShell(eshell,HM[iz-11]*eta);
525  } else if(63 > iz) {
526  res += f*18*LShell(eshell,HM[iz-11]*eta);
527  } else {
528  res += f*18*LShell(eshell,HM[52]*eta);
529  }
530  // Add N-shell
531  if(32 < iz) {
532  if(60 > iz) {
533  res += f*(iz - 28)*LShell(eshell,HN[iz-33]*eta);
534  } else if(63 > iz) {
535  res += 4*LShell(eshell,HN[iz-33]*eta);
536  } else {
537  res += 4*LShell(eshell,HN[30]*eta);
538  }
539  // Add O-P-shells
540  if(60 < iz) {
541  res += f*(iz - 60)*LShell(eshell,150*eta);
542  }
543  }
544  }
545  }
546  term += res*atomDensity[i]/Z;
547  }
548 
549  term /= material->GetTotNbOfAtomsPerVolume();
550  //G4cout << "# Shell Correction= " << term << G4endl;
551  return term;
552 }
553 
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
555 
556 G4double G4EmCorrections::DensityCorrection(const G4ParticleDefinition* p,
557  const G4Material* mat,
558  G4double e)
559 {
560  SetupKinematics(p, mat, e);
561 
562  G4double cden = material->GetIonisation()->GetCdensity();
563  G4double mden = material->GetIonisation()->GetMdensity();
564  G4double aden = material->GetIonisation()->GetAdensity();
565  G4double x0den = material->GetIonisation()->GetX0density();
566  G4double x1den = material->GetIonisation()->GetX1density();
567 
568  G4double dedx = 0.0;
569 
570  // density correction
571  G4double x = G4Log(bg2)/twoln10;
572  if ( x >= x0den ) {
573  dedx = twoln10*x - cden ;
574  if ( x < x1den ) dedx += aden*G4Exp(G4Log(x1den-x)*mden) ;
575  }
576 
577  return dedx;
578 }
579 
580 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
581 
582 G4double G4EmCorrections::BarkasCorrection(const G4ParticleDefinition* p,
583  const G4Material* mat,
584  G4double e)
585 {
586  // . Z^3 Barkas effect in the stopping power of matter for charged particles
587  // J.C Ashley and R.H.Ritchie
588  // Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397
589  // valid for kineticEnergy > 0.5 MeV
590 
591  SetupKinematics(p, mat, e);
592  G4double BarkasTerm = 0.0;
593 
594  for (G4int i = 0; i<numberOfElements; ++i) {
595 
596  G4double Z = (*theElementVector)[i]->GetZ();
597  G4int iz = G4lrint(Z);
598  if(iz == 47) {
599  BarkasTerm += atomDensity[i]*0.006812*G4Exp(-G4Log(beta)*0.9);
600  } else if(iz >= 64) {
601  BarkasTerm += atomDensity[i]*0.002833*G4Exp(-G4Log(beta)*1.2);
602  } else {
603 
604  G4double X = ba2 / Z;
605  G4double b = 1.3;
606  if(1 == iz) {
607  if(material->GetName() == "G4_lH2") { b = 0.6; }
608  else { b = 1.8; }
609  }
610  else if(2 == iz) { b = 0.6; }
611  else if(10 >= iz) { b = 1.8; }
612  else if(17 >= iz) { b = 1.4; }
613  else if(18 == iz) { b = 1.8; }
614  else if(25 >= iz) { b = 1.4; }
615  else if(50 >= iz) { b = 1.35;}
616 
617  G4double W = b/std::sqrt(X);
618 
619  G4double val = BarkasCorr->Value(W);
620  if(W > BarkasCorr->Energy(46)) {
621  val *= BarkasCorr->Energy(46)/W;
622  }
623  // G4cout << "i= " << i << " b= " << b << " W= " << W
624  // << " Z= " << Z << " X= " << X << " val= " << val<< G4endl;
625  BarkasTerm += val*atomDensity[i] / (std::sqrt(Z*X)*X);
626  }
627  }
628 
629  BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
630 
631  // temporary protection
632  //if(charge < -7.0 ) { BarkasTerm *= (-7.0/charge); }
633 
634  return BarkasTerm;
635 }
636 
637 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
638 
639 G4double G4EmCorrections::BlochCorrection(const G4ParticleDefinition* p,
640  const G4Material* mat,
641  G4double e)
642 {
643  SetupKinematics(p, mat, e);
644 
645  G4double y2 = q2/ba2;
646 
647  G4double term = 1.0/(1.0 + y2);
648  G4double del;
649  G4double j = 1.0;
650  do {
651  j += 1.0;
652  del = 1.0/(j* (j*j + y2));
653  term += del;
654  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
655  } while (del > 0.01*term);
656 
657  G4double res = -y2*term;
658  // temporary protection
659  //if(q2 > 49. && res < -0.2) { res = -0.2; }
660 
661  return res;
662 }
663 
664 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
665 
666 G4double G4EmCorrections::MottCorrection(const G4ParticleDefinition* p,
667  const G4Material* mat,
668  G4double e)
669 {
670  SetupKinematics(p, mat, e);
671  G4double mterm = CLHEP::pi*fine_structure_const*beta*charge;
672  return mterm;
673 }
674 
675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
676 
677 G4double G4EmCorrections::NuclearDEDX(const G4ParticleDefinition* p,
678  const G4Material* mat,
679  G4double e,
680  G4bool fluct)
681 {
682  G4double nloss = 0.0;
683  if(e <= 0.0) return nloss;
684  SetupKinematics(p, mat, e);
685 
686  lossFlucFlag = fluct;
687 
688  // Projectile nucleus
689  G4double z1 = std::abs(particle->GetPDGCharge()*inveplus);
690  G4double mass1 = mass/amu_c2;
691 
692  // loop for the elements in the material
693  for (G4int iel=0; iel<numberOfElements; iel++) {
694  const G4Element* element = (*theElementVector)[iel] ;
695  G4double z2 = element->GetZ();
696  G4double mass2 = element->GetN();
697  nloss += (NuclearStoppingPower(kinEnergy, z1, z2, mass1, mass2))
698  * atomDensity[iel] ;
699  }
700  nloss *= theZieglerFactor;
701  return nloss;
702 }
703 
704 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
705 
706 G4double G4EmCorrections::NuclearStoppingPower(G4double kineticEnergy,
707  G4double z1, G4double z2,
708  G4double mass1, G4double mass2)
709 {
710  G4double energy = kineticEnergy/keV ; // energy in keV
711  G4double nloss = 0.0;
712 
713  G4double rm;
714  if(z1 > 1.5) rm = (mass1 + mass2) * ( Z23[G4lrint(z1)] + Z23[G4lrint(z2)] ) ;
715  else rm = (mass1 + mass2) * nist->GetZ13(G4lrint(z2));
716 
717  G4double er = 32.536 * mass2 * energy / ( z1 * z2 * rm ) ; // reduced energy
718 
719  if (er >= ed[0]) { nloss = a[0]; }
720  else {
721  // the table is inverse in energy
722  for (G4int i=102; i>=0; --i) {
723  if (er <= ed[i]) {
724  nloss = (a[i] - a[i+1])*(er - ed[i+1])/(ed[i] - ed[i+1]) + a[i+1];
725  break;
726  }
727  }
728  }
729 
730  // Stragling
731  if(lossFlucFlag) {
732  G4double sig = 4.0 * mass1 * mass2 / ((mass1 + mass2)*(mass1 + mass2)*
733  (4.0 + 0.197/(er*er) + 6.584/er));
734 
735  nloss *= G4RandGauss::shoot(1.0,sig) ;
736  }
737 
738  nloss *= 8.462 * z1 * z2 * mass1 / rm ; // Return to [ev/(10^15 atoms/cm^2]
739 
740  if ( nloss < 0.0) nloss = 0.0 ;
741 
742  return nloss;
743 }
744 
745 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
746 
747 G4double
748 G4EmCorrections::EffectiveChargeCorrection(const G4ParticleDefinition* p,
749  const G4Material* mat,
750  G4double ekin)
751 {
752  G4double factor = 1.0;
753  if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || nIons <= 0) { return factor; }
754  /*
755  if(verbose > 1) {
756  G4cout << "EffectiveChargeCorrection: " << p->GetParticleName()
757  << " in " << mat->GetName()
758  << " ekin(MeV)= " << ekin/MeV << G4endl;
759  }
760  */
761  if(p != curParticle || mat != curMaterial) {
762  curParticle = p;
763  curMaterial = mat;
764  curVector = 0;
765  currentZ = p->GetAtomicNumber();
766  if(verbose > 1) {
767  G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= "
768  << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl;
769  }
770  massFactor = proton_mass_c2/p->GetPDGMass();
771  idx = -1;
772 
773  for(G4int i=0; i<nIons; ++i) {
774  if(materialList[i] == mat && currentZ == Zion[i]) {
775  idx = i;
776  break;
777  }
778  }
779  // G4cout << " idx= " << idx << " dz= " << G4endl;
780  if(idx >= 0) {
781  if(!ionList[idx]) { BuildCorrectionVector(); }
782  if(ionList[idx]) { curVector = stopData[idx]; }
783  } else { return factor; }
784  }
785  if(curVector) {
786  factor = curVector->Value(ekin*massFactor);
787  if(verbose > 1) {
788  G4cout << "E= " << ekin << " factor= " << factor << " massfactor= "
789  << massFactor << G4endl;
790  }
791  }
792  return factor;
793 }
794 
795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
796 
797 void G4EmCorrections::AddStoppingData(G4int Z, G4int A,
798  const G4String& mname,
799  G4PhysicsVector* dVector)
800 {
801  G4int i = 0;
802  for(; i<nIons; ++i) {
803  if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break;
804  }
805  if(i == nIons) {
806  Zion.push_back(Z);
807  Aion.push_back(A);
808  materialName.push_back(mname);
809  materialList.push_back(0);
810  ionList.push_back(0);
811  stopData.push_back(dVector);
812  nIons++;
813  if(verbose>1) {
814  G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname
815  << " idx= " << i << G4endl;
816  }
817  }
818 }
819 
820 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
821 
822 void G4EmCorrections::BuildCorrectionVector()
823 {
824  if(!ionLEModel || !ionHEModel) {
825  return;
826  }
827 
828  const G4ParticleDefinition* ion = curParticle;
829  G4int Z = Zion[idx];
830  if(currentZ != Z) {
831  ion = ionTable->GetIon(Z, Aion[idx], 0);
832  }
833  //G4cout << "BuildCorrectionVector: idx= " << idx << " Z= " << Z
834  // << " curZ= " << currentZ << G4endl;
835 
836  // G4double A = nist->GetAtomicMassAmu(Z);
837  G4double A = G4double(ion->GetBaryonNumber());
838  G4PhysicsVector* v = stopData[idx];
839 
840  const G4ParticleDefinition* p = G4GenericIon::GenericIon();
841  G4double massRatio = proton_mass_c2/ion->GetPDGMass();
842 
843  if(verbose>1) {
844  G4cout << "BuildCorrectionVector: Stopping for "
845  << curParticle->GetParticleName() << " in "
846  << materialName[idx] << " Ion Z= " << Z << " A= " << A
847  << " massRatio= " << massRatio << G4endl;
848  }
849 
850  G4PhysicsLogVector* vv =
851  new G4PhysicsLogVector(eCorrMin,eCorrMax,nbinCorr);
852  vv->SetSpline(true);
853  G4double e, eion, dedx, dedx1;
854  G4double eth0 = v->Energy(0);
855  G4double escal = eth/massRatio;
856  G4double qe =
857  effCharge.EffectiveChargeSquareRatio(ion, curMaterial, escal);
858  G4double dedxt =
859  ionLEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe;
860  G4double dedx1t =
861  ionHEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe
862  + ComputeIonCorrections(curParticle, curMaterial, escal);
863  G4double rest = escal*(dedxt - dedx1t);
864  //G4cout << "Escal(MeV)= "<<escal<<" dedxt0= " <<dedxt
865  // << " dedxt1= " << dedx1t << G4endl;
866 
867  for(G4int i=0; i<=nbinCorr; ++i) {
868  e = vv->Energy(i);
869  escal = e/massRatio;
870  eion = escal/A;
871  if(eion <= eth0) {
872  dedx = v->Value(eth0)*std::sqrt(eion/eth0);
873  } else {
874  dedx = v->Value(eion);
875  }
876  qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal);
877  if(e <= eth) {
878  dedx1 = ionLEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe;
879  } else {
880  dedx1 = ionHEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe +
881  ComputeIonCorrections(curParticle, curMaterial, escal) + rest/escal;
882  }
883  vv->PutValue(i, dedx/dedx1);
884  if(verbose>1) {
885  G4cout << " E(meV)= " << e/MeV << " Correction= " << dedx/dedx1
886  << " " << dedx << " " << dedx1
887  << " massF= " << massFactor << G4endl;
888  }
889  }
890  delete v;
891  ionList[idx] = ion;
892  stopData[idx] = vv;
893  if(verbose>1) { G4cout << "End data set " << G4endl; }
894 }
895 
896 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
897 
898 void G4EmCorrections::InitialiseForNewRun()
899 {
900  G4ProductionCutsTable* tb = G4ProductionCutsTable::GetProductionCutsTable();
901  ncouples = tb->GetTableSize();
902  if(currmat.size() != ncouples) {
903  currmat.resize(ncouples);
904  for(std::map< G4int, std::vector<G4double> >::iterator it =
905  thcorr.begin(); it != thcorr.end(); ++it){
906  (it->second).clear();
907  }
908  thcorr.clear();
909  for(size_t i=0; i<ncouples; ++i) {
910  currmat[i] = tb->GetMaterialCutsCouple(i)->GetMaterial();
911  G4String nam = currmat[i]->GetName();
912  for(G4int j=0; j<nIons; ++j) {
913  if(nam == materialName[j]) { materialList[j] = currmat[i]; }
914  }
915  }
916  }
917 }
918 
919 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
920 
921 void G4EmCorrections::Initialise()
922 {
923 // . Z^3 Barkas effect in the stopping power of matter for charged particles
924 // J.C Ashley and R.H.Ritchie
925 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
926  G4int i, j;
927  const G4double fTable[47][2] = {
928  { 0.02, 21.5},
929  { 0.03, 20.0},
930  { 0.04, 18.0},
931  { 0.05, 15.6},
932  { 0.06, 15.0},
933  { 0.07, 14.0},
934  { 0.08, 13.5},
935  { 0.09, 13.},
936  { 0.1, 12.2},
937  { 0.2, 9.25},
938  { 0.3, 7.0},
939  { 0.4, 6.0},
940  { 0.5, 4.5},
941  { 0.6, 3.5},
942  { 0.7, 3.0},
943  { 0.8, 2.5},
944  { 0.9, 2.0},
945  { 1.0, 1.7},
946  { 1.2, 1.2},
947  { 1.3, 1.0},
948  { 1.4, 0.86},
949  { 1.5, 0.7},
950  { 1.6, 0.61},
951  { 1.7, 0.52},
952  { 1.8, 0.5},
953  { 1.9, 0.43},
954  { 2.0, 0.42},
955  { 2.1, 0.3},
956  { 2.4, 0.2},
957  { 3.0, 0.13},
958  { 3.08, 0.1},
959  { 3.1, 0.09},
960  { 3.3, 0.08},
961  { 3.5, 0.07},
962  { 3.8, 0.06},
963  { 4.0, 0.051},
964  { 4.1, 0.04},
965  { 4.8, 0.03},
966  { 5.0, 0.024},
967  { 5.1, 0.02},
968  { 6.0, 0.013},
969  { 6.5, 0.01},
970  { 7.0, 0.009},
971  { 7.1, 0.008},
972  { 8.0, 0.006},
973  { 9.0, 0.0032},
974  { 10.0, 0.0025} };
975 
976  BarkasCorr = new G4LPhysicsFreeVector(47, 0.02, 10.);
977  for(i=0; i<47; ++i) { BarkasCorr->PutValues(i, fTable[i][0], fTable[i][1]); }
978  BarkasCorr->SetSpline(true);
979 
980  const G4double nuca[104][2] = {
981  { 1.0E+8, 5.831E-8},
982  { 8.0E+7, 7.288E-8},
983  { 6.0E+7, 9.719E-8},
984  { 5.0E+7, 1.166E-7},
985  { 4.0E+7, 1.457E-7},
986  { 3.0E+7, 1.942E-7},
987  { 2.0E+7, 2.916E-7},
988  { 1.5E+7, 3.887E-7},
989 
990  { 1.0E+7, 5.833E-7},
991  { 8.0E+6, 7.287E-7},
992  { 6.0E+6, 9.712E-7},
993  { 5.0E+6, 1.166E-6},
994  { 4.0E+6, 1.457E-6},
995  { 3.0E+6, 1.941E-6},
996  { 2.0E+6, 2.911E-6},
997  { 1.5E+6, 3.878E-6},
998 
999  { 1.0E+6, 5.810E-6},
1000  { 8.0E+5, 7.262E-6},
1001  { 6.0E+5, 9.663E-6},
1002  { 5.0E+5, 1.157E-5},
1003  { 4.0E+5, 1.442E-5},
1004  { 3.0E+5, 1.913E-5},
1005  { 2.0E+5, 2.845E-5},
1006  { 1.5E+5, 3.762E-5},
1007 
1008  { 1.0E+5, 5.554E-5},
1009  { 8.0E+4, 6.866E-5},
1010  { 6.0E+4, 9.020E-5},
1011  { 5.0E+4, 1.070E-4},
1012  { 4.0E+4, 1.319E-4},
1013  { 3.0E+4, 1.722E-4},
1014  { 2.0E+4, 2.499E-4},
1015  { 1.5E+4, 3.248E-4},
1016 
1017  { 1.0E+4, 4.688E-4},
1018  { 8.0E+3, 5.729E-4},
1019  { 6.0E+3, 7.411E-4},
1020  { 5.0E+3, 8.718E-4},
1021  { 4.0E+3, 1.063E-3},
1022  { 3.0E+3, 1.370E-3},
1023  { 2.0E+3, 1.955E-3},
1024  { 1.5E+3, 2.511E-3},
1025 
1026  { 1.0E+3, 3.563E-3},
1027  { 8.0E+2, 4.314E-3},
1028  { 6.0E+2, 5.511E-3},
1029  { 5.0E+2, 6.430E-3},
1030  { 4.0E+2, 7.756E-3},
1031  { 3.0E+2, 9.855E-3},
1032  { 2.0E+2, 1.375E-2},
1033  { 1.5E+2, 1.736E-2},
1034 
1035  { 1.0E+2, 2.395E-2},
1036  { 8.0E+1, 2.850E-2},
1037  { 6.0E+1, 3.552E-2},
1038  { 5.0E+1, 4.073E-2},
1039  { 4.0E+1, 4.802E-2},
1040  { 3.0E+1, 5.904E-2},
1041  { 1.5E+1, 9.426E-2},
1042 
1043  { 1.0E+1, 1.210E-1},
1044  { 8.0E+0, 1.377E-1},
1045  { 6.0E+0, 1.611E-1},
1046  { 5.0E+0, 1.768E-1},
1047  { 4.0E+0, 1.968E-1},
1048  { 3.0E+0, 2.235E-1},
1049  { 2.0E+0, 2.613E-1},
1050  { 1.5E+0, 2.871E-1},
1051 
1052  { 1.0E+0, 3.199E-1},
1053  { 8.0E-1, 3.354E-1},
1054  { 6.0E-1, 3.523E-1},
1055  { 5.0E-1, 3.609E-1},
1056  { 4.0E-1, 3.693E-1},
1057  { 3.0E-1, 3.766E-1},
1058  { 2.0E-1, 3.803E-1},
1059  { 1.5E-1, 3.788E-1},
1060 
1061  { 1.0E-1, 3.711E-1},
1062  { 8.0E-2, 3.644E-1},
1063  { 6.0E-2, 3.530E-1},
1064  { 5.0E-2, 3.444E-1},
1065  { 4.0E-2, 3.323E-1},
1066  { 3.0E-2, 3.144E-1},
1067  { 2.0E-2, 2.854E-1},
1068  { 1.5E-2, 2.629E-1},
1069 
1070  { 1.0E-2, 2.298E-1},
1071  { 8.0E-3, 2.115E-1},
1072  { 6.0E-3, 1.883E-1},
1073  { 5.0E-3, 1.741E-1},
1074  { 4.0E-3, 1.574E-1},
1075  { 3.0E-3, 1.372E-1},
1076  { 2.0E-3, 1.116E-1},
1077  { 1.5E-3, 9.559E-2},
1078 
1079  { 1.0E-3, 7.601E-2},
1080  { 8.0E-4, 6.668E-2},
1081  { 6.0E-4, 5.605E-2},
1082  { 5.0E-4, 5.008E-2},
1083  { 4.0E-4, 4.352E-2},
1084  { 3.0E-4, 3.617E-2},
1085  { 2.0E-4, 2.768E-2},
1086  { 1.5E-4, 2.279E-2},
1087 
1088  { 1.0E-4, 1.723E-2},
1089  { 8.0E-5, 1.473E-2},
1090  { 6.0E-5, 1.200E-2},
1091  { 5.0E-5, 1.052E-2},
1092  { 4.0E-5, 8.950E-3},
1093  { 3.0E-5, 7.246E-3},
1094  { 2.0E-5, 5.358E-3},
1095  { 1.5E-5, 4.313E-3},
1096  { 0.0, 3.166E-3}
1097  };
1098 
1099  for(i=0; i<104; ++i) {
1100  ed[i] = nuca[i][0];
1101  a[i] = nuca[i][1];
1102  }
1103 
1104  // Constants
1105  theZieglerFactor = eV*cm2*1.0e-15 ;
1106  alpha2 = fine_structure_const*fine_structure_const;
1107  lossFlucFlag = true;
1108 
1109  // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111.
1110  // "Shell corrections for K- and L- electrons
1111 
1112  nK = 20;
1113  nL = 26;
1114  nEtaK = 29;
1115  nEtaL = 28;
1116 
1117  const G4double d[11] =
1118  {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15};
1119  const G4double thek[20] =
1120  {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78,
1121  0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95};
1122  const G4double sk[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137,
1123  1.7754, 1.7396, 1.7223, 1.7063, 1.6752,
1124  1.6461, 1.6189, 1.5933, 1.5811, 1.5693,
1125  1.5467, 1.5254, 1.5053, 1.4863, 1.4772};
1126  const G4double tk[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608,
1127  2.4388, 2.4163, 2.4044, 2.3933, 2.3701,
1128  2.3466, 2.3229, 2.2992, 2.2872, 2.2753,
1129  2.2515, 2.2277, 2.2040, 2.1804, 2.1686};
1130  const G4double uk[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662,
1131  2.0817, 2.0945, 2.0999, 2.1049, 2.1132,
1132  2.1197, 2.1246, 2.1280, 2.1292, 2.1301,
1133  2.1310, 2.1310, 2.1300, 2.1283, 2.1271};
1134  const G4double vk[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247,
1135  8.3219, 8.3201, 8.3194, 8.3191, 8.3188,
1136  8.3191, 8.3199, 8.3211, 8.3218, 8.3226,
1137  8.3244, 8.3264, 8.3285, 8.3308, 8.3320};
1138 
1139  for(i=0; i<11; ++i) { ZD[i] = d[i];}
1140 
1141  for(i=0; i<nK; ++i) {
1142  TheK[i] = thek[i];
1143  SK[i] = sk[i];
1144  TK[i] = tk[i];
1145  UK[i] = uk[i];
1146  VK[i] = vk[i];
1147  }
1148 
1149  const G4double thel[26] =
1150  {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40,
1151  0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.52, 0.54, 0.55, 0.56,
1152  0.58, 0.60, 0.62, 0.64, 0.65, 0.66};
1153  const G4double sl[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283,
1154  10.3424, 10.0371, 9.7537, 9.2443, 8.8005,
1155  8.4114, 8.0683, 7.9117, 7.7641, 7.4931,
1156  7.2506, 7.0327, 6.8362, 6.7452, 6.6584,
1157  6.4969, 6.3498, 6.2154, 6.0923, 6.0345, 5.9792};
1158  const G4double tl[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226,
1159  28.6128, 28.1449, 27.6991, 26.8674, 26.1061,
1160  25.4058, 24.7587, 24.4531, 24.1583, 23.5992,
1161  23.0771, 22.5880, 22.1285, 21.9090, 21.6958,
1162  21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546};
1163  const G4double ul[26] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828,
1164  1.4379, 1.5032, 1.5617, 1.6608, 1.7401,
1165  1.8036, 1.8543, 1.8756, 1.8945, 1.9262,
1166  1.9508, 1.9696, 1.9836, 1.9890, 1.9935,
1167  2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028};
1168  for(i=0; i<nL; ++i) {
1169  TheL[i] = thel[i];
1170  SL[i] = sl[i];
1171  TL[i] = tl[i];
1172  UL[i] = ul[i];
1173  }
1174 
1175  const G4double eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02,
1176  0.03, 0.04, 0.05, 0.06, 0.08,
1177  0.1, 0.15, 0.2, 0.3, 0.4,
1178  0.5, 0.6, 0.7, 0.8, 1.0,
1179  1.2, 1.4, 1.5, 1.7, 2.0, 3.0, 5.0, 7.0, 10.};
1180 
1181  const G4double bk1[29][11] = {
1182  {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},
1183  {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},
1184  {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},
1185  {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},
1186  {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},
1187  {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},
1188  {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},
1189  {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},
1190  {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},
1191  {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},
1192  {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},
1193  {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},
1194  {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},
1195  {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},
1196  {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},
1197  {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},
1198  {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},
1199  {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},
1200  {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},
1201  {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.21317, 1.23921, 1.25257, 1.26616, 1.29408, 1.32303},
1202  {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.49875, 1.52842, 1.54364, 1.55913, 1.59095, 1.62396},
1203  {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.76244, 1.79526, 1.81210, 1.82925, 1.86448, 1.90104},
1204  {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.88617, 1.92042, 1.93800, 1.95590, 1.99269, 2.03087},
1205  {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.11827, 2.15555, 2.17409, 2.19337, 2.23302, 2.27419},
1206  {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.42395, 2.46429, 2.48401, 2.50612, 2.54955, 2.59468},
1207  {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.24748, 3.29664, 3.32192, 3.34770, 3.40081, 3.45611},
1208  {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.27010, 4.33056, 4.36172, 4.39353, 4.45918, 4.52772},
1209  {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.91402, 4.98200, 5.01707, 5.05290, 5.12695, 5.20436},
1210  {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359}
1211  };
1212 
1213  const G4double bk2[29][11] = {
1214  {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},
1215  {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},
1216  {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},
1217  {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},
1218  {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},
1219  {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},
1220  {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},
1221  {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},
1222  {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},
1223  {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},
1224  {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},
1225  {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},
1226  {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},
1227  {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},
1228  {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},
1229  {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},
1230  {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},
1231  {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},
1232  {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1.10667, 1.13615, 1.16692, 1.19907, 1.21570, 1.23272},
1233  {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1.45057, 1.48582, 1.52263, 1.56111, 1.58102, 1.60140},
1234  {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1.76954, 1.80983, 1.85192, 1.89596, 1.91876, 1.94211},
1235  {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2.06251, 2.10727, 2.15406, 2.20304, 2.22842, 2.25442},
1236  {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2.19963, 2.24644, 2.29539, 2.34666, 2.37323, 2.40045},
1237  {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2.45641, 2.50703, 2.56000, 2.61552, 2.64430, 2.67381},
1238  {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2.79481, 2.85052, 2.90887, 2.97009, 3.00185, 3.03442},
1239  {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3.70268, 3.77170, 3.84418, 3.92040, 3.96003, 4.00073},
1240  {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4.83543, 4.92217, 5.01353, 5.10992, 5.16014, 5.21181},
1241  {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5.55344, 5.65226, 5.79496, 5.90517, 5.96269, 6.02191},
1242  {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535}
1243  };
1244 
1245  const G4double bls1[28][10] = {
1246  {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},
1247  {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},
1248  {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},
1249  {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},
1250  {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},
1251  {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},
1252  {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},
1253  {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},
1254  {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},
1255  {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},
1256  {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},
1257  {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},
1258  {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.0037, 1.0372, 1.0723, 1.1092, 1.1284, 1.1480},
1259  {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.6438, 1.6899, 1.7382, 1.7632, 1.7889},
1260  {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.1659, 2.2199, 2.2765, 2.3059, 2.3360},
1261  {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.5587, 2.6190, 2.6821, 2.7148, 2.7484},
1262  {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.8884, 2.9538, 3.0222, 3.0577, 3.0941},
1263  {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.1652, 3.2349, 3.3079, 3.3457, 3.3845},
1264  {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.4232, 3.4965, 3.5734, 3.6133, 3.6542},
1265  {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.8866, 3.9661, 4.0495, 4.0928, 4.1371},
1266  {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.2130, 4.2974, 4.3861, 4.4321, 4.4794},
1267  {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.5170, 4.6056, 4.6988, 4.7471, 4.7968},
1268  {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.6521, 4.7426, 4.8378, 4.8872, 4.9379},
1269  {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.8949, 4.9888, 5.0876, 5.1388, 5.1915},
1270  {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.2053, 5.3036, 5.4070, 5.4607, 5.5159},
1271  {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.9489, 6.0581, 6.1730, 6.2328, 6.2943},
1272  {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.8003, 6.9231, 7.0525, 7.1199, 7.1892},
1273  {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384}
1274  };
1275 
1276  const G4double bls2[28][10] = {
1277  {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},
1278  {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},
1279  {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},
1280  {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},
1281  {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},
1282  {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},
1283  {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},
1284  {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},
1285  {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},
1286  {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},
1287  {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},
1288  {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},
1289  {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.3506, 1.3765, 1.4308, 1.4886, 1.5504},
1290  {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.0528, 2.0864, 2.1569, 2.2319, 2.3120},
1291  {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.6452, 2.6847, 2.7674, 2.8554, 2.9494},
1292  {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.0937, 3.1378, 3.2301, 3.3285, 3.4337},
1293  {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.4691, 3.5171, 3.6175, 3.7246, 3.8391},
1294  {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.7848, 3.8360, 3.9433, 4.0578, 4.1804},
1295  {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.0764, 4.1304, 4.2438, 4.3648, 4.4944},
1296  {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.5962, 4.6551, 4.7786, 4.9106, 5.0520},
1297  {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.9685, 5.0314, 5.1633, 5.3043, 5.4555},
1298  {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.3114, 5.3776, 5.5166, 5.6653, 5.8249},
1299  {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.4640, 5.5317, 5.6739, 5.8260, 5.9893},
1300  {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.7373, 5.8088, 5.9568, 6.1152, 6.2853},
1301  {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.0896, 6.1636, 6.3192, 6.4857, 6.6647},
1302  {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.9351, 7.0180, 7.1925, 7.3795, 7.5808},
1303  {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.9146, 8.0087, 8.2071, 8.4200, 8.6496},
1304  {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133}
1305  };
1306 
1307  const G4double bls3[28][9] = {
1308  {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},
1309  {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},
1310  {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},
1311  {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},
1312  {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},
1313  {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},
1314  {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},
1315  {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},
1316  {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},
1317  {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},
1318  {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},
1319  {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.3424, 1.4163, 1.4974, 1.5868},
1320  {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.8484, 1.9394, 2.0390, 2.1489},
1321  {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.6980, 2.8159, 2.9451, 3.0876},
1322  {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.4030, 3.5416, 3.6938, 3.8620},
1323  {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.9421, 4.0978, 4.2688, 4.4580},
1324  {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.3935, 4.5636, 4.7506, 4.9576},
1325  {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.7644, 4.9470, 5.1478, 5.3703},
1326  {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.1233, 5.3169, 5.5300, 5.7661},
1327  {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.7400, 5.9523, 6.1863, 6.4458},
1328  {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.1927, 6.4206, 6.6719, 6.9510},
1329  {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.6038, 6.8451, 7.1113, 7.4071},
1330  {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.7871, 7.0344, 7.3073, 7.6107},
1331  {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.1175, 7.3757, 7.6609, 7.9782},
1332  {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.5418, 7.8143, 8.1156, 8.4510},
1333  {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.5708, 8.8796, 9.2214, 9.6027},
1334  {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.7834, 10.1384, 10.5323, 10.9722},
1335  {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356}
1336  };
1337 
1338  const G4double bll1[28][10] = {
1339  {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},
1340  {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},
1341  {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},
1342  {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},
1343  {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},
1344  {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},
1345  {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},
1346  {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},
1347  {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},
1348  {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},
1349  {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},
1350  {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.1073, 1.1443, 1.1832, 1.2035, 1.2243},
1351  {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.5624, 1.6083, 1.6566, 1.6817, 1.7076},
1352  {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.3322, 2.3915, 2.4542, 2.4869, 2.5205},
1353  {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.9377, 3.0072, 3.0807, 3.1192, 3.1587},
1354  {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.4122, 3.4898, 3.5721, 3.6151, 3.6595},
1355  {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.8303, 3.9146, 4.0042, 4.0511, 4.0995},
1356  {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.1520, 4.2421, 4.3380, 4.3882, 4.4401},
1357  {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.4333, 4.5285, 4.6298, 4.6830, 4.7380},
1358  {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.9106, 5.0144, 5.1250, 5.1831, 5.2432},
1359  {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.2729, 5.3837, 5.5050, 5.5642, 5.6287},
1360  {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.5782, 5.6950, 5.8198, 5.8855, 5.9554},
1361  {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.7128, 5.8323, 5.9601, 6.0274, 6.0972},
1362  {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.9541, 6.0785, 6.2116, 6.2818, 6.3546},
1363  {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.2618, 6.3925, 6.5327, 6.6066, 6.6833},
1364  {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.0062, 7.1529, 7.3107, 7.3941, 7.4807},
1365  {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.8982, 8.0653, 8.2454, 8.3409, 8.4402},
1366  {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927}
1367  };
1368 
1369  const G4double bll2[28][10] = {
1370  {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},
1371  {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},
1372  {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},
1373  {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},
1374  {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},
1375  {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},
1376  {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},
1377  {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},
1378  {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},
1379  {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},
1380  {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},
1381  {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.4421, 1.4703, 1.5300, 1.5942, 1.6636},
1382  {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.9788, 2.0142, 2.0889, 2.1694, 2.2567},
1383  {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.8762, 2.9231, 3.0222, 3.1295, 3.2463},
1384  {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.5803, 3.6361, 3.7546, 3.8835, 4.0242},
1385  {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.1350, 4.1983, 4.3330, 4.4799, 4.6408},
1386  {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.6206, 4.6904, 4.8390, 5.0013, 5.1796},
1387  {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.0009, 5.0762, 5.2370, 5.4129, 5.6066},
1388  {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.3338, 5.4141, 5.5857, 5.7738, 5.9811},
1389  {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.8989, 5.9878, 6.1780, 6.3870, 6.6179},
1390  {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.3343, 6.4303, 6.6361, 6.8626, 7.1137},
1391  {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.7021, 6.8043, 7.0237, 7.2655, 7.5338},
1392  {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.8650, 6.9700, 7.1954, 7.4442, 7.7203},
1393  {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.1578, 7.2679, 7.5045, 7.7660, 8.0565},
1394  {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.5328, 7.6496, 7.9010, 8.1791, 8.4886},
1395  {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.4475, 8.5814, 8.8702, 9.1908, 9.5488},
1396  {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.5572, 9.7132, 10.0504, 10.4259, 10.8466},
1397  {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250}
1398  };
1399 
1400  const G4double bll3[28][9] = {
1401  {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},
1402  {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},
1403  {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},
1404  {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},
1405  {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},
1406  {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},
1407  {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},
1408  {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},
1409  {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},
1410  {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7836E-1, 9.3355E-1, 9.9462E-1, 1.0627, 1.1394},
1411  {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.2503, 1.3265, 1.4116, 1.5076},
1412  {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.0108, 2.1217, 2.2462, 2.3876},
1413  {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.6971, 2.8391, 2.9994, 3.1822},
1414  {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.8445, 4.0404, 4.2631, 4.5193},
1415  {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.7530, 4.9944, 5.2703, 5.5895},
1416  {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.4811, 5.7618, 6.0840, 6.4583},
1417  {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.1173, 6.4326, 6.7958, 7.2191},
1418  {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.6306, 6.9769, 7.3767, 7.8440},
1419  {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.0827, 7.4570, 7.8900, 8.3972},
1420  {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.8546, 8.2778, 8.7690, 9.3464},
1421  {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.4652, 8.9302, 9.4713, 10.1090},
1422  {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.9851, 9.4866, 10.0713, 10.7619},
1423  {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.2181, 9.7352, 10.3399, 11.0546},
1424  {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.6367, 10.1856, 10.8270, 11.5863},
1425  {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.1793, 10.7688, 11.4590, 12.2775},
1426  {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11.5229, 12.2172, 13.0332, 14.0048},
1427  {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 13.1898, 14.0213, 15.0024, 16.1752},
1428  {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420}
1429  };
1430 
1431  G4double b, bs;
1432  for(i=0; i<nEtaK; ++i) {
1433 
1434  G4double et = eta[i];
1435  G4double loget = G4Log(et);
1436  Eta[i] = et;
1437  //G4cout << "### eta["<<i<<"]="<<et<<" KShell: tet= "<<TheK[0]<<" - "<<TheK[nK-1]<<G4endl;
1438 
1439  for(j=0; j<nK; ++j) {
1440 
1441  if(j < 10) { b = bk2[i][10-j]; }
1442  else { b = bk1[i][20-j]; }
1443 
1444  CK[j][i] = SK[j]*loget + TK[j] - b;
1445 
1446  if(i == nEtaK-1) {
1447  ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]);
1448  //G4cout << "i= " << i << " j= " << j
1449  // << " CK[j][i]= " << CK[j][i]
1450  // << " ZK[j]= " << ZK[j] << " b= " << b << G4endl;
1451  }
1452  }
1453  //G4cout << G4endl;
1454  if(i < nEtaL) {
1455  //G4cout << " LShell:" <<G4endl;
1456  for(j=0; j<nL; ++j) {
1457 
1458  if(j < 8) {
1459  bs = bls3[i][8-j];
1460  b = bll3[i][8-j];
1461  } else if(j < 17) {
1462  bs = bls2[i][17-j];
1463  b = bll2[i][17-j];
1464  } else {
1465  bs = bls1[i][26-j];
1466  b = bll1[i][26-j];
1467  }
1468  G4double c = SL[j]*loget + TL[j];
1469  CL[j][i] = c - bs - 3.0*b;
1470  if(i == nEtaL-1) {
1471  VL[j] = et*(et*CL[j][i] - UL[j]);
1472  //G4cout << "i= " << i << " j= " << j
1473  // << " CL[j][i]= " << CL[j][i]
1474  // << " VL[j]= " << VL[j] << " b= " << b << " bs= " << bs
1475  // << " et= " << et << G4endl;
1476  //" UL= " << UL[j] << " TL= " << TL[j] << " SL= " << SL[j] <<G4endl;
1477  }
1478  }
1479  //G4cout << G4endl;
1480  }
1481  }
1482  const G4double hm[53] = {
1483  12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4,
1484  10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87,
1485  5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38,
1486  4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 3.95, 3.93, 3.91,
1487  3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89,
1488  3.90, 3.92, 3.93 };
1489  const G4double hn[31] = {
1490  75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8,
1491  24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7,
1492  19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4,
1493  18.2
1494  };
1495  for(i=0; i<53; ++i) {HM[i] = hm[i];}
1496  for(i=0; i<31; ++i) {HN[i] = hn[i];}
1497 
1498  const G4double xzk[34] = { 11.7711,
1499  13.3669, 15.5762, 17.1715, 18.7667, 20.8523, 23.0606, 24.901, 26.9861, 29.4394, 31.77,
1500  34.3457, 37.4119, 40.3555, 42.3177, 44.7705, 47.2234, 50.78, 53.8458, 56.4214, 58.3834,
1501  60.9586, 63.6567, 66.5998, 68.807, 71.8728, 74.5706, 77.3911, 81.8056, 85.7297, 89.8988,
1502  93.4549, 96.2753, 99.709};
1503  const G4double yzk[34] = { 0.70663,
1504  0.72033, 0.73651, 0.74647, 0.75518, 0.76388, 0.77258, 0.78129, 0.78625, 0.7937, 0.79991,
1505  0.80611, 0.8123, 0.8185, 0.82097, 0.82467, 0.82838, 0.83457, 0.83702, 0.84198, 0.8432,
1506  0.84565, 0.84936, 0.85181, 0.85303, 0.85548, 0.85794, 0.8604, 0.86283, 0.86527, 0.86646,
1507  0.86891, 0.87011, 0.87381};
1508 
1509  const G4double xzl[36] = { 15.5102,
1510  16.7347, 17.9592, 19.551, 21.0204, 22.6122, 24.9388, 27.3878, 29.5918, 31.3061, 32.898,
1511  34.4898, 36.2041, 38.4082, 40.3674, 42.5714, 44.898, 47.4694, 49.9184, 52.7347, 55.9184,
1512  59.3469, 61.9184, 64.6122, 67.4286, 71.4694, 75.2653, 78.3265, 81.2653, 85.551, 88.7347,
1513  91.551, 94.2449, 96.449, 98.4082, 99.7551};
1514  const G4double yzl[36] = { 0.29875,
1515  0.31746, 0.33368, 0.35239, 0.36985, 0.38732, 0.41102, 0.43472, 0.45343, 0.4659, 0.47713,
1516  0.4896, 0.50083, 0.51331, 0.52328, 0.53077, 0.54075, 0.54823, 0.55572, 0.56445, 0.57193,
1517  0.58191, 0.5869, 0.59189, 0.60062, 0.60686, 0.61435, 0.61809, 0.62183, 0.62931, 0.6343,
1518  0.6368, 0.64054, 0.64304, 0.64428, 0.64678};
1519 
1520  ThetaK = new G4LPhysicsFreeVector(34, xzk[0], xzk[33]);
1521  ThetaL = new G4LPhysicsFreeVector(36, xzl[0], xzl[35]);
1522  for(i=0; i<34; ++i) { ThetaK->PutValues(i, xzk[i], yzk[i]); }
1523  for(i=0; i<36; ++i) { ThetaL->PutValues(i, xzl[i], yzl[i]); }
1524  ThetaK->SetSpline(true);
1525  ThetaL->SetSpline(true);
1526 
1527  const G4double coseb[14] = {0.0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
1528  1.0,1.2,1.5,2.0};
1529  const G4double cosxi[14] = {1.0000, 0.9905, 0.9631, 0.9208, 0.8680,
1530  0.7478, 0.6303, 0.5290, 0.4471, 0.3323,
1531  0.2610, 0.2145, 0.1696, 0.1261};
1532  for(i=0; i<14; ++i) {
1533  COSEB[i] = coseb[i];
1534  COSXI[i] = cosxi[i];
1535  }
1536 
1537  G4Pow* g4pow = G4Pow::GetInstance();
1538  for(i=1; i<100; ++i) {
1539  Z23[i] = G4Exp(g4pow->logZ(i)*0.23);
1540  }
1541 }
1542 
1543 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1544 
G4Exp.hh
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4EmCorrections::COSXI
G4double COSXI[14]
Definition: G4EmCorrections.hh:215
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:68
G4Material::GetIonisation
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:226
G4EmCorrections::CK
G4double CK[20][29]
Definition: G4EmCorrections.hh:232
G4LPhysicsFreeVector::PutValues
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
G4EmCorrections::HighOrderCorrections
G4double HighOrderCorrections(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy, G4double cutEnergy)
Definition: G4EmCorrections.cc:119
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4GenericIon.hh
G4IonisParamMat::GetAdensity
G4double GetAdensity() const
Definition: G4IonisParamMat.hh:93
G4EmCorrections::ionLEModel
G4VEmModel * ionLEModel
Definition: G4EmCorrections.hh:276
cm2
static const double cm2
Definition: G4SIunits.hh:119
G4EmCorrections::DensityCorrection
G4double DensityCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:556
G4EmCorrections::ed
G4double ed[104]
Definition: G4EmCorrections.hh:201
G4EmCorrections::theZieglerFactor
G4double theZieglerFactor
Definition: G4EmCorrections.hh:203
G4EmCorrections::currmat
std::vector< const G4Material * > currmat
Definition: G4EmCorrections.hh:242
G4Pow
Definition: G4Pow.hh:56
G4EmCorrections::KShell
G4double KShell(G4double theta, G4double eta)
Definition: G4EmCorrections.cc:347
G4Element::GetN
G4double GetN() const
Definition: G4Element.hh:134
G4LPhysicsFreeVector.hh
G4MaterialCutsCouple.hh
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4EmCorrections::ThetaK
G4LPhysicsFreeVector * ThetaK
Definition: G4EmCorrections.hh:239
G4IonisParamMat::GetX1density
G4double GetX1density() const
Definition: G4IonisParamMat.hh:97
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4NistManager::GetZ13
G4double GetZ13(G4double Z)
Definition: G4NistManager.hh:548
G4VEmModel.hh
G4ionEffectiveCharge::EffectiveChargeSquareRatio
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *p, const G4Material *material, G4double kineticEnergy)
Definition: G4ionEffectiveCharge.hh:106
G4EmCorrections::material
const G4Material * material
Definition: G4EmCorrections.hh:248
G4EmCorrections::MottCorrection
G4double MottCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:666
G4IonisParamMat::GetTaul
G4double GetTaul() const
Definition: G4IonisParamMat.hh:81
G4EmCorrections::EffectiveChargeCorrection
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:748
G4PhysicsLogVector.hh
G4EmCorrections::Aion
std::vector< G4int > Aion
Definition: G4EmCorrections.hh:284
G4EmCorrections::nist
G4NistManager * nist
Definition: G4EmCorrections.hh:274
G4EmCorrections::BlochCorrection
G4double BlochCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:639
G4EmCorrections.hh
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4EmCorrections::SpinCorrection
G4double SpinCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:269
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4EmCorrections::ThetaL
G4LPhysicsFreeVector * ThetaL
Definition: G4EmCorrections.hh:240
bg2lim
static const G4double bg2lim
Definition: G4UrbanMscModel.cc:86
G4EmCorrections::ionHEModel
G4VEmModel * ionHEModel
Definition: G4EmCorrections.hh:277
G4IonTable.hh
G4ParticleDefinition::GetAtomicNumber
G4int GetAtomicNumber() const
G4EmCorrections::AddStoppingData
void AddStoppingData(G4int Z, G4int A, const G4String &materialName, G4PhysicsVector *dVector)
Definition: G4EmCorrections.cc:797
G4EmCorrections::TheL
G4double TheL[26]
Definition: G4EmCorrections.hh:225
G4PhysicsVector::SetSpline
void SetSpline(G4bool)
G4EmCorrections::G4EmCorrections
G4EmCorrections(G4int verb)
Definition: G4EmCorrections.cc:83
G4EmCorrections::thcorr
std::map< G4int, std::vector< G4double > > thcorr
Definition: G4EmCorrections.hh:243
G4EmCorrections::Value2
G4double Value2(G4double xv, G4double yv, G4double x1, G4double x2, G4double y1, G4double y2, G4double z11, G4double z21, G4double z12, G4double z22)
Definition: G4EmCorrections.hh:308
G4Pow::logZ
G4double logZ(G4int Z) const
Definition: G4Pow.hh:166
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:658
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:59
G4EmCorrections::NuclearStoppingPower
G4double NuclearStoppingPower(G4double e, G4double z1, G4double z2, G4double m1, G4double m2)
Definition: G4EmCorrections.cc:706
G4ParticleTable.hh
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4EmCorrections::Index
G4int Index(G4double x, G4double *y, G4int n)
Definition: G4EmCorrections.hh:294
G4EmCorrections::NuclearDEDX
G4double NuclearDEDX(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy, G4bool fluct=true)
Definition: G4EmCorrections.cc:677
G4bool
bool G4bool
Definition: G4Types.hh:79
G4EmCorrections::IonBarkasCorrection
G4double IonBarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:154
iz
G4double iz
Definition: TRTMaterials.hh:39
G4Proton.hh
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4EmCorrections::materialName
std::vector< G4String > materialName
Definition: G4EmCorrections.hh:285
G4EmCorrections::BarkasCorrection
G4double BarkasCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:582
G4EmCorrections::IonHighOrderCorrections
G4double IonHighOrderCorrections(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:206
nmax
const G4int nmax
Definition: G4UrQMD1_3Interface.hh:69
Randomize.hh
G4Log.hh
G4EmCorrections::particle
const G4ParticleDefinition * particle
Definition: G4EmCorrections.hh:246
G4PhysicsVector::Energy
G4double Energy(size_t index) const
G4PhysicsVector::Value
G4double Value(G4double theEnergy, size_t &lastidx) const
Definition: G4PhysicsVector.cc:498
G4EmCorrections::curParticle
const G4ParticleDefinition * curParticle
Definition: G4EmCorrections.hh:247
G4PhysicalConstants.hh
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4EmCorrections::materialList
std::vector< const G4Material * > materialList
Definition: G4EmCorrections.hh:289
G4Material::GetTotNbOfAtomsPerVolume
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4EmCorrections::SetupKinematics
void SetupKinematics(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.hh:348
G4EmCorrections::ShellCorrection
G4double ShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:467
G4EmCorrections::atomDensity
const G4double * atomDensity
Definition: G4EmCorrections.hh:251
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
eV
static const double eV
Definition: G4SIunits.hh:212
G4GenericIon::GenericIon
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
G4VEmModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.cc:249
factor
static const G4double factor
Definition: G4ContinuumGammaTransition.cc:63
G4IonisParamMat::GetX0density
G4double GetX0density() const
Definition: G4IonisParamMat.hh:95
G4ProductionCutsTable.hh
G4EmCorrections::TheK
G4double TheK[20]
Definition: G4EmCorrections.hh:218
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
pi
static const double pi
Definition: G4SIunits.hh:74
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4EmCorrections::BuildCorrectionVector
void BuildCorrectionVector()
Definition: G4EmCorrections.cc:822
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
G4IonisParamMat::GetCdensity
G4double GetCdensity() const
Definition: G4IonisParamMat.hh:89
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4IonisParamMat::GetShellCorrectionVector
G4double * GetShellCorrectionVector() const
Definition: G4IonisParamMat.hh:79
G4EmCorrections::ionList
std::vector< const G4ParticleDefinition * > ionList
Definition: G4EmCorrections.hh:287
x
const G4double x[NPOINTSGL]
Definition: G4PreCompoundFragment.cc:56
G4EmCorrections::inveplus
static const G4double inveplus
Definition: G4EmCorrections.hh:199
G4AtomicShells::GetNumberOfElectrons
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
Definition: G4AtomicShells.cc:768
G4EmCorrections::COSEB
G4double COSEB[14]
Definition: G4EmCorrections.hh:214
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:117
G4EmCorrections::ionTable
G4IonTable * ionTable
Definition: G4EmCorrections.hh:275
G4EmCorrections::KShellCorrection
G4double KShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:280
G4EmCorrections::LShell
G4double LShell(G4double theta, G4double eta)
Definition: G4EmCorrections.cc:390
G4EmCorrections::LShellCorrection
G4double LShellCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:309
G4EmCorrections::a
G4double a[104]
Definition: G4EmCorrections.hh:202
G4EmCorrections::stopData
std::vector< G4PhysicsVector * > stopData
Definition: G4EmCorrections.hh:290
G4IonisParamMat::GetMeanExcitationEnergy
G4double GetMeanExcitationEnergy() const
Definition: G4IonisParamMat.hh:71
G4EmCorrections::BarkasCorr
G4LPhysicsFreeVector * BarkasCorr
Definition: G4EmCorrections.hh:238
twoln10
static const G4double twoln10
Definition: G4EmCorrections.cc:81
G4EmCorrections::Z23
G4double Z23[100]
Definition: G4EmCorrections.hh:236
G4endl
#define G4endl
Definition: G4ios.hh:61
G4EmCorrections::effCharge
G4ionEffectiveCharge effCharge
Definition: G4EmCorrections.hh:272
keV
static const double keV
Definition: G4SIunits.hh:213
G4AtomicShells.hh
G4EmCorrections::ComputeIonCorrections
G4double ComputeIonCorrections(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:174
G4double
double G4double
Definition: G4Types.hh:76
G4EmCorrections::ShellCorrectionSTD
G4double ShellCorrectionSTD(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:433
G4SystemOfUnits.hh
G4EmCorrections::curMaterial
const G4Material * curMaterial
Definition: G4EmCorrections.hh:249
eplus
static const double eplus
Definition: G4SIunits.hh:196
G4EmCorrections::curVector
G4PhysicsVector * curVector
Definition: G4EmCorrections.hh:291
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4Pow.hh
G4IonisParamMat::GetMdensity
G4double GetMdensity() const
Definition: G4IonisParamMat.hh:91
G4EmCorrections::Value
G4double Value(G4double xv, G4double x1, G4double x2, G4double y1, G4double y2)
Definition: G4EmCorrections.hh:302
G4EmCorrections::~G4EmCorrections
virtual ~G4EmCorrections()
Definition: G4EmCorrections.cc:109
G4EmCorrections::Bethe
G4double Bethe(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:256
G4EmCorrections::CL
G4double CL[26][28]
Definition: G4EmCorrections.hh:233
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4EmCorrections::Zion
std::vector< G4int > Zion
Definition: G4EmCorrections.hh:283
G4EmCorrections::Eta
G4double Eta[29]
Definition: G4EmCorrections.hh:231
G4AtomicShells::GetNumberOfShells
static G4int GetNumberOfShells(G4int Z)
Definition: G4AtomicShells.cc:750