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