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