Geant4 Cross Reference |
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 // 26 // 23 // 27 // ------------------------------------------- 24 // ------------------------------------------------------------------- 28 // 25 // 29 // GEANT4 Class file 26 // GEANT4 Class file 30 // 27 // 31 // 28 // 32 // File name: G4hIonEffChargeSquare 29 // File name: G4hIonEffChargeSquare 33 // 30 // 34 // Author: V.Ivanchenko (Vladimir.Ivanc 31 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch) 35 // << 32 // 36 // Creation date: 20 July 2000 33 // Creation date: 20 July 2000 37 // 34 // 38 // Modifications: << 35 // Modifications: 39 // 20/07/2000 V.Ivanchenko First implementati 36 // 20/07/2000 V.Ivanchenko First implementation 40 // 18/06/2001 V.Ivanchenko Continuation for e 37 // 18/06/2001 V.Ivanchenko Continuation for eff.charge (small change of y) 41 // 08/10/2002 V.Ivanchenko The charge of the << 42 // DynamicParticle << 43 // 38 // 44 // Class Description: << 39 // Class Description: 45 // 40 // 46 // Ion effective charge model 41 // Ion effective charge model 47 // J.F.Ziegler and J.M.Manoyan, The stopping o 42 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 48 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988 43 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 49 // 44 // 50 // Class Description: End << 45 // Class Description: End 51 // 46 // 52 // ------------------------------------------- 47 // ------------------------------------------------------------------- 53 // 48 // 54 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 50 56 #include "G4hIonEffChargeSquare.hh" 51 #include "G4hIonEffChargeSquare.hh" 57 #include "G4PhysicalConstants.hh" << 58 #include "G4SystemOfUnits.hh" << 59 #include "G4DynamicParticle.hh" 52 #include "G4DynamicParticle.hh" 60 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh" 61 #include "G4Material.hh" 54 #include "G4Material.hh" 62 #include "G4Element.hh" 55 #include "G4Element.hh" 63 #include "G4Exp.hh" << 64 56 65 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 66 58 67 G4hIonEffChargeSquare::G4hIonEffChargeSquare(c 59 G4hIonEffChargeSquare::G4hIonEffChargeSquare(const G4String& name) 68 : G4VLowEnergyModel(name), << 60 : G4VLowEnergyModel(name), 69 theHeMassAMU(4.0026) 61 theHeMassAMU(4.0026) 70 {;} 62 {;} 71 63 72 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 73 65 74 G4hIonEffChargeSquare::~G4hIonEffChargeSquare( << 66 G4hIonEffChargeSquare::~G4hIonEffChargeSquare() 75 {;} 67 {;} 76 68 77 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 70 79 G4double G4hIonEffChargeSquare::TheValue(const 71 G4double G4hIonEffChargeSquare::TheValue(const G4DynamicParticle* particle, 80 const G4Mat << 72 const G4Material* material) 81 { 73 { 82 G4double energy = particle->GetKineticEnergy 74 G4double energy = particle->GetKineticEnergy() ; 83 G4double particleMass = particle->GetMass() 75 G4double particleMass = particle->GetMass() ; 84 G4double charge = (particle->GetDefinition() << 76 G4double charge = (particle->GetCharge())/eplus ; 85 77 86 G4double q = IonEffChargeSquare(material,ene 78 G4double q = IonEffChargeSquare(material,energy,particleMass,charge) ; 87 79 88 return q ; 80 return q ; 89 } 81 } 90 82 91 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 92 84 93 G4double G4hIonEffChargeSquare::TheValue(const 85 G4double G4hIonEffChargeSquare::TheValue(const G4ParticleDefinition* aParticle, 94 const G4Mat 86 const G4Material* material, 95 G4double kineticEnergy) << 87 G4double kineticEnergy) 96 { 88 { 97 // SetRateMass(aParticle) ; 89 // SetRateMass(aParticle) ; 98 G4double particleMass = aParticle->GetPDGMas 90 G4double particleMass = aParticle->GetPDGMass() ; 99 G4double charge = (aParticle->GetPDGCharge() 91 G4double charge = (aParticle->GetPDGCharge())/eplus ; 100 92 101 G4double q = IonEffChargeSquare(material,kin 93 G4double q = IonEffChargeSquare(material,kineticEnergy,particleMass,charge) ; 102 94 103 return q ; 95 return q ; 104 } 96 } 105 97 106 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 99 108 G4double G4hIonEffChargeSquare::HighEnergyLimi 100 G4double G4hIonEffChargeSquare::HighEnergyLimit( 109 const G4ParticleDefinition* , << 101 const G4ParticleDefinition* aParticle, 110 const G4Material* ) const << 102 const G4Material* material) const 111 { 103 { 112 return 1.0*TeV ; 104 return 1.0*TeV ; 113 } 105 } 114 106 115 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 108 117 G4double G4hIonEffChargeSquare::LowEnergyLimit 109 G4double G4hIonEffChargeSquare::LowEnergyLimit( 118 const G4ParticleDefinition* , << 110 const G4ParticleDefinition* aParticle, 119 const G4Material* ) const << 111 const G4Material* material) const 120 { 112 { 121 return 0.0 ; 113 return 0.0 ; 122 } 114 } 123 115 124 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 125 117 126 G4double G4hIonEffChargeSquare::HighEnergyLimi 118 G4double G4hIonEffChargeSquare::HighEnergyLimit( 127 const G4ParticleDefinition* ) cons << 119 const G4ParticleDefinition* aParticle) const 128 { 120 { 129 return 1.0*TeV ; 121 return 1.0*TeV ; 130 } 122 } 131 123 132 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 133 125 134 G4double G4hIonEffChargeSquare::LowEnergyLimit 126 G4double G4hIonEffChargeSquare::LowEnergyLimit( 135 const G4ParticleDefinition* ) << 127 const G4ParticleDefinition* aParticle) const 136 { 128 { 137 return 0.0 ; 129 return 0.0 ; 138 } 130 } 139 131 140 //....oooOO0OOooo........oooOO0OOooo........oo 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 141 << 133 142 G4bool G4hIonEffChargeSquare::IsInCharge(const << 134 G4bool G4hIonEffChargeSquare::IsInCharge(const G4DynamicParticle* particle, 143 const G4Material* << 135 const G4Material* material) const 144 { 136 { 145 return true ; 137 return true ; 146 } 138 } 147 139 148 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 << 141 150 G4bool G4hIonEffChargeSquare::IsInCharge(const << 142 G4bool G4hIonEffChargeSquare::IsInCharge(const G4ParticleDefinition* aParticle, 151 const G4Mat << 143 const G4Material* material) const 152 { 144 { 153 return true ; 145 return true ; 154 } 146 } 155 147 156 //....oooOO0OOooo........oooOO0OOooo........oo 148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 149 158 G4double G4hIonEffChargeSquare::IonEffChargeSq 150 G4double G4hIonEffChargeSquare::IonEffChargeSquare( 159 const G4Material* material, << 151 const G4Material* material, 160 G4double kineticEnergy, << 152 G4double kineticEnergy, 161 G4double particleMass, << 153 G4double particleMass, 162 G4double ionCharge) const << 154 G4double ionCharge) const 163 { 155 { 164 // The aproximation of ion effective charge << 156 // The aproximation of ion effective charge from: 165 // J.F.Ziegler, J.P. Biersack, U. Littmark 157 // J.F.Ziegler, J.P. Biersack, U. Littmark 166 // The Stopping and Range of Ions in Matter, 158 // The Stopping and Range of Ions in Matter, 167 // Vol.1, Pergamon Press, 1985 159 // Vol.1, Pergamon Press, 1985 168 160 169 // Fast ions or hadrons 161 // Fast ions or hadrons 170 G4double reducedEnergy = kineticEnergy * pro 162 G4double reducedEnergy = kineticEnergy * proton_mass_c2/particleMass ; 171 if(reducedEnergy < 1.0*keV) reducedEnergy = << 163 if( (reducedEnergy > ionCharge * 10.0 * MeV) || 172 if( (reducedEnergy > ionCharge * 10.0 * MeV) << 173 (ionCharge < 1.5) ) return ionCharge*ion 164 (ionCharge < 1.5) ) return ionCharge*ionCharge ; 174 165 175 static const G4double vFermi[92] = { << 166 static G4double vFermi[92] = { 176 1.0309, 0.15976, 0.59782, 1.0781, 1.0486 167 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 168 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 169 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 170 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 171 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 172 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 173 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, 174 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 175 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} ; 176 1.9578, 1.0257} ; 186 177 187 static const G4double c[6] = {0.2865, 0.126 << 178 static G4double lFactor[92] = { >> 179 1.0, 1.0, 1.1, 1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9, >> 180 0.82, 0.81, 0.83, 0.88, 1.0, 0.95, 0.97, 0.99, 0.98, 0.97, >> 181 0.98, 0.97, 0.96, 0.93, 0.91, 0.9, 0.88, 0.9, 0.9, 0.9, >> 182 0.9, 0.85, 0.9, 0.9, 0.91, 0.92, 0.9, 0.9, 0.9, 0.9, >> 183 0.9, 0.88, 0.9, 0.88, 0.88, 0.9, 0.9, 0.88, 0.9, 0.9, >> 184 0.9, 0.9, 0.96, 1.2, 0.9, 0.88, 0.88, 0.85, 0.9, 0.9, >> 185 0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1, 1.08, 1.08, >> 186 1.08, 1.08, 1.09, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, >> 187 1.17, 1.2, 1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16, >> 188 1.16, 1.16} ; >> 189 >> 190 static G4double c[6] = {0.2865, 0.1266, -0.001429, 188 0.02402,-0.01135, 0. 191 0.02402,-0.01135, 0.001475} ; 189 192 190 // get elements in the actual material, 193 // get elements in the actual material, 191 const G4ElementVector* theElementVector = ma 194 const G4ElementVector* theElementVector = material->GetElementVector() ; 192 const G4double* theAtomicNumDensityVector = << 195 const G4double* theAtomicNumDensityVector = 193 material->GetAtomicNu 196 material->GetAtomicNumDensityVector() ; 194 const G4int NumberOfElements = (G4int)materi << 197 const G4int NumberOfElements = material->GetNumberOfElements() ; 195 << 198 196 // loop for the elements in the material 199 // loop for the elements in the material 197 // to find out average values Z, vF, lF 200 // to find out average values Z, vF, lF 198 G4double z = 0.0, vF = 0.0, norm = 0.0 ; << 201 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ; 199 202 200 if( 1 == NumberOfElements ) { 203 if( 1 == NumberOfElements ) { 201 z = material->GetZ() ; 204 z = material->GetZ() ; 202 G4int iz = G4int(z) - 1 ; 205 G4int iz = G4int(z) - 1 ; 203 if(iz < 0) iz = 0 ; 206 if(iz < 0) iz = 0 ; 204 else if(iz > 91) iz = 91 ; 207 else if(iz > 91) iz = 91 ; 205 vF = vFermi[iz] ; 208 vF = vFermi[iz] ; >> 209 lF = lFactor[iz] ; 206 210 207 } else { 211 } else { 208 for (G4int iel=0; iel<NumberOfElements; ie 212 for (G4int iel=0; iel<NumberOfElements; iel++) 209 { 213 { 210 const G4Element* element = (*theElemen 214 const G4Element* element = (*theElementVector)[iel] ; 211 G4double z2 = element->GetZ() ; 215 G4double z2 = element->GetZ() ; 212 const G4double weight = theAtomicNumDe 216 const G4double weight = theAtomicNumDensityVector[iel] ; 213 norm += weight ; 217 norm += weight ; 214 z += z2 * weight ; 218 z += z2 * weight ; 215 G4int iz = G4int(z2) - 1 ; 219 G4int iz = G4int(z2) - 1 ; 216 if(iz < 0) iz = 0 ; 220 if(iz < 0) iz = 0 ; 217 else if(iz > 91) iz =91 ; 221 else if(iz > 91) iz =91 ; 218 vF += vFermi[iz] * weight ; 222 vF += vFermi[iz] * weight ; >> 223 lF += lFactor[iz] * weight ; 219 } 224 } 220 if (norm > 0.0) { << 225 z /= norm ; 221 z /= norm ; << 226 vF /= norm ; 222 vF /= norm ; << 227 lF /= norm ; 223 } << 224 } 228 } 225 229 226 // Helium ion case 230 // Helium ion case 227 if( ionCharge < 2.5 ) { 231 if( ionCharge < 2.5 ) { 228 232 229 G4double e = std::log(std::max(1.0, kineti << 233 G4double e = log(G4std::max(1.0, kineticEnergy / (keV*theHeMassAMU) )) ; 230 G4double x = c[0] ; 234 G4double x = c[0] ; 231 G4double y = 1.0 ; 235 G4double y = 1.0 ; 232 for (G4int i=1; i<6; i++) { 236 for (G4int i=1; i<6; i++) { 233 y *= e ; 237 y *= e ; 234 x += y * c[i] ; 238 x += y * c[i] ; 235 } 239 } 236 G4double q = 7.6 - e ; << 240 G4double q = 7.6 - e ; 237 q = 1.0 + ( 0.007 + 0.00005 * z ) * G4Exp( << 241 q = 1.0 + ( 0.007 + 0.00005 * z ) * exp( -q*q ) ; 238 return 4.0 * q * q * (1.0 - G4Exp(-x)) ; << 242 return 4.0 * q * q * (1.0 - exp(-x)) ; 239 243 240 // Heavy ion case 244 // Heavy ion case 241 } else { 245 } else { 242 246 243 // v1 is ion velocity in vF unit 247 // v1 is ion velocity in vF unit 244 G4double v1{0.0}, v2{0.0}; << 248 G4double v1 = 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 ; 249 G4double y ; 250 G4double z13 = std::pow(ionCharge, 0.3333) << 250 G4double z13 = pow(ionCharge, 0.3333) ; 251 251 252 // Faster than Fermi velocity 252 // Faster than Fermi velocity 253 if ( v1 > 1.0 ) { 253 if ( v1 > 1.0 ) { 254 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / 254 y = vF * v1 * ( 1.0 + 0.2 / (v1*v1) ) / (z13*z13) ; 255 255 256 // Slower than Fermi velocity 256 // Slower than Fermi velocity 257 } else { 257 } else { 258 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + 258 y = 0.6923 * vF * (1.0 + 2.0*v1*v1/3.0 + v1*v1*v1*v1/15.0) / (z13*z13) ; 259 } 259 } 260 260 261 G4double y3 = std::pow(y, 0.3) ; << 261 G4double y3 = pow(y, 0.3) ; 262 G4double q = 1.0 - G4Exp( 0.803*y3 - 1.316 << 262 G4double q = 1.0 - exp( 0.803*y3 - 1.3167*y3*y3 - 263 0.38157*y - 0.0089 << 263 0.38157*y - 0.008983*y*y ) ; 264 if( q < 0.0 ) q = 0.0 ; 264 if( q < 0.0 ) q = 0.0 ; 265 265 266 G4double sLocal = 7.6 - std::log(std::max << 266 G4double s = 7.6 - log(G4std::max(1.0, reducedEnergy/keV)) ; 267 sLocal = 1.0 + ( 0.18 + 0.0015 * z ) * G4E << 267 s = 1.0 + ( 0.18 + 0.0015 * z ) * exp( -s*s )/ (ionCharge*ionCharge) ; 268 268 269 // Screen length according to 269 // Screen length according to 270 // J.F.Ziegler and J.M.Manoyan, The stoppi 270 // J.F.Ziegler and J.M.Manoyan, The stopping of ions in compaunds, 271 // Nucl. Inst. & Meth. in Phys. Res. B35 ( 271 // Nucl. Inst. & Meth. in Phys. Res. B35 (1988) 215-228. 272 272 273 G4double lambda = 10.0 * vF * std::pow(1.0 << 273 G4double lambda = 10.0 * vF * pow(1.0-q, 0.6667) / (z13 * (6.0 + q)) ; 274 G4double qeff = ionCharge * sLocal * << 274 G4double qeff = ionCharge * s * 275 ( q + 0.5*(1.0-q) * std::log(1.0 + lambd << 275 ( q + 0.5*(1.0-q) * log(1.0 + lambda*lambda) / (vF*vF) ) ; 276 if( 0.1 > qeff ) qeff = 0.1 ; << 276 if( 1.0 > qeff ) qeff = 1.0 ; 277 return qeff*qeff ; << 277 return qeff*qeff ; 278 } 278 } 279 } 279 } >> 280 >> 281 280 282