Geant4  10.02.p02
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 95428 2016-02-10 15:00:35Z 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 
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
55 
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;
68 
69  // compute parameters
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 
82  : fMaterial(0), fShellCorrectionVector(0)
83 {
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 
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
115 
117 {
119  if (fDensityData) { delete fDensityData; }
120  fDensityData = 0;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
125 
127 {
128  // compute mean excitation energy and shell correction vector
129  fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
130 
132  fLogMeanExcEnergy = 0.;
133 
134  size_t nElements = fMaterial->GetNumberOfElements();
135  const G4ElementVector* elmVector = fMaterial->GetElementVector();
136  const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
137 
139 
140  if(ch != "") { fMeanExcitationEnergy = FindMeanExcitationEnergy(ch); }
141 
142  // Chemical formula defines mean excitation energy
143  if(fMeanExcitationEnergy > 0.0) {
145 
146  // Compute average
147  } else {
148  for (size_t i=0; i < nElements; i++) {
149  const G4Element* elm = (*elmVector)[i];
150  fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
152  }
155  }
156 
158 
159  for (G4int j=0; j<=2; j++)
160  {
161  fShellCorrectionVector[j] = 0.;
162 
163  for (size_t k=0; k<nElements; k++) {
164  fShellCorrectionVector[j] += nAtomsPerVolume[k]
165  *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
166  }
168  }
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
172 
174 {
175  return fDensityData;
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
179 
181 {
183 
184  // Check if density effect data exist in the table
185  // R.M. Sternheimer, Atomic Data and Nuclear Data Tables, 30: 261 (1984)
188  G4int Z0 = G4lrint((*(fMaterial->GetElementVector()))[0]->GetZ());
189 
190  // for simple non-NIST materials
191  G4double corr = 0.0;
192  if(idx < 0 && 1 == nelm) {
194 
195  // Correction for base material or for non-nominal density
196  // Except cases of very different density defined in user code
197  if(idx >= 0) {
198  const G4Material* bmat = fMaterial->GetBaseMaterial();
199  if(bmat) {
200  corr = G4Log(bmat->GetDensity()/fMaterial->GetDensity());
201  } else {
203  if(dens <= 0.0) { idx = -1; }
204  else {
205  corr = G4Log(dens/fMaterial->GetDensity());
206  }
207  }
208  // 1.0 is an arbitrary empirical limit
209  // parameterisation with very different density is not applicable
210  if(std::abs(corr) > 1.0) { idx = -1; }
211  }
212  }
213 
214  //G4cout<<"DensityEffect for "<<fMaterial->GetName()<<" "<< idx << G4endl;
215 
216  if(idx >= 0) {
217 
218  // Take parameters for the density effect correction from
219  // R.M. Sternheimer et al. Density Effect For The Ionization Loss
220  // of Charged Particles in Various Substances.
221  // Atom. Data Nucl. Data Tabl. 30 (1984) 261-271.
222 
224  //G4double Cdensity = fCdensity;
232 
233  // parameter C is computed and not taken from Sternheimer tables
234  //fCdensity = 1. + 2*G4Log(fMeanExcitationEnergy/fPlasmaEnergy);
235  //G4cout << "IonisParamMat: " << fMaterial->GetName()
236  // << " Cst= " << Cdensity << " C= " << fCdensity << G4endl;
237 
238  // correction on nominal density
239  fCdensity += corr;
240  fX0density += corr/twoln10;
241  fX1density += corr/twoln10;
242 
243  } else {
244 
245  static const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
247 
248  // Compute parameters for the density effect correction in DE/Dx formula.
249  // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
250  G4int icase;
251 
253  //
254  // condensed materials
255  //
256  if ((State == kStateSolid)||(State == kStateLiquid)) {
257 
258  static const G4double E100eV = 100.*eV;
259  static const G4double ClimiS[] = {3.681 , 5.215 };
260  static const G4double X0valS[] = {1.0 , 1.5 };
261  static const G4double X1valS[] = {2.0 , 3.0 };
262 
263  if(fMeanExcitationEnergy < E100eV) { icase = 0; }
264  else { icase = 1; }
265 
266  if(fCdensity < ClimiS[icase]) { fX0density = 0.2; }
267  else { fX0density = 0.326*fCdensity - X0valS[icase]; }
268 
269  fX1density = X1valS[icase]; fMdensity = 3.0;
270 
271  //special: Hydrogen
272  if (1 == nelm && 1 == Z0) {
273  fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
274  }
275  }
276  //
277  // gases
278  //
279  if (State == kStateGas) {
280 
281  fMdensity = 3.;
282  fX1density = 4.0;
283  //static const G4double ClimiG[] = {10.,10.5,11.,11.5,12.25,13.804};
284  //static const G4double X0valG[] = {1.6,1.7,1.8,1.9,2.0,2.0};
285  //static const G4double X1valG[] = {4.0,4.0,4.0,4.0,4.0,5.0};
286 
287  if(fCdensity < 10.) {
288  fX0density = 1.6;
289  } else if(fCdensity < 11.5) {
290  fX0density = 1.6 + 0.2*(fCdensity - 10.);
291  } else if(fCdensity < 12.25) {
292  fX0density = 1.9 + (fCdensity - 11.5)/7.5;
293  } else if(fCdensity < 13.804) {
294  fX0density = 2.0;
295  fX1density = 4.0 + (fCdensity - 12.25)/1.554;
296  } else {
297  fX0density = 0.326*fCdensity-2.5; fX1density = 5.0;
298  }
299 
300  //special: Hydrogen
301  if (1 == nelm && 1 == Z0) {
302  fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
303  }
304 
305  //special: Helium
306  if (1 == nelm && 2 == Z0) {
307  fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
308  }
309  }
310  }
311 
312  // change parameters if the gas is not in STP.
313  // For the correction the density(STP) is needed.
314  // Density(STP) is calculated here :
315 
316 
317  if (State == kStateGas) {
318  G4double Density = fMaterial->GetDensity();
319  G4double Pressure = fMaterial->GetPressure();
321 
322  G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*NTP_Temperature);
323 
324  G4double ParCorr = G4Log(Density/DensitySTP);
325 
326  fCdensity -= ParCorr;
327  fX0density -= ParCorr/twoln10;
328  fX1density -= ParCorr/twoln10;
329  }
330 
331  // fAdensity parameter can be fixed for not conductive materials
332  if(0.0 == fD0density) {
335  /std::pow((fX1density-fX0density),fMdensity);
336  }
337  /*
338  G4cout << "G4IonisParamMat: density effect data for <"
339  << fMaterial->GetName()
340  << "> " << G4endl;
341  G4cout << "Eplasma(eV)= " << fPlasmaEnergy/eV
342  << " rho= " << fAdjustmentFactor
343  << " -C= " << fCdensity
344  << " x0= " << fX0density
345  << " x1= " << fX1density
346  << " a= " << fAdensity
347  << " m= " << fMdensity
348  << G4endl;
349  */
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
353 
355 {
356  // compute parameters for the energy loss fluctuation model
357  // needs an 'effective Z'
358  G4double Zeff = 0.;
359  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
360  Zeff += (fMaterial->GetFractionVector())[i]
361  *((*(fMaterial->GetElementVector()))[i]->GetZ());
362  }
363  if (Zeff > 2.) { fF2fluct = 2./Zeff; }
364  else { fF2fluct = 0.; }
365 
366  fF1fluct = 1. - fF2fluct;
367  fEnergy2fluct = 10.*Zeff*Zeff*eV;
370  /fF1fluct;
372  fEnergy0fluct = 10.*eV;
373  fRateionexcfluct = 0.4;
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
377 
379 {
380  // get elements in the actual material,
381  const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
382  const G4double* theAtomicNumDensityVector =
384  const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
385 
386  // loop for the elements in the material
387  // to find out average values Z, vF, lF
388  G4double z(0.0), vF(0.0), lF(0.0), norm(0.0), a23(0.0);
389 
390  G4Pow* g4pow = G4Pow::GetInstance();
391  if( 1 == NumberOfElements ) {
392  const G4Element* element = (*theElementVector)[0];
393  z = element->GetZ();
394  vF= element->GetIonisation()->GetFermiVelocity();
395  lF= element->GetIonisation()->GetLFactor();
396  a23 = 1.0/g4pow->A23(element->GetN());
397 
398  } else {
399  for (G4int iel=0; iel<NumberOfElements; iel++)
400  {
401  const G4Element* element = (*theElementVector)[iel];
402  const G4double weight = theAtomicNumDensityVector[iel];
403  norm += weight ;
404  z += element->GetZ() * weight;
405  vF += element->GetIonisation()->GetFermiVelocity() * weight;
406  lF += element->GetIonisation()->GetLFactor() * weight;
407  a23 += weight/g4pow->A23(element->GetN());
408  }
409  z /= norm;
410  vF /= norm;
411  lF /= norm;
412  a23 /= norm;
413  }
414  fZeff = z;
415  fLfactor = lF;
416  fFermiEnergy = 25.*keV*vF*vF;
417  fInvA23 = a23;
418 }
419 
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
421 
423 {
424  if(value == fMeanExcitationEnergy || value <= 0.0) { return; }
425  if (G4NistManager::Instance()->GetVerbose() > 1) {
426  G4cout << "G4Material: Mean excitation energy is changed for "
427  << fMaterial->GetName()
428  << " Iold= " << fMeanExcitationEnergy/eV
429  << "eV; Inew= " << value/eV << " eV;"
430  << G4endl;
431  }
432 
433  fMeanExcitationEnergy = value;
434 
435  // add corrections to density effect
436  G4double newlog = G4Log(value);
437  G4double corr = 2*(newlog - fLogMeanExcEnergy);
438  fCdensity += corr;
439  fX0density += corr/twoln10;
440  fX1density += corr/twoln10;
441 
442  // recompute parameters of fluctuation model
443  fLogMeanExcEnergy = newlog;
445 }
446 
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
448 
450 {
451  // The data on mean excitation energy for compaunds
452  // from "Stopping Powers for Electrons and Positrons"
453  // ICRU Report N#37, 1984 (energy in eV)
454 
455  size_t numberOfMolecula = 54;
456  static const G4String name[54] = {
457  // gas 0 - 12
458  "NH_3", "C_4H_10", "CO_2", "C_2H_6", "C_7H_16-Gas",
459  // "G4_AMMONIA", "G4_BUTANE","G4_CARBON_DIOXIDE","G4_ETHANE", "G4_N-HEPTANE"
460  "C_6H_14-Gas", "CH_4", "NO", "N_2O", "C_8H_18-Gas",
461  // "G4_N-HEXANE" , "G4_METHANE", "x", "G4_NITROUS_OXIDE", "G4_OCTANE"
462  "C_5H_12-Gas", "C_3H_8", "H_2O-Gas",
463  // "G4_N-PENTANE", "G4_PROPANE", "G4_WATER_VAPOR"
464 
465  // liquid 13 - 39
466  "C_3H_6O", "C_6H_5NH_2", "C_6H_6", "C_4H_9OH", "CCl_4",
467  //"G4_ACETONE","G4_ANILINE","G4_BENZENE","G4_N-BUTYL_ALCOHOL","G4_CARBON_TETRACHLORIDE"
468  "C_6H_5Cl", "CHCl_3", "C_6H_12", "C_6H_4Cl_2", "C_4Cl_2H_8O",
469  //"G4_CHLOROBENZENE","G4_CHLOROFORM","G4_CYCLOHEXANE","G4_1,2-DICHLOROBENZENE","G4_DICHLORODIETHYL_ETHER"
470  "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH", "C_3H_5(OH)_3","C_7H_16",
471  //"G4_1,2-DICHLOROETHANE","G4_DIETHYL_ETHER","G4_ETHYL_ALCOHOL","G4_GLYCEROL","G4_N-HEPTANE"
472  "C_6H_14", "CH_3OH", "C_6H_5NO_2","C_5H_12", "C_3H_7OH",
473  //"G4_N-HEXANE","G4_METHANOL","G4_NITROBENZENE","G4_N-PENTANE","G4_N-PROPYL_ALCOHOL",
474  "C_5H_5N", "C_8H_8", "C_2Cl_4", "C_7H_8", "C_2Cl_3H",
475  //"G4_PYRIDINE","G4_POLYSTYRENE","G4_TETRACHLOROETHYLENE","G4_TOLUENE","G4_TRICHLOROETHYLENE"
476  "H_2O", "C_8H_10",
477  // "G4_WATER", "G4_XYLENE"
478 
479  // solid 40 - 53
480  "C_5H_5N_5", "C_5H_5N_5O", "(C_6H_11NO)-nylon", "C_25H_52",
481  // "G4_ADENINE", "G4_GUANINE", "G4_NYLON-6-6", "G4_PARAFFIN"
482  "(C_2H_4)-Polyethylene", "(C_5H_8O_2)-Polymethil_Methacrylate",
483  // "G4_ETHYLENE", "G4_PLEXIGLASS"
484  "(C_8H_8)-Polystyrene", "A-150-tissue", "Al_2O_3", "CaF_2",
485  // "G4_POLYSTYRENE", "G4_A-150_TISSUE", "G4_ALUMINUM_OXIDE", "G4_CALCIUM_FLUORIDE"
486  "LiF", "Photo_Emulsion", "(C_2F_4)-Teflon", "SiO_2"
487  // "G4_LITHIUM_FLUORIDE", "G4_PHOTO_EMULSION", "G4_TEFLON", "G4_SILICON_DIOXIDE"
488  } ;
489 
490  static const G4double meanExcitation[54] = {
491 
492  53.7, 48.3, 85.0, 45.4, 49.2,
493  49.1, 41.7, 87.8, 84.9, 49.5,
494  48.2, 47.1, 71.6,
495 
496  64.2, 66.2, 63.4, 59.9, 166.3,
497  89.1, 156.0, 56.4, 106.5, 103.3,
498  111.9, 60.0, 62.9, 72.6, 54.4,
499  54.0, 67.6, 75.8, 53.6, 61.1,
500  66.2, 64.0, 159.2, 62.5, 148.1,
501  75.0, 61.8,
502 
503  71.4, 75.0, 63.9, 48.3, 57.4,
504  74.0, 68.7, 65.1, 145.2, 166.,
505  94.0, 331.0, 99.1, 139.2
506  };
507 
509 
510  for(size_t i=0; i<numberOfMolecula; i++) {
511  if(chFormula == name[i]) {
512  x = meanExcitation[i]*eV;
513  break;
514  }
515  }
516  return x;
517 }
518 
519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
520 
G4double GetPressure() const
Definition: G4Material.hh:183
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double fLogEnergy2fluct
G4double GetTotNbOfElectPerVolume() const
Definition: G4Material.hh:212
std::vector< G4Element * > G4ElementVector
Definition: G4Pow.hh:56
G4double GetN() const
Definition: G4Element.hh:134
G4int GetElementIndex(G4int Z, G4State mState)
G4State
Definition: G4Material.hh:114
G4double z
Definition: TRTMaterials.hh:39
G4double FindMeanExcitationEnergy(const G4String &chFormula)
static const G4int numberOfMolecula
const G4String & GetChemicalFormula() const
Definition: G4Material.hh:179
void SetMeanExcitationEnergy(G4double value)
G4String name
Definition: TRTMaterials.hh:40
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetZ() const
Definition: G4Element.hh:131
#define State(X)
G4double fAdjustmentFactor
G4double GetDensity() const
Definition: G4Material.hh:180
G4double GetLFactor() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4double logZ(G4int Z) const
Definition: G4Pow.hh:166
G4double GetFermiVelocity() const
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:206
G4double GetMeanExcitationEnergy() const
G4double A23(G4double A) const
Definition: G4Pow.hh:160
G4double fRateionexcfluct
G4GLOB_DLL std::ostream G4cout
G4double GetX1density(G4int idx)
G4int GetIndex(const G4String &matName)
virtual ~G4IonisParamMat()
G4double GetPlasmaEnergy(G4int idx)
G4double fLogMeanExcEnergy
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:216
G4double fMeanEnergyPerIon
G4double fLogEnergy1fluct
static G4DensityEffectData * GetDensityEffectData()
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetAdensity(G4int idx)
G4double fMeanExcitationEnergy
static const double eV
Definition: G4SIunits.hh:212
G4double * fShellCorrectionVector
static const double pi
Definition: G4SIunits.hh:74
int G4lrint(double ad)
Definition: templates.hh:163
G4Material * fMaterial
G4IonisParamMat(G4Material *)
const G4double x[NPOINTSGL]
G4IonisParamElm * GetIonisation() const
Definition: G4Element.hh:198
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4double GetCdensity(G4int idx)
static const G4double NTP_Temperature
Definition: G4Material.hh:116
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
G4double GetMdensity(G4int idx)
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double keV
Definition: G4SIunits.hh:213
G4State GetState() const
Definition: G4Material.hh:181
G4double GetNominalDensity(G4int Z) const
double G4double
Definition: G4Types.hh:76
G4double GetDelta0density(G4int idx)
const G4double * GetFractionVector() const
Definition: G4Material.hh:194
G4double GetAdjustmentFactor(G4int idx)
G4double GetX0density(G4int idx)
static G4DensityEffectData * fDensityData
void ComputeMeanParameters()