source: trunk/source/materials/src/G4IonisParamMat.cc@ 822

Last change on this file since 822 was 822, checked in by garnier, 17 years ago

import all except CVS
File size: 15.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
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....
32
33// 09-07-98, data moved from G4Material, M.Maire
34// 18-07-98, bug corrected in ComputeDensityEffect() for gas
35// 16-01-01, bug corrected in ComputeDensityEffect() E100eV (L.Urban)
36// 08-02-01, fShellCorrectionVector correctly handled (mma)
37// 28-10-02, add setMeanExcitationEnergy (V.Ivanchenko)
38// 06-09-04, factor 2 to shell correction term (V.Ivanchenko)
39// 10-05-05, add a missing coma in FindMeanExcitationEnergy() - Bug#746 (mma)
40// 27-09-07, add computation of parameters for ions (V.Ivanchenko)
41
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
43
44#include "G4IonisParamMat.hh"
45#include "G4Material.hh"
46#include "G4NistManager.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
49
50G4IonisParamMat::G4IonisParamMat(G4Material* material)
51 : fMaterial(material)
52{
53 ComputeMeanParameters();
54 ComputeDensityEffect();
55 ComputeFluctModel();
56 ComputeIonParameters();
57}
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
60
61// Fake default constructor - sets only member data and allocates memory
62// for usage restricted to object persistency
63
64G4IonisParamMat::G4IonisParamMat(__void__&)
65 : fMaterial(0), fShellCorrectionVector(0)
66{
67}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
70
71void G4IonisParamMat::ComputeMeanParameters()
72{
73 // compute mean excitation energy and shell correction vector
74
75 fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
76
77 fMeanExcitationEnergy = 0.;
78 fLogMeanExcEnergy = 0.;
79
80
81 for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) {
82 fLogMeanExcEnergy +=
83 (fMaterial->GetVecNbOfAtomsPerVolume())[i]
84 *((*(fMaterial->GetElementVector()))[i]->GetZ())
85 *std::log((*(fMaterial->GetElementVector()))[i]->GetIonisation()
86 ->GetMeanExcitationEnergy());
87 }
88
89 fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
90 fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
91
92 fShellCorrectionVector = new G4double[3]; //[3]
93
94 for (G4int j=0; j<=2; j++)
95 {
96 fShellCorrectionVector[j] = 0.;
97
98 for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) {
99 fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k]
100 *((*(fMaterial->GetElementVector()))[k]->GetIonisation()
101 ->GetShellCorrectionVector()[j]);
102 }
103 fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
104 }
105}
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
108
109void G4IonisParamMat::ComputeDensityEffect()
110{
111 // Compute parameters for the density effect correction in DE/Dx formula.
112 // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
113
114 const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
115 const G4double twoln10 = 2.*std::log(10.);
116
117 G4int icase;
118
119 fCdensity = 1. + std::log(fMeanExcitationEnergy*fMeanExcitationEnergy
120 /(Cd2*fMaterial->GetTotNbOfElectPerVolume()));
121
122 //
123 // condensed materials
124 //
125 G4State State = fMaterial->GetState();
126
127 if ((State == kStateSolid)||(State == kStateLiquid)) {
128
129 const G4double E100eV = 100.*eV;
130 const G4double ClimiS[] = {3.681 , 5.215 };
131 const G4double X0valS[] = {1.0 , 1.5 };
132 const G4double X1valS[] = {2.0 , 3.0 };
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);
192}
193
194//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
195
196void G4IonisParamMat::ComputeFluctModel()
197{
198 // compute parameters for the energy loss fluctuation model
199
200 // need an 'effective Z' ?????
201 G4double Zeff = 0.;
202 for (size_t i=0;i<fMaterial->GetNumberOfElements();i++)
203 Zeff += (fMaterial->GetFractionVector())[i]
204 *((*(fMaterial->GetElementVector()))[i]->GetZ());
205
206 if (Zeff > 2.) fF2fluct = 2./Zeff ;
207 else fF2fluct = 0.;
208
209 fF1fluct = 1. - fF2fluct;
210 fEnergy2fluct = 10.*Zeff*Zeff*eV;
211 fLogEnergy2fluct = std::log(fEnergy2fluct);
212 fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
213 /fF1fluct;
214 fEnergy1fluct = std::exp(fLogEnergy1fluct);
215 fEnergy0fluct = 10.*eV;
216 fRateionexcfluct = 0.4;
217}
218
219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
220
221void G4IonisParamMat::ComputeIonParameters()
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
254 // get elements in the actual material,
255 const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
256 const G4double* theAtomicNumDensityVector =
257 fMaterial->GetAtomicNumDensityVector() ;
258 const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
259
260 // loop for the elements in the material
261 // to find out average values Z, vF, lF
262 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
263
264 if( 1 == NumberOfElements ) {
265 z = fMaterial->GetZ() ;
266 G4int iz = G4int(z) - 1 ;
267 if(iz < 0) iz = 0 ;
268 else if(iz > 91) iz = 91 ;
269 vF = vFermi[iz] ;
270 lF = lFactor[iz] ;
271
272 } else {
273 for (G4int iel=0; iel<NumberOfElements; iel++)
274 {
275 const G4Element* element = (*theElementVector)[iel] ;
276 G4double z2 = element->GetZ() ;
277 const G4double weight = theAtomicNumDensityVector[iel] ;
278 norm += weight ;
279 z += z2 * weight ;
280 G4int iz = G4int(z2) - 1 ;
281 if(iz < 0) iz = 0 ;
282 else if(iz > 91) iz =91 ;
283 vF += vFermi[iz] * weight ;
284 lF += lFactor[iz] * weight ;
285 }
286 z /= norm ;
287 vF /= norm ;
288 lF /= norm ;
289 }
290 fZeff = z;
291 fLfactor = lF;
292 fFermiEnergy = 25.*keV*vF*vF;
293}
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
296
297void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
298{
299 if(value == fMeanExcitationEnergy || value <= 0.0) return;
300
301 if (G4NistManager::Instance()->GetVerbose() > 0)
302 G4cout << "G4Material: Mean excitation energy is changed for "
303 << fMaterial->GetName()
304 << " Iold= " << fMeanExcitationEnergy/eV
305 << "eV; Inew= " << value/eV << " eV;"
306 << G4endl;
307
308 fMeanExcitationEnergy = value;
309 fLogMeanExcEnergy = std::log(value);
310 ComputeDensityEffect();
311 ComputeFluctModel();
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
315
316G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
317{
318
319 // The data on mean excitation energy for compaunds
320 // from "Stopping Powers for Electrons and Positrons"
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;
376}
377
378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
379
380G4IonisParamMat::~G4IonisParamMat()
381{
382 if (fShellCorrectionVector) delete [] fShellCorrectionVector;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
386
387G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
388{
389 *this = right;
390}
391
392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
393
394const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
395{
396 if (this != &right)
397 {
398 fMaterial = right.fMaterial;
399 fMeanExcitationEnergy = right.fMeanExcitationEnergy;
400 fLogMeanExcEnergy = right.fLogMeanExcEnergy;
401 if (fShellCorrectionVector) delete [] fShellCorrectionVector;
402 fShellCorrectionVector = new G4double[3];
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}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
425
426G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
427{
428 return (this == (G4IonisParamMat*) &right);
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
432
433G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
434{
435 return (this != (G4IonisParamMat*) &right);
436}
437
438//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
439
Note: See TracBrowser for help on using the repository browser.