source: trunk/source/materials/src/G4IonisParamElm.cc@ 1317

Last change on this file since 1317 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 9.0 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: G4IonisParamElm.cc,v 1.16 2010/04/29 11:11:56 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30//
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
32//
33// 09-07-98, data moved from G4Element. M.Maire
34// 22-11-00, tabulation of ionisation potential from
35// the ICRU Report N#37. V.Ivanchenko
36// 08-03-01, correct handling of fShellCorrectionVector. M.Maire
37// 17-10-02, Fix excitation energy interpolation. V.Ivanchenko
38// 06-09-04, Update calculated values after any change of ionisation
39// potential change. V.Ivanchenko
40// 29-04-10, Using G4Pow and mean ionisation energy from NIST V.Ivanchenko
41//
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
43
44#include "G4IonisParamElm.hh"
45#include "G4NistManager.hh"
46#include "G4Pow.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
49
50G4IonisParamElm::G4IonisParamElm(G4double AtomNumber)
51{
52 G4int Z = G4int(AtomNumber + 0.5);
53 if (Z < 1) {
54 G4Exception (" ERROR! It is not allowed to create an Element with Z < 1" );
55 }
56 G4Pow* g4pow = G4Pow::GetInstance();
57
58 // some basic functions of the atomic number
59 fZ = Z;
60 fZ3 = g4pow->Z13(Z);
61 fZZ3 = fZ3*g4pow->Z13(Z+1);
62 flogZ3 = g4pow->logZ(Z)/3.;
63
64 // Parameters for energy loss by ionisation
65 /*
66 // Mean excitation energy
67 // from "Stopping Powers for Electrons and Positrons"
68 // ICRU Report N#37, 1984 (energy in eV)
69 static double exc[100] = {
70 21.8, 20.9, 13.3, 15.9, 15.2, 13.0, 11.7, 11.9, 11.5, 13.7,
71 13.6, 13.0, 12.8, 12.4, 11.5, 11.3, 10.2, 10.4, 10.0, 9.6,
72 10.3, 10.6, 10.7, 10.7, 10.9, 11.0, 11.0, 11.1, 11.1, 11.0,
73 10.8, 10.4, 10.5, 10.2, 9.8, 9.8, 9.8, 9.6, 9.7, 9.8,
74 10.2, 10.1, 10.0, 10.0, 10.0, 10.2, 10.0, 9.8, 10.0, 9.8,
75 9.5, 9.3, 9.3, 8.9, 8.9, 8.8, 8.8, 8.8, 9.1, 9.1,
76 9.2, 9.3, 9.2, 9.2, 9.4, 9.5, 9.7, 9.7, 9.8, 9.8,
77 9.8, 9.8, 9.8, 9.8, 9.8, 9.8, 9.8, 10.1, 10.0, 10.0,
78 10.0, 10.0, 9.9, 9.9, 9.7, 9.2, 9.5, 9.4, 9.4, 9.4,
79 9.6, 9.7, 9.7, 9.8, 9.8, 9.8, 9.8, 9.9, 9.9, 9.9 };
80
81 fMeanExcitationEnergy = fZ * exc[iz] * eV ;
82 */
83
84 fMeanExcitationEnergy =
85 G4NistManager::Instance()->GetMeanIonisationEnergy(Z);
86
87 // compute parameters for ion transport
88 // The aproximation from:
89 // J.F.Ziegler, J.P. Biersack, U. Littmark
90 // The Stopping and Range of Ions in Matter,
91 // Vol.1, Pergamon Press, 1985
92 // Fast ions or hadrons
93
94 G4int iz = Z - 1;
95 if(91 < iz) { iz = 91; }
96
97 static G4double vFermi[92] = {
98 1.0309, 0.15976, 0.59782, 1.0781, 1.0486, 1.0, 1.058, 0.93942, 0.74562, 0.3424,
99 0.45259, 0.71074, 0.90519, 0.97411, 0.97184, 0.89852, 0.70827, 0.39816, 0.36552, 0.62712,
100 0.81707, 0.9943, 1.1423, 1.2381, 1.1222, 0.92705, 1.0047, 1.2, 1.0661, 0.97411,
101 0.84912, 0.95, 1.0903, 1.0429, 0.49715, 0.37755, 0.35211, 0.57801, 0.77773, 1.0207,
102 1.029, 1.2542, 1.122, 1.1241, 1.0882, 1.2709, 1.2542, 0.90094, 0.74093, 0.86054,
103 0.93155, 1.0047, 0.55379, 0.43289, 0.32636, 0.5131, 0.695, 0.72591, 0.71202, 0.67413,
104 0.71418, 0.71453, 0.5911, 0.70263, 0.68049, 0.68203, 0.68121, 0.68532, 0.68715, 0.61884,
105 0.71801, 0.83048, 1.1222, 1.2381, 1.045, 1.0733, 1.0953, 1.2381, 1.2879, 0.78654,
106 0.66401, 0.84912, 0.88433, 0.80746, 0.43357, 0.41923, 0.43638, 0.51464, 0.73087, 0.81065,
107 1.9578, 1.0257} ;
108
109 static G4double lFactor[92] = {
110 1.0, 1.0, 1.1, 1.06, 1.01, 1.03, 1.04, 0.99, 0.95, 0.9,
111 0.82, 0.81, 0.83, 0.88, 1.0, 0.95, 0.97, 0.99, 0.98, 0.97,
112 0.98, 0.97, 0.96, 0.93, 0.91, 0.9, 0.88, 0.9, 0.9, 0.9,
113 0.9, 0.85, 0.9, 0.9, 0.91, 0.92, 0.9, 0.9, 0.9, 0.9,
114 0.9, 0.88, 0.9, 0.88, 0.88, 0.9, 0.9, 0.88, 0.9, 0.9,
115 0.9, 0.9, 0.96, 1.2, 0.9, 0.88, 0.88, 0.85, 0.9, 0.9,
116 0.92, 0.95, 0.99, 1.03, 1.05, 1.07, 1.08, 1.1, 1.08, 1.08,
117 1.08, 1.08, 1.09, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15,
118 1.17, 1.2, 1.18, 1.17, 1.17, 1.16, 1.16, 1.16, 1.16, 1.16,
119 1.16, 1.16} ;
120
121 fVFermi = vFermi[iz];
122 fLFactor = lFactor[iz];
123
124 // obsolete parameters for ionisation
125 fTau0 = 0.1*fZ3*MeV/proton_mass_c2;
126 fTaul = 2.*MeV/proton_mass_c2;
127
128 // compute the Bethe-Bloch formula for energy = fTaul*particle mass
129 G4double rate = fMeanExcitationEnergy/electron_mass_c2 ;
130 G4double w = fTaul*(fTaul+2.) ;
131 fBetheBlochLow = (fTaul+1.)*(fTaul+1.)*std::log(2.*w/rate)/w - 1. ;
132 fBetheBlochLow = 2.*fZ*twopi_mc2_rcl2*fBetheBlochLow ;
133
134 fClow = std::sqrt(fTaul)*fBetheBlochLow;
135 fAlow = 6.458040 * fClow/fTau0;
136 G4double Taum = 0.035*fZ3*MeV/proton_mass_c2;
137 fBlow =-3.229020*fClow/(fTau0*std::sqrt(Taum));
138
139 // Shell correction parameterization
140 fShellCorrectionVector = new G4double[3]; //[3]
141 rate = 0.001*fMeanExcitationEnergy/eV;
142 G4double rate2 = rate*rate;
143 /*
144 fShellCorrectionVector[0] = ( 1.10289e5 + 5.14781e8*rate)*rate2 ;
145 fShellCorrectionVector[1] = ( 7.93805e3 - 2.22565e7*rate)*rate2 ;
146 fShellCorrectionVector[2] = (-9.92256e1 + 2.10823e5*rate)*rate2 ;
147 */
148 fShellCorrectionVector[0] = ( 0.422377 + 3.858019*rate)*rate2 ;
149 fShellCorrectionVector[1] = ( 0.0304043 - 0.1667989*rate)*rate2 ;
150 fShellCorrectionVector[2] = (-0.00038106 + 0.00157955*rate)*rate2 ;
151}
152
153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
154
155// Fake default constructor - sets only member data and allocates memory
156// for usage restricted to object persistency
157
158G4IonisParamElm::G4IonisParamElm(__void__&)
159 : fShellCorrectionVector(0)
160{}
161
162//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
163
164G4IonisParamElm::~G4IonisParamElm()
165{
166 if (fShellCorrectionVector) { delete [] fShellCorrectionVector; }
167}
168
169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
170
171G4IonisParamElm::G4IonisParamElm(G4IonisParamElm& right)
172{
173 *this = right;
174}
175
176//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
177
178const G4IonisParamElm& G4IonisParamElm::operator=(const G4IonisParamElm& right)
179{
180 if (this != &right)
181 {
182 fZ = right.fZ;
183 fZ3 = right.fZ3;
184 fZZ3 = right.fZZ3;
185 flogZ3 = right.flogZ3;
186 fTau0 = right.fTau0;
187 fTaul = right.fTaul;
188 fBetheBlochLow = right.fBetheBlochLow;
189 fAlow = right.fAlow;
190 fBlow = right.fBlow;
191 fClow = right.fClow;
192 fMeanExcitationEnergy = right.fMeanExcitationEnergy;
193 if (fShellCorrectionVector) { delete [] fShellCorrectionVector; }
194 fShellCorrectionVector = new G4double[3];
195 fShellCorrectionVector[0] = right.fShellCorrectionVector[0];
196 fShellCorrectionVector[1] = right.fShellCorrectionVector[1];
197 fShellCorrectionVector[2] = right.fShellCorrectionVector[2];
198 }
199 return *this;
200}
201
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
203
204G4int G4IonisParamElm::operator==(const G4IonisParamElm& right) const
205{
206 return (this == (G4IonisParamElm *) &right);
207}
208
209//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
210
211G4int G4IonisParamElm::operator!=(const G4IonisParamElm& right) const
212{
213 return (this != (G4IonisParamElm *) &right);
214}
215
216//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.