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

Last change on this file since 1201 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 17.2 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.33 2009/11/19 16:43:36 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
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// 30-10-09, add G4DensityEffectData class and density effect computation (VI)
43
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
45
46#include "G4IonisParamMat.hh"
47#include "G4Material.hh"
48#include "G4DensityEffectData.hh"
49
50G4DensityEffectData* G4IonisParamMat::fDensityData = 0;
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
53
54G4IonisParamMat::G4IonisParamMat(G4Material* material)
55  : fMaterial(material)
56{
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
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
94  size_t nElements = fMaterial->GetNumberOfElements();
95  const G4ElementVector* elmVector = fMaterial->GetElementVector();
96  const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
97
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);
115  }
116
117  fShellCorrectionVector = new G4double[3]; //[3]
118
119  for (G4int j=0; j<=2; j++)
120  {
121    fShellCorrectionVector[j] = 0.;
122
123    for (size_t k=0; k<nElements; k++) {
124      fShellCorrectionVector[j] += nAtomsPerVolume[k]
125        *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
126    }
127    fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
128  } 
129}
130
131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
132
133G4DensityEffectData* G4IonisParamMat::GetDensityEffectData()
134{
135  return fDensityData;
136}
137
138//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
139                   
140void G4IonisParamMat::ComputeDensityEffect()
141{
142  static const G4double twoln10 = 2.*std::log(10.);
143  G4State State = fMaterial->GetState();
144
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  }
151
152  if(idx >= 0) {
153
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    //
185 
186    if ((State == kStateSolid)||(State == kStateLiquid)) {
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      }
205    }
206
207    //
208    // gases
209    //
210    if (State == kStateGas) { 
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      }
230    }
231  }
232
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();
242     
243    G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
244
245    G4double ParCorr = std::log(Density/DensitySTP);
246 
247    fCdensity  -= ParCorr;
248    fX0density -= ParCorr/twoln10;
249    fX1density -= ParCorr/twoln10;
250  }
251
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  */
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.;
280  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
281     Zeff += (fMaterial->GetFractionVector())[i]
282             *((*(fMaterial->GetElementVector()))[i]->GetZ());
283  }
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 ) {
312    const G4Element* element = (*theElementVector)[0];
313    z = element->GetZ();
314    vF= element->GetIonisation()->GetFermiVelocity();
315    lF= element->GetIonisation()->GetLFactor();
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 ;
323        z    += element->GetZ() * weight ;
324        vF   += element->GetIonisation()->GetFermiVelocity() * weight ;
325        lF   += element->GetIonisation()->GetLFactor() * weight ;
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
342  /*
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;
349  */
350 
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{
425  if (fShellCorrectionVector) { delete [] fShellCorrectionVector; }
426  if (fDensityData) { delete fDensityData; }
427  fDensityData = 0;
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;
457      fD0density                = right.fD0density;
458      fPlasmaEnergy             = right.fPlasmaEnergy;
459      fAdjustmentFactor         = right.fAdjustmentFactor;
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;
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    } 
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.