Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4hIonEffChargeSquare 32 // File name: G4hIonEffChargeSquare 33 // 33 // 34 // Author: V.Ivanchenko (Vladimir.Ivanc 34 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch) 35 // 35 // 36 // Creation date: 20 July 2000 36 // Creation date: 20 July 2000 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 20/07/2000 V.Ivanchenko First implementati 39 // 20/07/2000 V.Ivanchenko First implementation 40 // 18/06/2001 V.Ivanchenko Continuation for e 40 // 18/06/2001 V.Ivanchenko Continuation for eff.charge (small change of y) 41 // 08/10/2002 V.Ivanchenko The charge of the 41 // 08/10/2002 V.Ivanchenko The charge of the nucleus is used not charge of 42 // DynamicParticle 42 // DynamicParticle 43 // 43 // 44 // Class Description: 44 // Class Description: 45 // 45 // 46 // Ion effective charge model 46 // Ion effective charge model 47 // J.F.Ziegler and J.M.Manoyan, The stopping o 47 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988 48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 49 // 49 // 50 // Class Description: End 50 // Class Description: End 51 // 51 // 52 // ------------------------------------------- 52 // ------------------------------------------------------------------- 53 // 53 // 54 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 55 56 #include "G4hIonEffChargeSquare.hh" 56 #include "G4hIonEffChargeSquare.hh" 57 #include "G4PhysicalConstants.hh" 57 #include "G4PhysicalConstants.hh" 58 #include "G4SystemOfUnits.hh" 58 #include "G4SystemOfUnits.hh" 59 #include "G4DynamicParticle.hh" 59 #include "G4DynamicParticle.hh" 60 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleDefinition.hh" 61 #include "G4Material.hh" 61 #include "G4Material.hh" 62 #include "G4Element.hh" 62 #include "G4Element.hh" 63 #include "G4Exp.hh" 63 #include "G4Exp.hh" 64 64 65 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 66 66 67 G4hIonEffChargeSquare::G4hIonEffChargeSquare(c 67 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name) 68 : G4VLowEnergyModel(name), 68 : G4VLowEnergyModel(name), 69 theHeMassAMU(4.0026) 69 theHeMassAMU(4.0026) 70 {;} 70 {;} 71 71 72 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 73 73 74 G4hIonEffChargeSquare::~G4hIonEffChargeSquare( 74 G4hIonEffChargeSquare::~G4hIonEffChargeSquare() 75 {;} 75 {;} 76 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 78 79 G4double G4hIonEffChargeSquare::TheValue(const 79 G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle, 80 const G4Mat 80 const G4Material* material) 81 { 81 { 82 G4double energy = particle->GetKineticEnergy 82 G4double energy = particle->GetKineticEnergy() ; 83 G4double particleMass = particle->GetMass() 83 G4double particleMass = particle->GetMass() ; 84 G4double charge = (particle->GetDefinition() 84 G4double charge = (particle->GetDefinition()->GetPDGCharge())/eplus ; 85 85 86 G4double q = IonEffChargeSquare(material,ene 86 G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ; 87 87 88 return q ; 88 return q ; 89 } 89 } 90 90 91 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 92 93 G4double G4hIonEffChargeSquare::TheValue(const 93 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle, 94 const G4Mat 94 const G4Material* material, 95 G4double kineticEnergy) 95 G4double kineticEnergy) 96 { 96 { 97 // SetRateMass(aParticle) ; 97 // SetRateMass(aParticle) ; 98 G4double particleMass = aParticle->GetPDGMas 98 G4double particleMass = aParticle->GetPDGMass() ; 99 G4double charge = (aParticle->GetPDGCharge() 99 G4double charge = (aParticle->GetPDGCharge())/eplus ; 100 100 101 G4double q = IonEffChargeSquare(material,kin 101 G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ; 102 102 103 return q ; 103 return q ; 104 } 104 } 105 105 106 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 107 108 G4double G4hIonEffChargeSquare::HighEnergyLimi 108 G4double G4hIonEffChargeSquare::HighEnergyLimit( 109 const G4ParticleDefinition* , 109 const G4ParticleDefinition* , 110 const G4Material* ) const 110 const G4Material* ) const 111 { 111 { 112 return 1.0*TeV ; 112 return 1.0*TeV ; 113 } 113 } 114 114 115 //....oooOO0OOooo........oooOO0OOooo........oo 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 116 117 G4double G4hIonEffChargeSquare::LowEnergyLimit 117 G4double G4hIonEffChargeSquare::LowEnergyLimit( 118 const G4ParticleDefinition* , 118 const G4ParticleDefinition* , 119 const G4Material* ) const 119 const G4Material* ) const 120 { 120 { 121 return 0.0 ; 121 return 0.0 ; 122 } 122 } 123 123 124 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 125 125 126 G4double G4hIonEffChargeSquare::HighEnergyLimi 126 G4double G4hIonEffChargeSquare::HighEnergyLimit( 127 const G4ParticleDefinition* ) cons 127 const G4ParticleDefinition* ) const 128 { 128 { 129 return 1.0*TeV ; 129 return 1.0*TeV ; 130 } 130 } 131 131 132 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 133 133 134 G4double G4hIonEffChargeSquare::LowEnergyLimit 134 G4double G4hIonEffChargeSquare::LowEnergyLimit( 135 const G4ParticleDefinition* ) 135 const G4ParticleDefinition* ) const 136 { 136 { 137 return 0.0 ; 137 return 0.0 ; 138 } 138 } 139 139 140 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 141 141 142 G4bool G4hIonEffChargeSquare::IsInCharge(const 142 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* , 143 const G4Material* 143 const G4Material* ) const 144 { 144 { 145 return true ; 145 return true ; 146 } 146 } 147 147 148 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 149 150 G4bool G4hIonEffChargeSquare::IsInCharge(const 150 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* , 151 const G4Mat 151 const G4Material* ) const 152 { 152 { 153 return true ; 153 return true ; 154 } 154 } 155 155 156 //....oooOO0OOooo........oooOO0OOooo........oo 156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 157 158 G4double G4hIonEffChargeSquare::IonEffChargeSq 158 G4double G4hIonEffChargeSquare::IonEffChargeSquare( 159 const G4Material* material, 159 const G4Material* material, 160 G4double kineticEnergy, 160 G4double kineticEnergy, 161 G4double particleMass, 161 G4double particleMass, 162 G4double ionCharge) const 162 G4double ionCharge) const 163 { 163 { 164 // The aproximation of ion effective charge 164 // The aproximation of ion effective charge from: 165 // J.F.Ziegler, J.P. Biersack, U. Littmark 165 // J.F.Ziegler, J.P. Biersack, U. Littmark 166 // The Stopping and Range of Ions in Matter, 166 // The Stopping and Range of Ions in Matter, 167 // Vol.1, Pergamon Press, 1985 167 // Vol.1, Pergamon Press, 1985 168 168 169 // Fast ions or hadrons 169 // Fast ions or hadrons 170 G4double reducedEnergy = kineticEnergy * pro 170 G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ; 171 if(reducedEnergy < 1.0*keV) reducedEnergy = 171 if(reducedEnergy < 1.0*keV) reducedEnergy = 1.0*keV; 172 if( (reducedEnergy > ionCharge * 10.0 * MeV) 172 if( (reducedEnergy > ionCharge * 10.0 * MeV) || 173 (ionCharge < 1.5) ) return ionCharge*ion 173 (ionCharge < 1.5) ) return ionCharge*ionCharge ; 174 174 175 static const G4double vFermi[92] = { 175 static const G4double vFermi[92] = { 176 1.0309, 0.15976, 0.59782, 1.0781, 1.0486 176 1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424, 177 0.45259, 0.71074, 0.90519, 0.97411, 0.9718 177 0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712, 178 0.81707, 0.9943, 1.1423, 1.2381, 1.1222 178 0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411, 179 0.84912, 0.95, 1.0903, 1.0429, 0.4971 179 0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207, 180 1.029, 1.2542, 1.122, 1.1241, 1.0882 180 1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054, 181 0.93155, 1.0047, 0.55379, 0.43289, 0.3263 181 0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413, 182 0.71418, 0.71453, 0.5911, 0.70263, 0.6804 182 0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884, 183 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 183 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654, 184 0.66401, 0.84912, 0.88433, 0.80746, 0.4335 184 0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065, 185 1.9578, 1.0257} ; 185 1.9578, 1.0257} ; 186 186 187 static const G4double c[6] = {0.2865, 0.126 187 static const G4double c[6] = {0.2865, 0.1266, -0.001429, 188 0.02402,-0.01135, 0. 188 0.02402,-0.01135, 0.001475} ; 189 189 190 // get elements in the actual material, 190 // get elements in the actual material, 191 const G4ElementVector* theElementVector = ma 191 const G4ElementVector* theElementVector = material->GetElementVector() ; 192 const G4double* theAtomicNumDensityVector = 192 const G4double* theAtomicNumDensityVector = 193 material->GetAtomicNu 193 material->GetAtomicNumDensityVector() ; 194 const G4int NumberOfElements = (G4int)materi << 194 const G4int NumberOfElements = material->GetNumberOfElements() ; 195 195 196 // loop for the elements in the material 196 // loop for the elements in the material 197 // to find out average values Z, vF, lF 197 // to find out average values Z, vF, lF 198 G4double z = 0.0, vF = 0.0, norm = 0.0 ; 198 G4double z = 0.0, vF = 0.0, norm = 0.0 ; 199 199 200 if( 1 == NumberOfElements ) { 200 if( 1 == NumberOfElements ) { 201 z = material->GetZ() ; 201 z = material->GetZ() ; 202 G4int iz = G4int(z) - 1 ; 202 G4int iz = G4int(z) - 1 ; 203 if(iz < 0) iz = 0 ; 203 if(iz < 0) iz = 0 ; 204 else if(iz > 91) iz = 91 ; 204 else if(iz > 91) iz = 91 ; 205 vF = vFermi[iz] ; 205 vF = vFermi[iz] ; 206 206 207 } else { 207 } else { 208 for (G4int iel=0; iel<NumberOfElements; ie 208 for (G4int iel=0; iel<NumberOfElements; iel++) 209 { 209 { 210 const G4Element* element = (*theElemen 210 const G4Element* element = (*theElementVector)[iel] ; 211 G4double z2 = element->GetZ() ; 211 G4double z2 = element->GetZ() ; 212 const G4double weight = theAtomicNumDe 212 const G4double weight = theAtomicNumDensityVector[iel] ; 213 norm += weight ; 213 norm += weight ; 214 z += z2 * weight ; 214 z += z2 * weight ; 215 G4int iz = G4int(z2) - 1 ; 215 G4int iz = G4int(z2) - 1 ; 216 if(iz < 0) iz = 0 ; 216 if(iz < 0) iz = 0 ; 217 else if(iz > 91) iz =91 ; 217 else if(iz > 91) iz =91 ; 218 vF += vFermi[iz] * weight ; 218 vF += vFermi[iz] * weight ; 219 } 219 } 220 if (norm > 0.0) { << 220 z /= norm ; 221 z /= norm ; << 221 vF /= norm ; 222 vF /= norm ; << 223 } << 224 } 222 } 225 223 226 // Helium ion case 224 // Helium ion case 227 if( ionCharge < 2.5 ) { 225 if( ionCharge < 2.5 ) { 228 226 229 G4double e = std::log(std::max(1.0, kineti 227 G4double e = std::log(std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ; 230 G4double x = c[0] ; 228 G4double x = c[0] ; 231 G4double y = 1.0 ; 229 G4double y = 1.0 ; 232 for (G4int i=1; i<6; i++) { 230 for (G4int i=1; i<6; i++) { 233 y *= e ; 231 y *= e ; 234 x += y * c[i] ; 232 x += y * c[i] ; 235 } 233 } 236 G4double q = 7.6 - e ; 234 G4double q = 7.6 - e ; 237 q = 1.0 + ( 0.007 + 0.00005 * z ) * G4Exp( 235 q = 1.0 + ( 0.007 + 0.00005 * z ) * G4Exp( -q*q ) ; 238 return 4.0 * q * q * (1.0 - G4Exp(-x)) ; 236 return 4.0 * q * q * (1.0 - G4Exp(-x)) ; 239 237 240 // Heavy ion case 238 // Heavy ion case 241 } else { 239 } else { 242 240 243 // v1 is ion velocity in vF unit 241 // v1 is ion velocity in vF unit 244 G4double v1{0.0}, v2{0.0}; << 242 G4double v1 = std::sqrt( reducedEnergy / (25.0 * keV) )/ vF ; 245 if (vF > 0.0) { << 246 v1 = std::sqrt( reducedEnergy / (25.0 * << 247 v2 = 1.0/ (vF*vF); << 248 } << 249 G4double y ; 243 G4double y ; 250 G4double z13 = std::pow(ionCharge, 0.3333) 244 G4double z13 = std::pow(ionCharge, 0.3333) ; 251 245 252 // Faster than Fermi velocity 246 // Faster than Fermi velocity 253 if ( v1 > 1.0 ) { 247 if ( v1 > 1.0 ) { 254 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / 248 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ; 255 249 256 // Slower than Fermi velocity 250 // Slower than Fermi velocity 257 } else { 251 } else { 258 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + 252 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ; 259 } 253 } 260 254 261 G4double y3 = std::pow(y, 0.3) ; 255 G4double y3 = std::pow(y, 0.3) ; 262 G4double q = 1.0 - G4Exp( 0.803*y3 - 1.316 256 G4double q = 1.0 - G4Exp( 0.803*y3 - 1.3167*y3*y3 - 263 0.38157*y - 0.0089 257 0.38157*y - 0.008983*y*y ) ; 264 if( q < 0.0 ) q = 0.0 ; 258 if( q < 0.0 ) q = 0.0 ; 265 259 266 G4double sLocal = 7.6 - std::log(std::max 260 G4double sLocal = 7.6 - std::log(std::max(1.0, reducedEnergy/keV)) ; 267 sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * G4E 261 sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * G4Exp( -sLocal*sLocal )/ (ionCharge*ionCharge) ; 268 262 269 // Screen length according to 263 // Screen length according to 270 // J.F.Ziegler and J.M.Manoyan, The stoppi 264 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 271 // Nucl. Inst. & Meth. in Phys. Res. B35 ( 265 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 272 266 273 G4double lambda = 10.0 * vF * std::pow(1.0 267 G4double lambda = 10.0 * vF * std::pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ; 274 G4double qeff = ionCharge * sLocal * 268 G4double qeff = ionCharge * sLocal * 275 ( q + 0.5*(1.0-q) * std::log(1.0 + lambd << 269 ( q + 0.5*(1.0-q) * std::log(1.0 + lambda*lambda) / (vF*vF) ) ; 276 if( 0.1 > qeff ) qeff = 0.1 ; 270 if( 0.1 > qeff ) qeff = 0.1 ; 277 return qeff*qeff ; 271 return qeff*qeff ; 278 } 272 } 279 } 273 } 280 274