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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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