| 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 //
>> 26 //
>> 27 // $Id: G4IonisParamMat.cc,v 1.20 2007/09/27 14:05:47 vnivanch Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
>> 29 //
>> 30 //
>> 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
25 32
26 // 09-07-98, data moved from G4Material, M.Mai 33 // 09-07-98, data moved from G4Material, M.Maire
27 // 18-07-98, bug corrected in ComputeDensityEf 34 // 18-07-98, bug corrected in ComputeDensityEffect() for gas
28 // 16-01-01, bug corrected in ComputeDensityEf 35 // 16-01-01, bug corrected in ComputeDensityEffect() E100eV (L.Urban)
29 // 08-02-01, fShellCorrectionVector correctly 36 // 08-02-01, fShellCorrectionVector correctly handled (mma)
30 // 28-10-02, add setMeanExcitationEnergy (V.Iv 37 // 28-10-02, add setMeanExcitationEnergy (V.Ivanchenko)
31 // 06-09-04, factor 2 to shell correction term << 38 // 06-09-04, factor 2 to shell correction term (V.Ivanchenko)
32 // 10-05-05, add a missing coma in FindMeanExc 39 // 10-05-05, add a missing coma in FindMeanExcitationEnergy() - Bug#746 (mma)
33 // 27-09-07, add computation of parameters for 40 // 27-09-07, add computation of parameters for ions (V.Ivanchenko)
34 // 04-03-08, remove reference to G4NistManager <<
35 // 30-10-09, add G4DensityEffectData class and <<
36 41
37 #include "G4IonisParamMat.hh" << 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
38 43
39 #include "G4AtomicShells.hh" << 44 #include "G4IonisParamMat.hh"
40 #include "G4AutoLock.hh" <<
41 #include "G4DensityEffectData.hh" <<
42 #include "G4Exp.hh" <<
43 #include "G4Log.hh" <<
44 #include "G4Material.hh" 45 #include "G4Material.hh"
45 #include "G4NistManager.hh" 46 #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 47
57 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
58 49
59 G4IonisParamMat::G4IonisParamMat(const G4Mater << 50 G4IonisParamMat::G4IonisParamMat(G4Material* material)
>> 51 : fMaterial(material)
60 { 52 {
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(); 53 ComputeMeanParameters();
76 ComputeDensityEffectParameters(); << 54 ComputeDensityEffect();
77 ComputeFluctModel(); 55 ComputeFluctModel();
78 ComputeIonParameters(); 56 ComputeIonParameters();
79 } 57 }
80 58
81 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
82 60
83 G4IonisParamMat::~G4IonisParamMat() << 61 // Fake default constructor - sets only member data and allocates memory
84 { << 62 // for usage restricted to object persistency
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 63
95 G4double G4IonisParamMat::GetDensityCorrection << 64 G4IonisParamMat::G4IonisParamMat(__void__&)
>> 65 : fMaterial(0), fShellCorrectionVector(0)
96 { 66 {
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 } 67 }
112 68
113 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
114 70
115 void G4IonisParamMat::ComputeMeanParameters() 71 void G4IonisParamMat::ComputeMeanParameters()
116 { 72 {
117 // compute mean excitation energy and shell 73 // compute mean excitation energy and shell correction vector
118 fTaul = (*(fMaterial->GetElementVector()))[0 <<
119 74
120 std::size_t nElements = fMaterial->GetNumber << 75 fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
121 const G4ElementVector* elmVector = fMaterial <<
122 const G4double* nAtomsPerVolume = fMaterial- <<
123 76
124 fMeanExcitationEnergy = FindMeanExcitationEn << 77 fMeanExcitationEnergy = 0.;
125 fLogMeanExcEnergy = 0.; 78 fLogMeanExcEnergy = 0.;
126 79
127 // Chemical formula defines mean excitation <<
128 if (fMeanExcitationEnergy > 0.0) { <<
129 fLogMeanExcEnergy = G4Log(fMeanExcitationE <<
130 80
131 // Compute average << 81 for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) {
132 } << 82 fLogMeanExcEnergy +=
133 else { << 83 (fMaterial->GetVecNbOfAtomsPerVolume())[i]
134 for (std::size_t i = 0; i < nElements; ++i << 84 *((*(fMaterial->GetElementVector()))[i]->GetZ())
135 const G4Element* elm = (*elmVector)[i]; << 85 *std::log((*(fMaterial->GetElementVector()))[i]->GetIonisation()
136 fLogMeanExcEnergy += << 86 ->GetMeanExcitationEnergy());
137 nAtomsPerVolume[i] * elm->GetZ() * G4L <<
138 } <<
139 fLogMeanExcEnergy /= fMaterial->GetTotNbOf <<
140 fMeanExcitationEnergy = G4Exp(fLogMeanExcE <<
141 } 87 }
142 88
143 fShellCorrectionVector = new G4double[3]; << 89 fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
>> 90 fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
144 91
145 for (G4int j = 0; j <= 2; ++j) { << 92 fShellCorrectionVector = new G4double[3]; //[3]
>> 93
>> 94 for (G4int j=0; j<=2; j++)
>> 95 {
146 fShellCorrectionVector[j] = 0.; 96 fShellCorrectionVector[j] = 0.;
147 97
148 for (std::size_t k = 0; k < nElements; ++k << 98 for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) {
149 fShellCorrectionVector[j] += << 99 fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k]
150 nAtomsPerVolume[k] * (((*elmVector)[k] << 100 *((*(fMaterial->GetElementVector()))[k]->GetIonisation()
>> 101 ->GetShellCorrectionVector()[j]);
151 } 102 }
152 fShellCorrectionVector[j] *= 2.0 / fMateri << 103 fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
153 } << 104 }
154 } <<
155 <<
156 //....oooOO0OOooo........oooOO0OOooo........oo <<
157 <<
158 G4DensityEffectData* G4IonisParamMat::GetDensi <<
159 return fDensityData; <<
160 } 105 }
161 106
162 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
163 << 108
164 void G4IonisParamMat::ComputeDensityEffectPara << 109 void G4IonisParamMat::ComputeDensityEffect()
165 { 110 {
166 G4State State = fMaterial->GetState(); << 111 // Compute parameters for the density effect correction in DE/Dx formula.
167 G4double density = fMaterial->GetDensity(); << 112 // 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 <<
216 // for compound non-NIST materials with one <<
217 if (idx < 0 && 1 < nelm) { <<
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 <<
238 if (idx >= 0) { <<
239 // Take parameters for the density effect <<
240 // R.M. Sternheimer et al. Density Effect <<
241 // of Charged Particles in Various Substan <<
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 <<
281 if (fMeanExcitationEnergy < E100eV) { <<
282 icase = 0; <<
283 } <<
284 else { <<
285 icase = 1; <<
286 } <<
287 <<
288 if (fCdensity < ClimiS[icase]) { <<
289 fX0density = 0.2; <<
290 } <<
291 else { <<
292 fX0density = 0.326 * fCdensity - X0val <<
293 } <<
294 <<
295 fX1density = X1valS[icase]; <<
296 fMdensity = 3.0; <<
297 <<
298 // special: Hydrogen <<
299 if (1 == nelm && 1 == Z0) { <<
300 fX0density = 0.425; <<
301 fX1density = 2.0; <<
302 fMdensity = 5.949; <<
303 } <<
304 } <<
305 else { <<
306 // <<
307 // gases <<
308 // <<
309 fMdensity = 3.; <<
310 fX1density = 4.0; <<
311 <<
312 if (fCdensity <= 10.) { <<
313 fX0density = 1.6; <<
314 } <<
315 else if (fCdensity <= 10.5) { <<
316 fX0density = 1.7; <<
317 } <<
318 else if (fCdensity <= 11.0) { <<
319 fX0density = 1.8; <<
320 } <<
321 else if (fCdensity <= 11.5) { <<
322 fX0density = 1.9; <<
323 } <<
324 else if (fCdensity <= 12.25) { <<
325 fX0density = 2.0; <<
326 } <<
327 else if (fCdensity <= 13.804) { <<
328 fX0density = 2.0; <<
329 fX1density = 5.0; <<
330 } <<
331 else { <<
332 fX0density = 0.326 * fCdensity - 2.5; <<
333 fX1density = 5.0; <<
334 } <<
335 <<
336 // special: Hydrogen <<
337 if (1 == nelm && 1 == Z0) { <<
338 fX0density = 1.837; <<
339 fX1density = 3.0; <<
340 fMdensity = 4.754; <<
341 } <<
342 <<
343 // special: Helium <<
344 if (1 == nelm && 2 == Z0) { <<
345 fX0density = 2.191; <<
346 fX1density = 3.0; <<
347 fMdensity = 3.297; <<
348 } <<
349 } <<
350 } <<
351 <<
352 // change parameters if the gas is not in ST <<
353 // For the correction the density(STP) is ne <<
354 // Density(STP) is calculated here : <<
355 113
356 if (State == kStateGas) { << 114 const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
357 G4double Pressure = fMaterial->GetPressure << 115 const G4double twoln10 = 2.*std::log(10.);
358 G4double Temp = fMaterial->GetTemperature( <<
359 116
360 G4double DensitySTP = density * STP_Pressu << 117 G4int icase;
>> 118
>> 119 fCdensity = 1. + std::log(fMeanExcitationEnergy*fMeanExcitationEnergy
>> 120 /(Cd2*fMaterial->GetTotNbOfElectPerVolume()));
361 121
362 G4double ParCorr = G4Log(density / Density << 122 //
363 << 123 // condensed materials
364 fCdensity -= ParCorr; << 124 //
365 fX0density -= ParCorr / twoln10; << 125 G4State State = fMaterial->GetState();
366 fX1density -= ParCorr / twoln10; << 126
367 } << 127 if ((State == kStateSolid)||(State == kStateLiquid)) {
368 128
369 // fAdensity parameter can be fixed for not << 129 const G4double E100eV = 100.*eV;
370 if (0.0 == fD0density) { << 130 const G4double ClimiS[] = {3.681 , 5.215 };
371 G4double Xa = fCdensity / twoln10; << 131 const G4double X0valS[] = {1.0 , 1.5 };
372 fAdensity = twoln10 * (Xa - fX0density) / << 132 const G4double X1valS[] = {2.0 , 3.0 };
373 } << 133
>> 134 if(fMeanExcitationEnergy < E100eV) icase = 0;
>> 135 else icase = 1;
>> 136
>> 137 if(fCdensity < ClimiS[icase]) fX0density = 0.2;
>> 138 else fX0density = 0.326*fCdensity-X0valS[icase];
>> 139
>> 140 fX1density = X1valS[icase] ; fMdensity = 3.0;
>> 141
>> 142 //special: Hydrogen
>> 143 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
>> 144 fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
>> 145 }
>> 146 }
>> 147
>> 148 //
>> 149 // gases
>> 150 //
>> 151 if (State == kStateGas) {
>> 152
>> 153 const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804};
>> 154 const G4double X0valG[] = { 1.6 , 1.7 , 1.8 , 1.9 , 2.0 , 2.0 };
>> 155 const G4double X1valG[] = { 4.0 , 4.0 , 4.0 , 4.0 , 4.0 , 5.0 };
>> 156
>> 157 icase = 5;
>> 158 fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ;
>> 159 while((icase > 0)&&(fCdensity < ClimiG[icase])) icase-- ;
>> 160 fX0density = X0valG[icase] ; fX1density = X1valG[icase] ;
>> 161
>> 162 //special: Hydrogen
>> 163 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
>> 164 fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
>> 165 }
>> 166
>> 167 //special: Helium
>> 168 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==2.)) {
>> 169 fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
>> 170 }
>> 171
>> 172 // change parameters if the gas is not in STP.
>> 173 // For the correction the density(STP) is needed.
>> 174 // Density(STP) is calculated here :
>> 175
>> 176 G4double Density = fMaterial->GetDensity();
>> 177 G4double Pressure = fMaterial->GetPressure();
>> 178 G4double Temp = fMaterial->GetTemperature();
>> 179
>> 180 G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
>> 181
>> 182 G4double ParCorr = std::log(Density/DensitySTP);
>> 183
>> 184 fCdensity -= ParCorr;
>> 185 fX0density -= ParCorr/twoln10;
>> 186 fX1density -= ParCorr/twoln10;
>> 187 }
>> 188
>> 189 G4double Xa = fCdensity/twoln10;
>> 190 fAdensity = twoln10*(Xa-fX0density)
>> 191 /std::pow((fX1density-fX0density),fMdensity);
374 } 192 }
375 193
376 //....oooOO0OOooo........oooOO0OOooo........oo 194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
377 195
378 void G4IonisParamMat::ComputeFluctModel() 196 void G4IonisParamMat::ComputeFluctModel()
379 { 197 {
380 // compute parameters for the energy loss fl 198 // compute parameters for the energy loss fluctuation model
381 // needs an 'effective Z' << 199
>> 200 // need an 'effective Z' ?????
382 G4double Zeff = 0.; 201 G4double Zeff = 0.;
383 for (std::size_t i = 0; i < fMaterial->GetNu << 202 for (size_t i=0;i<fMaterial->GetNumberOfElements();i++)
384 Zeff += (fMaterial->GetFractionVector())[i << 203 Zeff += (fMaterial->GetFractionVector())[i]
385 } << 204 *((*(fMaterial->GetElementVector()))[i]->GetZ());
386 if (Zeff > 2.1) { << 205
387 fF2fluct = 2.0 / Zeff; << 206 if (Zeff > 2.) fF2fluct = 2./Zeff ;
388 fF1fluct = 1. - fF2fluct; << 207 else fF2fluct = 0.;
389 fEnergy2fluct = 10. * Zeff * Zeff * CLHEP: << 208
390 fLogEnergy2fluct = G4Log(fEnergy2fluct); << 209 fF1fluct = 1. - fF2fluct;
391 fLogEnergy1fluct = (fLogMeanExcEnergy - fF << 210 fEnergy2fluct = 10.*Zeff*Zeff*eV;
392 } else if (Zeff > 1.1) { << 211 fLogEnergy2fluct = std::log(fEnergy2fluct);
393 fF2fluct = 0.0; << 212 fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
394 fF1fluct = 1.0; << 213 /fF1fluct;
395 fEnergy2fluct = 40. * CLHEP::eV; << 214 fEnergy1fluct = std::exp(fLogEnergy1fluct);
396 fLogEnergy2fluct = G4Log(fEnergy2fluct); << 215 fEnergy0fluct = 10.*eV;
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; 216 fRateionexcfluct = 0.4;
408 } 217 }
409 218
410 //....oooOO0OOooo........oooOO0OOooo........oo 219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
411 220
412 void G4IonisParamMat::ComputeIonParameters() 221 void G4IonisParamMat::ComputeIonParameters()
413 { 222 {
>> 223 // compute parameters for ion transport
>> 224 // The aproximation from:
>> 225 // J.F.Ziegler, J.P. Biersack, U. Littmark
>> 226 // The Stopping and Range of Ions in Matter,
>> 227 // Vol.1, Pergamon Press, 1985
>> 228 // Fast ions or hadrons
>> 229
>> 230 static G4double vFermi[92] = {
>> 231 1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424,
>> 232 0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
>> 233 0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411,
>> 234 0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
>> 235 1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054,
>> 236 0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
>> 237 0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
>> 238 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
>> 239 0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
>> 240 1.9578, 1.0257} ;
>> 241
>> 242 static G4double lFactor[92] = {
>> 243 1.0, 1.0, 1.1, 1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
>> 244 0.82, 0.81, 0.83, 0.88, 1.0, 0.95, 0.97, 0.99, 0.98, 0.97,
>> 245 0.98, 0.97, 0.96, 0.93, 0.91, 0.9, 0.88, 0.9, 0.9, 0.9,
>> 246 0.9, 0.85, 0.9, 0.9, 0.91, 0.92, 0.9, 0.9, 0.9, 0.9,
>> 247 0.9, 0.88, 0.9, 0.88, 0.88, 0.9, 0.9, 0.88, 0.9, 0.9,
>> 248 0.9, 0.9, 0.96, 1.2, 0.9, 0.88, 0.88, 0.85, 0.9, 0.9,
>> 249 0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1, 1.08, 1.08,
>> 250 1.08, 1.08, 1.09, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15,
>> 251 1.17, 1.2, 1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
>> 252 1.16, 1.16} ;
>> 253
414 // get elements in the actual material, 254 // get elements in the actual material,
415 const G4ElementVector* theElementVector = fM << 255 const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
416 const G4double* theAtomicNumDensityVector = << 256 const G4double* theAtomicNumDensityVector =
417 const auto NumberOfElements = (G4int)fMateri << 257 fMaterial->GetAtomicNumDensityVector() ;
>> 258 const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
418 259
419 // loop for the elements in the material 260 // loop for the elements in the material
420 // to find out average values Z, vF, lF 261 // to find out average values Z, vF, lF
421 G4double z(0.0), vF(0.0), lF(0.0), a23(0.0); << 262 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
422 263
423 G4Pow* g4pow = G4Pow::GetInstance(); << 264 if( 1 == NumberOfElements ) {
424 if (1 == NumberOfElements) { << 265 z = fMaterial->GetZ() ;
425 const G4Element* element = (*theElementVec << 266 G4int iz = G4int(z) - 1 ;
426 z = element->GetZ(); << 267 if(iz < 0) iz = 0 ;
427 vF = element->GetIonisation()->GetFermiVel << 268 else if(iz > 91) iz = 91 ;
428 lF = element->GetIonisation()->GetLFactor( << 269 vF = vFermi[iz] ;
429 a23 = 1.0 / g4pow->A23(element->GetN()); << 270 lF = lFactor[iz] ;
430 } << 271
431 else { << 272 } else {
432 G4double norm(0.0); << 273 for (G4int iel=0; iel<NumberOfElements; iel++)
433 for (G4int iel = 0; iel < NumberOfElements << 274 {
434 const G4Element* element = (*theElementV << 275 const G4Element* element = (*theElementVector)[iel] ;
435 const G4double weight = theAtomicNumDens << 276 G4double z2 = element->GetZ() ;
436 norm += weight; << 277 const G4double weight = theAtomicNumDensityVector[iel] ;
437 z += element->GetZ() * weight; << 278 norm += weight ;
438 vF += element->GetIonisation()->GetFermi << 279 z += z2 * weight ;
439 lF += element->GetIonisation()->GetLFact << 280 G4int iz = G4int(z2) - 1 ;
440 a23 += weight / g4pow->A23(element->GetN << 281 if(iz < 0) iz = 0 ;
441 } << 282 else if(iz > 91) iz =91 ;
442 if (norm > 0.0) { norm = 1.0/norm; } << 283 vF += vFermi[iz] * weight ;
443 z *= norm; << 284 lF += lFactor[iz] * weight ;
444 vF *= norm; << 285 }
445 lF *= norm; << 286 z /= norm ;
446 a23 *= norm; << 287 vF /= norm ;
447 } << 288 lF /= norm ;
448 fZeff = z; << 289 }
449 fLfactor = lF; << 290 fZeff = z;
450 fFermiEnergy = 25. * CLHEP::keV * vF * vF; << 291 fLfactor = lF;
451 fInvA23 = a23; << 292 fFermiEnergy = 25.*keV*vF*vF;
452 } 293 }
453 294
454 //....oooOO0OOooo........oooOO0OOooo........oo 295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
455 296
456 void G4IonisParamMat::SetMeanExcitationEnergy( 297 void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
457 { 298 {
458 if (value == fMeanExcitationEnergy || value << 299 if(value == fMeanExcitationEnergy || value <= 0.0) return;
459 return; << 300
460 } << 301 if (G4NistManager::Instance()->GetVerbose() > 0)
461 if (G4NistManager::Instance()->GetVerbose() << 302 G4cout << "G4Material: Mean excitation energy is changed for "
462 G4cout << "G4Material: Mean excitation ene << 303 << fMaterial->GetName()
463 << " Iold= " << fMeanExcitationEner << 304 << " Iold= " << fMeanExcitationEnergy/eV
>> 305 << "eV; Inew= " << value/eV << " eV;"
464 << G4endl; 306 << G4endl;
465 } <<
466 307
467 fMeanExcitationEnergy = value; 308 fMeanExcitationEnergy = value;
>> 309 fLogMeanExcEnergy = std::log(value);
>> 310 ComputeDensityEffect();
>> 311 ComputeFluctModel();
>> 312 }
>> 313
>> 314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
468 315
469 // add corrections to density effect << 316 G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
470 G4double newlog = G4Log(value); << 317 {
471 G4double corr = 2 * (newlog - fLogMeanExcEne <<
472 fCdensity += corr; <<
473 fX0density += corr / twoln10; <<
474 fX1density += corr / twoln10; <<
475 318
476 // recompute parameters of fluctuation model << 319 // The data on mean excitation energy for compaunds
477 fLogMeanExcEnergy = newlog; << 320 // from "Stopping Powers for Electrons and Positrons"
478 ComputeFluctModel(); << 321 // ICRU Report N#37, 1984 (energy in eV)
>> 322
>> 323 const size_t numberOfMolecula = 79 ;
>> 324
>> 325 static G4String name[numberOfMolecula] = {
>> 326
>> 327 // gas
>> 328 "NH_3", "C_4H_10", "CO_2", "C_2H_6", "C_7H_16",
>> 329 "C_6H_14", "CH_4", "NO", "N_2O", "C_8H_18",
>> 330 "C_5H_12", "C_3H_8", "H_2O-Gas",
>> 331
>> 332 // liquid
>> 333 "C_3H_6O", "C_6H_5NH_2", "C_6H_6", "C_4H_9OH", "CCl_4",
>> 334 "C_6H_5Cl", "CHCl_3", "C_6H_12", "C_6H_4Cl_2", "C_4Cl_2H_8O",
>> 335 "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", "C_3H_5(OH)_3","C_7H_16",
>> 336 "C_6H_14", "CH_3OH", "C_6H_5NO_2","C_5H_12", "C_3H_7OH",
>> 337 "C_5H_5N", "C_8H_8", "C_2Cl_4", "C_7H_8", "C_2Cl_3H",
>> 338 "H_2O", "C_8H_10",
>> 339
>> 340 //solid
>> 341 "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO)-nylon", "C_25H_52",
>> 342 "(C_2H_4)-Polyethylene", "(C_5H_8O-2)-Polymethil_Methacrylate",
>> 343 "(C_8H_8)-Polystyrene", "A-150-tissue", "Al_2O_3", "CaF_2",
>> 344 "LiF", "Photo_Emulsion", "(C_2F_4)-Teflon", "SiO_2"
>> 345
>> 346 } ;
>> 347
>> 348 static G4double meanExcitation[numberOfMolecula] = {
>> 349
>> 350 53.7, 48.3, 85.0, 45.4, 49.2,
>> 351 49.1, 41.7, 87.8, 84.9, 49.5,
>> 352 48.2, 47.1, 71.6,
>> 353
>> 354 64.2, 66.2, 63.4, 59.9, 166.3,
>> 355 89.1, 156.0, 56.4, 106.5, 103.3,
>> 356 111.9, 60.0, 62.9, 72.6, 54.4,
>> 357 54.0, 67.6, 75.8, 53.6, 61.1,
>> 358 66.2, 64.0, 159.2, 62.5, 148.1,
>> 359 75.0, 61.8,
>> 360
>> 361 71.4, 75.0, 63.9, 48.3, 57.4,
>> 362 74.0, 68.7, 65.1, 145.2, 166.,
>> 363 94.0, 331.0, 99.1, 139.2
>> 364
>> 365 } ;
>> 366
>> 367 G4double x = fMeanExcitationEnergy;
>> 368
>> 369 for(size_t i=0; i<numberOfMolecula; i++) {
>> 370 if(chFormula == name[i]) {
>> 371 x = meanExcitation[i]*eV;
>> 372 break;
>> 373 }
>> 374 }
>> 375 return x;
479 } 376 }
480 377
481 //....oooOO0OOooo........oooOO0OOooo........oo 378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
482 379
483 void G4IonisParamMat::SetDensityEffectParamete << 380 G4IonisParamMat::~G4IonisParamMat()
484 G4double cd, G4double md, G4double ad, G4dou <<
485 { 381 {
486 // no check on consistence of user parameter << 382 if (fShellCorrectionVector) delete [] fShellCorrectionVector;
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 } 383 }
496 384
497 //....oooOO0OOooo........oooOO0OOooo........oo 385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
498 386
499 void G4IonisParamMat::SetDensityEffectParamete << 387 G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
500 { << 388 {
501 G4AutoLock l(&ionisMutex); << 389 *this = 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 } 390 }
538 391
539 //....oooOO0OOooo........oooOO0OOooo........oo 392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
540 393
541 G4double G4IonisParamMat::FindMeanExcitationEn << 394 const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
542 { 395 {
543 G4double res = 0.0; << 396 if (this != &right)
544 // data from density effect data << 397 {
545 if (fDensityData != nullptr) { << 398 fMaterial = right.fMaterial;
546 G4int idx = fDensityData->GetIndex(mat->Ge << 399 fMeanExcitationEnergy = right.fMeanExcitationEnergy;
547 if (idx >= 0) { << 400 fLogMeanExcEnergy = right.fLogMeanExcEnergy;
548 res = fDensityData->GetMeanIonisationPot << 401 if (fShellCorrectionVector) delete [] fShellCorrectionVector;
549 } << 402 fShellCorrectionVector = new G4double[3];
550 } << 403 fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
>> 404 fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
>> 405 fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
>> 406 fTaul = right.fTaul;
>> 407 fCdensity = right.fCdensity;
>> 408 fMdensity = right.fMdensity;
>> 409 fAdensity = right.fAdensity;
>> 410 fX0density = right.fX0density;
>> 411 fX1density = right.fX1density;
>> 412 fF1fluct = right.fF1fluct;
>> 413 fF2fluct = right.fF2fluct;
>> 414 fEnergy1fluct = right.fEnergy1fluct;
>> 415 fLogEnergy1fluct = right.fLogEnergy1fluct;
>> 416 fEnergy2fluct = right.fEnergy2fluct;
>> 417 fLogEnergy2fluct = right.fLogEnergy2fluct;
>> 418 fEnergy0fluct = right.fEnergy0fluct;
>> 419 fRateionexcfluct = right.fRateionexcfluct;
>> 420 }
>> 421 return *this;
>> 422 }
551 423
552 // The data on mean excitation energy for co << 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
553 // from "Stopping Powers for Electrons and P << 425
554 // ICRU Report N#37, 1984 (energy in eV) << 426 G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
555 // this value overwrites Density effect data << 427 {
556 G4String chFormula = mat->GetChemicalFormula << 428 return (this == (G4IonisParamMat*) &right);
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 } 429 }
>> 430
>> 431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
>> 432
>> 433 G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
>> 434 {
>> 435 return (this != (G4IonisParamMat*) &right);
>> 436 }
>> 437
>> 438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
>> 439
623 440