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 // ******************************************************************** >> 22 // >> 23 // >> 24 // $Id: G4IonisParamMat.cc,v 1.12 2004/09/06 15:50:47 vnivanch Exp $ >> 25 // GEANT4 tag $Name: materials-V06-02-01 $ >> 26 // >> 27 // >> 28 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 25 29 26 // 09-07-98, data moved from G4Material, M.Mai << 30 // 06-09-04, Factor 2 to shell correction term (V.Ivanchenko) 27 // 18-07-98, bug corrected in ComputeDensityEf << 28 // 16-01-01, bug corrected in ComputeDensityEf << 29 // 08-02-01, fShellCorrectionVector correctly << 30 // 28-10-02, add setMeanExcitationEnergy (V.Iv 31 // 28-10-02, add setMeanExcitationEnergy (V.Ivanchenko) 31 // 06-09-04, factor 2 to shell correction term << 32 // 08-02-01, fShellCorrectionVector correctly handled (mma) 32 // 10-05-05, add a missing coma in FindMeanExc << 33 // 16-01-01, bug corrected in ComputeDensityEffect() E100eV (L.Urban) 33 // 27-09-07, add computation of parameters for << 34 // 18-07-98, bug corrected in ComputeDensityEffect() for gas 34 // 04-03-08, remove reference to G4NistManager << 35 // 09-07-98, data moved from G4Material, M.Maire 35 // 30-10-09, add G4DensityEffectData class and << 36 36 37 #include "G4IonisParamMat.hh" << 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 38 38 39 #include "G4AtomicShells.hh" << 39 #include "G4IonisParamMat.hh" 40 #include "G4AutoLock.hh" << 41 #include "G4DensityEffectData.hh" << 42 #include "G4Exp.hh" << 43 #include "G4Log.hh" << 44 #include "G4Material.hh" 40 #include "G4Material.hh" 45 #include "G4NistManager.hh" << 46 #include "G4PhysicalConstants.hh" << 47 #include "G4Pow.hh" << 48 #include "G4SystemOfUnits.hh" << 49 << 50 G4DensityEffectData* G4IonisParamMat::fDensity << 51 << 52 namespace << 53 { << 54 G4Mutex ionisMutex = G4MUTEX_INITIALIZER; << 55 } << 56 41 57 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 58 43 59 G4IonisParamMat::G4IonisParamMat(const G4Mater << 44 G4IonisParamMat::G4IonisParamMat(G4Material* material) >> 45 :fMaterial(material) 60 { 46 { 61 fBirks = 0.; << 62 fMeanEnergyPerIon = 0.0; << 63 twoln10 = 2. * G4Pow::GetInstance()->logZ(10 << 64 << 65 // minimal set of default parameters for den << 66 fCdensity = 0.0; << 67 fD0density = 0.0; << 68 fAdjustmentFactor = 1.0; << 69 if (fDensityData == nullptr) { << 70 fDensityData = new G4DensityEffectData(); << 71 } << 72 fDensityEffectCalc = nullptr; << 73 << 74 // compute parameters << 75 ComputeMeanParameters(); 47 ComputeMeanParameters(); 76 ComputeDensityEffectParameters(); << 48 ComputeDensityEffect(); 77 ComputeFluctModel(); 49 ComputeFluctModel(); 78 ComputeIonParameters(); << 79 } << 80 << 81 //....oooOO0OOooo........oooOO0OOooo........oo << 82 << 83 G4IonisParamMat::~G4IonisParamMat() << 84 { << 85 delete fDensityEffectCalc; << 86 delete[] fShellCorrectionVector; << 87 delete fDensityData; << 88 fDensityData = nullptr; << 89 fShellCorrectionVector = nullptr; << 90 fDensityEffectCalc = nullptr; << 91 } << 92 << 93 //....oooOO0OOooo........oooOO0OOooo........oo << 94 << 95 G4double G4IonisParamMat::GetDensityCorrection << 96 { << 97 // x = log10(beta*gamma) << 98 G4double y = 0.0; << 99 if (x < fX0density) { << 100 if (fD0density > 0.0) { << 101 y = fD0density * G4Exp(twoln10 * (x - fX << 102 } << 103 } << 104 else if (x >= fX1density) { << 105 y = twoln10 * x - fCdensity; << 106 } << 107 else { << 108 y = twoln10 * x - fCdensity + fAdensity * << 109 } << 110 return y; << 111 } 50 } 112 51 113 //....oooOO0OOooo........oooOO0OOooo........oo 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 114 53 115 void G4IonisParamMat::ComputeMeanParameters() 54 void G4IonisParamMat::ComputeMeanParameters() 116 { 55 { 117 // compute mean excitation energy and shell 56 // compute mean excitation energy and shell correction vector 118 fTaul = (*(fMaterial->GetElementVector()))[0 << 119 57 120 std::size_t nElements = fMaterial->GetNumber << 58 fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul(); 121 const G4ElementVector* elmVector = fMaterial << 122 const G4double* nAtomsPerVolume = fMaterial- << 123 59 124 fMeanExcitationEnergy = FindMeanExcitationEn << 60 fMeanExcitationEnergy = 0.; 125 fLogMeanExcEnergy = 0.; 61 fLogMeanExcEnergy = 0.; 126 62 127 // Chemical formula defines mean excitation << 128 if (fMeanExcitationEnergy > 0.0) { << 129 fLogMeanExcEnergy = G4Log(fMeanExcitationE << 130 63 131 // Compute average << 64 for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) { 132 } << 65 fLogMeanExcEnergy += (fMaterial->GetVecNbOfAtomsPerVolume())[i] 133 else { << 66 *((*(fMaterial->GetElementVector()))[i]->GetZ()) 134 for (std::size_t i = 0; i < nElements; ++i << 67 *log((*(fMaterial->GetElementVector()))[i]->GetIonisation() 135 const G4Element* elm = (*elmVector)[i]; << 68 ->GetMeanExcitationEnergy()); 136 fLogMeanExcEnergy += << 137 nAtomsPerVolume[i] * elm->GetZ() * G4L << 138 } << 139 fLogMeanExcEnergy /= fMaterial->GetTotNbOf << 140 fMeanExcitationEnergy = G4Exp(fLogMeanExcE << 141 } 69 } 142 70 >> 71 fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume(); >> 72 fMeanExcitationEnergy = exp(fLogMeanExcEnergy); >> 73 143 fShellCorrectionVector = new G4double[3]; 74 fShellCorrectionVector = new G4double[3]; 144 75 145 for (G4int j = 0; j <= 2; ++j) { << 76 for (G4int j=0; j<=2; j++) >> 77 { 146 fShellCorrectionVector[j] = 0.; 78 fShellCorrectionVector[j] = 0.; 147 79 148 for (std::size_t k = 0; k < nElements; ++k << 80 for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) { 149 fShellCorrectionVector[j] += << 81 fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k] 150 nAtomsPerVolume[k] * (((*elmVector)[k] << 82 *((*(fMaterial->GetElementVector()))[k]->GetIonisation() >> 83 ->GetShellCorrectionVector()[j]); 151 } 84 } 152 fShellCorrectionVector[j] *= 2.0 / fMateri << 85 fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume(); 153 } << 86 } 154 } 87 } 155 88 156 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 157 << 90 158 G4DensityEffectData* G4IonisParamMat::GetDensi << 91 void G4IonisParamMat::ComputeDensityEffect() 159 return fDensityData; << 160 } << 161 << 162 //....oooOO0OOooo........oooOO0OOooo........oo << 163 << 164 void G4IonisParamMat::ComputeDensityEffectPara << 165 { 92 { 166 G4State State = fMaterial->GetState(); << 93 // Compute parameters for the density effect correction in DE/Dx formula. 167 G4double density = fMaterial->GetDensity(); << 94 // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971) 168 << 169 // Check if density effect data exist in the << 170 // R.M. Sternheimer, Atomic Data and Nuclear << 171 // or is assign to one of data set in this t << 172 G4int idx = fDensityData->GetIndex(fMaterial << 173 auto nelm = (G4int)fMaterial->GetNumberOfEle << 174 G4int Z0 = ((*(fMaterial->GetElementVector() << 175 const G4Material* bmat = fMaterial->GetBaseM << 176 G4NistManager* nist = G4NistManager::Instanc << 177 << 178 // arbitrary empirical limits << 179 // parameterisation with very different dens << 180 static const G4double corrmax = 1.; << 181 static const G4double massfracmax = 0.9; << 182 << 183 // for simple non-NIST materials << 184 G4double corr = 0.0; << 185 << 186 if (idx < 0 && 1 == nelm) { << 187 G4int z = (1 == Z0 && State == kStateLiqui << 188 idx = fDensityData->GetElementIndex(z); << 189 << 190 // Correction for base material or for non << 191 // Except cases of very different density << 192 if (idx >= 0 && 0 < z) { << 193 G4double dens = nist->GetNominalDensity( << 194 if (dens <= 0.0) { << 195 idx = -1; << 196 } << 197 else { << 198 corr = G4Log(dens / density); << 199 if (std::abs(corr) > corrmax) { << 200 idx = -1; << 201 } << 202 } << 203 } << 204 } << 205 // for base material case << 206 if (idx < 0 && nullptr != bmat) { << 207 idx = fDensityData->GetIndex(bmat->GetName << 208 if (idx >= 0) { << 209 corr = G4Log(bmat->GetDensity() / densit << 210 if (std::abs(corr) > corrmax) { << 211 idx = -1; << 212 } << 213 } << 214 } << 215 95 216 // for compound non-NIST materials with one << 96 const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius; 217 if (idx < 0 && 1 < nelm) { << 97 const G4double twoln10 = 2.*log(10.); 218 const G4double tot = fMaterial->GetTotNbOf << 219 for (G4int i = 0; i < nelm; ++i) { << 220 const G4double frac = fMaterial->GetVecN << 221 if (frac > massfracmax) { << 222 Z0 = ((*(fMaterial->GetElementVector() << 223 idx = fDensityData->GetElementIndex(Z0 << 224 G4double dens = nist->GetNominalDensit << 225 if (idx >= 0 && dens > 0.0) { << 226 corr = G4Log(dens / density); << 227 if (std::abs(corr) > corrmax) { << 228 idx = -1; << 229 } << 230 else { << 231 break; << 232 } << 233 } << 234 } << 235 } << 236 } << 237 98 238 if (idx >= 0) { << 99 G4int icase; 239 // Take parameters for the density effect << 100 240 // R.M. Sternheimer et al. Density Effect << 101 fCdensity = 1. + log(fMeanExcitationEnergy*fMeanExcitationEnergy 241 // of Charged Particles in Various Substan << 102 /(Cd2*fMaterial->GetTotNbOfElectPerVolume())); 242 // Atom. Data Nucl. Data Tabl. 30 (1984) 2 << 243 << 244 fCdensity = fDensityData->GetCdensity(idx) << 245 fMdensity = fDensityData->GetMdensity(idx) << 246 fAdensity = fDensityData->GetAdensity(idx) << 247 fX0density = fDensityData->GetX0density(id << 248 fX1density = fDensityData->GetX1density(id << 249 fD0density = fDensityData->GetDelta0densit << 250 fPlasmaEnergy = fDensityData->GetPlasmaEne << 251 fAdjustmentFactor = fDensityData->GetAdjus << 252 << 253 // parameter C is computed and not taken f << 254 // fCdensity = 1. + 2*G4Log(fMeanExcitatio << 255 // G4cout << "IonisParamMat: " << fMateria << 256 // << " Cst= " << Cdensity << " C= " << 257 << 258 // correction on nominal density << 259 fCdensity += corr; << 260 fX0density += corr / twoln10; << 261 fX1density += corr / twoln10; << 262 } << 263 else { << 264 static const G4double Cd2 = 4 * CLHEP::pi << 265 fPlasmaEnergy = std::sqrt(Cd2 * fMaterial- << 266 << 267 // Compute parameters for the density effe << 268 // The parametrization is from R.M. Sternh << 269 G4int icase; << 270 << 271 fCdensity = 1. + 2 * G4Log(fMeanExcitation << 272 // << 273 // condensed materials << 274 // << 275 if ((State == kStateSolid) || (State == kS << 276 static const G4double E100eV = 100. * CL << 277 static const G4double ClimiS[] = {3.681, << 278 static const G4double X0valS[] = {1.0, 1 << 279 static const G4double X1valS[] = {2.0, 3 << 280 103 281 if (fMeanExcitationEnergy < E100eV) { << 104 // 282 icase = 0; << 105 // condensed materials 283 } << 106 // 284 else { << 107 G4State State = fMaterial->GetState(); 285 icase = 1; << 108 286 } << 109 if ((State == kStateSolid)||(State == kStateLiquid)) { 287 110 288 if (fCdensity < ClimiS[icase]) { << 111 const G4double E100eV = 100.*eV; 289 fX0density = 0.2; << 112 const G4double ClimiS[] = {3.681 , 5.215 }; 290 } << 113 const G4double X0valS[] = {1.0 , 1.5 }; 291 else { << 114 const G4double X1valS[] = {2.0 , 3.0 }; 292 fX0density = 0.326 * fCdensity - X0val << 115 293 } << 116 if(fMeanExcitationEnergy < E100eV) icase = 0; >> 117 else icase = 1; >> 118 >> 119 if(fCdensity < ClimiS[icase]) fX0density = 0.2; >> 120 else fX0density = 0.326*fCdensity-X0valS[icase]; >> 121 >> 122 fX1density = X1valS[icase] ; fMdensity = 3.0; >> 123 >> 124 //special: Hydrogen >> 125 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) { >> 126 fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949; >> 127 } >> 128 } >> 129 >> 130 // >> 131 // gases >> 132 // >> 133 if (State == kStateGas) { >> 134 >> 135 const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804}; >> 136 const G4double X0valG[] = { 1.6 , 1.7 , 1.8 , 1.9 , 2.0 , 2.0 }; >> 137 const G4double X1valG[] = { 4.0 , 4.0 , 4.0 , 4.0 , 4.0 , 5.0 }; >> 138 >> 139 icase = 5; >> 140 fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ; >> 141 while((icase > 0)&&(fCdensity < ClimiG[icase])) icase-- ; >> 142 fX0density = X0valG[icase] ; fX1density = X1valG[icase] ; >> 143 >> 144 //special: Hydrogen >> 145 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) { >> 146 fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754; >> 147 } >> 148 >> 149 //special: Helium >> 150 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==2.)) { >> 151 fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297; >> 152 } >> 153 >> 154 // change parameters if the gas is not in STP. >> 155 // For the correction the density(STP) is needed. Density(STP) is calculated here : >> 156 >> 157 G4double Density = fMaterial->GetDensity(); >> 158 G4double Pressure = fMaterial->GetPressure(); >> 159 G4double Temp = fMaterial->GetTemperature(); >> 160 >> 161 G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature); >> 162 >> 163 G4double ParCorr = log(Density/DensitySTP); >> 164 >> 165 fCdensity -= ParCorr; >> 166 fX0density -= ParCorr/twoln10; >> 167 fX1density -= ParCorr/twoln10; >> 168 } >> 169 >> 170 G4double Xa = fCdensity/twoln10; >> 171 fAdensity = twoln10*(Xa-fX0density) >> 172 /pow((fX1density-fX0density),fMdensity); >> 173 } 294 174 295 fX1density = X1valS[icase]; << 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 296 fMdensity = 3.0; << 297 176 298 // special: Hydrogen << 177 void G4IonisParamMat::ComputeFluctModel() 299 if (1 == nelm && 1 == Z0) { << 178 { 300 fX0density = 0.425; << 179 // compute parameters for the energy loss fluctuation model 301 fX1density = 2.0; << 302 fMdensity = 5.949; << 303 } << 304 } << 305 else { << 306 // << 307 // gases << 308 // << 309 fMdensity = 3.; << 310 fX1density = 4.0; << 311 180 312 if (fCdensity <= 10.) { << 181 // need an 'effective Z' ????? 313 fX0density = 1.6; << 182 G4double Zeff = 0.; 314 } << 183 for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) 315 else if (fCdensity <= 10.5) { << 184 Zeff += (fMaterial->GetFractionVector())[i] 316 fX0density = 1.7; << 185 *((*(fMaterial->GetElementVector()))[i]->GetZ()); 317 } << 186 318 else if (fCdensity <= 11.0) { << 187 if (Zeff > 2.) fF2fluct = 2./Zeff ; 319 fX0density = 1.8; << 188 else fF2fluct = 0.; 320 } << 189 321 else if (fCdensity <= 11.5) { << 190 fF1fluct = 1. - fF2fluct; 322 fX0density = 1.9; << 191 fEnergy2fluct = 10.*Zeff*Zeff*eV; 323 } << 192 fLogEnergy2fluct = log(fEnergy2fluct); 324 else if (fCdensity <= 12.25) { << 193 fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct) 325 fX0density = 2.0; << 194 /fF1fluct; 326 } << 195 fEnergy1fluct = exp(fLogEnergy1fluct); 327 else if (fCdensity <= 13.804) { << 196 fEnergy0fluct = 10.*eV; 328 fX0density = 2.0; << 197 fRateionexcfluct = 0.4; 329 fX1density = 5.0; << 198 } 330 } << 331 else { << 332 fX0density = 0.326 * fCdensity - 2.5; << 333 fX1density = 5.0; << 334 } << 335 199 336 // special: Hydrogen << 200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 337 if (1 == nelm && 1 == Z0) { << 338 fX0density = 1.837; << 339 fX1density = 3.0; << 340 fMdensity = 4.754; << 341 } << 342 201 343 // special: Helium << 202 void G4IonisParamMat::SetMeanExcitationEnergy(G4double value) 344 if (1 == nelm && 2 == Z0) { << 203 { 345 fX0density = 2.191; << 204 if(value == fMeanExcitationEnergy || value <= 0.0) return; 346 fX1density = 3.0; << 347 fMdensity = 3.297; << 348 } << 349 } << 350 } << 351 205 352 // change parameters if the gas is not in ST << 206 if(fMeanExcitationEnergy > 0.0) { 353 // For the correction the density(STP) is ne << 354 // Density(STP) is calculated here : << 355 207 356 if (State == kStateGas) { << 208 G4cout << "G4Material: Mean excitation energy is changed for " << fMaterial->GetName() 357 G4double Pressure = fMaterial->GetPressure << 209 << " Iold= " << fMeanExcitationEnergy/eV 358 G4double Temp = fMaterial->GetTemperature( << 210 << "eV; Inew= " << value/eV << " eV;" >> 211 << G4endl; >> 212 } >> 213 >> 214 fMeanExcitationEnergy = value; >> 215 fLogMeanExcEnergy = log(value); >> 216 ComputeDensityEffect(); >> 217 ComputeFluctModel(); >> 218 >> 219 } 359 220 360 G4double DensitySTP = density * STP_Pressu << 221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 361 222 362 G4double ParCorr = G4Log(density / Density << 223 G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula) >> 224 { 363 225 364 fCdensity -= ParCorr; << 226 // The data on mean excitation energy for compaunds 365 fX0density -= ParCorr / twoln10; << 227 // from "Stopping Powers for Electrons and Positrons" 366 fX1density -= ParCorr / twoln10; << 228 // ICRU Report N#37, 1984 (energy in eV) 367 } << 368 229 369 // fAdensity parameter can be fixed for not << 230 const size_t numberOfMolecula = 79 ; 370 if (0.0 == fD0density) { << 231 371 G4double Xa = fCdensity / twoln10; << 232 static G4String name[numberOfMolecula] = { 372 fAdensity = twoln10 * (Xa - fX0density) / << 233 >> 234 // gas >> 235 "NH_3", "C_4H_10", "CO_2", "C_2H_6", "C_7H_16", >> 236 "C_6H_14", "CH_4", "NO", "N_2O", "C_8H_18", >> 237 "C_5H_12", "C_3H_8", "H_2O-Gas", >> 238 >> 239 // liquid >> 240 "C_3H_6O", "C_6H_5NH_2", "C_6H_6", "C_4H_9OH", "CCl_4", >> 241 "C_6H_5Cl", "CHCl_3", "C_6H_12", "C_6H_4Cl_2", "C_4Cl_2H_8O", >> 242 "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", "C_3H_5(OH)_3","C_7H_16", >> 243 "C_6H_14", "CH_3OH", "C_6H_5NO_2","C_5H_12", "C_3H_7OH", >> 244 "C_5H_5N", "C_8H_8", "C_2Cl_4", "C_7H_8", "C_2Cl_3H", >> 245 "H_2O", "C_8H_10" >> 246 >> 247 //solid >> 248 "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO)-nylon", "C_25H_52", >> 249 "(C_2H_4)-Polyethylene", "(C_5H_8O-2)-Polymethil_Methacrylate", >> 250 "(C_8H_8)-Polystyrene", "A-150-tissue", "Al_2O_3", "CaF_2", >> 251 "LiF", "Photo_Emulsion", "(C_2F_4)-Teflon", "SiO_2" >> 252 >> 253 } ; >> 254 >> 255 static G4double meanExcitation[numberOfMolecula] = { >> 256 >> 257 53.7, 48.3, 85.0, 45.4, 49.2, >> 258 49.1, 41.7, 87.8, 84.9, 49.5, >> 259 48.2, 47.1, 71.6, >> 260 >> 261 64.2, 66.2, 63.4, 59.9, 166.3, >> 262 89.1, 156.0, 56.4, 106.5, 103.3, >> 263 111.9, 60.0, 62.9, 72.6, 54.4, >> 264 54.0, 67.6, 75.8, 53.6, 61.1, >> 265 66.2, 64.0, 159.2, 62.5, 148.1, >> 266 75.0, 61.8, >> 267 >> 268 71.4, 75.0, 63.9, 48.3, 57.4, >> 269 74.0, 68.7, 65.1, 145.2, 166., >> 270 94.0, 331.0, 99.1, 139.2 >> 271 >> 272 } ; >> 273 >> 274 G4double x = fMeanExcitationEnergy; >> 275 >> 276 for(size_t i=0; i<numberOfMolecula; i++) { >> 277 if(chFormula == name[i]) { >> 278 x = meanExcitation[i]*eV; >> 279 break; >> 280 } 373 } 281 } >> 282 return x; 374 } 283 } 375 284 376 //....oooOO0OOooo........oooOO0OOooo........oo 285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 377 286 378 void G4IonisParamMat::ComputeFluctModel() << 287 G4IonisParamMat::~G4IonisParamMat() 379 { 288 { 380 // compute parameters for the energy loss fl << 289 if (fShellCorrectionVector) delete [] fShellCorrectionVector; 381 // needs an 'effective Z' << 382 G4double Zeff = 0.; << 383 for (std::size_t i = 0; i < fMaterial->GetNu << 384 Zeff += (fMaterial->GetFractionVector())[i << 385 } << 386 if (Zeff > 2.1) { << 387 fF2fluct = 2.0 / Zeff; << 388 fF1fluct = 1. - fF2fluct; << 389 fEnergy2fluct = 10. * Zeff * Zeff * CLHEP: << 390 fLogEnergy2fluct = G4Log(fEnergy2fluct); << 391 fLogEnergy1fluct = (fLogMeanExcEnergy - fF << 392 } else if (Zeff > 1.1) { << 393 fF2fluct = 0.0; << 394 fF1fluct = 1.0; << 395 fEnergy2fluct = 40. * CLHEP::eV; << 396 fLogEnergy2fluct = G4Log(fEnergy2fluct); << 397 fLogEnergy1fluct = fLogMeanExcEnergy; << 398 } else { << 399 fF2fluct = 0.0; << 400 fF1fluct = 1.0; << 401 fEnergy2fluct = 10. * CLHEP::eV; << 402 fLogEnergy2fluct = G4Log(fEnergy2fluct); << 403 fLogEnergy1fluct = fLogMeanExcEnergy; << 404 } << 405 fEnergy1fluct = G4Exp(fLogEnergy1fluct); << 406 fEnergy0fluct = 10. * CLHEP::eV; << 407 fRateionexcfluct = 0.4; << 408 } 290 } 409 291 410 //....oooOO0OOooo........oooOO0OOooo........oo 292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 411 293 412 void G4IonisParamMat::ComputeIonParameters() << 294 G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right) 413 { 295 { 414 // get elements in the actual material, << 296 *this = right; 415 const G4ElementVector* theElementVector = fM << 416 const G4double* theAtomicNumDensityVector = << 417 const auto NumberOfElements = (G4int)fMateri << 418 << 419 // loop for the elements in the material << 420 // to find out average values Z, vF, lF << 421 G4double z(0.0), vF(0.0), lF(0.0), a23(0.0); << 422 << 423 G4Pow* g4pow = G4Pow::GetInstance(); << 424 if (1 == NumberOfElements) { << 425 const G4Element* element = (*theElementVec << 426 z = element->GetZ(); << 427 vF = element->GetIonisation()->GetFermiVel << 428 lF = element->GetIonisation()->GetLFactor( << 429 a23 = 1.0 / g4pow->A23(element->GetN()); << 430 } << 431 else { << 432 G4double norm(0.0); << 433 for (G4int iel = 0; iel < NumberOfElements << 434 const G4Element* element = (*theElementV << 435 const G4double weight = theAtomicNumDens << 436 norm += weight; << 437 z += element->GetZ() * weight; << 438 vF += element->GetIonisation()->GetFermi << 439 lF += element->GetIonisation()->GetLFact << 440 a23 += weight / g4pow->A23(element->GetN << 441 } << 442 if (norm > 0.0) { norm = 1.0/norm; } << 443 z *= norm; << 444 vF *= norm; << 445 lF *= norm; << 446 a23 *= norm; << 447 } << 448 fZeff = z; << 449 fLfactor = lF; << 450 fFermiEnergy = 25. * CLHEP::keV * vF * vF; << 451 fInvA23 = a23; << 452 } 297 } 453 298 454 //....oooOO0OOooo........oooOO0OOooo........oo 299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 455 300 456 void G4IonisParamMat::SetMeanExcitationEnergy( << 301 const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right) 457 { 302 { 458 if (value == fMeanExcitationEnergy || value << 303 if (this != &right) 459 return; << 304 { 460 } << 305 fMaterial = right.fMaterial; 461 if (G4NistManager::Instance()->GetVerbose() << 306 fMeanExcitationEnergy = right.fMeanExcitationEnergy; 462 G4cout << "G4Material: Mean excitation ene << 307 fLogMeanExcEnergy = right.fLogMeanExcEnergy; 463 << " Iold= " << fMeanExcitationEner << 308 if (fShellCorrectionVector) delete [] fShellCorrectionVector; 464 << G4endl; << 309 fShellCorrectionVector = new G4double[3]; 465 } << 310 fShellCorrectionVector[0] = right.fShellCorrectionVector[0]; 466 << 311 fShellCorrectionVector[1] = right.fShellCorrectionVector[1]; 467 fMeanExcitationEnergy = value; << 312 fShellCorrectionVector[2] = right.fShellCorrectionVector[2]; 468 << 313 fTaul = right.fTaul; 469 // add corrections to density effect << 314 fCdensity = right.fCdensity; 470 G4double newlog = G4Log(value); << 315 fMdensity = right.fMdensity; 471 G4double corr = 2 * (newlog - fLogMeanExcEne << 316 fAdensity = right.fAdensity; 472 fCdensity += corr; << 317 fX0density = right.fX0density; 473 fX0density += corr / twoln10; << 318 fX1density = right.fX1density; 474 fX1density += corr / twoln10; << 319 fF1fluct = right.fF1fluct; 475 << 320 fF2fluct = right.fF2fluct; 476 // recompute parameters of fluctuation model << 321 fEnergy1fluct = right.fEnergy1fluct; 477 fLogMeanExcEnergy = newlog; << 322 fLogEnergy1fluct = right.fLogEnergy1fluct; 478 ComputeFluctModel(); << 323 fEnergy2fluct = right.fEnergy2fluct; >> 324 fLogEnergy2fluct = right.fLogEnergy2fluct; >> 325 fEnergy0fluct = right.fEnergy0fluct; >> 326 fRateionexcfluct = right.fRateionexcfluct; >> 327 } >> 328 return *this; 479 } 329 } 480 330 481 //....oooOO0OOooo........oooOO0OOooo........oo 331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 482 332 483 void G4IonisParamMat::SetDensityEffectParamete << 333 G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const 484 G4double cd, G4double md, G4double ad, G4dou << 485 { 334 { 486 // no check on consistence of user parameter << 335 return (this == (G4IonisParamMat*) &right); 487 G4AutoLock l(&ionisMutex); << 488 fCdensity = cd; << 489 fMdensity = md; << 490 fAdensity = ad; << 491 fX0density = x0; << 492 fX1density = x1; << 493 fD0density = d0; << 494 l.unlock(); << 495 } 336 } 496 337 497 //....oooOO0OOooo........oooOO0OOooo........oo 338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 498 339 499 void G4IonisParamMat::SetDensityEffectParamete << 340 G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const 500 { << 341 { 501 G4AutoLock l(&ionisMutex); << 342 return (this != (G4IonisParamMat*) &right); 502 const G4IonisParamMat* ipm = bmat->GetIonisa << 503 fCdensity = ipm->GetCdensity(); << 504 fMdensity = ipm->GetMdensity(); << 505 fAdensity = ipm->GetAdensity(); << 506 fX0density = ipm->GetX0density(); << 507 fX1density = ipm->GetX1density(); << 508 fD0density = ipm->GetD0density(); << 509 << 510 G4double corr = G4Log(bmat->GetDensity() / f << 511 fCdensity += corr; << 512 fX0density += corr / twoln10; << 513 fX1density += corr / twoln10; << 514 l.unlock(); << 515 } << 516 << 517 //....oooOO0OOooo........oooOO0OOooo........oo << 518 << 519 void G4IonisParamMat::ComputeDensityEffectOnFl << 520 { << 521 if (val) { << 522 if (nullptr == fDensityEffectCalc) { << 523 G4int n = 0; << 524 for (std::size_t i = 0; i < fMaterial->G << 525 const G4int Z = fMaterial->GetElement( << 526 n += G4AtomicShells::GetNumberOfShells << 527 } << 528 // The last level is the conduction leve << 529 // make a dummy conductor level with zer << 530 fDensityEffectCalc = new G4DensityEffect << 531 } << 532 } << 533 else { << 534 delete fDensityEffectCalc; << 535 fDensityEffectCalc = nullptr; << 536 } << 537 } 343 } 538 344 539 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo.... 540 346 541 G4double G4IonisParamMat::FindMeanExcitationEn << 542 { << 543 G4double res = 0.0; << 544 // data from density effect data << 545 if (fDensityData != nullptr) { << 546 G4int idx = fDensityData->GetIndex(mat->Ge << 547 if (idx >= 0) { << 548 res = fDensityData->GetMeanIonisationPot << 549 } << 550 } << 551 << 552 // The data on mean excitation energy for co << 553 // from "Stopping Powers for Electrons and P << 554 // ICRU Report N#37, 1984 (energy in eV) << 555 // this value overwrites Density effect data << 556 G4String chFormula = mat->GetChemicalFormula << 557 if (! chFormula.empty()) { << 558 static const size_t numberOfMolecula = 54; << 559 // clang-format off << 560 static const G4String name[numberOfMolecul << 561 // gas 0 - 12 << 562 "NH_3", "C_4H_10", "CO_2", << 563 // "G4_AMMONIA", "G4_BUTANE","G4_CARBON_ << 564 "C_6H_14-Gas", "CH_4", "NO", << 565 // "G4_N-HEXANE" , "G4_METHANE", "x", "G << 566 "C_5H_12-Gas", "C_3H_8", "H_2O-Gas", << 567 // "G4_N-PENTANE", "G4_PROPANE", "G4_WAT << 568 << 569 // liquid 13 - 39 << 570 "C_3H_6O", "C_6H_5NH_2", "C_6H_6", << 571 //"G4_ACETONE","G4_ANILINE","G4_BENZENE" << 572 "C_6H_5Cl", "CHCl_3", "C_6H_12", << 573 //"G4_CHLOROBENZENE","G4_CHLOROFORM","G4 << 574 //"G4_DICHLORODIETHYL_ETHER" << 575 "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", << 576 //"G4_1,2-DICHLOROETHANE","G4_DIETHYL_ET << 577 "C_6H_14", "CH_3OH", "C_6H_5NO_2 << 578 //"G4_N-HEXANE","G4_METHANOL","G4_NITROB << 579 "C_5H_5N", "C_8H_8", "C_2Cl_4", << 580 //"G4_PYRIDINE","G4_POLYSTYRENE","G4_TET << 581 "H_2O", "C_8H_10", << 582 // "G4_WATER", "G4_XYLENE" << 583 << 584 // solid 40 - 53 << 585 "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO << 586 // "G4_ADENINE", "G4_GUANINE", "G4_NYLON << 587 "(C_2H_4)-Polyethylene", "(C_5H_8O_2 << 588 // "G4_ETHYLENE", "G4_PLEXIGLASS" << 589 "(C_8H_8)-Polystyrene", "A-150-tiss << 590 // "G4_POLYSTYRENE", "G4_A-150_TISSUE", << 591 "LiF", "Photo_Emulsion", "(C_2F_ << 592 // "G4_LITHIUM_FLUORIDE", "G4_PHOTO_EMUL << 593 } ; << 594 << 595 static const G4double meanExcitation[numbe << 596 << 597 53.7, 48.3, 85.0, 45.4, 49.2, << 598 49.1, 41.7, 87.8, 84.9, 49.5, << 599 48.2, 47.1, 71.6, << 600 << 601 64.2, 66.2, 63.4, 59.9, 166.3, << 602 89.1, 156.0, 56.4, 106.5, 103.3, << 603 111.9, 60.0, 62.9, 72.6, 54.4, << 604 54.0, 67.6, 75.8, 53.6, 61.1, << 605 66.2, 64.0, 159.2, 62.5, 148.1, << 606 75.0, 61.8, << 607 << 608 71.4, 75.0, 63.9, 48.3, 57.4, << 609 74.0, 68.7, 65.1, 145.2, 166., << 610 94.0, 331.0, 99.1, 139.2 << 611 }; << 612 // clang-format on << 613 << 614 for (std::size_t i = 0; i < numberOfMolecu << 615 if (chFormula == name[i]) { << 616 res = meanExcitation[i] * CLHEP::eV; << 617 break; << 618 } << 619 } << 620 } << 621 return res; << 622 } << 623 347