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

Last change on this file since 1058 was 986, checked in by garnier, 15 years ago

fichiers manquants
File size: 14.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.25 2008/07/08 10:34:56 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
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
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
44
45#include "G4IonisParamMat.hh"
46#include "G4Material.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
49
50G4IonisParamMat::G4IonisParamMat(G4Material* material)
51  : fMaterial(material)
52{
53  ComputeMeanParameters();
54  ComputeDensityEffect();
55  ComputeFluctModel();
56  ComputeIonParameters();
57 
58  fBirks = 0.;
59  fMeanEnergyPerIon = 0.;
60}
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
63
64// Fake default constructor - sets only member data and allocates memory
65//                            for usage restricted to object persistency
66
67G4IonisParamMat::G4IonisParamMat(__void__&)
68  : fMaterial(0), fShellCorrectionVector(0)
69{
70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
73
74void G4IonisParamMat::ComputeMeanParameters()
75{
76  // compute mean excitation energy and shell correction vector
77
78  fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
79
80  fMeanExcitationEnergy = 0.;
81  fLogMeanExcEnergy = 0.;
82
83  size_t nElements = fMaterial->GetNumberOfElements();
84  const G4ElementVector* elmVector = fMaterial->GetElementVector();
85  const G4double* nAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
86
87  const G4String ch = fMaterial->GetChemicalFormula();
88
89  if(ch != "") fMeanExcitationEnergy = FindMeanExcitationEnergy(ch);
90
91  // Chemical formula defines mean excitation energy
92  if(fMeanExcitationEnergy > 0.0) {
93    fLogMeanExcEnergy = std::log(fMeanExcitationEnergy);
94
95    // Compute average
96  } else {
97    for (size_t i=0; i < nElements; i++) {
98      const G4Element* elm = (*elmVector)[i];
99      fLogMeanExcEnergy += nAtomsPerVolume[i]*elm->GetZ()
100        *std::log(elm->GetIonisation()->GetMeanExcitationEnergy());
101    }
102    fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
103    fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
104  }
105
106  fShellCorrectionVector = new G4double[3]; //[3]
107
108  for (G4int j=0; j<=2; j++)
109  {
110    fShellCorrectionVector[j] = 0.;
111
112    for (size_t k=0; k<nElements; k++) {
113      fShellCorrectionVector[j] += nAtomsPerVolume[k]
114        *(((*elmVector)[k])->GetIonisation()->GetShellCorrectionVector())[j];
115    }
116    fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
117  } 
118}
119
120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
121                   
122void G4IonisParamMat::ComputeDensityEffect()
123{
124  // Compute parameters for the density effect correction in DE/Dx formula.
125  // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
126
127  const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
128  const G4double twoln10 = 2.*std::log(10.);
129
130  G4int icase;
131 
132  fCdensity = 1. + std::log(fMeanExcitationEnergy*fMeanExcitationEnergy
133              /(Cd2*fMaterial->GetTotNbOfElectPerVolume()));
134
135  //
136  // condensed materials
137  //
138  G4State State = fMaterial->GetState();
139 
140  if ((State == kStateSolid)||(State == kStateLiquid)) {
141
142      const G4double E100eV  = 100.*eV; 
143      const G4double ClimiS[] = {3.681 , 5.215 };
144      const G4double X0valS[] = {1.0   , 1.5   };
145      const G4double X1valS[] = {2.0   , 3.0   };
146                               
147      if(fMeanExcitationEnergy < E100eV) icase = 0;
148         else                            icase = 1;
149
150      if(fCdensity < ClimiS[icase]) fX0density = 0.2;
151         else                       fX0density = 0.326*fCdensity-X0valS[icase];
152
153      fX1density = X1valS[icase] ; fMdensity = 3.0;
154     
155      //special: Hydrogen
156      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
157         fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
158      }
159  }
160
161  //
162  // gases
163  //
164  if (State == kStateGas) { 
165
166      const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804};
167      const G4double X0valG[] = { 1.6 , 1.7 ,  1.8 ,  1.9 , 2.0   ,  2.0 };
168      const G4double X1valG[] = { 4.0 , 4.0 ,  4.0 ,  4.0 , 4.0   ,  5.0 };
169
170      icase = 5;
171      fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ; 
172      while((icase > 0)&&(fCdensity < ClimiG[icase])) icase-- ;
173      fX0density = X0valG[icase]  ; fX1density = X1valG[icase] ;
174     
175      //special: Hydrogen
176      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
177         fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
178      }
179     
180      //special: Helium
181      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==2.)) {
182         fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
183      }
184
185      // change parameters if the gas is not in STP.
186      // For the correction the density(STP) is needed.
187      // Density(STP) is calculated here :
188     
189      G4double Density  = fMaterial->GetDensity();
190      G4double Pressure = fMaterial->GetPressure();
191      G4double Temp     = fMaterial->GetTemperature();
192     
193     G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
194
195      G4double ParCorr = std::log(Density/DensitySTP);
196 
197      fCdensity  -= ParCorr;
198      fX0density -= ParCorr/twoln10;
199      fX1density -= ParCorr/twoln10;
200  }
201
202  G4double Xa = fCdensity/twoln10;
203  fAdensity = twoln10*(Xa-fX0density)
204              /std::pow((fX1density-fX0density),fMdensity);
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
208
209void G4IonisParamMat::ComputeFluctModel()
210{
211  // compute parameters for the energy loss fluctuation model
212
213  // need an 'effective Z' ?????
214  G4double Zeff = 0.;
215  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) {
216     Zeff += (fMaterial->GetFractionVector())[i]
217             *((*(fMaterial->GetElementVector()))[i]->GetZ());
218  }
219  if (Zeff > 2.) fF2fluct = 2./Zeff ;
220    else         fF2fluct = 0.;
221
222  fF1fluct         = 1. - fF2fluct;
223  fEnergy2fluct    = 10.*Zeff*Zeff*eV;
224  fLogEnergy2fluct = std::log(fEnergy2fluct);
225  fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
226                     /fF1fluct;
227  fEnergy1fluct    = std::exp(fLogEnergy1fluct);
228  fEnergy0fluct    = 10.*eV;
229  fRateionexcfluct = 0.4;
230}
231
232//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
233
234void G4IonisParamMat::ComputeIonParameters()
235{
236  // get elements in the actual material,
237  const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
238  const G4double* theAtomicNumDensityVector =
239                         fMaterial->GetAtomicNumDensityVector() ;
240  const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
241
242  //  loop for the elements in the material
243  //  to find out average values Z, vF, lF
244  G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
245
246  if( 1 == NumberOfElements ) {
247    const G4Element* element = (*theElementVector)[0];
248    z = element->GetZ();
249    vF= element->GetIonisation()->GetFermiVelocity();
250    lF= element->GetIonisation()->GetLFactor();
251
252  } else {
253    for (G4int iel=0; iel<NumberOfElements; iel++)
254      {
255        const G4Element* element = (*theElementVector)[iel] ;
256        const G4double weight = theAtomicNumDensityVector[iel] ;
257        norm += weight ;
258        z    += element->GetZ() * weight ;
259        vF   += element->GetIonisation()->GetFermiVelocity() * weight ;
260        lF   += element->GetIonisation()->GetLFactor() * weight ;
261      }
262    z  /= norm ;
263    vF /= norm ;
264    lF /= norm ;
265  } 
266  fZeff        = z;
267  fLfactor     = lF;
268  fFermiEnergy = 25.*keV*vF*vF;
269}
270
271//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
272
273void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
274{
275  if(value == fMeanExcitationEnergy || value <= 0.0) return;
276
277  /*
278  if (G4NistManager::Instance()->GetVerbose() > 0)
279    G4cout << "G4Material: Mean excitation energy is changed for "
280           << fMaterial->GetName()
281           << " Iold= " << fMeanExcitationEnergy/eV
282           << "eV; Inew= " << value/eV << " eV;"
283           << G4endl;
284  */
285 
286  fMeanExcitationEnergy = value;
287  fLogMeanExcEnergy = std::log(value);
288  ComputeDensityEffect();
289  ComputeFluctModel();
290}
291
292//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
293
294G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
295{
296
297  // The data on mean excitation energy for compaunds
298  // from "Stopping Powers for Electrons and Positrons"
299  // ICRU Report N#37, 1984  (energy in eV)
300
301  const size_t numberOfMolecula = 79 ;
302 
303  static G4String name[numberOfMolecula] = {
304
305    // gas
306    "NH_3",       "C_4H_10",    "CO_2",       "C_2H_6",      "C_7H_16",
307    "C_6H_14",    "CH_4",       "NO",         "N_2O",        "C_8H_18",
308    "C_5H_12",    "C_3H_8",     "H_2O-Gas", 
309
310    // liquid
311    "C_3H_6O",    "C_6H_5NH_2",  "C_6H_6",    "C_4H_9OH",    "CCl_4",   
312    "C_6H_5Cl",   "CHCl_3",      "C_6H_12",   "C_6H_4Cl_2",  "C_4Cl_2H_8O", 
313    "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH",  "C_3H_5(OH)_3","C_7H_16",     
314    "C_6H_14",    "CH_3OH",      "C_6H_5NO_2","C_5H_12",     "C_3H_7OH",   
315    "C_5H_5N",    "C_8H_8",      "C_2Cl_4",   "C_7H_8",      "C_2Cl_3H",   
316    "H_2O",       "C_8H_10",
317
318    //solid
319    "C_5H_5N_5",  "C_5H_5N_5O",  "(C_6H_11NO)-nylon",  "C_25H_52", 
320    "(C_2H_4)-Polyethylene",     "(C_5H_8O-2)-Polymethil_Methacrylate",   
321    "(C_8H_8)-Polystyrene",      "A-150-tissue",       "Al_2O_3",  "CaF_2", 
322    "LiF",        "Photo_Emulsion",  "(C_2F_4)-Teflon",  "SiO_2"     
323
324  } ;
325   
326  static G4double meanExcitation[numberOfMolecula] = {
327
328    53.7,   48.3,  85.0,  45.4,  49.2,
329    49.1,   41.7,  87.8,  84.9,  49.5,
330    48.2,   47.1,  71.6,
331
332    64.2,   66.2,  63.4,  59.9,  166.3,
333    89.1,  156.0,  56.4, 106.5,  103.3, 
334   111.9,   60.0,  62.9,  72.6,   54.4, 
335    54.0,  67.6,   75.8,  53.6,   61.1, 
336    66.2,  64.0,  159.2,  62.5,  148.1, 
337    75.0,  61.8,
338
339    71.4,  75.0,   63.9,  48.3,   57.4,
340    74.0,  68.7,   65.1, 145.2,  166.,
341    94.0, 331.0,   99.1, 139.2 
342
343  } ;
344
345  G4double x = fMeanExcitationEnergy;
346
347  for(size_t i=0; i<numberOfMolecula; i++) {
348    if(chFormula == name[i]) {
349      x = meanExcitation[i]*eV;
350      break;
351    }
352  }
353  return x;
354}
355
356//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
357
358G4IonisParamMat::~G4IonisParamMat()
359{
360  if (fShellCorrectionVector) delete [] fShellCorrectionVector;
361}
362
363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
364
365G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
366{
367  *this = right;
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
371
372const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
373{
374  if (this != &right)
375    {
376      fMaterial                 = right.fMaterial;
377      fMeanExcitationEnergy     = right.fMeanExcitationEnergy;
378      fLogMeanExcEnergy         = right.fLogMeanExcEnergy;
379      if (fShellCorrectionVector) delete [] fShellCorrectionVector;     
380      fShellCorrectionVector    = new G4double[3];             
381      fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
382      fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
383      fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
384      fTaul                     = right.fTaul;
385      fCdensity                 = right.fCdensity;
386      fMdensity                 = right.fMdensity;
387      fAdensity                 = right.fAdensity;
388      fX0density                = right.fX0density;
389      fX1density                = right.fX1density;
390      fF1fluct                  = right.fF1fluct;
391      fF2fluct                  = right.fF2fluct;
392      fEnergy1fluct             = right.fEnergy1fluct;
393      fLogEnergy1fluct          = right.fLogEnergy1fluct;     
394      fEnergy2fluct             = right.fEnergy2fluct;
395      fLogEnergy2fluct          = right.fLogEnergy2fluct;     
396      fEnergy0fluct             = right.fEnergy0fluct;
397      fRateionexcfluct          = right.fRateionexcfluct;
398     } 
399  return *this;
400}
401
402//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
403
404G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
405{
406  return (this == (G4IonisParamMat*) &right);
407}
408
409//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
410
411G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
412{
413  return (this != (G4IonisParamMat*) &right);
414}
415
416//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
417
Note: See TracBrowser for help on using the repository browser.