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