| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 26
>> 27 #include <cmath>
27 #include <iostream> 28 #include <iostream>
>> 29
28 #include "G4ecpssrBaseLixsModel.hh" 30 #include "G4ecpssrBaseLixsModel.hh"
>> 31
29 #include "globals.hh" 32 #include "globals.hh"
30 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
32 #include "G4AtomicTransitionManager.hh" 35 #include "G4AtomicTransitionManager.hh"
33 #include "G4NistManager.hh" 36 #include "G4NistManager.hh"
34 #include "G4Proton.hh" 37 #include "G4Proton.hh"
35 #include "G4Alpha.hh" 38 #include "G4Alpha.hh"
36 #include "G4LinLogInterpolation.hh" 39 #include "G4LinLogInterpolation.hh"
37 #include "G4Exp.hh" <<
38 40
39 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 42
41 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel() 43 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel()
42 { 44 {
43 verboseLevel=0; << 45 verboseLevel=0;
44 46
45 // Storing FLi data needed for 0.2 to 3.0 v << 47 // Storing FLi data needed for 0.2 to 3.0 velocities region
46 const char* path = G4FindDataDir("G4LEDATA") <<
47 48
48 if (!path) { << 49 char *path = getenv("G4LEDATA");
49 G4Exception("G4ecpssrLCrossSection::G4ecps <<
50 return; <<
51 } <<
52 std::ostringstream fileName1; <<
53 std::ostringstream fileName2; <<
54 50
55 fileName1 << path << "/pixe/uf/FL1.dat"; << 51 if (!path) {
56 fileName2 << path << "/pixe/uf/FL2.dat"; << 52 G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0006", FatalException ,"G4LEDATA environment variable not set");
>> 53 return;
>> 54 }
>> 55 std::ostringstream fileName1;
>> 56 std::ostringstream fileName2;
>> 57
>> 58 fileName1 << path << "/pixe/uf/FL1.dat";
>> 59 fileName2 << path << "/pixe/uf/FL2.dat";
57 60
58 // Reading of FL1.dat << 61 // Reading of FL1.dat
59 std::ifstream FL1(fileName1.str().c_str()); <<
60 if (!FL1) G4Exception("G4ecpssrLCrossSection <<
61 62
62 dummyVec1.push_back(0.); << 63 std::ifstream FL1(fileName1.str().c_str());
>> 64 if (!FL1) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003",FatalException, "error opening FL1 data file");
>> 65
>> 66 dummyVec1.push_back(0.);
63 67
64 while(!FL1.eof()) << 68 while(!FL1.eof())
65 { 69 {
66 double x1; << 70 double x1;
67 double y1; << 71 double y1;
68 72
69 FL1>>x1>>y1; << 73 FL1>>x1>>y1;
70 74
71 // Mandatory vector initialization << 75 // Mandatory vector initialization
72 if (x1 != dummyVec1.back()) << 76 if (x1 != dummyVec1.back())
73 { 77 {
74 dummyVec1.push_back(x1); 78 dummyVec1.push_back(x1);
75 aVecMap1[x1].push_back(-1.); 79 aVecMap1[x1].push_back(-1.);
76 } 80 }
77 81
78 FL1>>FL1Data[x1][y1]; << 82 FL1>>FL1Data[x1][y1];
79 83
80 if (y1 != aVecMap1[x1].back()) aVecMap1[ << 84 if (y1 != aVecMap1[x1].back()) aVecMap1[x1].push_back(y1);
81 } 85 }
82 86
83 // Reading of FL2.dat << 87 // Reading of FL2.dat
84 << 88
85 std::ifstream FL2(fileName2.str().c_str()); << 89 std::ifstream FL2(fileName2.str().c_str());
86 if (!FL2) G4Exception("G4ecpssrLCrossSection << 90 if (!FL2) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003", FatalException," error opening FL2 data file");
>> 91
>> 92 dummyVec2.push_back(0.);
87 93
88 dummyVec2.push_back(0.); << 94 while(!FL2.eof())
89 <<
90 while(!FL2.eof()) <<
91 { 95 {
92 double x2; << 96 double x2;
93 double y2; << 97 double y2;
94 98
95 FL2>>x2>>y2; << 99 FL2>>x2>>y2;
96 100
97 // Mandatory vector initialization << 101 // Mandatory vector initialization
98 if (x2 != dummyVec2.back()) << 102 if (x2 != dummyVec2.back())
99 { 103 {
100 dummyVec2.push_back(x2); 104 dummyVec2.push_back(x2);
101 aVecMap2[x2].push_back(-1.); 105 aVecMap2[x2].push_back(-1.);
102 } 106 }
103 107
104 FL2>>FL2Data[x2][y2]; << 108 FL2>>FL2Data[x2][y2];
105 109
106 if (y2 != aVecMap2[x2].back()) aVecMap2[ << 110 if (y2 != aVecMap2[x2].back()) aVecMap2[x2].push_back(y2);
107 } 111 }
>> 112
108 } 113 }
109 114
110 //....oooOO0OOooo........oooOO0OOooo........oo 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 116
112 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel( 117 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel()
113 {} 118 {}
114 119
115 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 121
117 G4double G4ecpssrBaseLixsModel::ExpIntFunction 122 G4double G4ecpssrBaseLixsModel::ExpIntFunction(G4int n,G4double x)
118 123
119 { 124 {
120 // this function allows fast evaluation of the 125 // this function allows fast evaluation of the n order exponential integral function En(x)
>> 126
121 G4int i; 127 G4int i;
122 G4int ii; 128 G4int ii;
123 G4int nm1; 129 G4int nm1;
124 G4double a; 130 G4double a;
125 G4double b; 131 G4double b;
126 G4double c; 132 G4double c;
127 G4double d; 133 G4double d;
128 G4double del; 134 G4double del;
129 G4double fact; 135 G4double fact;
130 G4double h; 136 G4double h;
131 G4double psi; 137 G4double psi;
132 G4double ans = 0; 138 G4double ans = 0;
133 static const G4double euler= 0.5772156649; << 139 const G4double euler= 0.5772156649;
134 static const G4int maxit= 100; << 140 const G4int maxit= 100;
135 static const G4double fpmin = 1.0e-30; << 141 const G4double fpmin = 1.0e-30;
136 static const G4double eps = 1.0e-7; << 142 const G4double eps = 1.0e-7;
137 nm1=n-1; 143 nm1=n-1;
138 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1 144 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
139 G4cout << "*** WARNING in G4ecpssrBaseLixsMo << 145 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction"
140 << G4endl; 146 << G4endl;
141 else { 147 else {
142 if (n==0) ans=G4Exp(-x)/x; << 148 if (n==0) ans=std::exp(-x)/x;
143 else { << 149 else {
144 if (x==0.0) ans=1.0/nm1; << 150 if (x==0.0) ans=1.0/nm1;
145 else { << 151 else {
146 if (x > 1.0) { << 152 if (x > 1.0) {
147 b=x+n; << 153 b=x+n;
148 c=1.0/fpmin; << 154 c=1.0/fpmin;
149 d=1.0/b; << 155 d=1.0/b;
150 h=d; << 156 h=d;
151 for (i=1;i<=maxit;i++) { << 157 for (i=1;i<=maxit;i++) {
152 a=-i*(nm1+i); << 158 a=-i*(nm1+i);
153 b +=2.0; << 159 b +=2.0;
154 d=1.0/(a*d+b); << 160 d=1.0/(a*d+b);
155 c=b+a/c; << 161 c=b+a/c;
156 del=c*d; << 162 del=c*d;
157 h *=del; << 163 h *=del;
158 if (std::fabs(del-1.0) < eps) { << 164 if (std::fabs(del-1.0) < eps) {
159 ans=h*G4Exp(-x); << 165 ans=h*std::exp(-x);
160 return ans; << 166 return ans;
161 } << 167 }
162 } << 168 }
163 } else { << 169 } else {
164 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-eul << 170 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
165 fact=1.0; << 171 fact=1.0;
166 for (i=1;i<=maxit;i++) { << 172 for (i=1;i<=maxit;i++) {
167 fact *=-x/i; << 173 fact *=-x/i;
168 if (i !=nm1) del = -fact/(i-nm1); << 174 if (i !=nm1) del = -fact/(i-nm1);
169 else { << 175 else {
170 psi = -euler; << 176 psi = -euler;
171 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii; << 177 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
172 del=fact*(-std::log(x)+psi); << 178 del=fact*(-std::log(x)+psi);
173 } << 179 }
174 ans += del; << 180 ans += del;
175 if (std::fabs(del) < std::fabs(ans)*eps) << 181 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
176 } << 182 }
>> 183 }
>> 184 }
177 } 185 }
178 } <<
179 } <<
180 } 186 }
181 return ans; << 187 return ans;
182 } 188 }
183 189
184 //....oooOO0OOooo........oooOO0OOooo........oo 190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 191
186 G4double G4ecpssrBaseLixsModel::CalculateL1Cro 192 G4double G4ecpssrBaseLixsModel::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
187 { 193 {
188 194
189 if (zTarget <=4) return 0.; 195 if (zTarget <=4) return 0.;
190 196
191 //this L1-CrossSection calculation method is << 197 //this L1-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
192 //and using data tables of O. Benka et al. A << 198 //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
193 199
194 G4NistManager* massManager = G4NistManager:: 200 G4NistManager* massManager = G4NistManager::Instance();
195 201
196 G4AtomicTransitionManager* transitionManage 202 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
197 203
198 G4double zIncident = 0; 204 G4double zIncident = 0;
199 G4Proton* aProtone = G4Proton::Proton(); 205 G4Proton* aProtone = G4Proton::Proton();
200 G4Alpha* aAlpha = G4Alpha::Alpha(); 206 G4Alpha* aAlpha = G4Alpha::Alpha();
201 207
202 if (massIncident == aProtone->GetPDGMass() ) 208 if (massIncident == aProtone->GetPDGMass() )
203 { << 209
204 zIncident = (aProtone->GetPDGCharge())/e << 210 zIncident = (aProtone->GetPDGCharge())/eplus;
205 } << 211
206 else 212 else
207 { 213 {
208 if (massIncident == aAlpha->GetPDGMass() << 214 if (massIncident == aAlpha->GetPDGMass())
209 { << 215
210 zIncident = (aAlpha->GetPDGCharge())/eplu 216 zIncident = (aAlpha->GetPDGCharge())/eplus;
211 } << 217
212 else 218 else
213 { 219 {
214 G4cout << "*** WARNING in G4ecpssrBaseLixs 220 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << G4endl;
215 G4cout << massIncident << ", " << aAlpha-> 221 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
216 return 0; 222 return 0;
217 } 223 }
218 } 224 }
219 225
220 G4double l1BindingEnergy = transitionManager 226 G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy(); //Observed binding energy of L1-subshell
>> 227
221 G4double massTarget = (massManager->GetAtomi 228 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
>> 229
222 G4double systemMass =((massIncident*massTarg 230 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
223 static const G4double zlshell= 4.15; << 231
>> 232 const G4double zlshell= 4.15;
224 // *** see Benka, ADANDT 22, p 223 233 // *** see Benka, ADANDT 22, p 223
>> 234
225 G4double screenedzTarget = zTarget-zlshell; 235 G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L1-sub shell
226 static const G4double rydbergMeV= 13.6056923 <<
227 236
228 static const G4double nl= 2.; << 237 const G4double rydbergMeV= 13.6056923e-6;
>> 238
>> 239 const G4double nl= 2.;
229 // *** see Benka, ADANDT 22, p 220, f3 240 // *** see Benka, ADANDT 22, p 220, f3
>> 241
230 G4double tetal1 = (l1BindingEnergy*nl*nl)/(( 242 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
231 // *** see Benka, ADANDT 22, p 220, f3 243 // *** see Benka, ADANDT 22, p 220, f3
232 244
233 if (verboseLevel>0) G4cout << " tetal1=" << << 245 if (verboseLevel>0) G4cout << " tetal1=" << tetal1<< G4endl;
234 246
235 G4double reducedEnergy = (energyIncident*ele 247 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
236 // *** also called etaS 248 // *** also called etaS
237 // *** see Benka, ADANDT 22, p 220, f3 249 // *** see Benka, ADANDT 22, p 220, f3
238 250
239 static const G4double bohrPow2Barn=(Bohr_rad << 251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
240 252
241 G4double sigma0 = 8.*pi*(zIncident*zIncident 253 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
242 // *** see Benka, ADANDT 22, p 220, f2, for 254 // *** see Benka, ADANDT 22, p 220, f2, for protons
243 // *** see Basbas, Phys Rev A7, p 1000 255 // *** see Basbas, Phys Rev A7, p 1000
244 256
245 G4double velocityl1 = CalculateVelocity(1, z 257 G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident); // Scaled velocity
246 258
247 if (verboseLevel>0) G4cout << " velocityl1= << 259 if (verboseLevel>0) G4cout << " velocityl1=" << velocityl1<< G4endl;
248 260
249 static const G4double l1AnalyticalApproximat << 261 const G4double l1AnalyticalApproximation= 1.5;
250 G4double x1 =(nl*l1AnalyticalApproximation)/ 262 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
251 // *** 1.5 is cK = cL1 (it is 1.25 for L2 & 263 // *** 1.5 is cK = cL1 (it is 1.25 for L2 & L3)
252 // *** see Brandt, Phys Rev A20, p 469, f16 264 // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h
253 << 265
254 if (verboseLevel>0) G4cout << " x1=" << x1< << 266 if (verboseLevel>0) G4cout << " x1=" << x1<< G4endl;
255 267
256 G4double electrIonizationEnergyl1=0.; 268 G4double electrIonizationEnergyl1=0.;
257 // *** see Basbas, Phys Rev A17, p1665, f27 269 // *** see Basbas, Phys Rev A17, p1665, f27
258 // *** see Brandt, Phys Rev A20, p469 270 // *** see Brandt, Phys Rev A20, p469
259 // *** see Liu, Comp Phys Comm 97, p325, f A 271 // *** see Liu, Comp Phys Comm 97, p325, f A5
260 272
261 if ( x1<=0.035) electrIonizationEnergyl1= 0 273 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*pi*(std::log(1./(x1*x1))-1.);
262 else 274 else
263 { 275 {
264 if ( x1<=3.) 276 if ( x1<=3.)
265 electrIonizationEnergyl1 =G4Exp(-2.*x1 << 277 electrIonizationEnergyl1 =std::exp(-2.*x1)/(0.031+(0.213*std::pow(x1,0.5))+(0.005*x1)-(0.069*std::pow(x1,3./2.))+(0.324*x1*x1));
266 else 278 else
267 {if ( x1<=11.) electrIonizationEnergyl1 =2.* << 279 {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
268 } 280 }
269 281
270 G4double hFunctionl1 =(electrIonizationEnerg 282 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3)); //takes into account the polarization effect
271 // *** see Brandt, Phys Rev A20, p 469, f16 283 // *** see Brandt, Phys Rev A20, p 469, f16
272 284
273 if (verboseLevel>0) G4cout << " hFunctionl1 << 285 if (verboseLevel>0) G4cout << " hFunctionl1=" << hFunctionl1<< G4endl;
274 286
275 G4double gFunctionl1 = (1.+(9.*velocityl1)+( 287 G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*std::pow(velocityl1,3.))+(162.*std::pow(velocityl1,4.))+(63.*std::pow(velocityl1,5.))+(18.*std::pow(velocityl1,6.))+(1.97*std::pow(velocityl1,7.)))/std::pow(1.+velocityl1,9.);//takes into account the reduced binding effect
276 // *** see Brandt, Phys Rev A20, p 469, f19 288 // *** see Brandt, Phys Rev A20, p 469, f19
277 289
278 if (verboseLevel>0) G4cout << " gFunctionl1 << 290 if (verboseLevel>0) G4cout << " gFunctionl1=" << gFunctionl1<< G4endl;
279 291
280 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/( 292 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //Binding-polarization factor
281 // *** also called dzeta 293 // *** also called dzeta
282 // *** also called epsilon 294 // *** also called epsilon
283 // *** see Basbas, Phys Rev A17, p1667, f45 295 // *** see Basbas, Phys Rev A17, p1667, f45
284 296
285 if (verboseLevel>0) G4cout << "sigmaPSS_l1 = << 297 if (verboseLevel>0) G4cout << "sigmaPSS_l1 =" << sigmaPSS_l1<< G4endl;
286 298
287 const G4double cNaturalUnit= 137.; 299 const G4double cNaturalUnit= 137.;
288 << 300
289 G4double yl1Formula=0.4*(screenedzTarget/cNa 301 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
290 // *** also called yS 302 // *** also called yS
291 // *** see Brandt, Phys Rev A20, p467, f6 303 // *** see Brandt, Phys Rev A20, p467, f6
292 // *** see Brandt, Phys Rev A23, p1728 304 // *** see Brandt, Phys Rev A23, p1728
293 << 305
294 G4double l1relativityCorrection = std::pow(( 306 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula; // Relativistic correction parameter
295 // *** also called mRS 307 // *** also called mRS
296 // *** see Brandt, Phys Rev A20, p467, f6 308 // *** see Brandt, Phys Rev A20, p467, f6
297 << 309
298 //G4double reducedVelocity_l1 = velocityl1*s 310 //G4double reducedVelocity_l1 = velocityl1*std::pow(l1relativityCorrection,0.5); //Reduced velocity parameter
299 311
300 G4double L1etaOverTheta2; 312 G4double L1etaOverTheta2;
301 313
302 G4double universalFunction_l1 = 0.; 314 G4double universalFunction_l1 = 0.;
303 315
304 G4double sigmaPSSR_l1; 316 G4double sigmaPSSR_l1;
305 317
306 // low velocity formula 318 // low velocity formula
307 // ***************** << 319 // *****************
308 if ( velocityl1 <20. ) 320 if ( velocityl1 <20. )
309 { 321 {
310 322
311 L1etaOverTheta2 =(reducedEnergy* l1relativ 323 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
312 // *** 1) RELATIVISTIC CORRECTION ADDED 324 // *** 1) RELATIVISTIC CORRECTION ADDED
313 // *** 2) sigma_PSS_l1 ADDED 325 // *** 2) sigma_PSS_l1 ADDED
314 // *** reducedEnergy is etaS, l1relativity 326 // *** reducedEnergy is etaS, l1relativityCorrection is mRS
315 // *** see Phys Rev A20, p468, top 327 // *** see Phys Rev A20, p468, top
316 << 328
317 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tet 329 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
>> 330
318 universalFunction_l1 = FunctionFL1((teta 331 universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
319 332
320 if (verboseLevel>0) G4cout << "at low velo << 333 if (verboseLevel>0) G4cout << "at low velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
321 334
322 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1 335 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
323 // *** see Benka, ADANDT 22, p220, f1 336 // *** see Benka, ADANDT 22, p220, f1
324 337
325 if (verboseLevel>0) G4cout << " at low ve << 338 if (verboseLevel>0) G4cout << " at low velocity range, sigma PWBA L1 CS = " << sigmaPSSR_l1<< G4endl;
326 } << 339
>> 340 }
>> 341
327 else 342 else
>> 343
328 { 344 {
>> 345
329 L1etaOverTheta2 = reducedEnergy/(tetal1*te 346 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
330 // Medium & high velocity 347 // Medium & high velocity
331 // *** 1) NO RELATIVISTIC CORRECTION 348 // *** 1) NO RELATIVISTIC CORRECTION
332 // *** 2) NO sigma_PSS_l1 349 // *** 2) NO sigma_PSS_l1
333 // *** see Benka, ADANDT 22, p223 350 // *** see Benka, ADANDT 22, p223
334 351
335 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) & << 352 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
>> 353
336 universalFunction_l1 = FunctionFL1(tetal 354 universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
337 355
338 if (verboseLevel>0) G4cout << "at medium a 356 if (verboseLevel>0) G4cout << "at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
339 357
340 sigmaPSSR_l1 = (sigma0/tetal1)*universalFu 358 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
341 // *** see Benka, ADANDT 22, p220, f1 359 // *** see Benka, ADANDT 22, p220, f1
342 360
343 if (verboseLevel>0) G4cout << " sigma PWB << 361 if (verboseLevel>0) G4cout << " sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<< G4endl;
344 } 362 }
345 << 363
346 G4double pssDeltal1 = (4./(systemMass *sigma 364 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
347 // *** also called dzeta*delta 365 // *** also called dzeta*delta
348 // *** see Brandt, Phys Rev A23, p1727, f B2 366 // *** see Brandt, Phys Rev A23, p1727, f B2
349 367
350 if (verboseLevel>0) G4cout << " pssDeltal1= << 368 if (verboseLevel>0) G4cout << " pssDeltal1=" << pssDeltal1<< G4endl;
351 369
352 if (pssDeltal1>1) return 0.; 370 if (pssDeltal1>1) return 0.;
353 371
354 G4double energyLossl1 = std::pow(1-pssDeltal 372 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
355 // *** also called z 373 // *** also called z
356 // *** see Brandt, Phys Rev A23, p1727, afte 374 // *** see Brandt, Phys Rev A23, p1727, after f B2
357 375
358 if (verboseLevel>0) G4cout << " energyLossl << 376 if (verboseLevel>0) G4cout << " energyLossl1=" << energyLossl1<< G4endl;
359 377
360 G4double coulombDeflectionl1 = << 378 G4double coulombDeflectionl1 =
361 (8.*pi*zIncident/systemMass)*std::pow(tetal 379 (8.*pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
362 // *** see Brandt, Phys Rev A20, v2s and f2 << 380 // *** see Brandt, Phys Rev A20, v2s and f2 and B2
363 // *** with factor n2 compared to Brandt, Ph 381 // *** with factor n2 compared to Brandt, Phys Rev A23, p1727, f B3
364 382
365 G4double cParameterl1 =2.* coulombDeflection 383 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
366 // *** see Brandt, Phys Rev A23, p1727, f B4 384 // *** see Brandt, Phys Rev A23, p1727, f B4
367 385
368 G4double coulombDeflectionFunction_l1 = 9.*E 386 G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1); //Coulomb-deflection effect correction
369 387
370 if (verboseLevel>0) G4cout << " coulombDefl << 388 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 << G4endl;
371 389
372 G4double crossSection_L1 = coulombDeflection 390 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
373 391
374 //ECPSSR L1 -subshell cross section is estim 392 //ECPSSR L1 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
375 //and reduced by the energy-loss(E),the Coul 393 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
376 394
377 if (verboseLevel>0) G4cout << " crossSectio << 395 if (verboseLevel>0) G4cout << " crossSection_L1 =" << crossSection_L1 << G4endl;
378 396
379 if (crossSection_L1 >= 0) << 397 if (crossSection_L1 >= 0)
380 { 398 {
381 return crossSection_L1 * barn; 399 return crossSection_L1 * barn;
382 } 400 }
383 << 401
384 else {return 0;} 402 else {return 0;}
385 } 403 }
386 404
387 //....oooOO0OOooo........oooOO0OOooo........oo 405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
388 406
389 G4double G4ecpssrBaseLixsModel::CalculateL2Cro 407 G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
390 408
391 { 409 {
392 if (zTarget <=13 ) return 0.; 410 if (zTarget <=13 ) return 0.;
393 411
394 // this L2-CrossSection calculation method i 412 // this L2-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
395 // and using data tables of O. Benka et al. 413 // and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
396 414
397 G4NistManager* massManager = G4NistManager:: 415 G4NistManager* massManager = G4NistManager::Instance();
398 416
399 G4AtomicTransitionManager* transitionManage 417 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
400 418
401 G4double zIncident = 0; 419 G4double zIncident = 0;
402 << 420
403 G4Proton* aProtone = G4Proton::Proton(); 421 G4Proton* aProtone = G4Proton::Proton();
404 G4Alpha* aAlpha = G4Alpha::Alpha(); 422 G4Alpha* aAlpha = G4Alpha::Alpha();
405 423
406 if (massIncident == aProtone->GetPDGMass() ) 424 if (massIncident == aProtone->GetPDGMass() )
>> 425
407 zIncident = (aProtone->GetPDGCharge())/eplu 426 zIncident = (aProtone->GetPDGCharge())/eplus;
408 427
409 else 428 else
410 { 429 {
411 if (massIncident == aAlpha->GetPDGMass() 430 if (massIncident == aAlpha->GetPDGMass())
412 zIncident = (aAlpha->GetPDGCharge())/eplus; << 431
>> 432 zIncident = (aAlpha->GetPDGCharge())/eplus;
413 433
414 else 434 else
415 { 435 {
416 G4cout << "*** WARNING in G4ecpssrBaseLixs 436 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << G4endl;
417 G4cout << massIncident << ", " << aAlpha-> 437 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
418 return 0; 438 return 0;
419 } 439 }
420 } 440 }
421 441
422 G4double l2BindingEnergy = transitionManager 442 G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy(); //Observed binding energy of L2-subshell
423 443
424 G4double massTarget = (massManager->GetAtomi 444 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
425 445
426 G4double systemMass =((massIncident*massTarg 446 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
427 447
428 const G4double zlshell= 4.15; 448 const G4double zlshell= 4.15;
429 449
430 G4double screenedzTarget = zTarget-zlshell; 450 G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L2-subshell
431 451
432 const G4double rydbergMeV= 13.6056923e-6; 452 const G4double rydbergMeV= 13.6056923e-6;
433 453
434 const G4double nl= 2.; 454 const G4double nl= 2.;
435 455
436 G4double tetal2 = (l2BindingEnergy*nl*nl)/(( 456 G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
437 457
438 if (verboseLevel>0) G4cout << " tetal2=" << << 458 if (verboseLevel>0) G4cout << " tetal2=" << tetal2<< G4endl;
439 459
440 G4double reducedEnergy = (energyIncident*ele << 460 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
441 461
442 const G4double bohrPow2Barn=(Bohr_radius*Boh 462 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
443 463
444 G4double sigma0 = 8.*pi*(zIncident*zIncident 464 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
445 465
446 G4double velocityl2 = CalculateVelocity(2, 466 G4double velocityl2 = CalculateVelocity(2, zTarget, massIncident, energyIncident); // Scaled velocity
447 467
448 if (verboseLevel>0) G4cout << " velocityl2= << 468 if (verboseLevel>0) G4cout << " velocityl2=" << velocityl2<< G4endl;
449 469
450 const G4double l23AnalyticalApproximation= 1 470 const G4double l23AnalyticalApproximation= 1.25;
451 471
452 G4double x2 = (nl*l23AnalyticalApproximation 472 G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
453 << 473
454 if (verboseLevel>0) G4cout << " x2=" << x2< << 474 if (verboseLevel>0) G4cout << " x2=" << x2<< G4endl;
455 475
456 G4double electrIonizationEnergyl2=0.; 476 G4double electrIonizationEnergyl2=0.;
457 477
458 if ( x2<=0.035) electrIonizationEnergyl2= 0 478 if ( x2<=0.035) electrIonizationEnergyl2= 0.75*pi*(std::log(1./(x2*x2))-1.);
459 else 479 else
460 { 480 {
461 if ( x2<=3.) 481 if ( x2<=3.)
462 electrIonizationEnergyl2 =G4Exp(-2.*x2 << 482 electrIonizationEnergyl2 =std::exp(-2.*x2)/(0.031+(0.213*std::pow(x2,0.5))+(0.005*x2)-(0.069*std::pow(x2,3./2.))+(0.324*x2*x2));
463 else 483 else
464 {if ( x2<=11.) electrIonizationEnergyl << 484 {if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6); }
465 } 485 }
466 486
467 G4double hFunctionl2 =(electrIonizationEnerg 487 G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3)); //takes into account the polarization effect
468 488
469 if (verboseLevel>0) G4cout << " hFunctionl2 << 489 if (verboseLevel>0) G4cout << " hFunctionl2=" << hFunctionl2<< G4endl;
470 490
471 G4double gFunctionl2 = (1.+(10.*velocityl2)+ 491 G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*std::pow(velocityl2,3.))+(331.*std::pow(velocityl2,4.))+(6.7*std::pow(velocityl2,5.))+(58.*std::pow(velocityl2,6.))+(7.8*std::pow(velocityl2,7.))+ (0.888*std::pow(velocityl2,8.)) )/std::pow(1.+velocityl2,10.);
472 //takes into account the reduced binding eff 492 //takes into account the reduced binding effect
473 493
474 if (verboseLevel>0) G4cout << " gFunctionl2 << 494 if (verboseLevel>0) G4cout << " gFunctionl2=" << gFunctionl2<< G4endl;
475 495
476 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/( 496 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); //Binding-polarization factor
477 497
478 if (verboseLevel>0) G4cout << " sigmaPSS_l2 << 498 if (verboseLevel>0) G4cout << " sigmaPSS_l2=" << sigmaPSS_l2<< G4endl;
479 499
480 const G4double cNaturalUnit= 137.; 500 const G4double cNaturalUnit= 137.;
481 << 501
482 G4double yl2Formula=0.15*(screenedzTarget/cN 502 G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
483 << 503
484 G4double l2relativityCorrection = std::pow(( 504 G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula; // Relativistic correction parameter
485 << 505
486 G4double L2etaOverTheta2; 506 G4double L2etaOverTheta2;
487 507
488 G4double universalFunction_l2 = 0.; 508 G4double universalFunction_l2 = 0.;
489 509
490 G4double sigmaPSSR_l2 ; 510 G4double sigmaPSSR_l2 ;
491 << 511
492 if ( velocityl2 < 20. ) 512 if ( velocityl2 < 20. )
493 { 513 {
494 514
495 L2etaOverTheta2 = (reducedEnergy*l2relativ 515 L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
496 516
497 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2* 517 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
>> 518
498 universalFunction_l2 = FunctionFL2((teta 519 universalFunction_l2 = FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
499 520
500 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2 521 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
501 522
502 if (verboseLevel>0) G4cout << " sigma PWB 523 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at low velocity range = " << sigmaPSSR_l2<< G4endl;
503 } << 524 }
504 else 525 else
505 { << 526 {
506 << 527
507 L2etaOverTheta2 = reducedEnergy /(tetal2*t 528 L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
508 529
509 if ( (tetal2>=0.2) && (tetal2<=2.6670) && 530 if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
510 universalFunction_l2 = FunctionFL2((teta <<
511 531
>> 532 universalFunction_l2 = FunctionFL2((tetal2),L2etaOverTheta2);
>> 533
512 sigmaPSSR_l2 = (sigma0/tetal2)*universalFu 534 sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
513 535
514 if (verboseLevel>0) G4cout << " sigma PWB << 536 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at medium and high velocity range = " << sigmaPSSR_l2<< G4endl;
515 537
516 } 538 }
517 << 539
518 G4double pssDeltal2 = (4./(systemMass*sigmaP 540 G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
519 541
520 if (pssDeltal2>1) return 0.; 542 if (pssDeltal2>1) return 0.;
521 543
522 G4double energyLossl2 = std::pow(1-pssDeltal 544 G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
523 << 545
524 if (verboseLevel>0) G4cout << " energyLos 546 if (verboseLevel>0) G4cout << " energyLossl2=" << energyLossl2<< G4endl;
525 547
526 G4double coulombDeflectionl2 << 548 G4double coulombDeflectionl2
527 =(8.*pi*zIncident/systemMass)*std::pow(tet 549 =(8.*pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
528 550
529 G4double cParameterl2 = 2.*coulombDeflection 551 G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
530 552
531 G4double coulombDeflectionFunction_l2 = 11.* 553 G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2); //Coulomb-deflection effect correction
532 // *** see Brandt, Phys Rev A10, p477, f25 554 // *** see Brandt, Phys Rev A10, p477, f25
533 << 555
534 if (verboseLevel>0) G4cout << " coulombDefl << 556 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l2 =" << coulombDeflectionFunction_l2 << G4endl;
535 557
536 G4double crossSection_L2 = coulombDeflection 558 G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
537 //ECPSSR L2 -subshell cross section is estim 559 //ECPSSR L2 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
538 //and reduced by the energy-loss(E),the Coul 560 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
539 561
540 if (verboseLevel>0) G4cout << " crossSectio << 562 if (verboseLevel>0) G4cout << " crossSection_L2 =" << crossSection_L2 << G4endl;
541 563
542 if (crossSection_L2 >= 0) << 564 if (crossSection_L2 >= 0)
543 { 565 {
544 return crossSection_L2 * barn; 566 return crossSection_L2 * barn;
545 } 567 }
546 else {return 0;} 568 else {return 0;}
547 } 569 }
548 570
549 //....oooOO0OOooo........oooOO0OOooo........oo 571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
550 572
551 573
552 G4double G4ecpssrBaseLixsModel::CalculateL3Cro 574 G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
553 575
554 { 576 {
555 if (zTarget <=13) return 0.; 577 if (zTarget <=13) return 0.;
556 578
557 //this L3-CrossSection calculation method is 579 //this L3-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
558 //and using data tables of O. Benka et al. A 580 //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
559 581
560 G4NistManager* massManager = G4NistManager:: 582 G4NistManager* massManager = G4NistManager::Instance();
561 583
562 G4AtomicTransitionManager* transitionManage 584 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
563 585
564 G4double zIncident = 0; 586 G4double zIncident = 0;
565 587
566 G4Proton* aProtone = G4Proton::Proton(); 588 G4Proton* aProtone = G4Proton::Proton();
567 G4Alpha* aAlpha = G4Alpha::Alpha(); 589 G4Alpha* aAlpha = G4Alpha::Alpha();
568 590
569 if (massIncident == aProtone->GetPDGMass() ) 591 if (massIncident == aProtone->GetPDGMass() )
570 592
571 zIncident = (aProtone->GetPDGCharge())/eplu 593 zIncident = (aProtone->GetPDGCharge())/eplus;
572 594
573 else 595 else
574 { 596 {
575 if (massIncident == aAlpha->GetPDGMass() 597 if (massIncident == aAlpha->GetPDGMass())
576 598
577 zIncident = (aAlpha->GetPDGCharge())/eplu 599 zIncident = (aAlpha->GetPDGCharge())/eplus;
578 600
579 else 601 else
580 { 602 {
581 G4cout << "*** WARNING in G4ecpssrBaseLixs 603 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << G4endl;
582 G4cout << massIncident << ", " << aAlpha-> 604 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
583 return 0; 605 return 0;
584 } 606 }
585 } 607 }
586 608
587 G4double l3BindingEnergy = transitionManager 609 G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
588 610
589 G4double massTarget = (massManager->GetAtomi 611 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
590 612
591 G4double systemMass =((massIncident*massTarg 613 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//Mass of the system (projectile, target)
592 614
593 const G4double zlshell= 4.15; 615 const G4double zlshell= 4.15;
594 616
595 G4double screenedzTarget = zTarget-zlshell;/ 617 G4double screenedzTarget = zTarget-zlshell;//Effective nuclear charge as seen by electrons in L3-subshell
596 618
597 const G4double rydbergMeV= 13.6056923e-6; 619 const G4double rydbergMeV= 13.6056923e-6;
598 620
599 const G4double nl= 2.; 621 const G4double nl= 2.;
600 622
601 G4double tetal3 = (l3BindingEnergy*nl*nl)/(( 623 G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);//Screening parameter
602 624
603 if (verboseLevel>0) G4cout << " tetal3=" 625 if (verboseLevel>0) G4cout << " tetal3=" << tetal3<< G4endl;
604 626
605 G4double reducedEnergy = (energyIncident*ele 627 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
606 628
607 const G4double bohrPow2Barn=(Bohr_radius*Boh 629 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;//Bohr radius of hydrogen
608 630
609 G4double sigma0 = 8.*pi*(zIncident*zIncident 631 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
610 632
611 G4double velocityl3 = CalculateVelocity(3, z 633 G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);// Scaled velocity
612 634
613 if (verboseLevel>0) G4cout << " velocityl 635 if (verboseLevel>0) G4cout << " velocityl3=" << velocityl3<< G4endl;
614 636
615 const G4double l23AnalyticalApproximation= 1 637 const G4double l23AnalyticalApproximation= 1.25;
616 638
617 G4double x3 = (nl*l23AnalyticalApproximation 639 G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
618 640
619 if (verboseLevel>0) G4cout << " x3=" << x 641 if (verboseLevel>0) G4cout << " x3=" << x3<< G4endl;
620 642
621 G4double electrIonizationEnergyl3=0.; 643 G4double electrIonizationEnergyl3=0.;
622 644
623 if ( x3<=0.035) electrIonizationEnergyl3= 0 645 if ( x3<=0.035) electrIonizationEnergyl3= 0.75*pi*(std::log(1./(x3*x3))-1.);
624 else 646 else
625 { 647 {
626 if ( x3<=3.) electrIonizationEnergyl3 =G << 648 if ( x3<=3.) electrIonizationEnergyl3 =std::exp(-2.*x3)/(0.031+(0.213*std::pow(x3,0.5))+(0.005*x3)-(0.069*std::pow(x3,3./2.))+(0.324*x3*x3));
627 else 649 else
628 { 650 {
629 if ( x3<=11.) electrIonizationEnergyl3 =2. << 651 if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
630 } 652 }
631 653
632 G4double hFunctionl3 =(electrIonizationEnerg 654 G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));//takes into account the polarization effect
633 655
634 if (verboseLevel>0) G4cout << " hFunction 656 if (verboseLevel>0) G4cout << " hFunctionl3=" << hFunctionl3<< G4endl;
635 657
636 G4double gFunctionl3 = (1.+(10.*velocityl3)+ 658 G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*std::pow(velocityl3,3.))+(331.*std::pow(velocityl3,4.))+(6.7*std::pow(velocityl3,5.))+(58.*std::pow(velocityl3,6.))+(7.8*std::pow(velocityl3,7.))+ (0.888*std::pow(velocityl3,8.)) )/std::pow(1.+velocityl3,10.);
637 //takes into account the reduced binding eff 659 //takes into account the reduced binding effect
638 660
639 if (verboseLevel>0) G4cout << " gFunction 661 if (verboseLevel>0) G4cout << " gFunctionl3=" << gFunctionl3<< G4endl;
640 662
641 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/( 663 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));//Binding-polarization factor
642 664
643 if (verboseLevel>0) G4cout << "sigmaPSS_l3 665 if (verboseLevel>0) G4cout << "sigmaPSS_l3 =" << sigmaPSS_l3<< G4endl;
644 666
645 const G4double cNaturalUnit= 137.; 667 const G4double cNaturalUnit= 137.;
646 << 668
647 G4double yl3Formula=0.15*(screenedzTarget/cN 669 G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
648 << 670
649 G4double l3relativityCorrection = std::pow(( 671 G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula; // Relativistic correction parameter
650 << 672
651 G4double L3etaOverTheta2; 673 G4double L3etaOverTheta2;
652 << 674
653 G4double universalFunction_l3 = 0.; 675 G4double universalFunction_l3 = 0.;
654 << 676
655 G4double sigmaPSSR_l3; 677 G4double sigmaPSSR_l3;
656 678
657 if ( velocityl3 < 20. ) 679 if ( velocityl3 < 20. )
658 { 680 {
659 681
660 L3etaOverTheta2 = (reducedEnergy* l3relati 682 L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
661 683
662 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3* 684 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
663 685
664 universalFunction_l3 = 2.*FunctionFL2((t 686 universalFunction_l3 = 2.*FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2 );
665 687
666 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3 688 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
667 689
668 if (verboseLevel>0) G4cout << " sigma PWB 690 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at low velocity range = " << sigmaPSSR_l3<< G4endl;
669 691
670 } 692 }
671 693
672 else << 694 else
673 695
674 { 696 {
675 697
676 L3etaOverTheta2 = reducedEnergy/(tetal3*te 698 L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
677 699
678 if ( (tetal3>=0.2) && (tetal3<=2.6670) && << 700 if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
679 << 701
680 universalFunction_l3 = 2.*FunctionFL2(te 702 universalFunction_l3 = 2.*FunctionFL2(tetal3, L3etaOverTheta2 );
681 703
682 sigmaPSSR_l3 = (sigma0/tetal3)*universalFu 704 sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
683 705
684 if (verboseLevel>0) G4cout << " sigma P 706 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at medium and high velocity range = " << sigmaPSSR_l3<< G4endl;
685 } 707 }
686 << 708
687 G4double pssDeltal3 = (4./(systemMass*sigmaP 709 G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
688 710
689 if (verboseLevel>0) G4cout << " pssDeltal 711 if (verboseLevel>0) G4cout << " pssDeltal3=" << pssDeltal3<< G4endl;
690 712
691 if (pssDeltal3>1) return 0.; 713 if (pssDeltal3>1) return 0.;
692 714
693 G4double energyLossl3 = std::pow(1-pssDeltal 715 G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
694 716
695 if (verboseLevel>0) G4cout << " energyLossl 717 if (verboseLevel>0) G4cout << " energyLossl3=" << energyLossl3<< G4endl;
696 718
697 G4double coulombDeflectionl3 = << 719 G4double coulombDeflectionl3 =
698 (8.*pi*zIncident/systemMass)*std::pow(teta 720 (8.*pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
699 721
700 G4double cParameterl3 = 2.*coulombDeflection 722 G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
701 723
702 G4double coulombDeflectionFunction_l3 = 11.* 724 G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);//Coulomb-deflection effect correction
703 // *** see Brandt, Phys Rev A10, p477, f25 725 // *** see Brandt, Phys Rev A10, p477, f25
704 726
705 if (verboseLevel>0) G4cout << " coulombDe 727 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l3 =" << coulombDeflectionFunction_l3 << G4endl;
706 728
707 G4double crossSection_L3 = coulombDeflectio 729 G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;
708 //ECPSSR L3 -subshell cross section is estim 730 //ECPSSR L3 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
709 //and reduced by the energy-loss(E),the Coul 731 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
710 732
711 if (verboseLevel>0) G4cout << " crossSect 733 if (verboseLevel>0) G4cout << " crossSection_L3 =" << crossSection_L3 << G4endl;
712 << 734
713 if (crossSection_L3 >= 0) << 735 if (crossSection_L3 >= 0)
714 { 736 {
715 return crossSection_L3 * barn; 737 return crossSection_L3 * barn;
716 } 738 }
717 else {return 0;} 739 else {return 0;}
718 } 740 }
719 741
720 //....oooOO0OOooo........oooOO0OOooo........oo 742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
721 743
722 G4double G4ecpssrBaseLixsModel::CalculateVeloc 744 G4double G4ecpssrBaseLixsModel::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
723 745
724 { 746 {
725 747
726 G4AtomicTransitionManager* transitionManage 748 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
727 749
728 G4double liBindingEnergy = transitionManager 750 G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
729 751
730 G4Proton* aProtone = G4Proton::Proton(); 752 G4Proton* aProtone = G4Proton::Proton();
731 G4Alpha* aAlpha = G4Alpha::Alpha(); 753 G4Alpha* aAlpha = G4Alpha::Alpha();
732 754
733 if (!((massIncident == aProtone->GetPDGMass( 755 if (!((massIncident == aProtone->GetPDGMass()) || (massIncident == aAlpha->GetPDGMass())))
734 { 756 {
735 G4cout << "*** WARNING in G4ecpssrBaseLi 757 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << G4endl;
736 G4cout << massIncident << ", " << aAlpha 758 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
737 return 0; 759 return 0;
738 } 760 }
739 761
740 constexpr G4double zlshell= 4.15; << 762 const G4double zlshell= 4.15;
741 763
742 G4double screenedzTarget = zTarget- zlshell; 764 G4double screenedzTarget = zTarget- zlshell;
743 765
744 constexpr G4double rydbergMeV= 13.6056923e-6 << 766 const G4double rydbergMeV= 13.6056923e-6;
745 767
746 constexpr G4double nl= 2.; << 768 const G4double nl= 2.;
747 769
748 G4double tetali = (liBindingEnergy*nl*nl)/(s 770 G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
749 771
750 G4double reducedEnergy = (energyIncident*ele 772 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
751 773
752 G4double velocity = 2.*nl*std::pow(reducedEn 774 G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
753 // *** see Brandt, Phys Rev A10, p10, f4 775 // *** see Brandt, Phys Rev A10, p10, f4
754 << 776
755 return velocity; 777 return velocity;
756 } 778 }
757 779
758 //....oooOO0OOooo........oooOO0OOooo........oo 780 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
759 781
760 G4double G4ecpssrBaseLixsModel::FunctionFL1(G4 782 G4double G4ecpssrBaseLixsModel::FunctionFL1(G4double k, G4double theta)
761 { 783 {
>> 784
762 G4double sigma = 0.; 785 G4double sigma = 0.;
763 G4double valueT1 = 0; 786 G4double valueT1 = 0;
764 G4double valueT2 = 0; 787 G4double valueT2 = 0;
765 G4double valueE21 = 0; 788 G4double valueE21 = 0;
766 G4double valueE22 = 0; 789 G4double valueE22 = 0;
767 G4double valueE12 = 0; 790 G4double valueE12 = 0;
768 G4double valueE11 = 0; 791 G4double valueE11 = 0;
769 G4double xs11 = 0; 792 G4double xs11 = 0;
770 G4double xs12 = 0; 793 G4double xs12 = 0;
771 G4double xs21 = 0; 794 G4double xs21 = 0;
772 G4double xs22 = 0; 795 G4double xs22 = 0;
773 796
774 // PROTECTION TO ALLOW INTERPOLATION AT MINI 797 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
775 798
776 if ( 799 if (
777 theta==8.66e-4 || 800 theta==8.66e-4 ||
778 theta==8.66e-3 || 801 theta==8.66e-3 ||
779 theta==8.66e-2 || 802 theta==8.66e-2 ||
780 theta==8.66e-1 || 803 theta==8.66e-1 ||
781 theta==8.66e+00 || 804 theta==8.66e+00 ||
782 theta==8.66e+01 805 theta==8.66e+01
783 ) theta=theta-1e-12; 806 ) theta=theta-1e-12;
784 807
785 if ( 808 if (
786 theta==1.e-4 || 809 theta==1.e-4 ||
787 theta==1.e-3 || 810 theta==1.e-3 ||
788 theta==1.e-2 || 811 theta==1.e-2 ||
789 theta==1.e-1 || 812 theta==1.e-1 ||
790 theta==1.e+00 || 813 theta==1.e+00 ||
791 theta==1.e+01 814 theta==1.e+01
792 ) theta=theta+1e-12; 815 ) theta=theta+1e-12;
793 816
794 // END PROTECTION 817 // END PROTECTION
795 818
796 auto t2 = std::upper_bound(dummyVec1.begin() << 819 std::vector<double>::iterator t2 = std::upper_bound(dummyVec1.begin(),dummyVec1.end(), k);
797 auto t1 = t2-1; << 820 std::vector<double>::iterator t1 = t2-1;
798 821
799 auto e12 = std::upper_bound(aVecMap1[(*t1)]. << 822 std::vector<double>::iterator e12 = std::upper_bound(aVecMap1[(*t1)].begin(),aVecMap1[(*t1)].end(), theta);
800 auto e11 = e12-1; << 823 std::vector<double>::iterator e11 = e12-1;
801 824
802 auto e22 = std::upper_bound(aVecMap1[(*t2)]. << 825 std::vector<double>::iterator e22 = std::upper_bound(aVecMap1[(*t2)].begin(),aVecMap1[(*t2)].end(), theta);
803 auto e21 = e22-1; << 826 std::vector<double>::iterator e21 = e22-1;
804 827
805 valueT1 =*t1; << 828 valueT1 =*t1;
806 valueT2 =*t2; << 829 valueT2 =*t2;
807 valueE21 =*e21; << 830 valueE21 =*e21;
808 valueE22 =*e22; << 831 valueE22 =*e22;
809 valueE12 =*e12; << 832 valueE12 =*e12;
810 valueE11 =*e11; << 833 valueE11 =*e11;
811 << 834
812 xs11 = FL1Data[valueT1][valueE11]; << 835 xs11 = FL1Data[valueT1][valueE11];
813 xs12 = FL1Data[valueT1][valueE12]; << 836 xs12 = FL1Data[valueT1][valueE12];
814 xs21 = FL1Data[valueT2][valueE21]; << 837 xs21 = FL1Data[valueT2][valueE21];
815 xs22 = FL1Data[valueT2][valueE22]; << 838 xs22 = FL1Data[valueT2][valueE22];
816 839
817 if (verboseLevel>0) << 840 if (verboseLevel>0)
818 G4cout << 841 G4cout
819 << valueT1 << " " 842 << valueT1 << " "
820 << valueT2 << " " 843 << valueT2 << " "
821 << valueE11 << " " 844 << valueE11 << " "
822 << valueE12 << " " 845 << valueE12 << " "
823 << valueE21 << " " 846 << valueE21 << " "
824 << valueE22 << " " 847 << valueE22 << " "
825 << xs11 << " " 848 << xs11 << " "
826 << xs12 << " " 849 << xs12 << " "
827 << xs21 << " " 850 << xs21 << " "
828 << xs22 << " " 851 << xs22 << " "
829 << G4endl; 852 << G4endl;
830 853
831 G4double xsProduct = xs11 * xs12 * xs21 * xs 854 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
832 855
833 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 856 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
834 857
835 if (xsProduct != 0.) 858 if (xsProduct != 0.)
836 { 859 {
837 sigma = QuadInterpolator( valueE11, value 860 sigma = QuadInterpolator( valueE11, valueE12,
838 valueE21, valueE22, 861 valueE21, valueE22,
839 xs11, xs12, 862 xs11, xs12,
840 xs21, xs22, 863 xs21, xs22,
841 valueT1, valueT2, 864 valueT1, valueT2,
842 k, theta ); 865 k, theta );
843 } 866 }
844 867
845 return sigma; 868 return sigma;
846 } 869 }
847 870
848 //....oooOO0OOooo........oooOO0OOooo........oo 871 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
849 872
850 G4double G4ecpssrBaseLixsModel::FunctionFL2(G4 873 G4double G4ecpssrBaseLixsModel::FunctionFL2(G4double k, G4double theta)
851 { 874 {
852 875
853 G4double sigma = 0.; 876 G4double sigma = 0.;
854 G4double valueT1 = 0; 877 G4double valueT1 = 0;
855 G4double valueT2 = 0; 878 G4double valueT2 = 0;
856 G4double valueE21 = 0; 879 G4double valueE21 = 0;
857 G4double valueE22 = 0; 880 G4double valueE22 = 0;
858 G4double valueE12 = 0; 881 G4double valueE12 = 0;
859 G4double valueE11 = 0; 882 G4double valueE11 = 0;
860 G4double xs11 = 0; 883 G4double xs11 = 0;
861 G4double xs12 = 0; 884 G4double xs12 = 0;
862 G4double xs21 = 0; 885 G4double xs21 = 0;
863 G4double xs22 = 0; 886 G4double xs22 = 0;
864 887
865 // PROTECTION TO ALLOW INTERPOLATION AT MINI 888 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
866 889
867 if ( 890 if (
868 theta==8.66e-4 || 891 theta==8.66e-4 ||
869 theta==8.66e-3 || 892 theta==8.66e-3 ||
870 theta==8.66e-2 || 893 theta==8.66e-2 ||
871 theta==8.66e-1 || 894 theta==8.66e-1 ||
872 theta==8.66e+00 || 895 theta==8.66e+00 ||
873 theta==8.66e+01 896 theta==8.66e+01
874 ) theta=theta-1e-12; 897 ) theta=theta-1e-12;
875 898
876 if ( 899 if (
877 theta==1.e-4 || 900 theta==1.e-4 ||
878 theta==1.e-3 || 901 theta==1.e-3 ||
879 theta==1.e-2 || 902 theta==1.e-2 ||
880 theta==1.e-1 || 903 theta==1.e-1 ||
881 theta==1.e+00 || 904 theta==1.e+00 ||
882 theta==1.e+01 905 theta==1.e+01
883 ) theta=theta+1e-12; 906 ) theta=theta+1e-12;
884 907
885 // END PROTECTION 908 // END PROTECTION
886 909
887 auto t2 = std::upper_bound(dummyVec2.begin() << 910 std::vector<double>::iterator t2 = std::upper_bound(dummyVec2.begin(),dummyVec2.end(), k);
888 auto t1 = t2-1; << 911 std::vector<double>::iterator t1 = t2-1;
889 auto e12 = std::upper_bound(aVecMap2[(*t1)]. <<
890 auto e11 = e12-1; <<
891 auto e22 = std::upper_bound(aVecMap2[(*t2)]. <<
892 auto e21 = e22-1; <<
893 <<
894 valueT1 =*t1; <<
895 valueT2 =*t2; <<
896 valueE21 =*e21; <<
897 valueE22 =*e22; <<
898 valueE12 =*e12; <<
899 valueE11 =*e11; <<
900 <<
901 xs11 = FL2Data[valueT1][valueE11]; <<
902 xs12 = FL2Data[valueT1][valueE12]; <<
903 xs21 = FL2Data[valueT2][valueE21]; <<
904 xs22 = FL2Data[valueT2][valueE22]; <<
905 912
906 if (verboseLevel>0) << 913 std::vector<double>::iterator e12 = std::upper_bound(aVecMap2[(*t1)].begin(),aVecMap2[(*t1)].end(), theta);
907 G4cout << 914 std::vector<double>::iterator e11 = e12-1;
>> 915
>> 916 std::vector<double>::iterator e22 = std::upper_bound(aVecMap2[(*t2)].begin(),aVecMap2[(*t2)].end(), theta);
>> 917 std::vector<double>::iterator e21 = e22-1;
>> 918
>> 919 valueT1 =*t1;
>> 920 valueT2 =*t2;
>> 921 valueE21 =*e21;
>> 922 valueE22 =*e22;
>> 923 valueE12 =*e12;
>> 924 valueE11 =*e11;
>> 925
>> 926 xs11 = FL2Data[valueT1][valueE11];
>> 927 xs12 = FL2Data[valueT1][valueE12];
>> 928 xs21 = FL2Data[valueT2][valueE21];
>> 929 xs22 = FL2Data[valueT2][valueE22];
>> 930
>> 931 if (verboseLevel>0)
>> 932 G4cout
908 << valueT1 << " " 933 << valueT1 << " "
909 << valueT2 << " " 934 << valueT2 << " "
910 << valueE11 << " " 935 << valueE11 << " "
911 << valueE12 << " " 936 << valueE12 << " "
912 << valueE21 << " " 937 << valueE21 << " "
913 << valueE22 << " " 938 << valueE22 << " "
914 << xs11 << " " 939 << xs11 << " "
915 << xs12 << " " 940 << xs12 << " "
916 << xs21 << " " 941 << xs21 << " "
917 << xs22 << " " 942 << xs22 << " "
918 << G4endl; 943 << G4endl;
919 944
920 G4double xsProduct = xs11 * xs12 * xs21 * xs 945 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
921 946
922 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 947 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
923 948
924 if (xsProduct != 0.) 949 if (xsProduct != 0.)
925 { 950 {
926 sigma = QuadInterpolator( valueE11, value 951 sigma = QuadInterpolator( valueE11, valueE12,
927 valueE21, valueE22, 952 valueE21, valueE22,
928 xs11, xs12, 953 xs11, xs12,
929 xs21, xs22, 954 xs21, xs22,
930 valueT1, valueT2, 955 valueT1, valueT2,
931 k, theta ); 956 k, theta );
932 } 957 }
933 958
934 return sigma; 959 return sigma;
935 } 960 }
936 961
937 //....oooOO0OOooo........oooOO0OOooo........oo 962 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
938 963
939 G4double G4ecpssrBaseLixsModel::LinLinInterpol 964 G4double G4ecpssrBaseLixsModel::LinLinInterpolate(G4double e1,
940 G4double e2, 965 G4double e2,
941 G4double e, 966 G4double e,
942 G4double xs1, 967 G4double xs1,
943 G4double xs2) 968 G4double xs2)
944 { 969 {
945 G4double value = xs1 + (xs2 - xs1)*(e - e1)/ 970 G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 - e1);
946 return value; 971 return value;
947 } 972 }
948 973
949 //....oooOO0OOooo........oooOO0OOooo........oo 974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
950 975
951 G4double G4ecpssrBaseLixsModel::LinLogInterpol 976 G4double G4ecpssrBaseLixsModel::LinLogInterpolate(G4double e1,
952 G4double e2, 977 G4double e2,
953 G4double e, 978 G4double e,
954 G4double xs1, 979 G4double xs1,
955 G4double xs2) 980 G4double xs2)
956 { 981 {
957 G4double d1 = std::log(xs1); 982 G4double d1 = std::log(xs1);
958 G4double d2 = std::log(xs2); 983 G4double d2 = std::log(xs2);
959 G4double value = G4Exp(d1 + (d2 - d1)*(e - e << 984 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
960 return value; 985 return value;
961 } 986 }
962 987
963 //....oooOO0OOooo........oooOO0OOooo........oo 988 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
964 989
965 G4double G4ecpssrBaseLixsModel::LogLogInterpol 990 G4double G4ecpssrBaseLixsModel::LogLogInterpolate(G4double e1,
966 G4double e2, 991 G4double e2,
967 G4double e, 992 G4double e,
968 G4double xs1, 993 G4double xs1,
969 G4double xs2) 994 G4double xs2)
970 { 995 {
971 G4double a = (std::log10(xs2)-std::log10(xs1 996 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
972 G4double b = std::log10(xs2) - a*std::log10( 997 G4double b = std::log10(xs2) - a*std::log10(e2);
973 G4double sigma = a*std::log10(e) + b; 998 G4double sigma = a*std::log10(e) + b;
974 G4double value = (std::pow(10.,sigma)); 999 G4double value = (std::pow(10.,sigma));
975 return value; 1000 return value;
976 } 1001 }
977 1002
978 //....oooOO0OOooo........oooOO0OOooo........oo 1003 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
979 1004
980 G4double G4ecpssrBaseLixsModel::QuadInterpolat 1005 G4double G4ecpssrBaseLixsModel::QuadInterpolator(G4double e11, G4double e12,
981 G4double e21, G4double e22, 1006 G4double e21, G4double e22,
982 G4double xs11, G4double xs1 1007 G4double xs11, G4double xs12,
983 G4double xs21, G4double xs2 1008 G4double xs21, G4double xs22,
984 G4double t1, G4double t2, 1009 G4double t1, G4double t2,
985 G4double t, G4double e) 1010 G4double t, G4double e)
986 { 1011 {
987 // Log-Log << 1012 // Log-Log
988 G4double interpolatedvalue1 = LogLogInterpol 1013 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
989 G4double interpolatedvalue2 = LogLogInterpol 1014 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
990 G4double value = LogLogInterpolate(t1, t2, t 1015 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
>> 1016
>> 1017 /*
>> 1018 // Lin-Log
>> 1019 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
>> 1020 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
>> 1021 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
>> 1022 */
>> 1023
>> 1024 /*
>> 1025 // Lin-Lin
>> 1026 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
>> 1027 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
>> 1028 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
>> 1029 */
991 return value; 1030 return value;
992 1031
993 } 1032 }
>> 1033
994 1034