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

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

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