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

Last change on this file since 828 was 822, checked in by garnier, 16 years ago

import all except CVS
File size: 15.9 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.20 2007/09/27 14:05:47 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-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
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
43
44#include "G4IonisParamMat.hh"
45#include "G4Material.hh"
46#include "G4NistManager.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
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
60
61// Fake default constructor - sets only member data and allocates memory
62//                            for usage restricted to object persistency
63
64G4IonisParamMat::G4IonisParamMat(__void__&)
65  : fMaterial(0), fShellCorrectionVector(0)
66{
67}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
70
71void G4IonisParamMat::ComputeMeanParameters()
72{
73  // compute mean excitation energy and shell correction vector
74
75  fTaul = (*(fMaterial->GetElementVector()))[0]->GetIonisation()->GetTaul();
76
77  fMeanExcitationEnergy = 0.;
78  fLogMeanExcEnergy = 0.;
79
80
81  for (size_t i=0; i < fMaterial->GetNumberOfElements(); i++) {
82    fLogMeanExcEnergy += 
83             (fMaterial->GetVecNbOfAtomsPerVolume())[i]
84            *((*(fMaterial->GetElementVector()))[i]->GetZ())
85            *std::log((*(fMaterial->GetElementVector()))[i]->GetIonisation()
86             ->GetMeanExcitationEnergy());
87  }
88
89  fLogMeanExcEnergy /= fMaterial->GetTotNbOfElectPerVolume();
90  fMeanExcitationEnergy = std::exp(fLogMeanExcEnergy);
91
92  fShellCorrectionVector = new G4double[3]; //[3]
93
94  for (G4int j=0; j<=2; j++)
95  {
96    fShellCorrectionVector[j] = 0.;
97
98    for (size_t k=0; k<fMaterial->GetNumberOfElements(); k++) {
99      fShellCorrectionVector[j] += (fMaterial->GetVecNbOfAtomsPerVolume())[k] 
100              *((*(fMaterial->GetElementVector()))[k]->GetIonisation()
101                                              ->GetShellCorrectionVector()[j]);
102    }
103    fShellCorrectionVector[j] *= 2.0/fMaterial->GetTotNbOfElectPerVolume();
104  } 
105}
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
108                   
109void G4IonisParamMat::ComputeDensityEffect()
110{
111  // Compute parameters for the density effect correction in DE/Dx formula.
112  // The parametrization is from R.M. Sternheimer, Phys. Rev.B,3:3681 (1971)
113
114  const G4double Cd2 = 4*pi*hbarc_squared*classic_electr_radius;
115  const G4double twoln10 = 2.*std::log(10.);
116
117  G4int icase;
118 
119  fCdensity = 1. + std::log(fMeanExcitationEnergy*fMeanExcitationEnergy
120              /(Cd2*fMaterial->GetTotNbOfElectPerVolume()));
121
122  //
123  // condensed materials
124  //
125  G4State State = fMaterial->GetState();
126 
127  if ((State == kStateSolid)||(State == kStateLiquid)) {
128
129      const G4double E100eV  = 100.*eV; 
130      const G4double ClimiS[] = {3.681 , 5.215 };
131      const G4double X0valS[] = {1.0   , 1.5   };
132      const G4double X1valS[] = {2.0   , 3.0   };
133                               
134      if(fMeanExcitationEnergy < E100eV) icase = 0;
135         else                            icase = 1;
136
137      if(fCdensity < ClimiS[icase]) fX0density = 0.2;
138         else                       fX0density = 0.326*fCdensity-X0valS[icase];
139
140      fX1density = X1valS[icase] ; fMdensity = 3.0;
141     
142      //special: Hydrogen
143      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
144         fX0density = 0.425; fX1density = 2.0; fMdensity = 5.949;
145      }
146  }
147
148  //
149  // gases
150  //
151  if (State == kStateGas) { 
152
153      const G4double ClimiG[] = { 10. , 10.5 , 11. , 11.5 , 12.25 , 13.804};
154      const G4double X0valG[] = { 1.6 , 1.7 ,  1.8 ,  1.9 , 2.0   ,  2.0 };
155      const G4double X1valG[] = { 4.0 , 4.0 ,  4.0 ,  4.0 , 4.0   ,  5.0 };
156
157      icase = 5;
158      fX0density = 0.326*fCdensity-2.5 ; fX1density = 5.0 ; fMdensity = 3. ; 
159      while((icase > 0)&&(fCdensity < ClimiG[icase])) icase-- ;
160      fX0density = X0valG[icase]  ; fX1density = X1valG[icase] ;
161     
162      //special: Hydrogen
163      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==1.)) {
164         fX0density = 1.837; fX1density = 3.0; fMdensity = 4.754;
165      }
166     
167      //special: Helium
168      if ((fMaterial->GetNumberOfElements()==1)&&(fMaterial->GetZ()==2.)) {
169         fX0density = 2.191; fX1density = 3.0; fMdensity = 3.297;
170      }
171
172      // change parameters if the gas is not in STP.
173      // For the correction the density(STP) is needed.
174      // Density(STP) is calculated here :
175     
176      G4double Density  = fMaterial->GetDensity();
177      G4double Pressure = fMaterial->GetPressure();
178      G4double Temp     = fMaterial->GetTemperature();
179     
180     G4double DensitySTP = Density*STP_Pressure*Temp/(Pressure*STP_Temperature);
181
182      G4double ParCorr = std::log(Density/DensitySTP);
183 
184      fCdensity  -= ParCorr;
185      fX0density -= ParCorr/twoln10;
186      fX1density -= ParCorr/twoln10;
187  }
188
189  G4double Xa = fCdensity/twoln10;
190  fAdensity = twoln10*(Xa-fX0density)
191              /std::pow((fX1density-fX0density),fMdensity);
192}
193
194//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
195
196void G4IonisParamMat::ComputeFluctModel()
197{
198  // compute parameters for the energy loss fluctuation model
199
200  // need an 'effective Z' ?????
201  G4double Zeff = 0.;
202  for (size_t i=0;i<fMaterial->GetNumberOfElements();i++) 
203     Zeff += (fMaterial->GetFractionVector())[i]
204             *((*(fMaterial->GetElementVector()))[i]->GetZ());
205
206  if (Zeff > 2.) fF2fluct = 2./Zeff ;
207    else         fF2fluct = 0.;
208
209  fF1fluct         = 1. - fF2fluct;
210  fEnergy2fluct    = 10.*Zeff*Zeff*eV;
211  fLogEnergy2fluct = std::log(fEnergy2fluct);
212  fLogEnergy1fluct = (fLogMeanExcEnergy - fF2fluct*fLogEnergy2fluct)
213                     /fF1fluct;
214  fEnergy1fluct    = std::exp(fLogEnergy1fluct);
215  fEnergy0fluct    = 10.*eV;
216  fRateionexcfluct = 0.4;
217}
218
219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
220
221void G4IonisParamMat::ComputeIonParameters()
222{
223  // compute parameters for ion transport
224  // The aproximation from:
225  // J.F.Ziegler, J.P. Biersack, U. Littmark
226  // The Stopping and Range of Ions in Matter,
227  // Vol.1, Pergamon Press, 1985
228  // Fast ions or hadrons
229
230  static G4double vFermi[92] = {
231    1.0309,  0.15976, 0.59782, 1.0781,  1.0486,  1.0,     1.058,   0.93942, 0.74562, 0.3424,
232    0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
233    0.81707, 0.9943,  1.1423,  1.2381,  1.1222,  0.92705, 1.0047,  1.2,     1.0661,  0.97411,
234    0.84912, 0.95,    1.0903,  1.0429,  0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
235    1.029,   1.2542,  1.122,   1.1241,  1.0882,  1.2709,  1.2542,  0.90094, 0.74093, 0.86054,
236    0.93155, 1.0047,  0.55379, 0.43289, 0.32636, 0.5131,  0.695,   0.72591, 0.71202, 0.67413,
237    0.71418, 0.71453, 0.5911,  0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
238    0.71801, 0.83048, 1.1222,  1.2381,  1.045,   1.0733,  1.0953,  1.2381,  1.2879,  0.78654,
239    0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
240    1.9578,  1.0257} ;
241
242  static G4double lFactor[92] = {
243    1.0,  1.0,  1.1,  1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
244    0.82, 0.81, 0.83, 0.88, 1.0,  0.95, 0.97, 0.99, 0.98, 0.97,
245    0.98, 0.97, 0.96, 0.93, 0.91, 0.9,  0.88, 0.9,  0.9,  0.9,
246    0.9,  0.85, 0.9,  0.9,  0.91, 0.92, 0.9,  0.9,  0.9,  0.9,
247    0.9,  0.88, 0.9,  0.88, 0.88, 0.9,  0.9,  0.88, 0.9,  0.9,
248    0.9,  0.9,  0.96, 1.2,  0.9,  0.88, 0.88, 0.85, 0.9,  0.9,
249    0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1,  1.08, 1.08,
250    1.08, 1.08, 1.09, 1.09, 1.1,  1.11, 1.12, 1.13, 1.14, 1.15,
251    1.17, 1.2,  1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
252    1.16, 1.16} ;
253
254  // get elements in the actual material,
255  const G4ElementVector* theElementVector = fMaterial->GetElementVector() ;
256  const G4double* theAtomicNumDensityVector =
257                         fMaterial->GetAtomicNumDensityVector() ;
258  const G4int NumberOfElements = fMaterial->GetNumberOfElements() ;
259
260  //  loop for the elements in the material
261  //  to find out average values Z, vF, lF
262  G4double z = 0.0, vF = 0.0, lF = 0.0, norm = 0.0 ;
263
264  if( 1 == NumberOfElements ) {
265    z = fMaterial->GetZ() ;
266    G4int iz = G4int(z) - 1 ;
267    if(iz < 0) iz = 0 ;
268    else if(iz > 91) iz = 91 ;
269    vF   = vFermi[iz] ;
270    lF   = lFactor[iz] ;
271
272  } else {
273    for (G4int iel=0; iel<NumberOfElements; iel++)
274      {
275        const G4Element* element = (*theElementVector)[iel] ;
276        G4double z2 = element->GetZ() ;
277        const G4double weight = theAtomicNumDensityVector[iel] ;
278        norm += weight ;
279        z    += z2 * weight ;
280        G4int iz = G4int(z2) - 1 ;
281        if(iz < 0) iz = 0 ;
282        else if(iz > 91) iz =91 ;
283        vF   += vFermi[iz] * weight ;
284        lF   += lFactor[iz] * weight ;
285      }
286    z  /= norm ;
287    vF /= norm ;
288    lF /= norm ;
289  } 
290  fZeff        = z;
291  fLfactor     = lF;
292  fFermiEnergy = 25.*keV*vF*vF;
293}
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
296
297void G4IonisParamMat::SetMeanExcitationEnergy(G4double value)
298{
299  if(value == fMeanExcitationEnergy || value <= 0.0) return;
300
301  if (G4NistManager::Instance()->GetVerbose() > 0) 
302    G4cout << "G4Material: Mean excitation energy is changed for "
303           << fMaterial->GetName()
304           << " Iold= " << fMeanExcitationEnergy/eV
305           << "eV; Inew= " << value/eV << " eV;"
306           << G4endl;
307
308  fMeanExcitationEnergy = value;
309  fLogMeanExcEnergy = std::log(value);
310  ComputeDensityEffect();
311  ComputeFluctModel();
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
315
316G4double G4IonisParamMat::FindMeanExcitationEnergy(const G4String& chFormula)
317{
318
319  // The data on mean excitation energy for compaunds
320  // from "Stopping Powers for Electrons and Positrons"
321  // ICRU Report N#37, 1984  (energy in eV)
322
323  const size_t numberOfMolecula = 79 ;
324 
325  static G4String name[numberOfMolecula] = {
326
327    // gas
328    "NH_3",       "C_4H_10",    "CO_2",       "C_2H_6",      "C_7H_16",
329    "C_6H_14",    "CH_4",       "NO",         "N_2O",        "C_8H_18",
330    "C_5H_12",    "C_3H_8",     "H_2O-Gas", 
331
332    // liquid
333    "C_3H_6O",    "C_6H_5NH_2",  "C_6H_6",    "C_4H_9OH",    "CCl_4",   
334    "C_6H_5Cl",   "CHCl_3",      "C_6H_12",   "C_6H_4Cl_2",  "C_4Cl_2H_8O", 
335    "C_2Cl_2H_4", "(C_2H_5)_2O", "C_2H_5OH",  "C_3H_5(OH)_3","C_7H_16",     
336    "C_6H_14",    "CH_3OH",      "C_6H_5NO_2","C_5H_12",     "C_3H_7OH",   
337    "C_5H_5N",    "C_8H_8",      "C_2Cl_4",   "C_7H_8",      "C_2Cl_3H",   
338    "H_2O",       "C_8H_10",
339
340    //solid
341    "C_5H_5N_5",  "C_5H_5N_5O",  "(C_6H_11NO)-nylon",  "C_25H_52", 
342    "(C_2H_4)-Polyethylene",     "(C_5H_8O-2)-Polymethil_Methacrylate",   
343    "(C_8H_8)-Polystyrene",      "A-150-tissue",       "Al_2O_3",  "CaF_2", 
344    "LiF",        "Photo_Emulsion",  "(C_2F_4)-Teflon",  "SiO_2"     
345
346  } ;
347   
348  static G4double meanExcitation[numberOfMolecula] = {
349
350    53.7,   48.3,  85.0,  45.4,  49.2,
351    49.1,   41.7,  87.8,  84.9,  49.5,
352    48.2,   47.1,  71.6,
353
354    64.2,   66.2,  63.4,  59.9,  166.3,
355    89.1,  156.0,  56.4, 106.5,  103.3, 
356   111.9,   60.0,  62.9,  72.6,   54.4, 
357    54.0,  67.6,   75.8,  53.6,   61.1, 
358    66.2,  64.0,  159.2,  62.5,  148.1, 
359    75.0,  61.8,
360
361    71.4,  75.0,   63.9,  48.3,   57.4,
362    74.0,  68.7,   65.1, 145.2,  166.,
363    94.0, 331.0,   99.1, 139.2 
364
365  } ;
366
367  G4double x = fMeanExcitationEnergy;
368
369  for(size_t i=0; i<numberOfMolecula; i++) {
370    if(chFormula == name[i]) {
371      x = meanExcitation[i]*eV;
372      break;
373    }
374  }
375  return x;
376}
377
378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
379
380G4IonisParamMat::~G4IonisParamMat()
381{
382  if (fShellCorrectionVector) delete [] fShellCorrectionVector;
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
386
387G4IonisParamMat::G4IonisParamMat(const G4IonisParamMat& right)
388{
389  *this = right;
390}
391
392//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
393
394const G4IonisParamMat& G4IonisParamMat::operator=(const G4IonisParamMat& right)
395{
396  if (this != &right)
397    {
398      fMaterial                 = right.fMaterial;
399      fMeanExcitationEnergy     = right.fMeanExcitationEnergy;
400      fLogMeanExcEnergy         = right.fLogMeanExcEnergy;
401      if (fShellCorrectionVector) delete [] fShellCorrectionVector;     
402      fShellCorrectionVector    = new G4double[3];             
403      fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
404      fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
405      fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
406      fTaul                     = right.fTaul;
407      fCdensity                 = right.fCdensity;
408      fMdensity                 = right.fMdensity;
409      fAdensity                 = right.fAdensity;
410      fX0density                = right.fX0density;
411      fX1density                = right.fX1density;
412      fF1fluct                  = right.fF1fluct;
413      fF2fluct                  = right.fF2fluct;
414      fEnergy1fluct             = right.fEnergy1fluct;
415      fLogEnergy1fluct          = right.fLogEnergy1fluct;     
416      fEnergy2fluct             = right.fEnergy2fluct;
417      fLogEnergy2fluct          = right.fLogEnergy2fluct;     
418      fEnergy0fluct             = right.fEnergy0fluct;
419      fRateionexcfluct          = right.fRateionexcfluct;
420     } 
421  return *this;
422}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
425
426G4int G4IonisParamMat::operator==(const G4IonisParamMat& right) const
427{
428  return (this == (G4IonisParamMat*) &right);
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
432
433G4int G4IonisParamMat::operator!=(const G4IonisParamMat& right) const
434{
435  return (this != (G4IonisParamMat*) &right);
436}
437
438//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
439
Note: See TracBrowser for help on using the repository browser.