Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmCorrections.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/utils/src/G4EmCorrections.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4EmCorrections.cc (Version 8.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
                                                   >>  23 // $Id: G4EmCorrections.cc,v 1.13 2005/11/26 16:59:26 vnivanch Exp $
                                                   >>  24 // GEANT4 tag $Name: geant4-08-00-patch-01 $
 26 //                                                 25 //
 27 // -------------------------------------------     26 // -------------------------------------------------------------------
 28 //                                                 27 //
 29 // GEANT4 Class                                    28 // GEANT4 Class 
 30 //                                                 29 //
 31 // File name:     G4EmCorrections                  30 // File name:     G4EmCorrections
 32 //                                                 31 //
 33 // Author:        Vladimir Ivanchenko              32 // Author:        Vladimir Ivanchenko
 34 //                                                 33 //
 35 // Creation date: 13.01.2005                       34 // Creation date: 13.01.2005
 36 //                                                 35 //
 37 // Modifications:                                  36 // Modifications:
 38 // 05.05.2005 V.Ivanchenko Fix misprint in Mot     37 // 05.05.2005 V.Ivanchenko Fix misprint in Mott term
 39 // 26.11.2005 V.Ivanchenko Fix effective charg     38 // 26.11.2005 V.Ivanchenko Fix effective charge for heavy ions using original paper
 40 // 28.04.2006 V.Ivanchenko General cleanup, ad << 
 41 // 13.05.2006 V.Ivanchenko Add corrections for << 
 42 // 08.05.2007 V.Ivanchenko Use G4IonTable for  << 
 43 //                         division by zero    << 
 44 // 29.02.2008 V.Ivanchenko use expantions for  << 
 45 // 21.04.2008 Updated computations for ions (V << 
 46 // 20.05.2008 Removed Finite Size correction ( << 
 47 // 19.04.2012 Fix reproducibility problem (A.R << 
 48 //                                                 39 //
 49 //                                                 40 //
 50 // Class Description:                              41 // Class Description:
 51 //                                                 42 //
 52 // This class provides calculation of EM corre <<  43 // This class provides calculation of EM corrections to ionisation 
 53 //                                                 44 //
 54                                                    45 
 55 // -------------------------------------------     46 // -------------------------------------------------------------------
 56 //                                                 47 //
 57                                                    48 
 58 #include "G4EmCorrections.hh"                      49 #include "G4EmCorrections.hh"
 59 #include "G4PhysicalConstants.hh"              <<  50 #include "G4Material.hh"
 60 #include "G4SystemOfUnits.hh"                  <<  51 #include "G4ParticleDefinition.hh"
 61 #include "G4ParticleTable.hh"                  <<  52 #include "Randomize.hh"
 62 #include "G4IonTable.hh"                       << 
 63 #include "G4VEmModel.hh"                       << 
 64 #include "G4Proton.hh"                         << 
 65 #include "G4GenericIon.hh"                     << 
 66 #include "G4PhysicsLogVector.hh"               << 
 67 #include "G4ProductionCutsTable.hh"            << 
 68 #include "G4MaterialCutsCouple.hh"             << 
 69 #include "G4AtomicShells.hh"                   << 
 70 #include "G4Log.hh"                            << 
 71 #include "G4Exp.hh"                            << 
 72 #include "G4Pow.hh"                            << 
 73                                                    53 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    55 
 76 namespace                                      <<  56 G4EmCorrections::G4EmCorrections()
 77 {                                                  57 {
 78   constexpr G4double inveplus = 1.0/CLHEP::epl <<  58   Initialise();
 79   constexpr G4double alpha2 = CLHEP::fine_stru << 
 80 }                                              << 
 81 const G4double G4EmCorrections::ZD[11] =       << 
 82     {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, << 
 83 const G4double G4EmCorrections::UK[20] = {1.99 << 
 84                            2.0817, 2.0945, 2.0 << 
 85                            2.1197, 2.1246, 2.1 << 
 86                            2.1310, 2.1310, 2.1 << 
 87 const G4double G4EmCorrections::VK[20] = {8.34 << 
 88                            8.3219, 8.3201, 8.3 << 
 89                            8.3191, 8.3199, 8.3 << 
 90                            8.3244, 8.3264, 8.3 << 
 91 G4double G4EmCorrections::ZK[] = {0.0};        << 
 92 const G4double G4EmCorrections::Eta[29] = {0.0 << 
 93                             0.03,  0.04,  0.05 << 
 94                             0.1,   0.15,  0.2, << 
 95                             0.5,   0.6,   0.7, << 
 96                             1.2,   1.4,   1.5, << 
 97 G4double G4EmCorrections::CK[20][29];          << 
 98 G4double G4EmCorrections::CL[26][28];          << 
 99 const G4double G4EmCorrections::UL[] = {0.1215 << 
100                            1.4379, 1.5032, 1.5 << 
101                            1.8036, 1.8543, 1.8 << 
102                            1.9508, 1.9696, 1.9 << 
103                            2.0001, 2.0039, 2.0 << 
104 G4double G4EmCorrections::VL[] = {0.0};        << 
105 G4double G4EmCorrections::sWmaxBarkas = 10.0;  << 
106                                                << 
107 G4PhysicsFreeVector* G4EmCorrections::sBarkasC << 
108 G4PhysicsFreeVector* G4EmCorrections::sThetaK  << 
109 G4PhysicsFreeVector* G4EmCorrections::sThetaL  << 
110                                                << 
111 G4EmCorrections::G4EmCorrections(G4int verb)   << 
112   : verbose(verb)                              << 
113 {                                              << 
114   eth = 2.0*CLHEP::MeV;                        << 
115   eCorrMin = 25.*CLHEP::keV;                   << 
116   eCorrMax = 1.*CLHEP::GeV;                    << 
117                                                << 
118   ionTable = G4ParticleTable::GetParticleTable << 
119   g4calc = G4Pow::GetInstance();               << 
120                                                << 
121   // fill vectors                              << 
122   if (nullptr == sBarkasCorr) {                << 
123     Initialise();                              << 
124     isInitializer = true;                      << 
125   }                                            << 
126 }                                                  59 }
127                                                    60 
128 //....oooOO0OOooo........oooOO0OOooo........oo     61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129                                                    62 
130 G4EmCorrections::~G4EmCorrections()                63 G4EmCorrections::~G4EmCorrections()
131 {                                              <<  64 {}
132   for (G4int i=0; i<nIons; ++i) { delete stopD << 
133   if (isInitializer) {                         << 
134     delete sBarkasCorr;                        << 
135     delete sThetaK;                            << 
136     delete sThetaL;                            << 
137     sBarkasCorr = sThetaK = sThetaL = nullptr; << 
138   }                                            << 
139 }                                              << 
140                                                << 
141 void G4EmCorrections::SetupKinematics(const G4 << 
142               const G4Material* mat,           << 
143               const G4double kineticEnergy)    << 
144 {                                              << 
145   if(kineticEnergy != kinEnergy || p != partic << 
146     particle = p;                              << 
147     kinEnergy = kineticEnergy;                 << 
148     mass  = p->GetPDGMass();                   << 
149     tau   = kineticEnergy / mass;              << 
150     gamma = 1.0 + tau;                         << 
151     bg2   = tau * (tau+2.0);                   << 
152     beta2 = bg2/(gamma*gamma);                 << 
153     beta  = std::sqrt(beta2);                  << 
154     ba2   = beta2/alpha2;                      << 
155     const G4double ratio = CLHEP::electron_mas << 
156     tmax  = 2.0*CLHEP::electron_mass_c2*bg2    << 
157       /(1. + 2.0*gamma*ratio + ratio*ratio);   << 
158     charge  = p->GetPDGCharge()*inveplus;      << 
159     if(charge > 1.5) { charge = effCharge.Effe << 
160     q2 = charge*charge;                        << 
161   }                                            << 
162   if(mat != material) {                        << 
163     material = mat;                            << 
164     theElementVector = material->GetElementVec << 
165     atomDensity  = material->GetAtomicNumDensi << 
166     numberOfElements = (G4int)material->GetNum << 
167   }                                            << 
168 }                                              << 
169                                                    65 
170 //....oooOO0OOooo........oooOO0OOooo........oo     66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171                                                    67 
172 G4double G4EmCorrections::HighOrderCorrections     68 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p,
173                                                <<  69                                                const G4Material* material,
174                                                <<  70                                                      G4double kineticEnergy)
175 {                                                  71 {
176   // . Z^3 Barkas effect in the stopping power <<  72 // . Z^3 Barkas effect in the stopping power of matter for charged particles
177   //   J.C Ashley and R.H.Ritchie              <<  73 //   J.C Ashley and R.H.Ritchie
178   //   Physical review B Vol.5 No.7 1 April 19 <<  74 //   Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
179   //   and ICRU49 report                       <<  75 //   and ICRU49 report
180   //   valid for kineticEnergy < 0.5 MeV       <<  76 //   valid for kineticEnergy < 0.5 MeV
181   //   Other corrections from S.P.Ahlen Rev. M <<  77 //   Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
182                                                <<  78 
183   SetupKinematics(p, mat, e);                  <<  79   G4double tau   = kineticEnergy / p->GetPDGMass();
184   if(tau <= 0.0) { return 0.0; }               <<  80   if(tau <= 0.0) return 0.0;
185                                                <<  81 
186   const G4double Barkas = BarkasCorrection(p,  <<  82   G4double gamma = 1.0 + tau;
187   const G4double Bloch  = BlochCorrection(p, m <<  83   G4double bg2   = tau * (tau+2.0);
188   const G4double Mott = MottCorrection(p, mat, <<  84   G4double beta2 = bg2/(gamma*gamma);
189                                                <<  85   G4double beta  = std::sqrt(beta2);
190   G4double sum = 2.0*(Barkas + Bloch) + Mott;  <<  86 
191                                                <<  87   G4double q  = p->GetPDGCharge()/eplus;
192   if(verbose > 1) {                            <<  88   if(q > 2.5)      q *= (1.0 - std::exp(-130.0*beta/std::pow(q,0.66666667)));
193     G4cout << "EmCorrections: E(MeV)= " << e/M <<  89   else if(q > 1.5) q  = effCharge.EffectiveCharge(p,material,kineticEnergy)/eplus;
194            << " Bloch= " << Bloch << " Mott= " << 
195            << " Sum= " << sum << " q2= " << q2 << 
196     G4cout << " ShellCorrection: " << ShellCor << 
197            << " Kshell= " << KShellCorrection( << 
198            << " Lshell= " << LShellCorrection( << 
199            << "   " << mat->GetName() << G4end << 
200   }                                            << 
201   sum *= material->GetElectronDensity()*q2*CLH << 
202   return sum;                                  << 
203 }                                              << 
204                                                    90 
205 //....oooOO0OOooo........oooOO0OOooo........oo <<  91   G4double q2 = q*q;
                                                   >>  92   G4double ba = beta2/alpha2;
                                                   >>  93   G4double BarkasTerm = 0.0;
                                                   >>  94   const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >>  95   const G4double* atomDensity  = material->GetAtomicNumDensityVector();
                                                   >>  96   G4int numberOfElements = material->GetNumberOfElements();
206                                                    97 
207 G4double G4EmCorrections::IonBarkasCorrection( <<  98   G4double Zeff = 0.0;
208                                                <<  99   G4double norm = 0.0;
209                                                << 
210 {                                              << 
211   SetupKinematics(p, mat, e);                  << 
212   return 2.0*BarkasCorrection(p, mat, e, true) << 
213     material->GetElectronDensity() * q2 * CLHE << 
214 }                                              << 
215                                                   100 
216 //....oooOO0OOooo........oooOO0OOooo........oo << 101   for (G4int i = 0; i<numberOfElements; i++) {
217                                                   102 
218 G4double G4EmCorrections::ComputeIonCorrection << 103     G4double Z = (*theElementVector)[i]->GetZ();
219                                                << 104     G4int   iz = G4int(Z);
220                                                << 105     G4double X = ba / Z;
221 {                                              << 
222   // . Z^3 Barkas effect in the stopping power << 
223   //   J.C Ashley and R.H.Ritchie              << 
224   //   Physical review B Vol.5 No.7 1 April 19 << 
225   //   and ICRU49 report                       << 
226   //   valid for kineticEnergy < 0.5 MeV       << 
227   //   Other corrections from S.P.Ahlen Rev. M << 
228   SetupKinematics(p, mat, e);                  << 
229   if(tau <= 0.0) { return 0.0; }               << 
230                                                << 
231   const G4double Barkas = BarkasCorrection (p, << 
232   const G4double Bloch  = BlochCorrection (p,  << 
233   const G4double Mott = MottCorrection (p, mat << 
234                                                << 
235   G4double sum = 2.0*(Barkas*(charge - 1.0)/ch << 
236                                                << 
237   if(verbose > 1) {                            << 
238     G4cout << "EmCorrections: E(MeV)= " << e/M << 
239            << " Bloch= " << Bloch << " Mott= " << 
240            << " Sum= " << sum << G4endl;       << 
241   }                                            << 
242   sum *= material->GetElectronDensity() * q2 * << 
243                                                   106 
244   if(verbose > 1) { G4cout << " Sum= " << sum  << 107     // Variables to compute L_1
245   return sum;                                  << 108     //    G4double Eta0Chi = 0.8;
246 }                                              << 109     //    G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( Z,-1.19 ) );
                                                   >> 110     //    G4double W = ( EtaChi * std::pow( Z,0.16666667 ) ) / std::sqrt(X);
247                                                   111 
248 //....oooOO0OOooo........oooOO0OOooo........oo << 112     G4double b = 1.3;
249                                                << 113     if(1 == iz) {
250 G4double G4EmCorrections::IonHighOrderCorrecti << 114       if(material->GetName() == "G4_lH2") b = 0.6;
251                                                << 115       else                                b = 1.8;
252                                                << 
253 {                                              << 
254   // . Z^3 Barkas effect in the stopping power << 
255   //   J.C Ashley and R.H.Ritchie              << 
256   //   Physical review B Vol.5 No.7 1 April 19 << 
257   //   and ICRU49 report                       << 
258   //   valid for kineticEnergy < 0.5 MeV       << 
259   //   Other corrections from S.P.Ahlen Rev. M << 
260                                                << 
261   G4double sum = 0.0;                          << 
262                                                << 
263   if (nullptr != ionHEModel) {                 << 
264     G4int Z = G4lrint(p->GetPDGCharge()*invepl << 
265     Z = std::max(std::min(Z, 99), 1);          << 
266                                                << 
267     const G4double ethscaled = eth*p->GetPDGMa << 
268     const G4int ionPDG = p->GetPDGEncoding();  << 
269     auto iter = thcorr.find(ionPDG);           << 
270     if (iter == thcorr.end()) {  // Not found: << 
271       std::vector<G4double> v;                 << 
272       for(std::size_t i=0; i<ncouples; ++i){   << 
273         v.push_back(ethscaled*ComputeIonCorrec << 
274       }                                        << 
275       thcorr.insert(std::pair< G4int, std::vec << 
276     }                                             116     }
277     G4double rest = 0.0;                       << 117     else if(2 == iz)  b = 0.6;
278     iter = thcorr.find(ionPDG);                << 118     else if(10 >= iz) b = 1.8;
279     if (iter != thcorr.end()) { rest = (iter-> << 119     else if(17 >= iz) b = 1.4;
280                                                << 120     else if(18 == iz) b = 1.8;
281     sum = ComputeIonCorrections(p,couple->GetM << 121     else if(25 >= iz) b = 1.4;
282                                                << 122     else if(50 >= iz) b = 1.35;
283     if(verbose > 1) {                          << 123 
284       G4cout << " Sum= " << sum << " dSum= " < << 124     G4double W = b/std::sqrt(X);
285     }                                          << 125 
                                                   >> 126     G4double val;
                                                   >> 127 
                                                   >> 128     if(W <= engBarkas[0])       val =  corBarkas[0];
                                                   >> 129     else if(W >= engBarkas[46]) val =  corBarkas[46]*engBarkas[46]/W;
                                                   >> 130     else {
                                                   >> 131       G4int iw = Index(W, engBarkas, 47);
                                                   >> 132       val = Value(W, engBarkas[iw], engBarkas[iw+1], corBarkas[iw], corBarkas[iw+1]);
                                                   >> 133     }
                                                   >> 134     BarkasTerm += 1.29*val*atomDensity[i] * std::sqrt(Z /X)/ X;
                                                   >> 135     Zeff += Z*atomDensity[i];
                                                   >> 136     norm += atomDensity[i];
                                                   >> 137   }
                                                   >> 138   Zeff /= norm; 
                                                   >> 139 
                                                   >> 140 
                                                   >> 141   BarkasTerm *= 2.0*q;
                                                   >> 142 
                                                   >> 143   // L2 correction
                                                   >> 144   G4double y2 = q2 / ba ;
                                                   >> 145   G4double BlochTerm = 1.202;
                                                   >> 146 
                                                   >> 147   if(y2 > 0.05) {
                                                   >> 148     BlochTerm = 1.0 / (1.0 + y2) ;
                                                   >> 149     G4double de = BlochTerm;
                                                   >> 150     G4double i = 1.0;
                                                   >> 151     do {
                                                   >> 152       i += 1.0;
                                                   >> 153       de = 1.0/( i * (i*i + y2));
                                                   >> 154       BlochTerm += de;
                                                   >> 155     } while (de > BlochTerm*0.01);
                                                   >> 156 
                                                   >> 157   }
                                                   >> 158   BlochTerm *= -2.0*y2;
                                                   >> 159 
                                                   >> 160   // Estimation of mean square root of the ionisation potential 
                                                   >> 161   G4double ze  = 2.0*Zeff;
                                                   >> 162   G4double ze1 = std::log(ze);
                                                   >> 163   G4double eexc= material->GetIonisation()->GetMeanExcitationEnergy()*ze1*ze1/ze;
                                                   >> 164 
                                                   >> 165   G4double invbeta = 1.0/beta;
                                                   >> 166   G4double invbeta2= invbeta*invbeta;
                                                   >> 167   G4double za  = q*fine_structure_const;
                                                   >> 168   G4double za2 = za*za;
                                                   >> 169   G4double za3 = za2*za;
                                                   >> 170   G4double x   = za*invbeta;
                                                   >> 171   G4double cosx;
                                                   >> 172   if(x < COSEB[13]) {
                                                   >> 173     G4int i = Index(x,COSEB,14);
                                                   >> 174     cosx    = Value(x,COSEB[i], COSEB[i+1],COSXI[i],COSXI[i+1]);
                                                   >> 175   } else {
                                                   >> 176     cosx = COSXI[13]*COSEB[13]/x;
286   }                                               177   }
287   return sum;                                  << 178 
                                                   >> 179   G4double mterm =
                                                   >> 180         za*beta*(1.725 + pi*cosx*(0.52 - 2.0*std::sqrt(eexc/(2.0*electron_mass_c2*bg2))))
                                                   >> 181       + za2*(3.246 - 0.451*beta2)
                                                   >> 182       + za3*(1.522*beta + 0.987*invbeta)
                                                   >> 183       + za2*za2*(4.569 - 0.494*beta2 - 2.696*invbeta2)
                                                   >> 184       + za3*za2*(1.254*beta + 0.222*invbeta - 1.17*invbeta*invbeta2);
                                                   >> 185 
                                                   >> 186   G4double eloss = (BarkasTerm + (BlochTerm + mterm)*material->GetElectronDensity()) * 
                                                   >> 187     q2 *  twopi_mc2_rcl2 *invbeta2;
                                                   >> 188   return eloss;
288 }                                                 189 }
289                                                   190 
290 //....oooOO0OOooo........oooOO0OOooo........oo    191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291                                                   192 
292 G4double G4EmCorrections::Bethe(const G4Partic    193 G4double G4EmCorrections::Bethe(const G4ParticleDefinition* p,
293                                 const G4Materi << 194                                 const G4Material* material,
294                                 const G4double << 195                                       G4double kineticEnergy)
295 {                                                 196 {
296   SetupKinematics(p, mat, e);                  << 197   G4double mass  = p->GetPDGMass();
297   const G4double eexc  = material->GetIonisati << 198   G4double ratio = electron_mass_c2/mass;
298   const G4double eexc2 = eexc*eexc;            << 199   G4double tau   = kineticEnergy / mass;
299   return 0.5*G4Log(2.0*electron_mass_c2*bg2*tm << 200   G4double gamma = 1.0 + tau;
                                                   >> 201   G4double bg2   = tau * (tau+2.0);
                                                   >> 202   G4double beta2 = bg2/(gamma*gamma);
                                                   >> 203   G4double tmax  = 2.0*electron_mass_c2*bg2 /(1. + 2.0*gamma*ratio + ratio*ratio);
                                                   >> 204   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
                                                   >> 205   G4double eexc2 = eexc*eexc;
                                                   >> 206   G4double dedx = 0.5*std::log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2;
                                                   >> 207   return dedx;
300 }                                                 208 }
301                                                   209 
302 //....oooOO0OOooo........oooOO0OOooo........oo    210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303                                                   211 
304 G4double G4EmCorrections::SpinCorrection(const    212 G4double G4EmCorrections::SpinCorrection(const G4ParticleDefinition* p,
305                                          const << 213                                          const G4Material*,
306                                          const << 214                                                G4double kineticEnergy)
307 {                                                 215 {
308   SetupKinematics(p, mat, e);                  << 216   G4double mass  = p->GetPDGMass();
309   const G4double dedx  = 0.5*tmax/(kinEnergy + << 217   G4double ratio = electron_mass_c2/mass;
                                                   >> 218   G4double tau   = kineticEnergy / mass;
                                                   >> 219   G4double gamma = 1.0 + tau;
                                                   >> 220   G4double bg2   = tau * (tau+2.0);
                                                   >> 221   G4double tmax  = 2.0*electron_mass_c2*bg2 /(1. + 2.0*gamma*ratio + ratio*ratio);
                                                   >> 222   G4double dedx  = 0.5*tmax/(kineticEnergy + mass);
310   return 0.5*dedx*dedx;                           223   return 0.5*dedx*dedx;
311 }                                                 224 }
312                                                   225 
313 //....oooOO0OOooo........oooOO0OOooo........oo    226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
314                                                   227 
315 G4double G4EmCorrections:: KShellCorrection(co    228 G4double G4EmCorrections:: KShellCorrection(const G4ParticleDefinition* p,
316                                             co << 229                                             const G4Material* material,
317                                             co << 230                                                   G4double kineticEnergy)
318 {                                                 231 {
319   SetupKinematics(p, mat, e);                  << 232   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 233   G4double gamma = 1.0 + tau;
                                                   >> 234   G4double beta2 = tau*(tau + 2.0)/(gamma*gamma);
                                                   >> 235   G4double ba2   = beta2/alpha2;
                                                   >> 236 
320   G4double term = 0.0;                            237   G4double term = 0.0;
321   for (G4int i = 0; i<numberOfElements; ++i) { << 238   const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 239   const G4double* atomDensity  = material->GetAtomicNumDensityVector(); 
                                                   >> 240   G4int numberOfElements = material->GetNumberOfElements();
                                                   >> 241 
                                                   >> 242   for (G4int i = 0; i<numberOfElements; i++) {
322                                                   243 
323     G4double Z = (*theElementVector)[i]->GetZ(    244     G4double Z = (*theElementVector)[i]->GetZ();
324     G4int   iz = (*theElementVector)[i]->GetZa << 245     G4int   iz = G4int(Z);
325     G4double f = 1.0;                             246     G4double f = 1.0;
326     G4double Z2= (Z-0.3)*(Z-0.3);                 247     G4double Z2= (Z-0.3)*(Z-0.3);
327     if(1 == iz) {                                 248     if(1 == iz) {
328       f  = 0.5;                                   249       f  = 0.5;
329       Z2 = 1.0;                                   250       Z2 = 1.0;
330     }                                             251     }
331     const G4double eta = ba2/Z2;               << 252     G4double e0= 13.6*eV*Z2;
332     const G4double tet = (11 < iz) ? sThetaK-> << 253     term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
333     term += f*atomDensity[i]*KShell(tet,eta)/Z << 
334   }                                               254   }
335                                                   255 
336   term /= material->GetTotNbOfAtomsPerVolume()    256   term /= material->GetTotNbOfAtomsPerVolume();
337                                                   257 
338   return term;                                    258   return term;
339 }                                                 259 }
340                                                   260 
341 //....oooOO0OOooo........oooOO0OOooo........oo    261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342                                                   262 
343 G4double G4EmCorrections:: LShellCorrection(co    263 G4double G4EmCorrections:: LShellCorrection(const G4ParticleDefinition* p,
344                                             co << 264                                             const G4Material* material,
345                                             co << 265                                                   G4double kineticEnergy)
346 {                                                 266 {
347   SetupKinematics(p, mat, e);                  << 267   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 268   G4double gamma = 1.0 + tau;
                                                   >> 269   G4double beta2 = tau*(tau + 2.0)/(gamma*gamma);
                                                   >> 270   G4double ba2   = beta2/alpha2;
                                                   >> 271 
348   G4double term = 0.0;                            272   G4double term = 0.0;
349   for (G4int i = 0; i<numberOfElements; ++i) { << 273   const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 274   const G4double* atomDensity  = material->GetAtomicNumDensityVector(); 
                                                   >> 275   G4int numberOfElements = material->GetNumberOfElements();
                                                   >> 276 
                                                   >> 277   for (G4int i = 0; i<numberOfElements; i++) {
350                                                   278 
351     const G4double Z = (*theElementVector)[i]- << 279     G4double Z = (*theElementVector)[i]->GetZ();
352     const G4int iz = (*theElementVector)[i]->G << 280     G4int   iz = G4int(Z);
353     if(2 < iz) {                                  281     if(2 < iz) {
354       const G4double Zeff = (iz < 10) ? Z - ZD << 282       G4double Zeff = Z - ZD[10];
355       const G4double Z2= Zeff*Zeff;            << 283       if(iz < 10) Zeff = Z - ZD[iz];
356       const G4double eta = ba2/Z2;             << 284       G4double Z2= Zeff*Zeff;
357       G4double tet = sThetaL->Value(Z);        << 285       G4double e0= 13.6*eV*Z2*0.25;
358       G4int nmax = std::min(4, G4AtomicShells: << 286       G4double f = 0.125;
359       for (G4int j=1; j<nmax; ++j) {           << 287       G4int nmax = std::min(4,shells.GetNumberOfShells(iz));
360         G4int ne = G4AtomicShells::GetNumberOf << 288       for(G4int j=1; j<nmax; j++) {
361         if (15 >= iz) {                        << 289         G4double ne = G4double(shells.GetNumberOfElectrons(iz,j));
362           tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 290         G4double e1 = shells.GetBindingEnergy(iz,j);
363             0.25*Z2*(1.0 + Z2*alpha2/16.);     << 291   //        G4cout << "LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV << " e0(eV)= " << e0/eV << G4endl; 
364         }                                      << 292         term += f*ne*atomDensity[i]*LShell(e1/e0,ba2/Z2)/Z;
365         //G4cout << " LShell: j= " << j << " n << 
366         //       << " ThetaL= " << tet << G4en << 
367         term += 0.125*ne*atomDensity[i]*LShell << 
368       }                                           293       }
369     }                                             294     }
370   }                                               295   }
371                                                   296 
372   term /= material->GetTotNbOfAtomsPerVolume()    297   term /= material->GetTotNbOfAtomsPerVolume();
373                                                   298 
374   return term;                                    299   return term;
375 }                                                 300 }
376                                                   301 
377 //....oooOO0OOooo........oooOO0OOooo........oo    302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378                                                   303 
379 G4double G4EmCorrections::KShell(const G4doubl << 304 G4double G4EmCorrections::KShell(G4double tet, G4double eta)
380 {                                                 305 {
381   G4double corr = 0.0;                            306   G4double corr = 0.0;
382                                                   307 
383   static const G4double TheK[20] =             << 
384     {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, << 
385      0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, << 
386                                                << 
387                                                << 
388   G4double x = tet;                               308   G4double x = tet;
389   G4int itet = 0;                              << 309   if(tet < TheK[0]) x =  TheK[0]; 
390   G4int ieta = 0;                              << 310   G4int itet = Index(x, TheK, nK);
391   if(tet < TheK[0]) {                          << 311 
392     x =  TheK[0];                              << 312   // assimptotic case
393   } else if(tet > TheK[nK-1]) {                << 
394     x =  TheK[nK-1];                           << 
395     itet = nK-2;                               << 
396   } else {                                     << 
397     itet = Index(x, TheK, nK);                 << 
398   }                                            << 
399   // asymptotic case                           << 
400   if(eta >= Eta[nEtaK-1]) {                       313   if(eta >= Eta[nEtaK-1]) {
401     corr =                                     << 314     corr = (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) + 
402       (Value(x, TheK[itet], TheK[itet+1], UK[i << 315             Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta +
403        Value(x, TheK[itet], TheK[itet+1], VK[i << 316             Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta;
404        Value(x, TheK[itet], TheK[itet+1], ZK[i << 
405   } else {                                        317   } else {
406     G4double y = eta;                             318     G4double y = eta;
407     if(eta < Eta[0]) {                         << 319     if(eta < Eta[0])  y =  Eta[0];
408       y = Eta[0];                              << 320     G4int ieta = Index(y, Eta, nEtaK);
409     } else {                                   << 
410       ieta = Index(y, Eta, nEtaK);             << 
411     }                                          << 
412     corr = Value2(x, y, TheK[itet], TheK[itet+    321     corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
413                   CK[itet][ieta], CK[itet+1][i << 322                   CK[itet][ieta], CK[itet+1][ieta], CK[itet][ieta+1], CK[itet+1][ieta+1]);
414                   CK[itet][ieta+1], CK[itet+1] << 323     //G4cout << "   x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]<<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
415     //G4cout << "   x= " <<x<<" y= "<<y<<" tet << 324     //       <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta]<<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl;
416     //           <<" "<< TheK[itet+1]<<" eta=  << 
417     //           <<" CK= " << CK[itet][ieta]<< << 
418     //           <<" "<< CK[itet][ieta+1]<<" " << 
419   }                                               325   }
420   //G4cout << "Kshell:  tet= " << tet << " eta << 326   //G4cout << "Kshell:  tet= " << tet << " eta= " << eta << "  C= " << corr << " itet,ieta= " << itet <<G4endl;
421   //         << " itet= " << itet << "  ieta=  << 
422   return corr;                                    327   return corr;
423 }                                                 328 }
424                                                   329 
425 //....oooOO0OOooo........oooOO0OOooo........oo    330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426                                                   331 
427 G4double G4EmCorrections::LShell(const G4doubl << 332 G4double G4EmCorrections::LShell(G4double tet, G4double eta) 
428 {                                                 333 {
429   G4double corr = 0.0;                            334   G4double corr = 0.0;
430                                                   335 
431   static const G4double TheL[26] =             << 
432     {0.24, 0.26, 0.28, 0.30, 0.32,   0.34, 0.3 << 
433      0.42, 0.44, 0.45, 0.46, 0.48,   0.50, 0.5 << 
434      0.58, 0.60, 0.62, 0.64, 0.65,   0.66};    << 
435                                                << 
436   G4double x = tet;                               336   G4double x = tet;
437   G4int itet = 0;                              << 337   if(tet < TheL[0]) x =  TheL[0]; 
438   G4int ieta = 0;                              << 338   G4int itet = Index(x, TheL, nL);
439   if(tet < TheL[0]) {                          << 
440     x =  TheL[0];                              << 
441   } else if(tet > TheL[nL-1]) {                << 
442     x =  TheL[nL-1];                           << 
443     itet = nL-2;                               << 
444   } else {                                     << 
445     itet = Index(x, TheL, nL);                 << 
446   }                                            << 
447                                                   339 
448   // asymptotic case                           << 340   // assimptotic case
449   if(eta >= Eta[nEtaL-1]) {                       341   if(eta >= Eta[nEtaL-1]) {
450     corr = (Value(x, TheL[itet], TheL[itet+1], << 342     corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1])  
451               + Value(x, TheL[itet], TheL[itet << 343     + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta
452             )/eta;                             << 344            )/eta;
453   } else {                                        345   } else {
454     G4double y = eta;                             346     G4double y = eta;
455     if(eta < Eta[0]) {                         << 347     if(eta < Eta[0])  y =  Eta[0];
456       y =  Eta[0];                             << 348     G4int ieta = Index(y, Eta, nEtaL);
457     } else {                                   << 
458       ieta = Index(y, Eta, nEtaL);             << 
459     }                                          << 
460     corr = Value2(x, y, TheL[itet], TheL[itet+    349     corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1],
461                   CL[itet][ieta], CL[itet+1][i << 350                CL[itet][ieta], CL[itet+1][ieta], CL[itet][ieta+1], CL[itet+1][ieta+1]);
462                   CL[itet][ieta+1], CL[itet+1] << 351     //G4cout << "   x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet]<<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
463     //G4cout << "   x= " <<x<<" y= "<<y<<" tet << 352     //       <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta]<<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl;
464     //           <<" "<< TheL[itet+1]<<" eta=  << 
465     //           <<" CL= " << CL[itet][ieta]<< << 
466     //           <<" "<< CL[itet][ieta+1]<<" " << 
467   }                                               353   }
468   //G4cout<<"Lshell:  tet= "<<tet<<" eta= "<<e << 354   //  G4cout << "Lshell:  tet= " << tet << " eta= " << eta << "  C= " << corr << " itet= " << itet <<G4endl;
469   //        <<" ieta= "<<ieta<<"  Corr= "<<cor << 
470   return corr;                                    355   return corr;
471 }                                                 356 }
472                                                   357 
473 //....oooOO0OOooo........oooOO0OOooo........oo    358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
474                                                   359 
475 G4double G4EmCorrections::ShellCorrectionSTD(c    360 G4double G4EmCorrections::ShellCorrectionSTD(const G4ParticleDefinition* p,
476                                              c << 361                const G4Material* material,
477                                              c << 362                      G4double kineticEnergy)
478 {                                                 363 {
479   SetupKinematics(p, mat, e);                  << 364   G4double tau   = kineticEnergy / p->GetPDGMass();
480   G4double taulim= 8.0*MeV/mass;               << 365   G4double bg2   = tau * (tau+2.0);
                                                   >> 366   G4double taulim= 8.0*MeV/p->GetPDGMass();
481   G4double bg2lim= taulim * (taulim+2.0);         367   G4double bg2lim= taulim * (taulim+2.0);
482                                                   368 
483   G4double* shellCorrectionVector =               369   G4double* shellCorrectionVector =
484             material->GetIonisation()->GetShel    370             material->GetIonisation()->GetShellCorrectionVector();
485   G4double sh = 0.0;                              371   G4double sh = 0.0;
486   G4double x  = 1.0;                              372   G4double x  = 1.0;
487   G4double taul  = material->GetIonisation()->    373   G4double taul  = material->GetIonisation()->GetTaul();
488                                                   374 
489   if ( bg2 >= bg2lim ) {                          375   if ( bg2 >= bg2lim ) {
490     for (G4int k=0; k<3; ++k) {                << 376     for (G4int k=0; k<3; k++) {
491         x *= bg2 ;                                377         x *= bg2 ;
492         sh += shellCorrectionVector[k]/x;         378         sh += shellCorrectionVector[k]/x;
493     }                                             379     }
494                                                   380 
495   } else {                                        381   } else {
496     for (G4int k=0; k<3; ++k) {                << 382     for (G4int k=0; k<3; k++) {
497         x *= bg2lim ;                             383         x *= bg2lim ;
498         sh += shellCorrectionVector[k]/x;         384         sh += shellCorrectionVector[k]/x;
499     }                                             385     }
500     sh *= G4Log(tau/taul)/G4Log(taulim/taul);  << 386     sh *= std::log(tau/taul)/std::log(taulim/taul);
501   }                                               387   }
502   sh *= 0.5;                                      388   sh *= 0.5;
503   return sh;                                      389   return sh;
504 }                                                 390 }
505                                                   391 
506                                                   392 
507 //....oooOO0OOooo........oooOO0OOooo........oo    393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
508                                                   394 
509 G4double G4EmCorrections::ShellCorrection(cons    395 G4double G4EmCorrections::ShellCorrection(const G4ParticleDefinition* p,
510                                           cons << 396                                            const G4Material* material,
511                                           cons << 397                                                  G4double kineticEnergy)
512 {                                                 398 {
513   SetupKinematics(p, mat, ekin);               << 399   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 400   G4double gamma = 1.0 + tau;
                                                   >> 401   G4double beta2 = tau*(tau + 2.0)/(gamma*gamma);
                                                   >> 402   G4double ba2   = beta2/alpha2;
                                                   >> 403 
514   G4double term = 0.0;                            404   G4double term = 0.0;
515   //G4cout << "### G4EmCorrections::ShellCorre << 405   const G4ElementVector* theElementVector = material->GetElementVector();
516   //       << "   " << ekin/MeV << " MeV " <<  << 406   const G4double* atomDensity  = material->GetAtomicNumDensityVector();
517   for (G4int i = 0; i<numberOfElements; ++i) { << 407   G4int numberOfElements = material->GetNumberOfElements();
518                                                << 408 
519     G4double res = 0.0;                        << 409   for (G4int i = 0; i<numberOfElements; i++) {
520     G4double res0 = 0.0;                       << 410 
521     const G4double Z = (*theElementVector)[i]- << 411     G4double Z = (*theElementVector)[i]->GetZ();
522     const G4int iz = (*theElementVector)[i]->G << 412     G4int   iz = G4int(Z);
523     G4double Z2= (Z-0.3)*(Z-0.3);                 413     G4double Z2= (Z-0.3)*(Z-0.3);
524     G4double f = 1.0;                             414     G4double f = 1.0;
525     if(1 == iz) {                                 415     if(1 == iz) {
526       f  = 0.5;                                   416       f  = 0.5;
527       Z2 = 1.0;                                   417       Z2 = 1.0;
528     }                                             418     }
529     G4double eta = ba2/Z2;                     << 419     G4double e0= 13.6*eV*Z2;
530     G4double tet = (11 < iz) ? sThetaK->Value( << 420     term += f*atomDensity[i]*KShell(shells.GetBindingEnergy(iz,0)/e0,ba2/Z2)/Z;
531     res0 = f*KShell(tet,eta);                  << 
532     res += res0;                               << 
533     //G4cout << " Z= " << iz << " Shell 0" <<  << 
534     //       << " eta= " << eta << "  resK= "  << 
535                                                << 
536     if(2 < iz) {                                  421     if(2 < iz) {
537       const G4double Zeff = (iz < 10) ? Z - ZD << 422       G4double Zeff = Z - ZD[10];
                                                   >> 423       if(iz < 10) Zeff = Z - ZD[iz];
538       Z2= Zeff*Zeff;                              424       Z2= Zeff*Zeff;
539       eta = ba2/Z2;                            << 425       e0= 13.6*eV*Z2*0.25;
540       tet = sThetaL->Value(Z);                 << 
541       f = 0.125;                                  426       f = 0.125;
542       const G4int ntot = G4AtomicShells::GetNu << 427       G4double eta = ba2/Z2;
543       const G4int nmax = std::min(4, ntot);    << 428       G4int ntot = shells.GetNumberOfShells(iz);
                                                   >> 429       G4int nmax = std::min(4, ntot);
544       G4double norm   = 0.0;                      430       G4double norm   = 0.0;
545       G4double eshell = 0.0;                      431       G4double eshell = 0.0;
546       for(G4int j=1; j<nmax; ++j) {            << 432       for(G4int j=1; j<nmax; j++) {
547         G4int ne = G4AtomicShells::GetNumberOf << 433         G4double x = G4double(shells.GetNumberOfElectrons(iz,j));
548         if(15 >= iz) {                         << 434         G4double e1 = shells.GetBindingEnergy(iz,j);
549           tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 435   norm   += x;
550       0.25*Z2*(1.0 + Z2*alpha2/16.);           << 436   eshell += e1*x;
551         }                                      << 437         term += f*x*atomDensity[i]*LShell(e1/e0,eta)/Z;
552         norm   += ne;                          << 
553         eshell += tet*ne;                      << 
554         res0 = f*ne*LShell(tet,eta);           << 
555         res += res0;                           << 
556         //G4cout << " Zeff= " << Zeff << " She << 
557         //       << " tet= " << tet << " eta=  << 
558         //       << "  resL= " << res0 << G4en << 
559       }                                           438       }
560       if(ntot > nmax) {                        << 439       if(10 < iz) {
561   if (norm > 0.0) { norm = 1.0/norm; }         << 440   eshell /= norm;
562         eshell *= norm;                        << 441   G4double eeff = eshell*eta;
563                                                << 442   for(G4int k=nmax; k<ntot; k++) {
564   static const G4double HM[53] = {             << 443           G4double x = G4double(shells.GetNumberOfElectrons(iz,k));
565     12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5,  << 444           G4double e1 = shells.GetBindingEnergy(iz,k);
566     10.0,  9.51, 8.97, 8.52, 8.03, 7.46, 6.95, << 445           term += f*x*atomDensity[i]*LShell(e1/e0,eeff/e1)/Z;
567     5.61,  5.39, 5.19, 5.01, 4.86, 4.72, 4.62, << 446     //          term += f*x*atomDensity[i]*LShell(eshell/e0,eeff/e1)/Z;
568     4.32,  4.26, 4.20, 4.15, 4.1,  4.04, 4.00, << 447   }
569     3.90,  3.89, 3.89, 3.88, 3.88, 3.88, 3.88, << 448   /*
570     3.90, 3.92, 3.93 };                        << 449         if(28 >= iz) {
571   static const G4double HN[31] = {             << 450           term += f*(Z - 10.)*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 
572     75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9,  << 451         } else if(32 >= iz) {
573     24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5,  << 452           term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 
574     19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7,  << 453         } else if(60 >= iz) {
575     18.2};                                     << 454           term += f*18.0*atomDensity[i]*LShell(eshell,HM[iz-11]*eta)/Z; 
576                                                << 455           term += f*(Z - 28.)*atomDensity[i]*LShell(eshell,HN[iz-33]*eta)/Z; 
577         // Add M-shell                         << 
578         if(28 > iz) {                          << 
579           res += f*(iz - 10)*LShell(eshell,HM[ << 
580         } else if(63 > iz) {                   << 
581           res += f*18*LShell(eshell,HM[iz-11]* << 
582         } else {                                  456         } else {
583           res += f*18*LShell(eshell,HM[52]*eta << 457           term += f*18.0*atomDensity[i]*LShell(eshell,HM[53]*eta)/Z;
584         }                                      << 458           term += f*32.0*atomDensity[i]*LShell(eshell,HN[30]*eta)/Z; 
585         // Add N-shell                         << 459           term += f*(Z - 60.)*atomDensity[i]*LShell(eshell,150.*eta)/Z; 
586         if(32 < iz) {                          << 460   }
587           if(60 > iz) {                        << 461   */
588             res += f*(iz - 28)*LShell(eshell,H << 
589           } else if(63 > iz) {                 << 
590             res += 4*LShell(eshell,HN[iz-33]*e << 
591           } else {                             << 
592             res += 4*LShell(eshell,HN[30]*eta) << 
593           }                                    << 
594           // Add O-P-shells                    << 
595           if(60 < iz) {                        << 
596             res += f*(iz - 60)*LShell(eshell,1 << 
597           }                                    << 
598         }                                      << 
599       }                                           462       }
600     }                                             463     }
601     term += res*atomDensity[i]/Z;              << 464     //term += atomDensity[i]*MSH[iz]/(ba2*ba2);
602   }                                               465   }
603                                                   466 
604   term /= material->GetTotNbOfAtomsPerVolume()    467   term /= material->GetTotNbOfAtomsPerVolume();
605   //G4cout << "##Shell Correction=" << term << << 
606   return term;                                    468   return term;
607 }                                                 469 }
608                                                   470 
609 //....oooOO0OOooo........oooOO0OOooo........oo    471 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
610                                                   472 
611 G4double G4EmCorrections::DensityCorrection(co    473 G4double G4EmCorrections::DensityCorrection(const G4ParticleDefinition* p,
612                                             co << 474                                             const G4Material* material,
613                                             co << 475                                                   G4double kineticEnergy)
614 {                                                 476 {
615   SetupKinematics(p, mat, e);                  << 477   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 478   G4double bg2   = tau * (tau+2.0);
616                                                   479 
617   G4double cden  = material->GetIonisation()->    480   G4double cden  = material->GetIonisation()->GetCdensity();
618   G4double mden  = material->GetIonisation()->    481   G4double mden  = material->GetIonisation()->GetMdensity();
619   G4double aden  = material->GetIonisation()->    482   G4double aden  = material->GetIonisation()->GetAdensity();
620   G4double x0den = material->GetIonisation()->    483   G4double x0den = material->GetIonisation()->GetX0density();
621   G4double x1den = material->GetIonisation()->    484   G4double x1den = material->GetIonisation()->GetX1density();
622                                                   485 
                                                   >> 486   G4double twoln10 = 2.0*std::log(10.0);
623   G4double dedx = 0.0;                            487   G4double dedx = 0.0;
624                                                   488 
625   // density correction                           489   // density correction
626   static const G4double twoln10 = 2.0*G4Log(10 << 490   G4double x = std::log(bg2)/twoln10;
627   G4double x = G4Log(bg2)/twoln10;             << 
628   if ( x >= x0den ) {                             491   if ( x >= x0den ) {
629     dedx = twoln10*x - cden ;                     492     dedx = twoln10*x - cden ;
630     if ( x < x1den ) dedx += aden*G4Exp(G4Log( << 493     if ( x < x1den ) dedx += aden*std::pow((x1den-x),mden) ;
631   }                                               494   }
632                                                   495 
633   return dedx;                                    496   return dedx;
634 }                                                 497 }
635                                                   498 
636 //....oooOO0OOooo........oooOO0OOooo........oo    499 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
637                                                   500 
638 G4double G4EmCorrections::BarkasCorrection(con    501 G4double G4EmCorrections::BarkasCorrection(const G4ParticleDefinition* p,
639                                            con << 502                                            const G4Material* material,
640                                            con << 503                                                  G4double kineticEnergy)
641                                            con << 
642 {                                                 504 {
643   // . Z^3 Barkas effect in the stopping power << 505 // . Z^3 Barkas effect in the stopping power of matter for charged particles
644   //   J.C Ashley and R.H.Ritchie              << 506 //   J.C Ashley and R.H.Ritchie
645   //   Physical review B Vol.5 No.7 1 April 19 << 507 //   Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397
646   //   valid for kineticEnergy > 0.5 MeV       << 508 //   valid for kineticEnergy > 0.5 MeV
                                                   >> 509 
                                                   >> 510   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 511   G4double gamma = 1.0 + tau;
                                                   >> 512   G4double beta2 = tau*(tau + 2.0)/(gamma*gamma);
                                                   >> 513 
                                                   >> 514   G4double q  = effCharge.EffectiveCharge(p,material,kineticEnergy)/eplus;
                                                   >> 515   G4double ba2 = beta2/alpha2;
647                                                   516 
648   if (!isInitialized) { SetupKinematics(p, mat << 
649   G4double BarkasTerm = 0.0;                      517   G4double BarkasTerm = 0.0;
                                                   >> 518   const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 519   const G4double* atomDensity  = material->GetAtomicNumDensityVector();
                                                   >> 520   G4int numberOfElements = material->GetNumberOfElements();
                                                   >> 521 
                                                   >> 522   for (G4int i = 0; i<numberOfElements; i++) {
                                                   >> 523 
                                                   >> 524     G4double Z = (*theElementVector)[i]->GetZ();
                                                   >> 525     G4int iz = G4int(Z);
                                                   >> 526 
                                                   >> 527     G4double X = ba2 / Z;
650                                                   528 
651   for (G4int i = 0; i<numberOfElements; ++i) { << 529     // Variables to compute L_1
                                                   >> 530     //    G4double Eta0Chi = 0.8;
                                                   >> 531     //    G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( Z,-1.19 ) );
                                                   >> 532     //    G4double W = ( EtaChi * std::pow( Z,0.16666667 ) ) / std::sqrt(X);
652                                                   533 
653     const G4int iz = (*theElementVector)[i]->G << 534     G4double b = 1.3;
654     if(iz == 47) {                             << 535     if(1 == iz) {
655       BarkasTerm += atomDensity[i]*0.006812*G4 << 536       if(material->GetName() == "G4_lH2") b = 0.6;
656     } else if(iz >= 64) {                      << 537       else                                b = 1.8;
657       BarkasTerm += atomDensity[i]*0.002833*G4 << 538     } 
658     } else {                                   << 539     else if(2 == iz)  b = 0.6;
659                                                << 540     else if(10 >= iz) b = 1.8;
660       const G4double Z = (*theElementVector)[i << 541     else if(17 >= iz) b = 1.4;
661       const G4double X = ba2 / Z;              << 542     else if(18 == iz) b = 1.8;
662       G4double b = 1.3;                        << 543     else if(25 >= iz) b = 1.4;
663       if(1 == iz) { b = (material->GetName() = << 544     else if(50 >= iz) b = 1.35;
664       else if(2 == iz)  { b = 0.6; }           << 545 
665       else if(10 >= iz) { b = 1.8; }           << 546     G4double W = b/std::sqrt(X);
666       else if(17 >= iz) { b = 1.4; }           << 547 
667       else if(18 == iz) { b = 1.8; }           << 548     G4double val;
668       else if(25 >= iz) { b = 1.4; }           << 549     if(W <= engBarkas[0])       val =  corBarkas[0];
669       else if(50 >= iz) { b = 1.35;}           << 550     else if(W >= engBarkas[46]) val =  corBarkas[46]*engBarkas[46]/W;
670                                                << 551     else {
671       const G4double W = b/std::sqrt(X);       << 552       G4int iw = Index(W, engBarkas, 47);
672                                                << 553       val = Value(W, engBarkas[iw], engBarkas[iw+1], corBarkas[iw], corBarkas[iw+1]);
673       G4double val = sBarkasCorr->Value(W, idx << 
674       if (W > sWmaxBarkas) { val *= (sWmaxBark << 
675       //    G4cout << "i= " << i << " b= " <<  << 
676       // << " Z= " << Z << " X= " << X << " va << 
677       BarkasTerm += val*atomDensity[i] / (std: << 
678     }                                             554     }
                                                   >> 555     BarkasTerm += val*atomDensity[i] * std::sqrt(1.0/(Z*X))/ X;
679   }                                               556   }
680                                                   557 
681   BarkasTerm *= 1.29*charge/material->GetTotNb << 558   BarkasTerm *= q/material->GetTotNbOfAtomsPerVolume();
682                                                   559 
683   return BarkasTerm;                              560   return BarkasTerm;
684 }                                                 561 }
685                                                << 
686 //....oooOO0OOooo........oooOO0OOooo........oo    562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
687                                                   563 
688 G4double G4EmCorrections::BlochCorrection(cons    564 G4double G4EmCorrections::BlochCorrection(const G4ParticleDefinition* p,
689                                           cons << 565                                           const G4Material* material,
690                                           cons << 566                                                 G4double kineticEnergy)
691                                           cons << 
692 {                                                 567 {
693   if (!isInitialized) { SetupKinematics(p, mat << 568   G4double tau   = kineticEnergy / p->GetPDGMass();
                                                   >> 569   G4double gamma = 1.0 + tau;
                                                   >> 570   G4double beta2 = tau*(tau + 2.0)/(gamma*gamma);
694                                                   571 
695   G4double y2 = q2/ba2;                        << 572   G4double q  = effCharge.EffectiveCharge(p,material,kineticEnergy)/eplus;
                                                   >> 573   G4double q2 = q*q;
696                                                   574 
                                                   >> 575   G4double y2 = q2*alpha2 / beta2;
                                                   >> 576   
697   G4double term = 1.0/(1.0 + y2);                 577   G4double term = 1.0/(1.0 + y2);
698   G4double del;                                   578   G4double del;
699   G4double j = 1.0;                               579   G4double j = 1.0;
700   do {                                            580   do {
701     j += 1.0;                                     581     j += 1.0;
702     del = 1.0/(j* (j*j + y2));                    582     del = 1.0/(j* (j*j + y2));
703     term += del;                                  583     term += del;
704     // Loop checking, 03-Aug-2015, Vladimir Iv << 
705   } while (del > 0.01*term);                      584   } while (del > 0.01*term);
706                                                   585 
707   return -y2*term;                                586   return -y2*term;
708 }                                                 587 }
709                                                   588 
710 //....oooOO0OOooo........oooOO0OOooo........oo    589 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
711                                                   590 
712 G4double G4EmCorrections::MottCorrection(const    591 G4double G4EmCorrections::MottCorrection(const G4ParticleDefinition* p,
713                                          const << 592                                          const G4Material* material,
714                                          const << 593                                                G4double kineticEnergy)
715                                          const << 
716 {                                                 594 {
717   if (!isInitialized) { SetupKinematics(p, mat << 595   G4double tau   = kineticEnergy / p->GetPDGMass();
718   return CLHEP::pi*CLHEP::fine_structure_const << 596   G4double gamma = 1.0 + tau;
719 }                                              << 597   G4double bg2   = tau*(tau + 2.0);
                                                   >> 598   G4double beta2 = bg2/(gamma*gamma);
                                                   >> 599 
                                                   >> 600   G4double q  = effCharge.EffectiveCharge(p,material,kineticEnergy)/eplus;
                                                   >> 601 
                                                   >> 602   G4double beta = std::sqrt(beta2);
                                                   >> 603   G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
                                                   >> 604 
                                                   >> 605   G4double mterm = 0.0;
                                                   >> 606 
                                                   >> 607   if(beta > 0.0) {
                                                   >> 608 
                                                   >> 609     G4double invbeta = 1.0/beta;
                                                   >> 610     G4double invbeta2= invbeta*invbeta;
                                                   >> 611     G4double za  = q*fine_structure_const;
                                                   >> 612     G4double za2 = za*za;
                                                   >> 613     G4double za3 = za2*za;
                                                   >> 614     G4double x   = za*invbeta;
                                                   >> 615     G4double cosx;
                                                   >> 616     if(x < COSEB[13]) {
                                                   >> 617       G4int i = Index(x,COSEB,14);
                                                   >> 618       cosx    = Value(x,COSEB[i], COSEB[i+1],COSXI[i],COSXI[i+1]);
                                                   >> 619     } else {
                                                   >> 620       cosx = COSXI[13]*COSEB[13]/x;
                                                   >> 621     }
720                                                   622 
721 //....oooOO0OOooo........oooOO0OOooo........oo << 623     mterm =
                                                   >> 624         za*beta*(1.725 + pi*cosx*(0.52 - 2.0*std::sqrt(eexc/(2.0*electron_mass_c2*bg2))))
                                                   >> 625       + za2*(3.246 - 0.451*beta2)
                                                   >> 626       + za3*(1.522*beta + 0.987*invbeta)
                                                   >> 627       + za2*za2*(4.569 - 0.494*beta2 - 2.696*invbeta2)
                                                   >> 628       + za3*za2*(1.254*beta + 0.222*invbeta - 1.17*invbeta*invbeta2);
722                                                   629 
723 G4double                                       << 
724 G4EmCorrections::EffectiveChargeCorrection(con << 
725                                            con << 
726                                            con << 
727 {                                              << 
728   G4double factor = 1.0;                       << 
729   if(p->GetPDGCharge() <= 2.5*CLHEP::eplus ||  << 
730                                                << 
731   if(verbose > 1) {                            << 
732     G4cout << "EffectiveChargeCorrection: " << << 
733            << " in " << mat->GetName()         << 
734            << " ekin(MeV)= " << ekin/MeV << G4 << 
735   }                                               630   }
736                                                   631 
737   if(p != curParticle || mat != curMaterial) { << 632   return mterm;
738     curParticle = p;                           << 
739     curMaterial = mat;                         << 
740     curVector   = nullptr;                     << 
741     currentZ = p->GetAtomicNumber();           << 
742     if(verbose > 1) {                          << 
743       G4cout << "G4EmCorrections::EffectiveCha << 
744              << currentZ << " Aion= " << p->Ge << 
745     }                                          << 
746     massFactor = CLHEP::proton_mass_c2/p->GetP << 
747     idx = -1;                                  << 
748                                                << 
749     for(G4int i=0; i<nIons; ++i) {             << 
750       if(materialList[i] == mat && currentZ == << 
751         idx = i;                               << 
752         break;                                 << 
753       }                                        << 
754     }                                          << 
755     //G4cout << " idx= " << idx << " dz= " <<  << 
756     if(idx >= 0) {                             << 
757       if(nullptr == ionList[idx]) { BuildCorre << 
758       curVector = stopData[idx];               << 
759     } else {                                   << 
760       return factor;                           << 
761     }                                          << 
762   }                                            << 
763   if(nullptr != curVector) {                   << 
764     factor = curVector->Value(ekin*massFactor) << 
765     if(verbose > 1) {                          << 
766       G4cout << "E= " << ekin << " factor= " < << 
767              << massFactor << G4endl;          << 
768     }                                          << 
769   }                                            << 
770   return factor;                               << 
771 }                                                 633 }
772                                                   634 
773 //....oooOO0OOooo........oooOO0OOooo........oo    635 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
774                                                   636 
775 void G4EmCorrections::AddStoppingData(const G4 << 637 G4double G4EmCorrections::NuclearDEDX(const G4ParticleDefinition* p,
776                                       const G4 << 638                                       const G4Material* material,
777                                       G4Physic << 639                     G4double kineticEnergy,
                                                   >> 640                                             G4bool fluct)
778 {                                                 641 {
779   G4int i = 0;                                 << 642   G4double nloss = 0.0;
780   for(; i<nIons; ++i) {                        << 643   if(kineticEnergy <= 0.0) return nloss; 
781     if(Z == Zion[i] && A == Aion[i] && mname = << 
782   }                                            << 
783   if(i == nIons) {                             << 
784     Zion.push_back(Z);                         << 
785     Aion.push_back(A);                         << 
786     materialName.push_back(mname);             << 
787     materialList.push_back(nullptr);           << 
788     ionList.push_back(nullptr);                << 
789     stopData.push_back(dVector);               << 
790     nIons++;                                   << 
791     if(verbose > 1) {                          << 
792       G4cout << "AddStoppingData Z= " << Z <<  << 
793              << "  idx= " << i << G4endl;      << 
794     }                                          << 
795   }                                            << 
796 }                                              << 
797                                                   644 
798 //....oooOO0OOooo........oooOO0OOooo........oo << 645   lossFlucFlag = fluct;
799                                                   646 
800 void G4EmCorrections::BuildCorrectionVector()  << 647   // Projectile nucleus
801 {                                              << 648   G4double z1 = std::abs((p->GetPDGCharge())/eplus) ;
802   if(nullptr == ionLEModel || nullptr == ionHE << 649   G4double m1 = p->GetPDGMass()/amu_c2;
803     return;                                    << 
804   }                                            << 
805                                                   650 
806   const G4ParticleDefinition* ion = curParticl << 651   G4int NumberOfElements = material->GetNumberOfElements() ;
807   const G4ParticleDefinition* gion = G4Generic << 652   const G4ElementVector* theElementVector = 
808   G4int Z = Zion[idx];                         << 653                                  material->GetElementVector() ;
809   G4double A = Aion[idx];                      << 654   const G4double* theAtomicNumDensityVector = 
810   G4PhysicsVector* v = stopData[idx];          << 655                                  material->GetAtomicNumDensityVector() ;
811                                                << 
812   if(verbose > 1) {                            << 
813     G4cout << "BuildCorrectionVector: Stopping << 
814            << curParticle->GetParticleName() < << 
815            << materialName[idx] << " Ion Z= "  << 
816            << " massFactor= " << massFactor << << 
817     G4cout << "    Nbins=" << nbinCorr << " Em << 
818      << " Emax(MeV)=" << eCorrMax << " ion: "  << 
819            << ion->GetParticleName() << G4endl << 
820   }                                            << 
821                                                   656 
822   auto vv = new G4PhysicsLogVector(eCorrMin,eC << 657   //  loop for the elements in the material
823   const G4double eth0 = v->Energy(0);          << 658 
824   const G4double escal = eth/massFactor;       << 659   for (G4int iel=0; iel<NumberOfElements; iel++) {
825   G4double qe =                                << 660     const G4Element* element = (*theElementVector)[iel] ;
826     effCharge.EffectiveChargeSquareRatio(curPa << 661     G4double z2 = element->GetZ();
827   const G4double dedxt =                       << 662     G4double m2 = element->GetA()*mole/g ;
828     ionLEModel->ComputeDEDXPerVolume(curMateri << 663     nloss += (NuclearStoppingPower(kineticEnergy, z1, z2, m1, m2))
829   const G4double dedx1t =                      << 664            * theAtomicNumDensityVector[iel] ;    
830     ionHEModel->ComputeDEDXPerVolume(curMateri << 
831     + ComputeIonCorrections(curParticle, curMa << 
832   const G4double rest = escal*(dedxt - dedx1t) << 
833   if(verbose > 1) {                            << 
834     G4cout << "Escal(MeV)= " << escal << " qe= << 
835            << " dedxt= " << dedxt << " dedx1t= << 
836   }                                            << 
837   for(G4int i=0; i<=nbinCorr; ++i) {           << 
838     // energy in the local table (GenericIon)  << 
839     G4double e = vv->Energy(i);                << 
840     // energy of the real ion                  << 
841     G4double eion = e/massFactor;              << 
842     // energy in the imput stopping data vecto << 
843     G4double e1 = eion/A;                      << 
844     G4double dedx = (e1 < eth0)                << 
845       ? v->Value(eth0)*std::sqrt(e1/eth0) : v- << 
846     qe = effCharge.EffectiveChargeSquareRatio( << 
847     G4double dedx1 = (e <= eth)                << 
848       ? ionLEModel->ComputeDEDXPerVolume(curMa << 
849       : ionHEModel->ComputeDEDXPerVolume(curMa << 
850       ComputeIonCorrections(curParticle, curMa << 
851     vv->PutValue(i, dedx/dedx1);               << 
852     if(verbose > 1) {                          << 
853       G4cout << "E(MeV)=" << e/CLHEP::MeV << " << 
854        << " e1=" << e1 << " dedxRatio= " << de << 
855              << " dedx="  << dedx << " dedx1=" << 
856              << " qe=" << qe << " rest/eion="  << 
857     }                                          << 
858   }                                               665   }
859   delete v;                                    << 666   nloss *= theZieglerFactor;
860   ionList[idx]  = ion;                         << 667   return nloss;
861   stopData[idx] = vv;                          << 
862   if(verbose > 1) { G4cout << "End data set "  << 
863 }                                                 668 }
864                                                   669 
865 //....oooOO0OOooo........oooOO0OOooo........oo    670 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
866                                                   671 
867 void G4EmCorrections::InitialiseForNewRun()    << 672 G4double G4EmCorrections::NuclearStoppingPower(G4double kineticEnergy, 
                                                   >> 673                                                G4double z1, G4double z2, 
                                                   >> 674                                                G4double m1, G4double m2)
868 {                                                 675 {
869   G4ProductionCutsTable* tb = G4ProductionCuts << 676   G4double energy = kineticEnergy/keV ;  // energy in keV
870   ncouples = tb->GetTableSize();               << 677   G4double nloss = 0.0;
871   if(currmat.size() != ncouples) {             << 678   
872     currmat.resize(ncouples);                  << 679   G4double rm;
873     for(auto it = thcorr.begin(); it != thcorr << 680   if(z1 > 1.5) rm = (m1 + m2) * ( std::pow(z1, .23) + std::pow(z2, .23) ) ;
874       (it->second).clear();                    << 681   else         rm = (m1 + m2) * std::pow(z2, 0.333333);  
875     }                                          << 682 
876     thcorr.clear();                            << 683   G4double er = 32.536 * m2 * energy / ( z1 * z2 * rm ) ;  // reduced energy
877     for(std::size_t i=0; i<ncouples; ++i) {    << 684     
878       currmat[i] = tb->GetMaterialCutsCouple(( << 685   for (G4int i=1; i<104; i++)
879       G4String nam = currmat[i]->GetName();    << 686     {
880       for(G4int j=0; j<nIons; ++j) {           << 687       if (er > e[i]) {
881         if(nam == materialName[j]) { materialL << 688   nloss = (a[i] - a[i-1])*(er-e[i-1])/(e[i] - e[i-1]) + a[i-1];   
                                                   >> 689   break;
882       }                                           690       }
883     }                                             691     }
884   }                                            << 692 
                                                   >> 693   // Stragling
                                                   >> 694   if(lossFlucFlag) {
                                                   >> 695     G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)*
                                                   >> 696                   (4.0 + 0.197*std::pow(er,-1.6991)+6.584*std::pow(er,-1.0494))) ;
                                                   >> 697 
                                                   >> 698     nloss *= G4RandGauss::shoot(1.0,sig) ;
                                                   >> 699   } 
                                                   >> 700    
                                                   >> 701   nloss *= 8.462 * z1 * z2 * m1 / rm ; // Return to [ev/(10^15 atoms/cm^2]
                                                   >> 702 
                                                   >> 703   if ( nloss < 0.0) nloss = 0.0 ;
                                                   >> 704   
                                                   >> 705   return nloss;
885 }                                                 706 }
886                                                   707 
887 //....oooOO0OOooo........oooOO0OOooo........oo    708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888                                                   709 
889 void G4EmCorrections::Initialise()                710 void G4EmCorrections::Initialise()
890 {                                                 711 {
891   // Z^3 Barkas effect in the stopping power o << 712 // . Z^3 Barkas effect in the stopping power of matter for charged particles
892   // J.C Ashley and R.H.Ritchie                << 713 //   J.C Ashley and R.H.Ritchie
893   // Physical review B Vol.5 No.7 1 April 1972 << 714 //   Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
894   // G.S. Khandelwal Nucl. Phys. A116(1968)97  << 
895   // "Shell corrections for K- and L- electron << 
896                                                << 
897   G4int i, j;                                     715   G4int i, j;
898   static const G4double fTable[47][2] = {      << 716   const G4double fTable[47][2] = {
899    { 0.02, 21.5},                                 717    { 0.02, 21.5},
900    { 0.03, 20.0},                                 718    { 0.03, 20.0},
901    { 0.04, 18.0},                                 719    { 0.04, 18.0},
902    { 0.05, 15.6},                                 720    { 0.05, 15.6},
903    { 0.06, 15.0},                                 721    { 0.06, 15.0},
904    { 0.07, 14.0},                                 722    { 0.07, 14.0},
905    { 0.08, 13.5},                                 723    { 0.08, 13.5},
906    { 0.09, 13.},                                  724    { 0.09, 13.},
907    { 0.1,  12.2},                                 725    { 0.1,  12.2},
908    { 0.2,  9.25},                                 726    { 0.2,  9.25},
909    { 0.3,  7.0},                                  727    { 0.3,  7.0},
910    { 0.4,  6.0},                                  728    { 0.4,  6.0},
911    { 0.5,  4.5},                                  729    { 0.5,  4.5},
912    { 0.6,  3.5},                                  730    { 0.6,  3.5},
913    { 0.7,  3.0},                                  731    { 0.7,  3.0},
914    { 0.8,  2.5},                                  732    { 0.8,  2.5},
915    { 0.9,  2.0},                                  733    { 0.9,  2.0},
916    { 1.0,  1.7},                                  734    { 1.0,  1.7},
917    { 1.2,  1.2},                                  735    { 1.2,  1.2},
918    { 1.3,  1.0},                                  736    { 1.3,  1.0},
919    { 1.4,  0.86},                                 737    { 1.4,  0.86},
920    { 1.5,  0.7},                                  738    { 1.5,  0.7},
921    { 1.6,  0.61},                                 739    { 1.6,  0.61},
922    { 1.7,  0.52},                                 740    { 1.7,  0.52},
923    { 1.8,  0.5},                                  741    { 1.8,  0.5},
924    { 1.9,  0.43},                                 742    { 1.9,  0.43},
925    { 2.0,  0.42},                                 743    { 2.0,  0.42},
926    { 2.1,  0.3},                                  744    { 2.1,  0.3},
927    { 2.4,  0.2},                                  745    { 2.4,  0.2},
928    { 3.0,  0.13},                                 746    { 3.0,  0.13},
929    { 3.08, 0.1},                                  747    { 3.08, 0.1},
930    { 3.1,  0.09},                                 748    { 3.1,  0.09},
931    { 3.3,  0.08},                                 749    { 3.3,  0.08},
932    { 3.5,  0.07},                                 750    { 3.5,  0.07},
933    { 3.8,  0.06},                                 751    { 3.8,  0.06},
934    { 4.0,  0.051},                                752    { 4.0,  0.051},
935    { 4.1,  0.04},                                 753    { 4.1,  0.04},
936    { 4.8,  0.03},                                 754    { 4.8,  0.03},
937    { 5.0,  0.024},                                755    { 5.0,  0.024},
938    { 5.1,  0.02},                                 756    { 5.1,  0.02},
939    { 6.0,  0.013},                                757    { 6.0,  0.013},
940    { 6.5,  0.01},                                 758    { 6.5,  0.01},
941    { 7.0,  0.009},                                759    { 7.0,  0.009},
942    { 7.1,  0.008},                                760    { 7.1,  0.008},
943    { 8.0,  0.006},                                761    { 8.0,  0.006},
944    { 9.0,  0.0032},                               762    { 9.0,  0.0032},
945    { 10.0, 0.0025} };                             763    { 10.0, 0.0025} };
946                                                   764 
947   sBarkasCorr = new G4PhysicsFreeVector(47, fa << 765   for(i=0; i<47; i++) {
948   for(i=0; i<47; ++i) { sBarkasCorr->PutValues << 766     engBarkas[i] = fTable[i][0];
                                                   >> 767     corBarkas[i] = fTable[i][1];
                                                   >> 768   }
                                                   >> 769 
                                                   >> 770   const G4double nuca[104][2] = {
                                                   >> 771   { 1.0E+8, 5.831E-8},
                                                   >> 772   { 8.0E+7, 7.288E-8},
                                                   >> 773   { 6.0E+7, 9.719E-8},
                                                   >> 774   { 5.0E+7, 1.166E-7},
                                                   >> 775   { 4.0E+7, 1.457E-7},
                                                   >> 776   { 3.0E+7, 1.942E-7},
                                                   >> 777   { 2.0E+7, 2.916E-7},
                                                   >> 778   { 1.5E+7, 3.887E-7},
                                                   >> 779 
                                                   >> 780   { 1.0E+7, 5.833E-7},
                                                   >> 781   { 8.0E+6, 7.287E-7},
                                                   >> 782   { 6.0E+6, 9.712E-7},
                                                   >> 783   { 5.0E+6, 1.166E-6},
                                                   >> 784   { 4.0E+6, 1.457E-6},
                                                   >> 785   { 3.0E+6, 1.941E-6},
                                                   >> 786   { 2.0E+6, 2.911E-6},
                                                   >> 787   { 1.5E+6, 3.878E-6},
                                                   >> 788 
                                                   >> 789   { 1.0E+6, 5.810E-6},
                                                   >> 790   { 8.0E+5, 7.262E-6},
                                                   >> 791   { 6.0E+5, 9.663E-6},
                                                   >> 792   { 5.0E+5, 1.157E-5},
                                                   >> 793   { 4.0E+5, 1.442E-5},
                                                   >> 794   { 3.0E+5, 1.913E-5},
                                                   >> 795   { 2.0E+5, 2.845E-5},
                                                   >> 796   { 1.5E+5, 3.762E-5},
                                                   >> 797 
                                                   >> 798   { 1.0E+5, 5.554E-5},
                                                   >> 799   { 8.0E+4, 6.866E-5},
                                                   >> 800   { 6.0E+4, 9.020E-5},
                                                   >> 801   { 5.0E+4, 1.070E-4},
                                                   >> 802   { 4.0E+4, 1.319E-4},
                                                   >> 803   { 3.0E+4, 1.722E-4},
                                                   >> 804   { 2.0E+4, 2.499E-4},
                                                   >> 805   { 1.5E+4, 3.248E-4},
                                                   >> 806 
                                                   >> 807   { 1.0E+4, 4.688E-4},
                                                   >> 808   { 8.0E+3, 5.729E-4},
                                                   >> 809   { 6.0E+3, 7.411E-4},
                                                   >> 810   { 5.0E+3, 8.718E-4},
                                                   >> 811   { 4.0E+3, 1.063E-3},
                                                   >> 812   { 3.0E+3, 1.370E-3},
                                                   >> 813   { 2.0E+3, 1.955E-3},
                                                   >> 814   { 1.5E+3, 2.511E-3},
                                                   >> 815 
                                                   >> 816   { 1.0E+3, 3.563E-3},
                                                   >> 817   { 8.0E+2, 4.314E-3},
                                                   >> 818   { 6.0E+2, 5.511E-3},
                                                   >> 819   { 5.0E+2, 6.430E-3},
                                                   >> 820   { 4.0E+2, 7.756E-3},
                                                   >> 821   { 3.0E+2, 9.855E-3},
                                                   >> 822   { 2.0E+2, 1.375E-2},
                                                   >> 823   { 1.5E+2, 1.736E-2},
                                                   >> 824 
                                                   >> 825   { 1.0E+2, 2.395E-2},
                                                   >> 826   { 8.0E+1, 2.850E-2},
                                                   >> 827   { 6.0E+1, 3.552E-2},
                                                   >> 828   { 5.0E+1, 4.073E-2},
                                                   >> 829   { 4.0E+1, 4.802E-2},
                                                   >> 830   { 3.0E+1, 5.904E-2},
                                                   >> 831   { 1.5E+1, 9.426E-2},
                                                   >> 832 
                                                   >> 833   { 1.0E+1, 1.210E-1},
                                                   >> 834   { 8.0E+0, 1.377E-1},
                                                   >> 835   { 6.0E+0, 1.611E-1},
                                                   >> 836   { 5.0E+0, 1.768E-1},
                                                   >> 837   { 4.0E+0, 1.968E-1},
                                                   >> 838   { 3.0E+0, 2.235E-1},
                                                   >> 839   { 2.0E+0, 2.613E-1},
                                                   >> 840   { 1.5E+0, 2.871E-1},
                                                   >> 841 
                                                   >> 842   { 1.0E+0, 3.199E-1},
                                                   >> 843   { 8.0E-1, 3.354E-1},
                                                   >> 844   { 6.0E-1, 3.523E-1},
                                                   >> 845   { 5.0E-1, 3.609E-1},
                                                   >> 846   { 4.0E-1, 3.693E-1},
                                                   >> 847   { 3.0E-1, 3.766E-1},
                                                   >> 848   { 2.0E-1, 3.803E-1},
                                                   >> 849   { 1.5E-1, 3.788E-1},
                                                   >> 850 
                                                   >> 851   { 1.0E-1, 3.711E-1},
                                                   >> 852   { 8.0E-2, 3.644E-1},
                                                   >> 853   { 6.0E-2, 3.530E-1},
                                                   >> 854   { 5.0E-2, 3.444E-1},
                                                   >> 855   { 4.0E-2, 3.323E-1},
                                                   >> 856   { 3.0E-2, 3.144E-1},
                                                   >> 857   { 2.0E-2, 2.854E-1},
                                                   >> 858   { 1.5E-2, 2.629E-1},
                                                   >> 859 
                                                   >> 860   { 1.0E-2, 2.298E-1},
                                                   >> 861   { 8.0E-3, 2.115E-1},
                                                   >> 862   { 6.0E-3, 1.883E-1},
                                                   >> 863   { 5.0E-3, 1.741E-1},
                                                   >> 864   { 4.0E-3, 1.574E-1},
                                                   >> 865   { 3.0E-3, 1.372E-1},
                                                   >> 866   { 2.0E-3, 1.116E-1},
                                                   >> 867   { 1.5E-3, 9.559E-2},
                                                   >> 868 
                                                   >> 869   { 1.0E-3, 7.601E-2},
                                                   >> 870   { 8.0E-4, 6.668E-2},
                                                   >> 871   { 6.0E-4, 5.605E-2},
                                                   >> 872   { 5.0E-4, 5.008E-2},
                                                   >> 873   { 4.0E-4, 4.352E-2},
                                                   >> 874   { 3.0E-4, 3.617E-2},
                                                   >> 875   { 2.0E-4, 2.768E-2},
                                                   >> 876   { 1.5E-4, 2.279E-2},
                                                   >> 877 
                                                   >> 878   { 1.0E-4, 1.723E-2},
                                                   >> 879   { 8.0E-5, 1.473E-2},
                                                   >> 880   { 6.0E-5, 1.200E-2},
                                                   >> 881   { 5.0E-5, 1.052E-2},
                                                   >> 882   { 4.0E-5, 8.950E-3},
                                                   >> 883   { 3.0E-5, 7.246E-3},
                                                   >> 884   { 2.0E-5, 5.358E-3},
                                                   >> 885   { 1.5E-5, 4.313E-3},
                                                   >> 886   { 0.0, 3.166E-3}
                                                   >> 887   };
949                                                   888 
950   const G4double SK[20] = {1.9477, 1.9232, 1.8 << 889   for(i=0; i<104; i++) {
                                                   >> 890     e[i] = nuca[i][0];
                                                   >> 891     a[i] = nuca[i][1];
                                                   >> 892   }
                                                   >> 893 
                                                   >> 894   // Constants
                                                   >> 895   theZieglerFactor = eV*cm2*1.0e-15 ; 
                                                   >> 896   alpha2           = fine_structure_const*fine_structure_const;
                                                   >> 897   lossFlucFlag     = true;
                                                   >> 898 
                                                   >> 899   // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111.
                                                   >> 900   // "Shell corrections for K- and L- electrons
                                                   >> 901   nK = 20;
                                                   >> 902   nL = 26;
                                                   >> 903   nEtaK = 29;
                                                   >> 904   nEtaL = 28;
                                                   >> 905 
                                                   >> 906   const G4double d[11] = {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15};
                                                   >> 907   const G4double thek[20] = {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78,
                                                   >> 908                              0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95};
                                                   >> 909   const G4double sk[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137,
951                            1.7754, 1.7396, 1.7    910                            1.7754, 1.7396, 1.7223, 1.7063, 1.6752,
952                            1.6461, 1.6189, 1.5    911                            1.6461, 1.6189, 1.5933, 1.5811, 1.5693,
953                            1.5467, 1.5254, 1.5    912                            1.5467, 1.5254, 1.5053, 1.4863, 1.4772};
954   const G4double TK[20] = {2.5222, 2.5125, 2.5 << 913   const G4double tk[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608,
955                            2.4388, 2.4163, 2.4    914                            2.4388, 2.4163, 2.4044, 2.3933, 2.3701,
956                            2.3466, 2.3229, 2.2    915                            2.3466, 2.3229, 2.2992, 2.2872, 2.2753,
957                            2.2515, 2.2277, 2.2    916                            2.2515, 2.2277, 2.2040, 2.1804, 2.1686};
                                                   >> 917   const G4double uk[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662,
                                                   >> 918                            2.0817, 2.0945, 2.0999, 2.1049, 2.1132,
                                                   >> 919                            2.1197, 2.1246, 2.1280, 2.1292, 2.1301,
                                                   >> 920                            2.1310, 2.1310, 2.1300, 2.1283, 2.1271};
                                                   >> 921   const G4double vk[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247,
                                                   >> 922                            8.3219, 8.3201, 8.3194, 8.3191, 8.3188,
                                                   >> 923                            8.3191, 8.3199, 8.3211, 8.3218, 8.3226,
                                                   >> 924                            8.3244, 8.3264, 8.3285, 8.3308, 8.3320};
                                                   >> 925 
                                                   >> 926   for(i=0; i<11; i++) { ZD[i] = d[i];}
                                                   >> 927 
                                                   >> 928   for(i=0; i<nK; i++) {
                                                   >> 929     TheK[i] = thek[i];
                                                   >> 930     SK[i]   = sk[i];
                                                   >> 931     TK[i]   = tk[i];
                                                   >> 932     UK[i]   = uk[i];
                                                   >> 933     VK[i]   = vk[i];
                                                   >> 934   }
958                                                   935 
959   const G4double SL[26] = {15.3343, 13.9389, 1 << 936   const G4double thel[26] = {0.24, 0.26, 0.28, 0.30, 0.32,   0.34, 0.35, 0.36, 0.38, 0.40,
                                                   >> 937                              0.42, 0.44, 0.45, 0.46, 0.48,   0.50, 0.52, 0.54, 0.55, 0.56,
                                                   >> 938                              0.58, 0.60, 0.62, 0.64, 0.65,   0.66};
                                                   >> 939   const G4double sl[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283,
960                            10.3424, 10.0371,      940                            10.3424, 10.0371,  9.7537,  9.2443,  8.8005,
961                            8.4114,  8.0683,    << 941          8.4114,  8.0683,   7.9117, 7.7641, 7.4931,
962                            7.2506,  7.0327,       942                            7.2506,  7.0327,   6.8362, 6.7452, 6.6584,
963                            6.4969,  6.3498,       943                            6.4969,  6.3498,   6.2154, 6.0923, 6.0345, 5.9792};
964   const G4double TL[26] = {35.0669, 33.4344, 3 << 944   const G4double tl[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226,
965                            28.6128, 28.1449, 2 << 945                            28.6128, 28.4149, 27.6991, 26.8674, 26.1061,
966                            25.4058, 24.7587, 2    946                            25.4058, 24.7587, 24.4531, 24.1583, 23.5992,
967                            23.0771, 22.5880, 2    947                            23.0771, 22.5880, 22.1285, 21.9090, 21.6958,
968                            21.2872, 20.9006, 2    948                            21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546};
                                                   >> 949   const G4double ul[26] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828,
                                                   >> 950                            1.4379, 1.5032, 1.5617, 1.6608, 1.7401,
                                                   >> 951                            1.8036, 1.8543, 1.8756, 1.8945, 1.9262,
                                                   >> 952                            1.9508, 1.9696, 1.9836, 1.9890, 1.9935,
                                                   >> 953                            2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028};  
                                                   >> 954   for(i=0; i<nL; i++) {
                                                   >> 955     TheL[i] = thel[i];
                                                   >> 956     SL[i]   = sl[i];
                                                   >> 957     TL[i]   = tl[i];
                                                   >> 958     UL[i]   = ul[i];
                                                   >> 959   }
                                                   >> 960 
                                                   >> 961   const G4double eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02,
                                                   >> 962                               0.03,  0.04,  0.05, 0.06,  0.08,
                                                   >> 963                               0.1,   0.15,  0.2,  0.3,   0.4,
                                                   >> 964                               0.5,   0.6,   0.7,  0.8,   1.0,
                                                   >> 965                               1.2,   1.4,   1.5,  1.7,   2.0, 3.0, 5.0, 7.0, 10.};
969                                                   966 
970   const G4double bk1[29][11] = {                  967   const G4double bk1[29][11] = { 
971   {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8,     968   {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}, 
972   {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8,     969   {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}, 
973   {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5    970   {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}, 
974   {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6,     971   {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}, 
975   {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1    972   {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}, 
976   {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7    973   {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}, 
977   {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2    974   {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}, 
978   {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6    975   {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}, 
979   {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1    976   {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}, 
980   {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3    977   {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}, 
981   {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7.    978   {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}, 
982   {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2    979   {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}, 
983   {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5.    980   {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}, 
984   {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1.    981   {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}, 
985   {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2.    982   {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}, 
986   {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4.    983   {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}, 
987   {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5.    984   {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}, 
988   {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7.    985   {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}, 
989   {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8.    986   {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}, 
990   {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.    987   {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.21317, 1.23921, 1.25257, 1.26616, 1.29408, 1.32303}, 
991   {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.    988   {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.49875, 1.52842, 1.54364, 1.55913, 1.59095, 1.62396}, 
992   {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.    989   {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.76244, 1.79526, 1.81210, 1.82925, 1.86448, 1.90104}, 
993   {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.    990   {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.88617, 1.92042, 1.93800, 1.95590, 1.99269, 2.03087}, 
994   {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.    991   {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.11827, 2.15555, 2.17409, 2.19337, 2.23302, 2.27419}, 
995   {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.    992   {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.42395, 2.46429, 2.48401, 2.50612, 2.54955, 2.59468}, 
996   {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.    993   {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.24748, 3.29664, 3.32192, 3.34770, 3.40081, 3.45611}, 
997   {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.    994   {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.27010, 4.33056, 4.36172, 4.39353, 4.45918, 4.52772}, 
998   {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.    995   {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.91402, 4.98200, 5.01707, 5.05290, 5.12695, 5.20436}, 
999   {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5    996   {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359}
1000   };                                             997   }; 
1001                                                  998 
1002   const G4double bk2[29][11] = {                 999   const G4double bk2[29][11] = { 
1003   {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8,    1000   {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}, 
1004   {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7,    1001   {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}, 
1005   {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6,     1002   {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}, 
1006   {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5,    1003   {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}, 
1007   {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5,     1004   {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}, 
1008   {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4,     1005   {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}, 
1009   {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4,     1006   {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}, 
1010   {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3,     1007   {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}, 
1011   {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3,     1008   {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}, 
1012   {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3,     1009   {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}, 
1013   {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1    1010   {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}, 
1014   {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2,     1011   {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}, 
1015   {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9    1012   {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}, 
1016   {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2    1013   {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}, 
1017   {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3    1014   {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}, 
1018   {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5    1015   {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}, 
1019   {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7    1016   {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}, 
1020   {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9    1017   {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}, 
1021   {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1    1018   {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1.10667, 1.13615, 1.16692, 1.19907, 1.21570, 1.23272}, 
1022   {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1    1019   {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1.45057, 1.48582, 1.52263, 1.56111, 1.58102, 1.60140}, 
1023   {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1    1020   {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1.76954, 1.80983, 1.85192, 1.89596, 1.91876, 1.94211}, 
1024   {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2    1021   {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2.06251, 2.10727, 2.15406, 2.20304, 2.22842, 2.25442}, 
1025   {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2    1022   {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2.19963, 2.24644, 2.29539, 2.34666, 2.37323, 2.40045}, 
1026   {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2    1023   {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2.45641, 2.50703, 2.56000, 2.61552, 2.64430, 2.67381}, 
1027   {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2    1024   {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2.79481, 2.85052, 2.90887, 2.97009, 3.00185, 3.03442}, 
1028   {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3    1025   {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3.70268, 3.77170, 3.84418, 3.92040, 3.96003, 4.00073}, 
1029   {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4    1026   {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4.83543, 4.92217, 5.01353, 5.10992, 5.16014, 5.21181}, 
1030   {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5    1027   {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5.55344, 5.65226, 5.79496, 5.90517, 5.96269, 6.02191}, 
1031   {10.0, 5.98474, 6.08046, 6.13015, 6.18112,     1028   {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535}
1032   };                                             1029   }; 
1033                                                  1030 
1034   const G4double bls1[28][10] = {                1031   const G4double bls1[28][10] = { 
1035   {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3.    1032   {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}, 
1036   {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7.    1033   {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}, 
1037   {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7    1034   {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}, 
1038   {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4.    1035   {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}, 
1039   {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0    1036   {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}, 
1040   {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0    1037   {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}, 
1041   {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5    1038   {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}, 
1042   {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8    1039   {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}, 
1043   {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1    1040   {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}, 
1044   {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2    1041   {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}, 
1045   {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.37    1042   {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}, 
1046   {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6    1043   {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}, 
1047   {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.00    1044   {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.0037, 1.0372, 1.0723, 1.1092, 1.1284, 1.1480}, 
1048   {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.643    1045   {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.6438, 1.6899, 1.7382, 1.7632, 1.7889}, 
1049   {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.165    1046   {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.1659, 2.2199, 2.2765, 2.3059, 2.3360}, 
1050   {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.558    1047   {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.5587, 2.6190, 2.6821, 2.7148, 2.7484}, 
1051   {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.888    1048   {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.8884, 2.9538, 3.0222, 3.0577, 3.0941}, 
1052   {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.165    1049   {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.1652, 3.2349, 3.3079, 3.3457, 3.3845}, 
1053   {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.423    1050   {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.4232, 3.4965, 3.5734, 3.6133, 3.6542}, 
1054   {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.886    1051   {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.8866, 3.9661, 4.0495, 4.0928, 4.1371}, 
1055   {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.213    1052   {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.2130, 4.2974, 4.3861, 4.4321, 4.4794}, 
1056   {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.517    1053   {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.5170, 4.6056, 4.6988, 4.7471, 4.7968}, 
1057   {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.652    1054   {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.6521, 4.7426, 4.8378, 4.8872, 4.9379}, 
1058   {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.894    1055   {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.8949, 4.9888, 5.0876, 5.1388, 5.1915}, 
1059   {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.205    1056   {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.2053, 5.3036, 5.4070, 5.4607, 5.5159}, 
1060   {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.948    1057   {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.9489, 6.0581, 6.1730, 6.2328, 6.2943}, 
1061   {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.800    1058   {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.8003, 6.9231, 7.0525, 7.1199, 7.1892}, 
1062   {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.220    1059   {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384}
1063   };                                             1060   };
1064                                                  1061  
                                                   >> 1062 
1065   const G4double bls2[28][10] = {                1063   const G4double bls2[28][10] = { 
1066   {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7.    1064   {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}, 
1067   {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1.    1065   {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}, 
1068   {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4    1066   {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}, 
1069   {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8.    1067   {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}, 
1070   {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7    1068   {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}, 
1071   {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6    1069   {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}, 
1072   {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1    1070   {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}, 
1073   {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4    1071   {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}, 
1074   {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0    1072   {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}, 
1075   {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5    1073   {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}, 
1076   {0.1,  4.3613E-1, 4.5956E-1,  4.852E-1, 5.1    1074   {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}, 
1077   {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1    1075   {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}, 
1078   {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.350    1076   {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.3506, 1.3765, 1.4308, 1.4886, 1.5504}, 
1079   {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.052    1077   {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.0528, 2.0864, 2.1569, 2.2319, 2.3120}, 
1080   {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.645    1078   {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.6452, 2.6847, 2.7674, 2.8554, 2.9494}, 
1081   {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.093    1079   {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.0937, 3.1378, 3.2301, 3.3285, 3.4337}, 
1082   {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.469    1080   {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.4691, 3.5171, 3.6175, 3.7246, 3.8391}, 
1083   {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.784    1081   {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.7848, 3.8360, 3.9433, 4.0578, 4.1804}, 
1084   {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.076    1082   {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.0764, 4.1304, 4.2438, 4.3648, 4.4944}, 
1085   {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.596    1083   {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.5962, 4.6551, 4.7786, 4.9106, 5.0520}, 
1086   {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.968    1084   {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.9685, 5.0314, 5.1633, 5.3043, 5.4555}, 
1087   {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.311    1085   {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.3114, 5.3776, 5.5166, 5.6653, 5.8249}, 
1088   {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.464    1086   {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.4640, 5.5317, 5.6739, 5.8260, 5.9893}, 
1089   {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.737    1087   {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.7373, 5.8088, 5.9568, 6.1152, 6.2853}, 
1090   {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.089    1088   {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.0896, 6.1636, 6.3192, 6.4857, 6.6647}, 
1091   {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.935    1089   {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.9351, 7.0180, 7.1925, 7.3795, 7.5808}, 
1092   {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.914    1090   {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.9146, 8.0087, 8.2071, 8.4200, 8.6496}, 
1093   {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.419    1091   {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133}
1094   };                                             1092   }; 
1095                                                  1093  
1096   const G4double bls3[28][9] = {                 1094   const G4double bls3[28][9] = { 
1097   {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1.    1095   {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}, 
1098   {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3.    1096   {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}, 
1099   {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3    1097   {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}, 
1100   {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2    1098   {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}, 
1101   {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2    1099   {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}, 
1102   {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0    1100   {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}, 
1103   {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7    1101   {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}, 
1104   {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6    1102   {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}, 
1105   {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5    1103   {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}, 
1106   {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5    1104   {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}, 
1107   {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.61    1105   {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}, 
1108   {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.34    1106   {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.3424, 1.4163, 1.4974, 1.5868}, 
1109   {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.848    1107   {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.8484, 1.9394, 2.0390, 2.1489}, 
1110   {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.698    1108   {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.6980, 2.8159, 2.9451, 3.0876}, 
1111   {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.403    1109   {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.4030, 3.5416, 3.6938, 3.8620}, 
1112   {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.942    1110   {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.9421, 4.0978, 4.2688, 4.4580}, 
1113   {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.393    1111   {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.3935, 4.5636, 4.7506, 4.9576}, 
1114   {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.764    1112   {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.7644, 4.9470, 5.1478, 5.3703}, 
1115   {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.123    1113   {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.1233, 5.3169, 5.5300, 5.7661}, 
1116   {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.740    1114   {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.7400, 5.9523, 6.1863, 6.4458}, 
1117   {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.192    1115   {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.1927, 6.4206, 6.6719, 6.9510}, 
1118   {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.603    1116   {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.6038, 6.8451, 7.1113, 7.4071}, 
1119   {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.787    1117   {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.7871, 7.0344, 7.3073, 7.6107}, 
1120   {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.117    1118   {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.1175, 7.3757, 7.6609, 7.9782}, 
1121   {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.541    1119   {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.5418, 7.8143, 8.1156, 8.4510}, 
1122   {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.570    1120   {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.5708, 8.8796, 9.2214, 9.6027}, 
1123   {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.783    1121   {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.7834, 10.1384, 10.5323, 10.9722}, 
1124   {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4    1122   {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356}
1125   };                                             1123   }; 
1126                                                  1124 
1127   const G4double bll1[28][10] = {                1125   const G4double bll1[28][10] = { 
1128   {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5.    1126   {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}, 
1129   {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2.    1127   {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}, 
1130   {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2    1128   {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}, 
1131   {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5.    1129   {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}, 
1132   {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4    1130   {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}, 
1133   {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7    1131   {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}, 
1134   {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7    1132   {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}, 
1135   {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6    1133   {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}, 
1136   {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3    1134   {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}, 
1137   {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0    1135   {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}, 
1138   {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.96    1136   {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}, 
1139   {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.10    1137   {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.1073, 1.1443, 1.1832, 1.2035, 1.2243}, 
1140   {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.562    1138   {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.5624, 1.6083, 1.6566, 1.6817, 1.7076}, 
1141   {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.332    1139   {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.3322, 2.3915, 2.4542, 2.4869, 2.5205}, 
1142   {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.937    1140   {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.9377, 3.0072, 3.0807, 3.1192, 3.1587}, 
1143   {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.412    1141   {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.4122, 3.4898, 3.5721, 3.6151, 3.6595}, 
1144   {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.830    1142   {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.8303, 3.9146, 4.0042, 4.0511, 4.0995}, 
1145   {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.152    1143   {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.1520, 4.2421, 4.3380, 4.3882, 4.4401}, 
1146   {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.433    1144   {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.4333, 4.5285, 4.6298, 4.6830, 4.7380}, 
1147   {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.910    1145   {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.9106, 5.0144, 5.1250, 5.1831, 5.2432}, 
1148   {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.272    1146   {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.2729, 5.3837, 5.5050, 5.5642, 5.6287}, 
1149   {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.578    1147   {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.5782, 5.6950, 5.8198, 5.8855, 5.9554}, 
1150   {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.712    1148   {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.7128, 5.8323, 5.9601, 6.0274, 6.0972}, 
1151   {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.954    1149   {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.9541, 6.0785, 6.2116, 6.2818, 6.3546}, 
1152   {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.261    1150   {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.2618, 6.3925, 6.5327, 6.6066, 6.6833}, 
1153   {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.006    1151   {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.0062, 7.1529, 7.3107, 7.3941, 7.4807}, 
1154   {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.898    1152   {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.8982, 8.0653, 8.2454, 8.3409, 8.4402}, 
1155   {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.206    1153   {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927}
1156   };                                             1154   }; 
1157                                                  1155 
1158   const G4double bll2[28][10] = {                1156   const G4double bll2[28][10] = { 
1159   {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3.    1157   {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}, 
1160   {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1.    1158   {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}, 
1161   {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8    1159   {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}, 
1162   {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1.    1160   {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}, 
1163   {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6    1161   {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}, 
1164   {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5    1162   {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}, 
1165   {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7    1163   {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}, 
1166   {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7    1164   {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}, 
1167   {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7    1165   {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}, 
1168   {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0    1166   {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}, 
1169   {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.35    1167   {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}, 
1170   {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.44    1168   {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.4421, 1.4703, 1.5300, 1.5942, 1.6636}, 
1171   {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.978    1169   {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.9788, 2.0142, 2.0889, 2.1694, 2.2567}, 
1172   {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.876    1170   {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.8762, 2.9231, 3.0222, 3.1295, 3.2463}, 
1173   {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.580    1171   {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.5803, 3.6361, 3.7546, 3.8835, 4.0242}, 
1174   {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.135    1172   {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.1350, 4.1983, 4.3330, 4.4799, 4.6408}, 
1175   {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.620    1173   {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.6206, 4.6904, 4.8390, 5.0013, 5.1796}, 
1176   {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.000    1174   {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.0009, 5.0762, 5.2370, 5.4129, 5.6066}, 
1177   {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.333    1175   {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.3338, 5.4141, 5.5857, 5.7738, 5.9811}, 
1178   {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.898    1176   {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.8989, 5.9878, 6.1780, 6.3870, 6.6179}, 
1179   {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.334    1177   {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.3343, 6.4303, 6.6361, 6.8626, 7.1137}, 
1180   {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.702    1178   {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.7021, 6.8043, 7.0237, 7.2655, 7.5338}, 
1181   {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.865    1179   {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.8650, 6.9700, 7.1954, 7.4442, 7.7203}, 
1182   {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.157    1180   {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.1578, 7.2679, 7.5045, 7.7660, 8.0565}, 
1183   {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.532    1181   {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.5328, 7.6496, 7.9010, 8.1791, 8.4886}, 
1184   {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.447    1182   {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.4475, 8.5814, 8.8702, 9.1908, 9.5488}, 
1185   {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.557    1183   {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.5572, 9.7132, 10.0504, 10.4259, 10.8466}, 
1186   {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.00    1184   {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250}
1187   };                                             1185   }; 
1188                                                  1186 
1189   const G4double bll3[28][9] = {                 1187   const G4double bll3[28][9] = { 
1190   {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2.    1188   {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}, 
1191   {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6.    1189   {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}, 
1192   {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6    1190   {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}, 
1193   {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4.    1191   {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}, 
1194   {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1    1192   {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}, 
1195   {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7    1193   {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}, 
1196   {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0    1194   {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}, 
1197   {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3    1195   {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}, 
1198   {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7    1196   {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}, 
1199   {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7    1197   {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7836E-1, 9.3355E-1, 9.9462E-1, 1.0627, 1.1394}, 
1200   {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.250    1198   {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.2503, 1.3265, 1.4116, 1.5076}, 
1201   {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.01    1199   {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.0108, 2.1217, 2.2462, 2.3876}, 
1202   {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.697    1200   {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.6971, 2.8391, 2.9994, 3.1822}, 
1203   {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.844    1201   {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.8445, 4.0404, 4.2631, 4.5193}, 
1204   {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.753    1202   {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.7530, 4.9944, 5.2703, 5.5895}, 
1205   {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.481    1203   {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.4811, 5.7618, 6.0840, 6.4583}, 
1206   {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.117    1204   {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.1173, 6.4326, 6.7958, 7.2191}, 
1207   {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.630    1205   {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.6306, 6.9769, 7.3767, 7.8440}, 
1208   {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.082    1206   {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.0827, 7.4570, 7.8900, 8.3972}, 
1209   {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.854    1207   {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.8546, 8.2778, 8.7690, 9.3464}, 
1210   {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.465    1208   {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.4652, 8.9302, 9.4713, 10.1090}, 
1211   {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.985    1209   {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.9851, 9.4866, 10.0713, 10.7619}, 
1212   {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.218    1210   {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.2181, 9.7352, 10.3399, 11.0546}, 
1213   {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.636    1211   {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.6367, 10.1856, 10.8270, 11.5863}, 
1214   {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.17    1212   {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.1793, 10.7688, 11.4590, 12.2775}, 
1215   {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11    1213   {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11.5229, 12.2172, 13.0332, 14.0048}, 
1216   {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 1    1214   {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 13.1898, 14.0213, 15.0024, 16.1752}, 
1217   {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 1    1215   {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420}
1218   };                                             1216   };
                                                   >> 1217 
1219                                                  1218                              
1220   G4double b, bs;                                1219   G4double b, bs; 
1221   for(i=0; i<nEtaK; ++i) {                    << 1220   for(i=0; i<nEtaK; i++) {
1222                                                  1221 
1223     G4double et = Eta[i];                     << 1222     G4double et = eta[i];
1224     G4double loget = G4Log(et);               << 1223     G4double loget = std::log(et);
                                                   >> 1224     Eta[i] = et;
                                                   >> 1225     // G4cout << "### eta[" << i << "]= " << et << "  KShell: tet= " << TheK[0]<<" - " <<TheK[nK-1]<< G4endl;
1225                                                  1226 
1226     for(j=0; j<nK; ++j) {                     << 1227     for(j=0; j<nK; j++) {
1227                                                  1228 
1228       if(j < 10) { b = bk2[i][10-j]; }        << 1229       if(j < 10) b = bk2[i][10-j];
1229       else       { b = bk1[i][20-j]; }        << 1230       else       b = bk1[i][20-j];
1230                                                  1231 
1231       CK[j][i] = SK[j]*loget + TK[j] - b;        1232       CK[j][i] = SK[j]*loget + TK[j] - b;
1232                                               << 1233       //G4cout << " " << CK[j][i];
1233       if(i == nEtaK-1) {                      << 1234       if(i == nEtaK-1) ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]); 
1234         ZK[j] = et*(et*et*CK[j][i] - et*UK[j] << 
1235         //G4cout << "i= " << i << " j= " << j << 
1236         //       << " CK[j][i]= " <<  CK[j][i << 
1237         //       << " ZK[j]= " << ZK[j] << "  << 
1238       }                                       << 
1239     }                                            1235     }
1240     //G4cout << G4endl;                          1236     //G4cout << G4endl;
1241     if(i < nEtaL) {                              1237     if(i < nEtaL) {
1242       //G4cout << "  LShell:" <<G4endl;          1238       //G4cout << "  LShell:" <<G4endl;
1243       for(j=0; j<nL; ++j) {                   << 1239       for(j=0; j<nL; j++) {
1244                                                  1240 
1245         if(j < 8) {                              1241         if(j < 8) {
1246           bs = bls3[i][8-j];                     1242           bs = bls3[i][8-j];
1247           b  = bll3[i][8-j];                     1243           b  = bll3[i][8-j];
1248         } else if(j < 17) {                      1244         } else if(j < 17) {
1249           bs = bls2[i][17-j];                    1245           bs = bls2[i][17-j];
1250           b  = bll2[i][17-j];                    1246           b  = bll2[i][17-j];
1251         } else {                                 1247         } else {
1252           bs = bls1[i][26-j];                    1248           bs = bls1[i][26-j];
1253           b  = bll1[i][26-j];                    1249           b  = bll1[i][26-j];
1254         }                                     << 1250   }
1255         G4double c = SL[j]*loget + TL[j];     << 1251   G4double c = SL[j]*loget + TL[j]; 
1256         CL[j][i] = c - bs - 3.0*b;               1252         CL[j][i] = c - bs - 3.0*b;
1257         if(i == nEtaL-1) {                    << 1253         //G4cout << " " << CL[j][i];
1258           VL[j] = et*(et*CL[j][i] - UL[j]);   << 1254         if(i == nEtaL-1) VL[j] = et*(et*CL[j][i] - UL[j]); 
1259           //G4cout << "i= " << i << " j= " << << 
1260           //         << " CL[j][i]= " <<  CL[ << 
1261           //         << " VL[j]= " << VL[j] < << 
1262           //         << " et= " << et << G4en << 
1263           //" UL= " << UL[j] << " TL= " << TL << 
1264         }                                     << 
1265       }                                          1255       }
1266       //G4cout << G4endl;                        1256       //G4cout << G4endl;
1267     }                                            1257     }
1268   }                                              1258   }
                                                   >> 1259   const G4double hm[53] = {
                                                   >> 1260     12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4,
                                                   >> 1261     10.0,  9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87, 
                                                   >> 1262     5.61,  5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38, 
                                                   >> 1263     4.32,  4.26, 4.20, 4.15, 4.1,  4.04, 4.00, 3.95, 3.93, 3.91, 
                                                   >> 1264     3.90,  3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89, 
                                                   >> 1265     3.90, 3.92, 3.93 };
                                                   >> 1266   const G4double hn[31] = {
                                                   >> 1267     75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8, 
                                                   >> 1268     24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7, 
                                                   >> 1269     19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4, 
                                                   >> 1270     18.2
                                                   >> 1271   };
                                                   >> 1272   for(i=0; i<53; i++) {HM[i] = hm[i];}
                                                   >> 1273   for(i=0; i<31; i++) {HN[i] = hn[i];}
                                                   >> 1274 
                                                   >> 1275   const G4double mm[93] = {
                                                   >> 1276      0.0,
                                                   >> 1277      /*
                                                   >> 1278  -0.0001577, 5.396e-05, 0.0004194, -0.0001969, -0.0003447, -0.0001904, 9.612e-05, 4.16e-05, 0.0001899, 0.0001322,
                                                   >> 1279  0.0001485, 4.698e-05, -7.726e-05, -6.448e-05, 3.269e-05, -0.0001293, 0.0001483, -8.559e-05, 4.018e-05, 6.239e-05,
                                                   >> 1280  4.447e-05, 2.212e-05, -5.568e-06, 3.782e-05, 4.398e-05, 8.942e-05, 0.0002202, 0.0001357, 0.0001195, 0.0001812,
                                                   >> 1281  0.0002365, 4.508e-05, 0.0001629, 7.75e-05, 0.0001886, 0.0001601, 0.00015, 5.679e-05, 5.956e-06, -7.309e-05,
                                                   >> 1282  -0.0001144, -8.585e-05, -0.0001051, -8.209e-05, -5.871e-05, -3.744e-05, -6.707e-05, -9.199e-06, 1.181e-05, 2.252e-05,
                                                   >> 1283  1.051e-05, 7.63e-05, 1.71e-05, 1.705e-05, 1.561e-05, 8.987e-06, 4.957e-06, -1.473e-05, -1.902e-05, -1.491e-05,
                                                   >> 1284  -9.891e-08, 2.584e-05, 3.696e-05, 8.651e-06, -3.601e-06, 1.647e-05, 4.536e-05, 9.954e-05, 8.922e-05, 0.0001026,
                                                   >> 1285  4.094e-05, 3.788e-05, 1.578e-05, 2.731e-05, -6.035e-05, -9.112e-05, -8.846e-05, -0.0001361, -0.0001959, -0.0001883,
                                                   >> 1286  -0.0002854, -0.0004129, -0.0004611, -0.0005523, -0.0005869, -0.0005008, -0.000659, -0.0007045, -0.0007322, -0.0008374,
                                                   >> 1287  -0.0008731, -0.0009327,
                                                   >> 1288      
                                                   >> 1289  -0.0003427, -4.685e-05, 0.0007304, -0.0003158, -0.0004104, -0.001133, -0.0007987, -0.0001528, 8.976e-05, 8.63e-05,
                                                   >> 1290  1.764e-05, -0.000136, -0.0002895, -0.0002618, -0.0002878, -0.0003813, 1.417e-05, -0.0004216, -0.0003859, -0.0001028,
                                                   >> 1291  -0.0001617, -0.0002141, -0.0002178, 5.59e-05, 9.964e-07, 7.528e-05, 0.0002261, 4.241e-05, 0.0005255, 0.0002139,
                                                   >> 1292  -5.821e-05, -2.661e-05, 0.0001039, 3.359e-05, 0.0004229, -0.0004872, 0.0001374, -2.79e-05, -1.319e-05, -0.0001077,
                                                   >> 1293  -0.0001923, -0.0001189, -0.0001091, -8.486e-05, 7.506e-05, 0.0001053, 0.0004596, 0.000177, 4.805e-05, 3.97e-05,
                                                   >> 1294  0.0004008, 0.0002194, 0.0001639, 0.0003285, 0.000219, 0.0001339, 6.323e-05, 2.013e-05, 2.568e-05, 4.346e-05,
                                                   >> 1295  8.455e-05, 8.785e-05, 0.0001393, 9.907e-05, 0.0001003, 0.0001508, 0.0002189, 0.0003885, 0.0003603, 0.0003839,
                                                   >> 1296  0.0003321, 0.0002207, 0.0001627, 0.0002282, 0.0002177, 0.0002156, 0.0002855, 0.0002361, 0.0007735, 0.0001157,
                                                   >> 1297  -0.0001214, -0.0002944, -0.0003193, -0.0004025, -0.0002436, -7.573e-05, -0.0003741, -0.0003678, -0.0004839, -0.0003934,
                                                   >> 1298  -0.0005817, -0.000726,
                                                   >> 1299 
                                                   >> 1300  -118, 40.37, 313.8, -147.3, -257.9, -142.4, 71.91, 31.12, 142.1, 98.92,
                                                   >> 1301  112.3, 41.36, -47.15, -17.84, 46.46, -66.91, 141.3, -19.96, 84.66, 110.2,
                                                   >> 1302  113.2, 108.7, 92.89, 127.3, 139.3, 187, 288.5, 233.4, 225.3, 284.6,
                                                   >> 1303  344.5, 246.2, 359.7, 324.1, 455, 457.1, 481, 459.8, 447.9, 407.5,
                                                   >> 1304  399.7, 457.7, 470.4, 511.8, 545.9, 581, 574.1, 638.6, 685.6, 719.6,
                                                   >> 1305  740.8, 817.2, 808.3, 839.7, 866.2, 873.2, 879.3, 886.3, 902, 912.7,
                                                   >> 1306  927.7, 942.1, 967.1, 962.2, 961.8, 993.4, 1006, 1060, 1046, 1061,
                                                   >> 1307  1036, 1060, 1062, 1084, 1048, 1063, 1098, 1093, 1089, 1117,
                                                   >> 1308  1085, 1044, 1033, 1027, 1051, 1150, 1064, 1078, 1090, 1140,
                                                   >> 1309  1119, 1097,
                                                   >> 1310 */
                                                   >> 1311 
                                                   >> 1312  -511.3, -69.91, 1090, -471.2, -612.4, -1690, -1192, -228, 133.9, 128.8,
                                                   >> 1313  28.71, -190.6, -408.7, -326.6, -378.5, -507.3, 93.63, -531, -459.2, -21.83,
                                                   >> 1314  -74.93, -121.7, -122.8, 285.2, 217.4, 366.3, 598.7, 341.8, 1078, 634.2,
                                                   >> 1315  280.5, 430.7, 663.5, 596.4, 1307, 6.814, 992.6, 887, 967.6, 868,
                                                   >> 1316  775.9, 1034, 1050, 1204, 1458, 1541, 2186, 1806, 1676, 1718,
                                                   >> 1317  2416, 2204, 2197, 2496, 2371, 2370, 2330, 2318, 2379, 2434,
                                                   >> 1318  2513, 2533, 2650, 2650, 2684, 2788, 2868, 3157, 3150, 3196,
                                                   >> 1319  3162, 3055, 3008, 3133, 3176, 3212, 3356, 3324, 4170, 3320,
                                                   >> 1320  3109, 2941, 2973, 2858, 3088, 3447, 3134, 3121, 3086, 3199,
                                                   >> 1321  3115, 2979,
                                                   >> 1322 
                                                   >> 1323   };
                                                   >> 1324 
                                                   >> 1325   const G4double tau[93] = {
                                                   >> 1326      0.0,
                                                   >> 1327  -251.1, 512.8, 1724, 649, 658.8, 594.1, 946.6, 969.2, 1154, 997.8,
                                                   >> 1328  939.3, 775.5, 619.9, 848.3, 973.4, 1035, 1252, 698.8, 1043, 1080,
                                                   >> 1329  957, 885.6, 772.1, 847.1, 799.7, 901.7, 1042, 880.8, 857.7, 968,
                                                   >> 1330  850.4, 482.8, 781.4, 685, 1093, 847.4, 824.4, 751.5, 638, 489.2,
                                                   >> 1331  266.5, 371, 319.2, 356.6, 424.9, 405, 390.8, 610.6, 404.3, 444.6,
                                                   >> 1332  439, 579, 466.4, 520.5, 551.1, 441.6, 328.8, 264.3, 260.9, 274,
                                                   >> 1333  264.8, 235.1, 290, 196, 178.2, 248.4, 224.3, 268.2, 239.5, 253.9,
                                                   >> 1334  219.2, 172.7, 169.5, 199, 53.35, 65.78, 120.2, 77.13, 32.76, 45.53,
                                                   >> 1335  -148, -331.4, -386.4, -430.8, -361.1, -108.3, -382.4, -400.1, -458.4, -430.8,
                                                   >> 1336  -536.9, -627.6,
                                                   >> 1337   };
1269                                                  1338 
1270   const G4double xzk[34] = { 11.7711,         << 1339   for(i=0; i<93; i++) {
1271     13.3669, 15.5762, 17.1715, 18.7667, 20.85 << 1340     MSH[i] = mm[i];
1272     34.3457, 37.4119, 40.3555, 42.3177, 44.77 << 1341     TAU[i] = tau[i];
1273     60.9586, 63.6567, 66.5998, 68.807, 71.872 << 1342   }
1274                              93.4549, 96.2753 << 1343   const G4double coseb[14] = {0.0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
1275   const G4double yzk[34] = { 0.70663,         << 1344                               1.0,1.2,1.5,2.0};
1276     0.72033, 0.73651, 0.74647, 0.75518, 0.763 << 1345   const G4double cosxi[14] = {1.0000, 0.9905, 0.9631, 0.9208, 0.8680,
1277     0.80611, 0.8123, 0.8185, 0.82097, 0.82467 << 1346                               0.7478, 0.6303, 0.5290, 0.4471, 0.3323,
1278     0.84565, 0.84936, 0.85181, 0.85303, 0.855 << 1347                               0.2610, 0.2145, 0.1696, 0.1261};
1279                              0.86891, 0.87011 << 1348   for(i=0; i<14; i++) {
1280                                               << 1349     COSEB[i] = coseb[i];
1281   const G4double xzl[36] = { 15.5102,         << 1350     COSXI[i] = cosxi[i];
1282     16.7347, 17.9592, 19.551, 21.0204, 22.612 << 1351   }
1283     34.4898, 36.2041, 38.4082, 40.3674, 42.57 << 
1284     59.3469, 61.9184, 64.6122, 67.4286, 71.46 << 
1285                              91.551, 94.2449, << 
1286   const G4double yzl[36] = { 0.29875,         << 
1287     0.31746, 0.33368, 0.35239, 0.36985, 0.387 << 
1288     0.4896, 0.50083, 0.51331, 0.52328, 0.5307 << 
1289     0.58191, 0.5869, 0.59189, 0.60062, 0.6068 << 
1290                               0.6368, 0.64054 << 
1291                                               << 
1292   sThetaK = new G4PhysicsFreeVector(34, false << 
1293   for(i=0; i<34; ++i) { sThetaK->PutValues(i, << 
1294   sThetaL = new G4PhysicsFreeVector(36, false << 
1295   for(i=0; i<36; ++i) { sThetaL->PutValues(i, << 
1296 }                                                1352 }
1297                                                  1353 
1298 //....oooOO0OOooo........oooOO0OOooo........o    1354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1299                                                  1355 
1300                                                  1356