00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 #include "G4hICRU49p.hh"
00056
00057 #include "globals.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4UnitsTable.hh"
00060 #include "G4Material.hh"
00061
00062
00063
00064 G4hICRU49p::G4hICRU49p():G4VhElectronicStoppingPower(),
00065 iMolecula(0),
00066 protonMassAMU(1.007276)
00067 {;}
00068
00069
00070
00071 G4hICRU49p::~G4hICRU49p()
00072 {;}
00073
00074
00075
00076 G4bool G4hICRU49p::HasMaterial(const G4Material* material)
00077 {
00078 G4String chFormula = material->GetChemicalFormula() ;
00079 G4String myFormula = G4String(" ") ;
00080
00081 if (myFormula == chFormula ) {
00082 if(1 == (material->GetNumberOfElements())) return true;
00083 return false ;
00084 }
00085
00086
00087 const size_t numberOfMolecula = 11 ;
00088 static G4String name[numberOfMolecula] = {
00089 "Al_2O_3", "CO_2", "CH_4",
00090 "(C_2H_4)_N-Polyethylene", "(C_2H_4)_N-Polypropylene", "(C_8H_8)_N",
00091 "C_3H_8", "SiO_2", "H_2O",
00092 "H_2O-Gas", "Graphite" } ;
00093
00094
00095 const G4State theState = material->GetState() ;
00096
00097 myFormula = G4String("H_2O");
00098 if( theState == kStateGas && myFormula == chFormula) {
00099 chFormula = G4String("H_2O-Gas");
00100 }
00101
00102
00103 for (size_t i=0; i<numberOfMolecula; i++) {
00104 if (chFormula == name[i]) {
00105 SetMoleculaNumber(i) ;
00106 return true ;
00107 }
00108 }
00109 return false ;
00110 }
00111
00112
00113
00114 G4double G4hICRU49p::StoppingPower(const G4Material* material,
00115 G4double kineticEnergy)
00116 {
00117 G4double ionloss = 0.0 ;
00118
00119
00120 if(1 == (material->GetNumberOfElements())) {
00121 G4double z = material->GetZ() ;
00122 ionloss = ElectronicStoppingPower( z, kineticEnergy ) ;
00123
00124 } else if (iMolecula < 11) {
00125
00126
00127
00128
00129
00130 G4double T = kineticEnergy/(keV*protonMassAMU) ;
00131
00132 static G4double a[11][5] = {
00133 {1.187E+1, 1.343E+1, 1.069E+4, 7.723E+2, 2.153E-2},
00134 {7.802E+0, 8.814E+0, 8.303E+3, 7.446E+2, 7.966E-3},
00135 {7.294E+0, 8.284E+0, 5.010E+3, 4.544E+2, 8.153E-3},
00136 {8.646E+0, 9.800E+0, 7.066E+3, 4.581E+2, 9.383E-3},
00137 {1.286E+1, 1.462E+1, 5.625E+3, 2.621E+3, 3.512E-2},
00138 {3.229E+1, 3.696E+1, 8.918E+3, 3.244E+3, 1.273E-1},
00139 {1.604E+1, 1.825E+1, 6.967E+3, 2.307E+3, 3.775E-2},
00140 {8.049E+0, 9.099E+0, 9.257E+3, 3.846E+2, 1.007E-2},
00141 {4.015E+0, 4.542E+0, 3.955E+3, 4.847E+2, 7.904E-3},
00142 {4.571E+0, 5.173E+0, 4.346E+3, 4.779E+2, 8.572E-3},
00143 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2} };
00144
00145
00146 if ( T < 10.0 ) {
00147 ionloss = a[iMolecula][0] * std::sqrt(T) ;
00148
00149 } else if ( T < 10000.0 ) {
00150 G4double slow = a[iMolecula][1] * std::pow(T, 0.45) ;
00151 G4double shigh = std::log( 1.0 + a[iMolecula][3]/T
00152 + a[iMolecula][4]*T ) * a[iMolecula][2]/T ;
00153 ionloss = slow*shigh / (slow + shigh) ;
00154 }
00155
00156 if ( ionloss < 0.0) ionloss = 0.0 ;
00158
00159
00160
00161
00162
00163
00164
00165 if ( 10 == iMolecula ) {
00166 if (T < 100.0) {
00167 ionloss *= (1.0+0.023+0.0066*std::log10(T));
00168 }
00169 else if (T < 700.0) {
00170 ionloss *=(1.0+0.089-0.0248*std::log10(T-99.));
00171 }
00172 else if (T < 10000.0) {
00173 ionloss *=(1.0+0.089-0.0248*std::log10(700.-99.));
00174 }
00175 }
00176 }
00177
00178 return ionloss;
00179 }
00180
00181
00182
00183 G4double G4hICRU49p::ElectronicStoppingPower(G4double z,
00184 G4double kineticEnergy) const
00185 {
00186 G4double ionloss ;
00187 G4int i = G4int(z)-1 ;
00188 if(i < 0) i = 0 ;
00189 if(i > 91) i = 91 ;
00190
00191
00192
00193
00194
00195 G4double T = kineticEnergy/(keV*protonMassAMU) ;
00196
00197 static G4double a[92][5] = {
00198 {1.254E+0, 1.440E+0, 2.426E+2, 1.200E+4, 1.159E-1},
00199 {1.229E+0, 1.397E+0, 4.845E+2, 5.873E+3, 5.225E-2},
00200 {1.411E+0, 1.600E+0, 7.256E+2, 3.013E+3, 4.578E-2},
00201 {2.248E+0, 2.590E+0, 9.660E+2, 1.538E+2, 3.475E-2},
00202 {2.474E+0, 2.815E+0, 1.206E+3, 1.060E+3, 2.855E-2},
00203 {2.631E+0, 2.601E+0, 1.701E+3, 1.279E+3, 1.638E-2},
00204 {2.954E+0, 3.350E+0, 1.683E+3, 1.900E+3, 2.513E-2},
00205 {2.652E+0, 3.000E+0, 1.920E+3, 2.000E+3, 2.230E-2},
00206 {2.085E+0, 2.352E+0, 2.157E+3, 2.634E+3, 1.816E-2},
00207 {1.951E+0, 2.199E+0, 2.393E+3, 2.699E+3, 1.568E-2},
00208 {2.542E+0, 2.869E+0, 2.628E+3, 1.854E+3, 1.472E-2},
00209 {3.791E+0, 4.293E+0, 2.862E+3, 1.009E+3, 1.397E-2},
00210 {4.154E+0, 4.739E+0, 2.766E+3, 1.645E+2, 2.023E-2},
00211 {4.914E+0, 5.598E+0, 3.193E+3, 2.327E+2, 1.419E-2},
00212 {3.232E+0, 3.647E+0, 3.561E+3, 1.560E+3, 1.267E-2},
00213 {3.447E+0, 3.891E+0, 3.792E+3, 1.219E+3, 1.211E-2},
00214 {5.301E+0, 6.008E+0, 3.969E+3, 6.451E+2, 1.183E-2},
00215 {5.731E+0, 6.500E+0, 4.253E+3, 5.300E+2, 1.123E-2},
00216 {5.152E+0, 5.833E+0, 4.482E+3, 5.457E+2, 1.129E-2},
00217 {5.521E+0, 6.252E+0, 4.710E+3, 5.533E+2, 1.112E-2},
00218 {5.201E+0, 5.884E+0, 4.938E+3, 5.609E+2, 9.995E-3},
00219 {4.858E+0, 5.489E+0, 5.260E+3, 6.511E+2, 8.930E-3},
00220 {4.479E+0, 5.055E+0, 5.391E+3, 9.523E+2, 9.117E-3},
00221 {3.983E+0, 4.489E+0, 5.616E+3, 1.336E+3, 8.413E-3},
00222 {3.469E+0, 3.907E+0, 5.725E+3, 1.461E+3, 8.829E-3},
00223 {3.519E+0, 3.963E+0, 6.065E+3, 1.243E+3, 7.782E-3},
00224 {3.140E+0, 3.535E+0, 6.288E+3, 1.372E+3, 7.361E-3},
00225 {3.553E+0, 4.004E+0, 6.205E+3, 5.551E+2, 8.763E-3},
00226 {3.696E+0, 4.194E+0, 4.649E+3, 8.113E+1, 2.242E-2},
00227 {4.210E+0, 4.750E+0, 6.953E+3, 2.952E+2, 6.809E-3},
00228 {5.041E+0, 5.697E+0, 7.173E+3, 2.026E+2, 6.725E-3},
00229 {5.554E+0, 6.300E+0, 6.496E+3, 1.100E+2, 9.689E-3},
00230 {5.323E+0, 6.012E+0, 7.611E+3, 2.925E+2, 6.447E-3},
00231 {5.874E+0, 6.656E+0, 7.395E+3, 1.175E+2, 7.684E-3},
00232 {6.658E+0, 7.536E+0, 7.694E+3, 2.223E+2, 6.509E-3},
00233 {6.413E+0, 7.240E+0, 1.185E+4, 1.537E+2, 2.880E-3},
00234 {5.694E+0, 6.429E+0, 8.478E+3, 2.929E+2, 6.087E-3},
00235 {6.339E+0, 7.159E+0, 8.693E+3, 3.303E+2, 6.003E-3},
00236 {6.407E+0, 7.234E+0, 8.907E+3, 3.678E+2, 5.889E-3},
00237 {6.734E+0, 7.603E+0, 9.120E+3, 4.052E+2, 5.765E-3},
00238 {6.901E+0, 7.791E+0, 9.333E+3, 4.427E+2, 5.587E-3},
00239 {6.424E+0, 7.248E+0, 9.545E+3, 4.802E+2, 5.376E-3},
00240 {6.799E+0, 7.671E+0, 9.756E+3, 5.176E+2, 5.315E-3},
00241 {6.109E+0, 6.887E+0, 9.966E+3, 5.551E+2, 5.151E-3},
00242 {5.924E+0, 6.677E+0, 1.018E+4, 5.925E+2, 4.919E-3},
00243 {5.238E+0, 5.900E+0, 1.038E+4, 6.300E+2, 4.758E-3},
00244 {5.345E+0, 6.038E+0, 6.790E+3, 3.978E+2, 1.676E-2},
00245 {5.814E+0, 6.554E+0, 1.080E+4, 3.555E+2, 4.626E-3},
00246 {6.229E+0, 7.024E+0, 1.101E+4, 3.709E+2, 4.540E-3},
00247 {6.409E+0, 7.227E+0, 1.121E+4, 3.864E+2, 4.474E-3},
00248 {7.500E+0, 8.480E+0, 8.608E+3, 3.480E+2, 9.074E-3},
00249 {6.979E+0, 7.871E+0, 1.162E+4, 3.924E+2, 4.402E-3},
00250 {7.725E+0, 8.716E+0, 1.183E+4, 3.948E+2, 4.376E-3},
00251 {8.337E+0, 9.425E+0, 1.051E+4, 2.696E+2, 6.206E-3},
00252 {7.287E+0, 8.218E+0, 1.223E+4, 3.997E+2, 4.447E-3},
00253 {7.899E+0, 8.911E+0, 1.243E+4, 4.021E+2, 4.511E-3},
00254 {8.041E+0, 9.071E+0, 1.263E+4, 4.045E+2, 4.540E-3},
00255 {7.488E+0, 8.444E+0, 1.283E+4, 4.069E+2, 4.420E-3},
00256 {7.291E+0, 8.219E+0, 1.303E+4, 4.093E+2, 4.298E-3},
00257 {7.098E+0, 8.000E+0, 1.323E+4, 4.118E+2, 4.182E-3},
00258 {6.909E+0, 7.786E+0, 1.343E+4, 4.142E+2, 4.058E-3},
00259 {6.728E+0, 7.580E+0, 1.362E+4, 4.166E+2, 3.976E-3},
00260 {6.551E+0, 7.380E+0, 1.382E+4, 4.190E+2, 3.877E-3},
00261 {6.739E+0, 7.592E+0, 1.402E+4, 4.214E+2, 3.863E-3},
00262 {6.212E+0, 6.996E+0, 1.421E+4, 4.239E+2, 3.725E-3},
00263 {5.517E+0, 6.210E+0, 1.440E+4, 4.263E+2, 3.632E-3},
00264 {5.220E+0, 5.874E+0, 1.460E+4, 4.287E+2, 3.498E-3},
00265 {5.071E+0, 5.706E+0, 1.479E+4, 4.330E+2, 3.405E-3},
00266 {4.926E+0, 5.542E+0, 1.498E+4, 4.335E+2, 3.342E-3},
00267 {4.788E+0, 5.386E+0, 1.517E+4, 4.359E+2, 3.292E-3},
00268 {4.893E+0, 5.505E+0, 1.536E+4, 4.384E+2, 3.243E-3},
00269 {5.028E+0, 5.657E+0, 1.555E+4, 4.408E+2, 3.195E-3},
00270 {4.738E+0, 5.329E+0, 1.574E+4, 4.432E+2, 3.186E-3},
00271 {4.587E+0, 5.160E+0, 1.541E+4, 4.153E+2, 3.406E-3},
00272 {5.201E+0, 5.851E+0, 1.612E+4, 4.416E+2, 3.122E-3},
00273 {5.071E+0, 5.704E+0, 1.630E+4, 4.409E+2, 3.082E-3},
00274 {4.946E+0, 5.563E+0, 1.649E+4, 4.401E+2, 2.965E-3},
00275 {4.477E+0, 5.034E+0, 1.667E+4, 4.393E+2, 2.871E-3},
00276 {4.844E+0, 5.458E+0, 7.852E+3, 9.758E+2, 2.077E-2},
00277 {4.307E+0, 4.843E+0, 1.704E+4, 4.878E+2, 2.882E-3},
00278 {4.723E+0, 5.311E+0, 1.722E+4, 5.370E+2, 2.913E-3},
00279 {5.319E+0, 5.982E+0, 1.740E+4, 5.863E+2, 2.871E-3},
00280 {5.956E+0, 6.700E+0, 1.780E+4, 6.770E+2, 2.660E-3},
00281 {6.158E+0, 6.928E+0, 1.777E+4, 5.863E+2, 2.812E-3},
00282 {6.203E+0, 6.979E+0, 1.795E+4, 5.863E+2, 2.776E-3},
00283 {6.181E+0, 6.954E+0, 1.812E+4, 5.863E+2, 2.748E-3},
00284 {6.949E+0, 7.820E+0, 1.830E+4, 5.863E+2, 2.737E-3},
00285 {7.506E+0, 8.448E+0, 1.848E+4, 5.863E+2, 2.727E-3},
00286 {7.648E+0, 8.609E+0, 1.866E+4, 5.863E+2, 2.697E-3},
00287 {7.711E+0, 8.679E+0, 1.883E+4, 5.863E+2, 2.641E-3},
00288 {7.407E+0, 8.336E+0, 1.901E+4, 5.863E+2, 2.603E-3},
00289 {7.290E+0, 8.204E+0, 1.918E+4, 5.863E+2, 2.673E-3}
00290 };
00291
00292 G4double fac = 1.0 ;
00293
00294
00295 if ( T < 40.0 && 5 == i) {
00296 fac = std::sqrt(T/40.0) ;
00297 T = 40.0 ;
00298
00299
00300 } else if ( T < 10.0 ) {
00301 fac = std::sqrt(T*0.1) ;
00302 T =10.0 ;
00303 }
00304
00305
00306 G4double slow = a[i][1] * std::pow(T, 0.45) ;
00307 G4double shigh = std::log( 1.0 + a[i][3]/T + a[i][4]*T ) * a[i][2]/T ;
00308 ionloss = slow*shigh*fac / (slow + shigh) ;
00309
00310 if ( ionloss < 0.0) ionloss = 0.0 ;
00311
00312 return ionloss;
00313 }
00314
00315
00316