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

Last change on this file was 1337, checked in by garnier, 14 years ago

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