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