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

Last change on this file since 1240 was 1228, checked in by garnier, 14 years ago

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