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

Last change on this file since 1233 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 17.3 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.34 2009/11/30 15:48:04 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03 $
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// 04-03-08, remove reference to G4NistManager. Add fBirks constant (mma)
42// 30-10-09, add G4DensityEffectData class and density effect computation (VI)
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
45
46#include "G4IonisParamMat.hh"
47#include "G4Material.hh"
48#include "G4DensityEffectData.hh"
49
50G4DensityEffectData* G4IonisParamMat::fDensityData = 0;
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
53
54G4IonisParamMat::G4IonisParamMat(G4Material* material)
55 : fMaterial(material)
56{
57 fBirks = 0.;
58 fMeanEnergyPerIon = 0.0;
59
60 // minimal set of default parameters for density effect
61 fCdensity = 0.0;
62 fD0density = 0.0;
63 fAdjustmentFactor = 1.0;
64 if(!fDensityData) fDensityData = new G4DensityEffectData();
65
66 // compute parameters
67 ComputeMeanParameters();
68 ComputeDensityEffect();
69 ComputeFluctModel();
70 ComputeIonParameters();
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
74
75// Fake default constructor - sets only member data and allocates memory
76// for usage restricted to object persistency
77
78G4IonisParamMat::G4IonisParamMat(__void__&)
79 : fMaterial(0), fShellCorrectionVector(0)
80{
81}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
84
85void G4IonisParamMat::ComputeMeanParameters()
86{
87 // compute mean excitation energy and shell correction vector
88
89 fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
90
91 fMeanExcitationEnergy = 0.;
92 fLogMeanExcEnergy = 0.;
93
94 size_t nElements = fMaterial->GetNumberOfElements();
95 const G4ElementVector* elmVector = fMaterial->GetElementVector();
96 const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
97
98 const G4String ch = fMaterial->GetChemicalFormula();
99
100 if(ch != "") fMeanExcitationEnergy = FindMeanExcitationEnergy(ch);
101
102 // Chemical formula defines mean excitation energy
103 if(fMeanExcitationEnergy > 0.0) {
104 fLogMeanExcEnergy = std::log(fMeanExcitationEnergy);
105
106 // Compute average
107 } else {
108 for (size_t i=0; i < nElements; i++) {
109 const G4Element* elm = (*elmVector)[i];
110 fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
111 *std::log(elm->GetIonisation()->GetMeanExcitationEnergy());
112 }
113 fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
114 fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
115 }
116
117 fShellCorrectionVector = new G4double[3]; //[3]
118
119 for (G4int j=0; j<=2; j++)
120 {
121 fShellCorrectionVector[j] = 0.;
122
123 for (size_t k=0; k<nElements; k++) {
124 fShellCorrectionVector[j] += nAtomsPerVolume[k]
125 *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
126 }
127 fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
128 }
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
132
133G4DensityEffectData* G4IonisParamMat::GetDensityEffectData()
134{
135 return fDensityData;
136}
137
138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
139
140void G4IonisParamMat::ComputeDensityEffect()
141{
142 static const G4double twoln10 = 2.*std::log(10.);
143 G4State State = fMaterial->GetState();
144
145 // Check if density effect data exist in the table
146 // R.M. Sternheimer, Atomic Data and Nuclear Data Tables, 30: 261 (1984)
147 G4int idx = fDensityData->GetIndex(fMaterial->GetName());
148 // if(idx < 0 && fMaterial->GetNumberOfElements() == 1) {
149 // idx = fDensityData->GetIndex(G4int(fMaterial->GetZ()));
150 // }
151
152 //G4cout << "DensityEffect for " << fMaterial->GetName() << " " << idx << G4endl;
153
154 if(idx >= 0) {
155
156 // Take parameters for the density effect correction from
157 // R.M. Sternheimer et al. Density Effect For The Ionization Loss
158 // of Charged Particles in Various Substances.
159 // Atom. Data Nucl. Data Tabl. 30 (1984) 261-271.
160
161 fCdensity = fDensityData->GetCdensity(idx);
162 fMdensity = fDensityData->GetMdensity(idx);
163 fAdensity = fDensityData->GetAdensity(idx);
164 fX0density = fDensityData->GetX0density(idx);
165 fX1density = fDensityData->GetX1density(idx);
166 fD0density = fDensityData->GetDelta0density(idx);
167 fPlasmaEnergy = fDensityData->GetPlasmaEnergy(idx);
168 fAdjustmentFactor = fDensityData->GetAdjustmentFactor(idx);
169
170 } else {
171
172 const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
173 fPlasmaEnergy = std::sqrt(Cd2*fMaterial->GetTotNbOfElectPerVolume());
174
175 // Compute parameters for the density effect correction in DE/Dx formula.
176 // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
177 G4int icase;
178
179 fCdensity = 1. + 2*std::log(fMeanExcitationEnergy/fPlasmaEnergy);
180
181 //fCdensity = 1. + std::log(fMeanExcitationEnergy*fMeanExcitationEnergy
182 // /(Cd2*fMaterial->GetTotNbOfElectPerVolume()));
183
184 //
185 // condensed materials
186 //
187
188 if ((State == kStateSolid)||(State == kStateLiquid)) {
189
190 const G4double E100eV = 100.*eV;
191 const G4double ClimiS[] = {3.681 , 5.215 };
192 const G4double X0valS[] = {1.0 , 1.5 };
193 const G4double X1valS[] = {2.0 , 3.0 };
194
195 if(fMeanExcitationEnergy < E100eV) icase = 0;
196 else icase = 1;
197
198 if(fCdensity < ClimiS[icase]) fX0density = 0.2;
199 else fX0density = 0.326*fCdensity-X0valS[icase];
200
201 fX1density = X1valS[icase] ; fMdensity = 3.0;
202
203 //special: Hydrogen
204 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
205 fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
206 }
207 }
208
209 //
210 // gases
211 //
212 if (State == kStateGas) {
213
214 const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804};
215 const G4double X0valG[] = { 1.6 , 1.7 , 1.8 , 1.9 , 2.0 , 2.0 };
216 const G4double X1valG[] = { 4.0 , 4.0 , 4.0 , 4.0 , 4.0 , 5.0 };
217
218 icase = 5;
219 fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ;
220 while((icase > 0)&&(fCdensity < ClimiG[icase])) icase-- ;
221 fX0density = X0valG[icase] ; fX1density = X1valG[icase] ;
222
223 //special: Hydrogen
224 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
225 fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
226 }
227
228 //special: Helium
229 if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==2.)) {
230 fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
231 }
232 }
233 }
234
235 // change parameters if the gas is not in STP.
236 // For the correction the density(STP) is needed.
237 // Density(STP) is calculated here :
238
239
240 if (State == kStateGas) {
241 G4double Density = fMaterial->GetDensity();
242 G4double Pressure = fMaterial->GetPressure();
243 G4double Temp = fMaterial->GetTemperature();
244
245 G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
246
247 G4double ParCorr = std::log(Density/DensitySTP);
248
249 fCdensity -= ParCorr;
250 fX0density -= ParCorr/twoln10;
251 fX1density -= ParCorr/twoln10;
252 }
253
254 // fAdensity parameter can be fixed for not conductive materials
255 if(0.0 == fD0density) {
256 G4double Xa = fCdensity/twoln10;
257 fAdensity = twoln10*(Xa-fX0density)
258 /std::pow((fX1density-fX0density),fMdensity);
259 }
260 /*
261 G4cout << "G4IonisParamMat: density effect data for <" << fMaterial->GetName()
262 << "> " << G4endl;
263 G4cout << "Eplasma(eV)= " << fPlasmaEnergy/eV
264 << " rho= " << fAdjustmentFactor
265 << " -C= " << fCdensity
266 << " x0= " << fX0density
267 << " x1= " << fX1density
268 << " a= " << fAdensity
269 << " m= " << fMdensity
270 << G4endl;
271 */
272}
273
274//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
275
276void G4IonisParamMat::ComputeFluctModel()
277{
278 // compute parameters for the energy loss fluctuation model
279
280 // need an 'effective Z' ?????
281 G4double Zeff = 0.;
282 for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
283 Zeff += (fMaterial->GetFractionVector())[i]
284 *((*(fMaterial->GetElementVector()))[i]->GetZ());
285 }
286 if (Zeff > 2.) fF2fluct = 2./Zeff ;
287 else fF2fluct = 0.;
288
289 fF1fluct = 1. - fF2fluct;
290 fEnergy2fluct = 10.*Zeff*Zeff*eV;
291 fLogEnergy2fluct = std::log(fEnergy2fluct);
292 fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
293 /fF1fluct;
294 fEnergy1fluct = std::exp(fLogEnergy1fluct);
295 fEnergy0fluct = 10.*eV;
296 fRateionexcfluct = 0.4;
297}
298
299//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
300
301void G4IonisParamMat::ComputeIonParameters()
302{
303 // get elements in the actual material,
304 const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
305 const G4double* theAtomicNumDensityVector =
306 fMaterial->GetAtomicNumDensityVector() ;
307 const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
308
309 // loop for the elements in the material
310 // to find out average values Z, vF, lF
311 G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
312
313 if( 1 == NumberOfElements ) {
314 const G4Element* element = (*theElementVector)[0];
315 z = element->GetZ();
316 vF= element->GetIonisation()->GetFermiVelocity();
317 lF= element->GetIonisation()->GetLFactor();
318
319 } else {
320 for (G4int iel=0; iel<NumberOfElements; iel++)
321 {
322 const G4Element* element = (*theElementVector)[iel] ;
323 const G4double weight = theAtomicNumDensityVector[iel] ;
324 norm += weight ;
325 z += element->GetZ() * weight ;
326 vF += element->GetIonisation()->GetFermiVelocity() * weight ;
327 lF += element->GetIonisation()->GetLFactor() * weight ;
328 }
329 z /= norm ;
330 vF /= norm ;
331 lF /= norm ;
332 }
333 fZeff = z;
334 fLfactor = lF;
335 fFermiEnergy = 25.*keV*vF*vF;
336}
337
338//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
339
340void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
341{
342 if(value == fMeanExcitationEnergy || value <= 0.0) return;
343
344 /*
345 if (G4NistManager::Instance()->GetVerbose() > 0)
346 G4cout << "G4Material: Mean excitation energy is changed for "
347 << fMaterial->GetName()
348 << " Iold= " << fMeanExcitationEnergy/eV
349 << "eV; Inew= " << value/eV << " eV;"
350 << G4endl;
351 */
352
353 fMeanExcitationEnergy = value;
354 fLogMeanExcEnergy = std::log(value);
355 ComputeDensityEffect();
356 ComputeFluctModel();
357}
358
359//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
360
361G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
362{
363
364 // The data on mean excitation energy for compaunds
365 // from "Stopping Powers for Electrons and Positrons"
366 // ICRU Report N#37, 1984 (energy in eV)
367
368 const size_t numberOfMolecula = 79 ;
369
370 static G4String name[numberOfMolecula] = {
371
372 // gas
373 "NH_3", "C_4H_10", "CO_2", "C_2H_6", "C_7H_16",
374 "C_6H_14", "CH_4", "NO", "N_2O", "C_8H_18",
375 "C_5H_12", "C_3H_8", "H_2O-Gas",
376
377 // liquid
378 "C_3H_6O", "C_6H_5NH_2", "C_6H_6", "C_4H_9OH", "CCl_4",
379 "C_6H_5Cl", "CHCl_3", "C_6H_12", "C_6H_4Cl_2", "C_4Cl_2H_8O",
380 "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", "C_3H_5(OH)_3","C_7H_16",
381 "C_6H_14", "CH_3OH", "C_6H_5NO_2","C_5H_12", "C_3H_7OH",
382 "C_5H_5N", "C_8H_8", "C_2Cl_4", "C_7H_8", "C_2Cl_3H",
383 "H_2O", "C_8H_10",
384
385 //solid
386 "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO)-nylon", "C_25H_52",
387 "(C_2H_4)-Polyethylene", "(C_5H_8O-2)-Polymethil_Methacrylate",
388 "(C_8H_8)-Polystyrene", "A-150-tissue", "Al_2O_3", "CaF_2",
389 "LiF", "Photo_Emulsion", "(C_2F_4)-Teflon", "SiO_2"
390
391 } ;
392
393 static G4double meanExcitation[numberOfMolecula] = {
394
395 53.7, 48.3, 85.0, 45.4, 49.2,
396 49.1, 41.7, 87.8, 84.9, 49.5,
397 48.2, 47.1, 71.6,
398
399 64.2, 66.2, 63.4, 59.9, 166.3,
400 89.1, 156.0, 56.4, 106.5, 103.3,
401 111.9, 60.0, 62.9, 72.6, 54.4,
402 54.0, 67.6, 75.8, 53.6, 61.1,
403 66.2, 64.0, 159.2, 62.5, 148.1,
404 75.0, 61.8,
405
406 71.4, 75.0, 63.9, 48.3, 57.4,
407 74.0, 68.7, 65.1, 145.2, 166.,
408 94.0, 331.0, 99.1, 139.2
409
410 } ;
411
412 G4double x = fMeanExcitationEnergy;
413
414 for(size_t i=0; i<numberOfMolecula; i++) {
415 if(chFormula == name[i]) {
416 x = meanExcitation[i]*eV;
417 break;
418 }
419 }
420 return x;
421}
422
423//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
424
425G4IonisParamMat::~G4IonisParamMat()
426{
427 if (fShellCorrectionVector) { delete [] fShellCorrectionVector; }
428 if (fDensityData) { delete fDensityData; }
429 fDensityData = 0;
430}
431
432//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
433
434G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
435{
436 *this = right;
437}
438
439//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
440
441const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
442{
443 if (this != &right)
444 {
445 fMaterial = right.fMaterial;
446 fMeanExcitationEnergy = right.fMeanExcitationEnergy;
447 fLogMeanExcEnergy = right.fLogMeanExcEnergy;
448 if (fShellCorrectionVector) delete [] fShellCorrectionVector;
449 fShellCorrectionVector = new G4double[3];
450 fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
451 fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
452 fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
453 fTaul = right.fTaul;
454 fCdensity = right.fCdensity;
455 fMdensity = right.fMdensity;
456 fAdensity = right.fAdensity;
457 fX0density = right.fX0density;
458 fX1density = right.fX1density;
459 fD0density = right.fD0density;
460 fPlasmaEnergy = right.fPlasmaEnergy;
461 fAdjustmentFactor = right.fAdjustmentFactor;
462 fF1fluct = right.fF1fluct;
463 fF2fluct = right.fF2fluct;
464 fEnergy1fluct = right.fEnergy1fluct;
465 fLogEnergy1fluct = right.fLogEnergy1fluct;
466 fEnergy2fluct = right.fEnergy2fluct;
467 fLogEnergy2fluct = right.fLogEnergy2fluct;
468 fEnergy0fluct = right.fEnergy0fluct;
469 fRateionexcfluct = right.fRateionexcfluct;
470 fZeff = right.fZeff;
471 fFermiEnergy = right.fFermiEnergy;
472 fLfactor = right.fLfactor;
473 fBirks = right.fBirks;
474 fMeanEnergyPerIon = right.fMeanEnergyPerIon;
475 fDensityData = right.fDensityData;
476 }
477 return *this;
478}
479
480//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
481
482G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
483{
484 return (this == (G4IonisParamMat*) &right);
485}
486
487//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
488
489G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
490{
491 return (this != (G4IonisParamMat*) &right);
492}
493
494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
495
Note: See TracBrowser for help on using the repository browser.