Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4hICRU49p.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/lowenergy/src/G4hICRU49p.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4hICRU49p.cc (Version 3.0)


                                                   >>   1 // This code implementation is the intellectual property of
                                                   >>   2 // the GEANT4 collaboration.
  1 //                                                  3 //
  2 // ******************************************* <<   4 // By copying, distributing or modifying the Program (or any work
  3 // * License and Disclaimer                    <<   5 // based on the Program) you indicate your acceptance of this statement,
  4 // *                                           <<   6 // and all its terms.
  5 // * The  Geant4 software  is  copyright of th << 
  6 // * the Geant4 Collaboration.  It is provided << 
  7 // * conditions of the Geant4 Software License << 
  8 // * LICENSE and available at  http://cern.ch/ << 
  9 // * include a list of copyright holders.      << 
 10 // *                                           << 
 11 // * Neither the authors of this software syst << 
 12 // * institutes,nor the agencies providing fin << 
 13 // * work  make  any representation or  warran << 
 14 // * regarding  this  software system or assum << 
 15 // * use.  Please see the license in the file  << 
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                           << 
 18 // * This  code  implementation is the result  << 
 19 // * technical work of the GEANT4 collaboratio << 
 20 // * By using,  copying,  modifying or  distri << 
 21 // * any work based  on the software)  you  ag << 
 22 // * use  in  resulting  scientific  publicati << 
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // ******************************************* << 
 25 //                                             << 
 26 //                                                  7 //
 27 // -------------------------------------------      8 // -------------------------------------------------------------------
 28 //                                                  9 //
 29 // GEANT4 Class file                               10 // GEANT4 Class file
 30 //                                                 11 //
                                                   >>  12 // For information related to this code contact:
                                                   >>  13 // Geant4 Collaboration
 31 //                                                 14 //
 32 // File name:     G4hICRU49p                       15 // File name:     G4hICRU49p
 33 //                                                 16 //
 34 // Author:        V.Ivanchenko (Vladimir.Ivanc     17 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
 35 //                                                 18 // 
 36 // Creation date: 20 July 2000                     19 // Creation date: 20 July 2000
 37 //                                                 20 //
 38 // Modifications:                                  21 // Modifications: 
 39 // 20/07/2000  V.Ivanchenko First implementati     22 // 20/07/2000  V.Ivanchenko First implementation
 40 // 18/09/2000  V.Ivanchenko clean up - all var     23 // 18/09/2000  V.Ivanchenko clean up - all variable are the same as in ICRU
 41 // 03/10/2000  V.Ivanchenko clean up accoding      24 // 03/10/2000  V.Ivanchenko clean up accoding to CodeWizard
 42 // 10/05/2001  V.Ivanchenko Clean up againist  << 
 43 //                                                 25 //
 44 // Class Description:                              26 // Class Description: 
 45 //                                                 27 //
 46 // Electronic stopping power parametrised acco     28 // Electronic stopping power parametrised according to
 47 // ICRU Report N49, 1993, for protons.             29 // ICRU Report N49, 1993, for protons.
 48 //                                                 30 //
 49 // Class Description: End                          31 // Class Description: End 
 50 //                                                 32 //
 51 // -------------------------------------------     33 // -------------------------------------------------------------------
 52 //                                                 34 //
 53 //....oooOO0OOooo........oooOO0OOooo........oo     35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 54                                                    36 
 55 #include "G4hICRU49p.hh"                           37 #include "G4hICRU49p.hh"
 56                                                << 
 57 #include "globals.hh"                          << 
 58 #include "G4SystemOfUnits.hh"                  << 
 59 #include "G4UnitsTable.hh"                         38 #include "G4UnitsTable.hh"
                                                   >>  39 #include "globals.hh"
 60 #include "G4Material.hh"                           40 #include "G4Material.hh"
 61                                                    41 
 62 //....oooOO0OOooo........oooOO0OOooo........oo     42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63                                                    43 
 64 G4hICRU49p::G4hICRU49p():                      <<  44 G4hICRU49p::G4hICRU49p():G4VhElectronicStoppingPower(), 
 65   G4VhElectronicStoppingPower(),               << 
 66   protonMassAMU(1.007276),                         45   protonMassAMU(1.007276),
 67   iMolecula(0)                                     46   iMolecula(0)
 68 {;}                                                47 {;}
 69                                                    48 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 71                                                    50 
 72 G4hICRU49p::~G4hICRU49p()                          51 G4hICRU49p::~G4hICRU49p() 
 73 {;}                                                52 {;}
 74                                                    53 
 75 //....oooOO0OOooo........oooOO0OOooo........oo     54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 76                                                    55 
 77 G4bool G4hICRU49p::HasMaterial(const G4Materia     56 G4bool G4hICRU49p::HasMaterial(const G4Material* material) 
 78 {                                                  57 {
 79   G4String chFormula = material->GetChemicalFo     58   G4String chFormula = material->GetChemicalFormula() ;
 80   G4String myFormula = G4String(" ") ;             59   G4String myFormula = G4String(" ") ;
 81                                                    60 
 82   if (myFormula ==  chFormula ) {                  61   if (myFormula ==  chFormula ) {
 83     if(1 == (material->GetNumberOfElements()))     62     if(1 == (material->GetNumberOfElements())) return true;
 84     return false ;                                 63     return false ;
 85   }                                                64   }
 86                                                    65 
 87   // ICRU Report N49, 1993. Power's model for      66   // ICRU Report N49, 1993. Power's model for He.
 88   const G4int numberOfMolecula = 11 ;          <<  67   const size_t numberOfMolecula = 11 ;    
 89   static const G4String name[numberOfMolecula] <<  68   static G4String name[numberOfMolecula] = {
 90     "Al_2O_3",                 "CO_2",         <<  69     "Al_2O_3",                 "CO_2",                      "CH_4",  
 91     "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Pol     70     "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene",  "(C_8H_8)_N",  
 92     "C_3H_8",          "SiO_2",                <<  71     "C_3H_8",                  "SiO_2",                     "H_2O",
 93     "H_2O-Gas",                "Graphite"};    <<  72     "H_2O-Gas",                "Graphite" } ;      
 94                                                    73   
 95   // Special treatment for water in gas state      74   // Special treatment for water in gas state
 96   const G4State theState = material->GetState(     75   const G4State theState = material->GetState() ;
 97                                                    76 
 98   myFormula = G4String("H_2O");                    77   myFormula = G4String("H_2O");
 99   if( theState == kStateGas && myFormula == ch     78   if( theState == kStateGas && myFormula == chFormula) {
100     chFormula = G4String("H_2O-Gas");              79     chFormula = G4String("H_2O-Gas");
101   }                                                80   }
                                                   >>  81 
                                                   >>  82   G4int i;
102                                                    83   
103   // Search for the material in the table          84   // Search for the material in the table
104   for (G4int i=0; i<numberOfMolecula; ++i) {   <<  85   for (i=0; i<numberOfMolecula; i++) {
105       if (chFormula == name[i]) {                  86       if (chFormula == name[i]) {
106         SetMoleculaNumber(i) ;                     87         SetMoleculaNumber(i) ;    
107   return true ;                                    88   return true ;
108       }                                            89       }
109   }                                                90   }
110   return false ;                                   91   return false ;
111 }                                                  92 }
112                                                    93 
113 //....oooOO0OOooo........oooOO0OOooo........oo     94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114                                                    95 
115 G4double G4hICRU49p::StoppingPower(const G4Mat     96 G4double G4hICRU49p::StoppingPower(const G4Material* material,
116                                          G4dou     97                                          G4double kineticEnergy) 
117 {                                                  98 {
118   G4double ionloss = 0.0 ;                         99   G4double ionloss = 0.0 ;
119                                                   100 
120   // pure material (normally not the case for     101   // pure material (normally not the case for this function)
121   if(1 == (material->GetNumberOfElements())) {    102   if(1 == (material->GetNumberOfElements())) {
122     G4double z = material->GetZ() ;               103     G4double z = material->GetZ() ;
123     ionloss = ElectronicStoppingPower( z, kine << 104     G4double ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;  
124                                                   105 
125   } else if (iMolecula < 11) {                    106   } else if (iMolecula < 11) {
126                                                   107   
127     // The data and the fit from:                 108     // The data and the fit from: 
128     // ICRU Report N49, 1993. Ziegler's model     109     // ICRU Report N49, 1993. Ziegler's model for protons.
129     // Proton kinetic energy for parametrisati    110     // Proton kinetic energy for parametrisation (keV/amu)  
130                                                   111 
131     G4double T = kineticEnergy/(keV*protonMass    112     G4double T = kineticEnergy/(keV*protonMassAMU) ; 
132                                                   113 
133      static const G4double a[11][5] = {        << 114      static G4double a[11][5] = {
134    {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2. << 115     1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2, 
135    {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7. << 116     7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3, 
136    {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8. << 117     7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3, 
137    {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9. << 118     8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3, 
138    {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3. << 119     1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2, 
139    {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1. << 120     3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1, 
140    {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3. << 121     1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2, 
141    {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1. << 122     8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2, 
142    {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7. << 123     4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3, 
143    {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8. << 124     4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3, 
144    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1. << 125     2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2 };
145                                                   126       
146                                                   127 
147     if ( T < 10.0 ) {                             128     if ( T < 10.0 ) {
148       ionloss = a[iMolecula][0] * std::sqrt(T) << 129       ionloss = a[iMolecula][0] * sqrt(T) ;
149                                                   130     
150     } else if ( T < 10000.0 ) {                   131     } else if ( T < 10000.0 ) {
151       G4double slow  = a[iMolecula][1] * std:: << 132       G4double slow  = a[iMolecula][1] * pow(T, 0.45) ;
152       G4double shigh = std::log( 1.0 + a[iMole << 133       G4double shigh = log( 1.0 + a[iMolecula][3]/T  
153                      + a[iMolecula][4]*T ) * a    134                      + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
154       ionloss = slow*shigh / (slow + shigh) ;     135       ionloss = slow*shigh / (slow + shigh) ;     
155     }                                             136     } 
156                                                   137 
157     if ( ionloss < 0.0) ionloss = 0.0 ;           138     if ( ionloss < 0.0) ionloss = 0.0 ;
158      /////////////////////////////////////////    139      /////////////////////////////////////////////////////////////////
159     // Graphite may be implemented in a very a    140     // Graphite may be implemented in a very approximate way (scaling 
160     // amorphous results according to rough fi    141     // amorphous results according to rough fits to ICRU tables of results:
161     // 1-100 keV: *(1+0.023+0.0066*std::log10( << 142     // 1-100 keV: *(1+0.023+0.0066*log10(E))
162     // 100-700 keV: *(1+0.089-0.0248*std::log1 << 143     // 100-700 keV: *(1+0.089-0.0248*log10(E-99.))
163     // 700-10000 keV: *(1+0.089-0.0248*std::lo << 144     // 700-10000 keV: *(1+0.089-0.0248*log10(700.-99.))
164     // continuity is (should!) be garanteed, b    145     // continuity is (should!) be garanteed, but not continuity of the
165     // first derivative. A better fit is in or    146     // first derivative. A better fit is in order!       
166     if ( 10 == iMolecula ) {                      147     if ( 10 == iMolecula ) { 
167       if (T < 100.0) {                            148       if (T < 100.0) {    
168   ionloss *= (1.0+0.023+0.0066*std::log10(T)); << 149   ionloss *= (1.0+0.023+0.0066*log10(T));  
169       }                                           150       }
170       else if (T < 700.0) {                       151       else if (T < 700.0) {   
171   ionloss *=(1.0+0.089-0.0248*std::log10(T-99. << 152   ionloss *=(1.0+0.089-0.0248*log10(T-99.));
172       }                                           153       } 
173       else if (T < 10000.0) {                     154       else if (T < 10000.0) {    
174   ionloss *=(1.0+0.089-0.0248*std::log10(700.- << 155   ionloss *=(1.0+0.089-0.0248*log10(700.-99.));
175       }                                           156       }
176     }                                             157     }
177   }                                               158   }
178                                                   159   
179   return ionloss;                                 160   return ionloss;
180 }                                                 161 }  
181                                                   162 
182 //....oooOO0OOooo........oooOO0OOooo........oo    163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183                                                   164 
184 G4double G4hICRU49p::ElectronicStoppingPower(G    165 G4double G4hICRU49p::ElectronicStoppingPower(G4double z,
185                                              G    166                                              G4double kineticEnergy) const
186 {                                                 167 {
187   G4double ionloss ;                              168   G4double ionloss ;
188   G4int i = G4int(z)-1 ;  // index of atom        169   G4int i = G4int(z)-1 ;  // index of atom
189   if(i < 0)  i = 0 ;                              170   if(i < 0)  i = 0 ;
190   if(i > 91) i = 91 ;                             171   if(i > 91) i = 91 ;
191                                                   172   
192   // The data and the fit from:                   173   // The data and the fit from: 
193   // ICRU Report 49, 1993. Ziegler's type of p    174   // ICRU Report 49, 1993. Ziegler's type of parametrisations.
194   // Proton kinetic energy for parametrisation    175   // Proton kinetic energy for parametrisation (keV/amu)  
195                                                   176 
196   G4double T = kineticEnergy/(keV*protonMassAM    177   G4double T = kineticEnergy/(keV*protonMassAMU) ; 
197                                                   178   
198   static const G4double a[92][5] = {           << 179   static G4double a[92][5] = {
199    {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1. << 180     1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1,
200    {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5. << 181     1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2,
201    {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4. << 182     1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2,
202    {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3. << 183     2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2,
203    {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2. << 184     2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2,
204    {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1. << 185     2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2,
205    {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2. << 186     2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2,
206    {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2. << 187     2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2,
207    {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1. << 188     2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2,
208    {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1. << 189     1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2,
209    {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1. << 190     2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2,
210    {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1. << 191     3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2,
211    {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2. << 192     4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2,
212    {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1. << 193     4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2,
213    {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1. << 194     3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2,
214    {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1. << 195     3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2,
215    {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1. << 196     5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2,
216    {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1. << 197     5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2,
217    {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1. << 198     5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2,
218    {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1. << 199     5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2,
219    {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9. << 200     5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3,
220    {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8. << 201     4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3,
221    {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9. << 202     4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3,
222    {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8. << 203     3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3,
223    {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8. << 204     3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3,
224    {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7. << 205     3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3,
225    {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7. << 206     3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3,
226    {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8. << 207     3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3,
227    {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2. << 208     3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2,
228    {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6. << 209     4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3,
229    {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6. << 210     5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3,
230    {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9. << 211     5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3,
231    {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6. << 212     5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3,
232    {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7. << 213     5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3,
233    {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6. << 214     6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3,
234    {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2. << 215     6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3,
235    {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6. << 216     5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3,
236    {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6. << 217     6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3,
237    {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5. << 218     6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3,
238    {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5. << 219     6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3,
239    {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5. << 220     6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3,
240    {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5. << 221     6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3,
241    {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5. << 222     6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3,
242    {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5. << 223     6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3,
243    {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4. << 224     5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3,
244    {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4. << 225     5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3,
245    {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1. << 226     5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2,
246    {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4. << 227     5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3,
247    {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4. << 228     6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3,
248    {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4. << 229     6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3,
249    {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9. << 230     7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3,
250    {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4. << 231     6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3,
251    {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4. << 232     7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3,
252    {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6. << 233     8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3,
253    {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4. << 234     7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3,
254    {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4. << 235     7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3,
255    {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4. << 236     8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3,
256    {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4. << 237     7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3,
257    {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4. << 238     7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3,
258    {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4. << 239     7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3,
259    {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4. << 240     6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3,
260    {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3. << 241     6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3,
261    {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3. << 242     6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3,
262    {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3. << 243     6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3,
263    {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3. << 244     6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3,
264    {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3. << 245     5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3,
265    {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3. << 246     5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3,
266    {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3. << 247     5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3,
267    {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3. << 248     4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3,
268    {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3. << 249     4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3,
269    {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3. << 250     4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3,
270    {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3. << 251     5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3,
271    {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3. << 252     4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3,
272    {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3. << 253     4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3,
273    {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3. << 254     5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3,
274    {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3. << 255     5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3,
275    {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2. << 256     4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3,
276    {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2. << 257     4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3,
277    {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2. << 258     4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2,
278    {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2. << 259     4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3,
279    {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2. << 260     4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3,
280    {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2. << 261     5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3,
281    {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2. << 262     5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3,
282    {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2. << 263     6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3,
283    {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2. << 264     6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3,
284    {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2. << 265     6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3,
285    {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2. << 266     6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3,
286    {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2. << 267     7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3,
287    {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2. << 268     7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3,
288    {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2. << 269     7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3,
289    {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2. << 270     7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3,
290    {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2. << 271     7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3
291   };                                              272   };
292                                                   273 
293   G4double fac = 1.0 ;                            274   G4double fac = 1.0 ;
294                                                   275 
295     // Carbon specific case for E < 40 keV        276     // Carbon specific case for E < 40 keV
296   if ( T < 40.0 && 5 == i) {                      277   if ( T < 40.0 && 5 == i) {
297     fac = std::sqrt(T/40.0) ;                  << 278     fac = sqrt(T/40.0) ;
298     T = 40.0 ;                                    279     T = 40.0 ;  
299                                                   280 
300     // Free electron gas model                    281     // Free electron gas model
301   } else if ( T < 10.0 ) {                        282   } else if ( T < 10.0 ) { 
302     fac = std::sqrt(T*0.1) ;                   << 283     fac = sqrt(T*0.1) ;
303     T =10.0 ;                                     284     T =10.0 ;  
304   }                                               285   }
305                                                   286 
306   // Main parametrisation                         287   // Main parametrisation
307   G4double slow  = a[i][1] * std::pow(T, 0.45) << 288   G4double slow  = a[i][1] * pow(T, 0.45) ;
308   G4double shigh = std::log( 1.0 + a[i][3]/T + << 289   G4double shigh = log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
309   ionloss = slow*shigh*fac / (slow + shigh) ;     290   ionloss = slow*shigh*fac / (slow + shigh) ;     
310                                                   291   
311   if ( ionloss < 0.0) ionloss = 0.0 ;             292   if ( ionloss < 0.0) ionloss = 0.0 ;
312                                                   293   
313   return ionloss;                                 294   return ionloss;
314 }                                                 295 }
315                                                   296 
316                                                   297 
317                                                   298 
318                                                   299