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