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

Last change on this file since 1317 was 1315, checked in by garnier, 15 years ago

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