| 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 //....oooOO0OOooo........oooOO0OOooo........oo 27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
28 28
29 #include <cmath> 29 #include <cmath>
30 #include <iostream> 30 #include <iostream>
31 #include "G4ecpssrBaseKxsModel.hh" 31 #include "G4ecpssrBaseKxsModel.hh"
32 #include "globals.hh" 32 #include "globals.hh"
33 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh"
35 #include "G4AtomicTransitionManager.hh" 35 #include "G4AtomicTransitionManager.hh"
36 #include "G4NistManager.hh" 36 #include "G4NistManager.hh"
37 #include "G4Proton.hh" 37 #include "G4Proton.hh"
38 #include "G4Alpha.hh" 38 #include "G4Alpha.hh"
39 #include "G4SemiLogInterpolation.hh" 39 #include "G4SemiLogInterpolation.hh"
40 #include "G4Exp.hh" 40 #include "G4Exp.hh"
41 41
42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43 43
44 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel() 44 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()
45 { 45 {
46 verboseLevel=0; 46 verboseLevel=0;
47 47
48 // Storing C coefficients for high velocit 48 // Storing C coefficients for high velocity formula
49 G4String fileC1("pixe/uf/c1"); 49 G4String fileC1("pixe/uf/c1");
50 tableC1 = new G4CrossSectionDataSet(new G4 50 tableC1 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
51 51
52 G4String fileC2("pixe/uf/c2"); 52 G4String fileC2("pixe/uf/c2");
53 tableC2 = new G4CrossSectionDataSet(new G4 53 tableC2 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
54 54
55 G4String fileC3("pixe/uf/c3"); 55 G4String fileC3("pixe/uf/c3");
56 tableC3 = new G4CrossSectionDataSet(new G4 56 tableC3 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
57 57
58 // Storing FK data needed for medium veloc 58 // Storing FK data needed for medium velocities region
59 const char* path = G4FindDataDir("G4LEDATA 59 const char* path = G4FindDataDir("G4LEDATA");
60 60
61 if (!path) { 61 if (!path) {
62 G4Exception("G4ecpssrBaseKxsModel::G4ecp 62 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0006", FatalException,"G4LEDATA environment variable not set" );
63 return; 63 return;
64 } 64 }
65 65
66 std::ostringstream fileName; 66 std::ostringstream fileName;
67 fileName << path << "/pixe/uf/FK.dat"; 67 fileName << path << "/pixe/uf/FK.dat";
68 std::ifstream FK(fileName.str().c_str()); 68 std::ifstream FK(fileName.str().c_str());
69 69
70 if (!FK) 70 if (!FK)
71 G4Exception("G4ecpssrBaseKxsModel::G4ecp 71 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0003", FatalException,"error opening FK data file" );
72 72
73 dummyVec.push_back(0.); 73 dummyVec.push_back(0.);
74 74
75 while(!FK.eof()) 75 while(!FK.eof())
76 { 76 {
77 double x; 77 double x;
78 double y; 78 double y;
79 79
80 FK>>x>>y; 80 FK>>x>>y;
81 81
82 // Mandatory vector initialization 82 // Mandatory vector initialization
83 if (x != dummyVec.back()) 83 if (x != dummyVec.back())
84 { 84 {
85 dummyVec.push_back(x); 85 dummyVec.push_back(x);
86 aVecMap[x].push_back(-1.); 86 aVecMap[x].push_back(-1.);
87 } 87 }
88 88
89 FK>>FKData[x][y]; 89 FK>>FKData[x][y];
90 90
91 if (y != aVecMap[x].back()) aVecMap[x] 91 if (y != aVecMap[x].back()) aVecMap[x].push_back(y);
92 92
93 } 93 }
94 94
95 tableC1->LoadData(fileC1); 95 tableC1->LoadData(fileC1);
96 tableC2->LoadData(fileC2); 96 tableC2->LoadData(fileC2);
97 tableC3->LoadData(fileC3); 97 tableC3->LoadData(fileC3);
98 98
99 } 99 }
100 100
101 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 102
103 void print (G4double elem) 103 void print (G4double elem)
104 { 104 {
105 G4cout << elem << " "; 105 G4cout << elem << " ";
106 } 106 }
107 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 108
109 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel() 109 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel()
110 { 110 {
111 delete tableC1; 111 delete tableC1;
112 delete tableC2; 112 delete tableC2;
113 delete tableC3; 113 delete tableC3;
114 } 114 }
115 115
116 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 117
118 G4double G4ecpssrBaseKxsModel::ExpIntFunction( 118 G4double G4ecpssrBaseKxsModel::ExpIntFunction(G4int n,G4double x)
119 119
120 { 120 {
121 // this "ExpIntFunction" function allows fast 121 // this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)
122 G4int i; 122 G4int i;
123 G4int ii; 123 G4int ii;
124 G4int nm1; 124 G4int nm1;
125 G4double a; 125 G4double a;
126 G4double b; 126 G4double b;
127 G4double c; 127 G4double c;
128 G4double d; 128 G4double d;
129 G4double del; 129 G4double del;
130 G4double fact; 130 G4double fact;
131 G4double h; 131 G4double h;
132 G4double psi; 132 G4double psi;
133 G4double ans = 0; 133 G4double ans = 0;
134 const G4double euler= 0.5772156649; 134 const G4double euler= 0.5772156649;
135 const G4int maxit= 100; 135 const G4int maxit= 100;
136 const G4double fpmin = 1.0e-30; 136 const G4double fpmin = 1.0e-30;
137 const G4double eps = 1.0e-7; 137 const G4double eps = 1.0e-7;
138 nm1=n-1; 138 nm1=n-1;
139 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1 139 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1))) {
140 G4cout << "*** WARNING in G4ecpssrBaseKxsM 140 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::ExpIntFunction: bad arguments in ExpIntFunction" << G4endl;
141 G4cout << n << ", " << x << G4endl; 141 G4cout << n << ", " << x << G4endl;
142 } 142 }
143 else { 143 else {
144 if (n==0) ans=G4Exp(-x)/x; 144 if (n==0) ans=G4Exp(-x)/x;
145 else { 145 else {
146 if (x==0.0) ans=1.0/nm1; 146 if (x==0.0) ans=1.0/nm1;
147 else { 147 else {
148 if (x > 1.0) { 148 if (x > 1.0) {
149 b=x+n; 149 b=x+n;
150 c=1.0/fpmin; 150 c=1.0/fpmin;
151 d=1.0/b; 151 d=1.0/b;
152 h=d; 152 h=d;
153 for (i=1;i<=maxit;i++) { 153 for (i=1;i<=maxit;i++) {
154 a=-i*(nm1+i); 154 a=-i*(nm1+i);
155 b +=2.0; 155 b +=2.0;
156 d=1.0/(a*d+b); 156 d=1.0/(a*d+b);
157 c=b+a/c; 157 c=b+a/c;
158 del=c*d; 158 del=c*d;
159 h *=del; 159 h *=del;
160 if (std::fabs(del-1.0) < eps) { 160 if (std::fabs(del-1.0) < eps) {
161 ans=h*G4Exp(-x); 161 ans=h*G4Exp(-x);
162 return ans; 162 return ans;
163 } 163 }
164 } 164 }
165 } else { 165 } else {
166 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x 166 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
167 fact=1.0; 167 fact=1.0;
168 for (i=1;i<=maxit;i++) { 168 for (i=1;i<=maxit;i++) {
169 fact *=-x/i; 169 fact *=-x/i;
170 if (i !=nm1) del = -fact/(i-nm1); 170 if (i !=nm1) del = -fact/(i-nm1);
171 else { 171 else {
172 psi = -euler; 172 psi = -euler;
173 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii; 173 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
174 del=fact*(-std::log(x)+psi); 174 del=fact*(-std::log(x)+psi);
175 } 175 }
176 ans += del; 176 ans += del;
177 if (std::fabs(del) < std::fabs(ans) 177 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
178 } 178 }
179 } 179 }
180 } 180 }
181 } 181 }
182 } 182 }
183 return ans; 183 return ans;
184 } 184 }
185 185
186 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 187
188 G4double G4ecpssrBaseKxsModel::CalculateCrossS 188 G4double G4ecpssrBaseKxsModel::CalculateCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
189 189
190 { 190 {
191 // this K-CrossSection calculation method is 191 // this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//
192 G4NistManager* massManager = G4NistManager:: 192 G4NistManager* massManager = G4NistManager::Instance();
193 193
194 G4AtomicTransitionManager* transitionManage 194 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
195 195
196 G4double zIncident = 0; 196 G4double zIncident = 0;
197 G4Proton* aProtone = G4Proton::Proton(); 197 G4Proton* aProtone = G4Proton::Proton();
198 G4Alpha* aAlpha = G4Alpha::Alpha(); 198 G4Alpha* aAlpha = G4Alpha::Alpha();
199 199
200 if (massIncident == aProtone->GetPDGMass() ) 200 if (massIncident == aProtone->GetPDGMass() )
201 { 201 {
202 zIncident = (aProtone->GetPDGCharge())/eplu 202 zIncident = (aProtone->GetPDGCharge())/eplus;
203 } 203 }
204 else 204 else
205 { 205 {
206 if (massIncident == aAlpha->GetPDGMass() 206 if (massIncident == aAlpha->GetPDGMass())
207 { 207 {
208 zIncident = (aAlpha->GetPDGCharge())/eplu 208 zIncident = (aAlpha->GetPDGCharge())/eplus;
209 } 209 }
210 else 210 else
211 { 211 {
212 G4cout << "*** WARNING in G4ecpssrBaseKxsM 212 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " << G4endl;
213 return 0; 213 return 0;
214 } 214 }
215 } 215 }
216 216
217 if (verboseLevel>0) G4cout << " massIncid 217 if (verboseLevel>0) G4cout << " massIncident=" << massIncident<< G4endl;
218 218
219 G4double kBindingEnergy = transitionManager- 219 G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy();
220 220
221 if (verboseLevel>0) G4cout << " kBindingE 221 if (verboseLevel>0) G4cout << " kBindingEnergy=" << kBindingEnergy/eV<< G4endl;
222 222
223 G4double massTarget = (massManager->GetAtomi 223 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
224 224
225 if (verboseLevel>0) G4cout << " massTarge 225 if (verboseLevel>0) G4cout << " massTarget=" << massTarget<< G4endl;
226 226
227 G4double systemMass =((massIncident*massTarg 227 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //the mass of the system (projectile, target)
228 228
229 if (verboseLevel>0) G4cout << " systemMas 229 if (verboseLevel>0) G4cout << " systemMass=" << systemMass<< G4endl;
230 230
231 constexpr G4double zkshell= 0.3; 231 constexpr G4double zkshell= 0.3;
232 // *** see Brandt, Phys Rev A23, p 1727 232 // *** see Brandt, Phys Rev A23, p 1727
233 233
234 G4double screenedzTarget = zTarget-zkshell; 234 G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target
235 // *** see Brandt, Phys Rev A23, p 1727 235 // *** see Brandt, Phys Rev A23, p 1727
236 236
237 constexpr G4double rydbergMeV= 13.6056923e-6 237 constexpr G4double rydbergMeV= 13.6056923e-6;
238 238
239 G4double tetaK = kBindingEnergy/((screenedzT 239 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV); //tetaK denotes the reduced binding energy of the electron
240 // *** see Rice, ADANDT 20, p 504, f 2 240 // *** see Rice, ADANDT 20, p 504, f 2
241 241
242 if (verboseLevel>0) G4cout << " tetaK=" < 242 if (verboseLevel>0) G4cout << " tetaK=" << tetaK<< G4endl;
243 243
244 G4double velocity =(2./(tetaK*screenedzTarge 244 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
245 // *** also called xiK 245 // *** also called xiK
246 // *** see Brandt, Phys Rev A23, p 1727 246 // *** see Brandt, Phys Rev A23, p 1727
247 // *** see Basbas, Phys Rev A17, p 1656, f4 247 // *** see Basbas, Phys Rev A17, p 1656, f4
248 248
249 if (verboseLevel>0) G4cout << " velocity= 249 if (verboseLevel>0) G4cout << " velocity=" << velocity<< G4endl;
250 250
251 const G4double bohrPow2Barn=(Bohr_radius*Boh 251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
252 252
253 if (verboseLevel>0) G4cout << " bohrPow2B 253 if (verboseLevel>0) G4cout << " bohrPow2Barn=" << bohrPow2Barn<< G4endl;
254 254
255 G4double sigma0 = 8.*pi*(zIncident*zIncident 255 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.); //sigma0 is the initial cross section of K shell at stable state
256 // *** see Benka, ADANDT 22, p 220, f2, for 256 // *** see Benka, ADANDT 22, p 220, f2, for protons
257 // *** see Basbas, Phys Rev A7, p 1000 257 // *** see Basbas, Phys Rev A7, p 1000
258 258
259 if (verboseLevel>0) G4cout << " sigma0=" << 259 if (verboseLevel>0) G4cout << " sigma0=" << sigma0<< G4endl;
260 260
261 const G4double kAnalyticalApproximation= 1.5 261 const G4double kAnalyticalApproximation= 1.5;
262 G4double x = kAnalyticalApproximation/veloci 262 G4double x = kAnalyticalApproximation/velocity;
263 // *** see Brandt, Phys Rev A23, p 1727 263 // *** see Brandt, Phys Rev A23, p 1727
264 // *** see Brandt, Phys Rev A20, p 469, f16 264 // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h
265 265
266 if (verboseLevel>0) G4cout << " x=" << x< 266 if (verboseLevel>0) G4cout << " x=" << x<< G4endl;
267 267
268 G4double electrIonizationEnergy; 268 G4double electrIonizationEnergy;
269 // *** see Basbas, Phys Rev A17, p1665, f27 269 // *** see Basbas, Phys Rev A17, p1665, f27
270 // *** see Brandt, Phys Rev A20, p469 270 // *** see Brandt, Phys Rev A20, p469
271 // *** see Liu, Comp Phys Comm 97, p325, f A 271 // *** see Liu, Comp Phys Comm 97, p325, f A5
272 272
273 if ((0.< x) && (x <= 0.035)) 273 if ((0.< x) && (x <= 0.035))
274 { 274 {
275 electrIonizationEnergy= 0.75*pi*(std::lo 275 electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.);
276 } 276 }
277 else 277 else
278 { 278 {
279 if ( (0.035 < x) && (x <=3.)) 279 if ( (0.035 < x) && (x <=3.))
280 { 280 {
281 electrIonizationEnergy =G4Exp(-2.*x)/(0.03 281 electrIonizationEnergy =G4Exp(-2.*x)/(0.031+(0.213*std::pow(x,0.5))+(0.005*x)-(0.069*std::pow(x,3./2.))+(0.324*x*x));
282 } 282 }
283 283
284 else 284 else
285 { 285 {
286 if ( (3.< x) && (x<=11.)) 286 if ( (3.< x) && (x<=11.))
287 { 287 {
288 electrIonizationEnergy =2.*G4Exp(-2.*x) 288 electrIonizationEnergy =2.*G4Exp(-2.*x)/std::pow(x,1.6);
289 } 289 }
290 290
291 else electrIonizationEnergy =0.; 291 else electrIonizationEnergy =0.;
292 } 292 }
293 } 293 }
294 294
295 if (verboseLevel>0) G4cout << " electrIon 295 if (verboseLevel>0) G4cout << " electrIonizationEnergy=" << electrIonizationEnergy<< G4endl;
296 296
297 G4double hFunction =(electrIonizationEnergy* 297 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet
298 // *** see Brandt, Phys Rev A20, p 469, f16 298 // *** see Brandt, Phys Rev A20, p 469, f16
299 299
300 if (verboseLevel>0) G4cout << " hFunction 300 if (verboseLevel>0) G4cout << " hFunction=" << hFunction<< G4endl;
301 301
302 G4double gFunction = (1.+(9.*velocity)+(31.* 302 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
303 +(4.2*std::pow(velocity,6.))+(0.515*std: 303 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet
304 // *** see Brandt, Phys Rev A20, p 469, f19 304 // *** see Brandt, Phys Rev A20, p 469, f19
305 305
306 if (verboseLevel>0) G4cout << " gFunction 306 if (verboseLevel>0) G4cout << " gFunction=" << gFunction<< G4endl;
307 307
308 //------------------------------------------ 308 //-----------------------------------------------------------------------------------------------------------------------------
309 309
310 G4double sigmaPSS = 1.+(((2.*zIncident)/(scr 310 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon
311 // *** also called dzeta 311 // *** also called dzeta
312 // *** also called epsilon 312 // *** also called epsilon
313 // *** see Basbas, Phys Rev A17, p1667, f45 313 // *** see Basbas, Phys Rev A17, p1667, f45
314 314
315 if (verboseLevel>0) G4cout << " sigmaPSS= 315 if (verboseLevel>0) G4cout << " sigmaPSS=" << sigmaPSS<< G4endl;
316 316
317 if (verboseLevel>0) G4cout << " sigmaPSS* 317 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK<< G4endl;
318 318
319 //------------------------------------------ 319 //----------------------------------------------------------------------------------------------------------------------------
320 320
321 const G4double cNaturalUnit= 1/fine_structur 321 const G4double cNaturalUnit= 1/fine_structure_const; // it's the speed of light according to Atomic-Unit-System
322 322
323 if (verboseLevel>0) G4cout << " cNaturalU 323 if (verboseLevel>0) G4cout << " cNaturalUnit=" << cNaturalUnit<< G4endl;
324 324
325 G4double ykFormula=0.4*(screenedzTarget/cNat 325 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
326 // *** also called yS 326 // *** also called yS
327 // *** see Brandt, Phys Rev A20, p467, f6 327 // *** see Brandt, Phys Rev A20, p467, f6
328 // *** see Brandt, Phys Rev A23, p1728 328 // *** see Brandt, Phys Rev A23, p1728
329 329
330 if (verboseLevel>0) G4cout << " ykFormula 330 if (verboseLevel>0) G4cout << " ykFormula=" << ykFormula<< G4endl;
331 331
332 G4double relativityCorrection = std::pow((1. 332 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter
333 // *** also called mRS 333 // *** also called mRS
334 // *** see Brandt, Phys Rev A20, p467, f6 334 // *** see Brandt, Phys Rev A20, p467, f6
335 335
336 if (verboseLevel>0) G4cout << " relativit 336 if (verboseLevel>0) G4cout << " relativityCorrection=" << relativityCorrection<< G4endl;
337 337
338 G4double reducedVelocity = velocity*std::pow 338 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5); // presents the reduced collision velocity parameter
339 // *** also called xiR 339 // *** also called xiR
340 // *** see Brandt, Phys Rev A20, p468, f7 340 // *** see Brandt, Phys Rev A20, p468, f7
341 // *** see Brandt, Phys Rev A23, p1728 341 // *** see Brandt, Phys Rev A23, p1728
342 342
343 if (verboseLevel>0) G4cout << " reducedVe 343 if (verboseLevel>0) G4cout << " reducedVelocity=" << reducedVelocity<< G4endl;
344 344
345 G4double etaOverTheta2 = (energyIncident*ele 345 G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
346 /(sigmaPSS*tetaK)/( 346 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
347 // *** see Benka, ADANDT 22, p220, f4 for et 347 // *** see Benka, ADANDT 22, p220, f4 for eta
348 // then we use sigmaPSS*tetaK == epsilon*tet 348 // then we use sigmaPSS*tetaK == epsilon*tetaK
349 349
350 if (verboseLevel>0) G4cout << " etaOverTh 350 if (verboseLevel>0) G4cout << " etaOverTheta2=" << etaOverTheta2<< G4endl;
351 351
352 G4double universalFunction = 0; 352 G4double universalFunction = 0;
353 353
354 // low velocity formula 354 // low velocity formula
355 // ***************** 355 // *****************
356 if ( velocity < 1. ) 356 if ( velocity < 1. )
357 // OR 357 // OR
358 //if ( reducedVelocity/sigmaPSS < 1.) 358 //if ( reducedVelocity/sigmaPSS < 1.)
359 // *** see Brandt, Phys Rev A23, p1727 359 // *** see Brandt, Phys Rev A23, p1727
360 // *** reducedVelocity/sigmaPSS is also call 360 // *** reducedVelocity/sigmaPSS is also called xiR/dzeta
361 // ***************** 361 // *****************
362 { 362 {
363 if (verboseLevel>0) G4cout << " Notice 363 if (verboseLevel>0) G4cout << " Notice : FK is computed from low velocity formula" << G4endl;
364 364
365 universalFunction = (std::pow(2.,9.)/45. 365 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);// is the reduced universal cross section
366 // *** see Brandt, Phys Rev A23, p1728 366 // *** see Brandt, Phys Rev A23, p1728
367 367
368 if (verboseLevel>0) G4cout << " univers 368 if (verboseLevel>0) G4cout << " universalFunction by Brandt 1981 =" << universalFunction<< G4endl;
369 369
370 } 370 }
371 else 371 else
372 { 372 {
373 373
374 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tet 374 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 )
375 { 375 {
376 // High and medium energies. Method from 376 // High and medium energies. Method from Rice ADANDT 20, p506, 1977 on tables from Benka 1978
377 377
378 if (verboseLevel>0) G4cout << " Notice 378 if (verboseLevel>0) G4cout << " Notice : FK is computed from high velocity formula" << G4endl;
379 379
380 if (verboseLevel>0) G4cout << " sigmaPS 380 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK << G4endl;
381 381
382 G4double C1= tableC1->FindValue(sigmaPSS 382 G4double C1= tableC1->FindValue(sigmaPSS*tetaK);
383 G4double C2= tableC2->FindValue(sigmaPSS 383 G4double C2= tableC2->FindValue(sigmaPSS*tetaK);
384 G4double C3= tableC3->FindValue(sigmaPSS 384 G4double C3= tableC3->FindValue(sigmaPSS*tetaK);
385 385
386 if (verboseLevel>0) G4cout << " C1=" 386 if (verboseLevel>0) G4cout << " C1=" << C1 << G4endl;
387 if (verboseLevel>0) G4cout << " C2=" 387 if (verboseLevel>0) G4cout << " C2=" << C2 << G4endl;
388 if (verboseLevel>0) G4cout << " C3=" 388 if (verboseLevel>0) G4cout << " C3=" << C3 << G4endl;
389 389
390 G4double etaK = (energyIncident*electron 390 G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
391 // *** see Benka, ADANDT 22, p220, f4 fo 391 // *** see Benka, ADANDT 22, p220, f4 for eta
392 392
393 if (verboseLevel>0) G4cout << " etaK= 393 if (verboseLevel>0) G4cout << " etaK=" << etaK << G4endl;
394 394
395 G4double etaT = (sigmaPSS*tetaK)*(sigmaP 395 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6); // at any theta, the largest tabulated etaOverTheta2 is 86.6
396 // *** see Rice, ADANDT 20, p506 396 // *** see Rice, ADANDT 20, p506
397 397
398 if (verboseLevel>0) G4cout << " etaT= 398 if (verboseLevel>0) G4cout << " etaT=" << etaT << G4endl;
399 399
400 G4double fKT = FunctionFK((sigmaPSS*teta 400 G4double fKT = FunctionFK((sigmaPSS*tetaK),86.6)*(etaT/(sigmaPSS*tetaK));
401 // *** see Rice, ADANDT 20, p506 401 // *** see Rice, ADANDT 20, p506
402 402
403 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0 403 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0.)
404 { 404 {
405 G4cout << 405 G4cout <<
406 "*** WARNING in G4ecpssrBaseKxsModel:: 406 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" << G4endl;
407 return 0; 407 return 0;
408 } 408 }
409 409
410 if (verboseLevel>0) G4cout << " Funct 410 if (verboseLevel>0) G4cout << " FunctionFK=" << FunctionFK((sigmaPSS*tetaK),86.6) << G4endl;
411 411
412 if (verboseLevel>0) G4cout << " fKT=" 412 if (verboseLevel>0) G4cout << " fKT=" << fKT << G4endl;
413 413
414 G4double GK = C2/(4*etaK) + C3/(32*etaK* 414 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK);
415 415
416 if (verboseLevel>0) G4cout << " GK=" << GK 416 if (verboseLevel>0) G4cout << " GK=" << GK << G4endl;
417 417
418 G4double GT = C2/(4*etaT) + C3/(32*etaT* 418 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT);
419 419
420 if (verboseLevel>0) G4cout << " GT=" 420 if (verboseLevel>0) G4cout << " GT=" << GT << G4endl;
421 421
422 G4double DT = fKT - C1*std::log(etaT) + 422 G4double DT = fKT - C1*std::log(etaT) + GT;
423 423
424 if (verboseLevel>0) G4cout << " DT=" 424 if (verboseLevel>0) G4cout << " DT=" << DT << G4endl;
425 425
426 G4double fKK = C1*std::log(etaK) + DT - 426 G4double fKK = C1*std::log(etaK) + DT - GK;
427 427
428 if (verboseLevel>0) G4cout << " fKK=" 428 if (verboseLevel>0) G4cout << " fKK=" << fKK << G4endl;
429 429
430 G4double universalFunction3= fKK/(etaK/t 430 G4double universalFunction3= fKK/(etaK/tetaK);
431 // *** see Rice, ADANDT 20, p505, f7 431 // *** see Rice, ADANDT 20, p505, f7
432 432
433 if (verboseLevel>0) G4cout << " unive 433 if (verboseLevel>0) G4cout << " universalFunction3=" << universalFunction3 << G4endl;
434 434
435 universalFunction=universalFunction3; 435 universalFunction=universalFunction3;
436 436
437 } 437 }
438 else if ( etaOverTheta2 >= 1.e-3 && etaOve 438 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 )
439 { 439 {
440 // From Benka 1978 440 // From Benka 1978
441 441
442 if (verboseLevel>0) G4cout << " Notice 442 if (verboseLevel>0) G4cout << " Notice : FK is computed from INTERPOLATED data" << G4endl;
443 443
444 G4double universalFunction2 = FunctionFK 444 G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2);
445 445
446 if (universalFunction2<=0) 446 if (universalFunction2<=0)
447 { 447 {
448 G4cout << 448 G4cout <<
449 "*** WARNING : G4ecpssrBaseKxsModel::C 449 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" << G4endl;
450 return 0; 450 return 0;
451 } 451 }
452 452
453 if (verboseLevel>0) G4cout << " univers 453 if (verboseLevel>0) G4cout << " universalFunction2=" << universalFunction2 << " for theta=" << sigmaPSS*tetaK << " and etaOverTheta2=" << etaOverTheta2 << G4endl;
454 454
455 universalFunction=universalFunction2; 455 universalFunction=universalFunction2;
456 } 456 }
457 457
458 } 458 }
459 459
460 //------------------------------------------ 460 //----------------------------------------------------------------------------------------------------------------------
461 461
462 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK 462 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section
463 // *** see Benka, ADANDT 22, p220, f1 463 // *** see Benka, ADANDT 22, p220, f1
464 464
465 if (verboseLevel>0) G4cout << " sigmaPSSR 465 if (verboseLevel>0) G4cout << " sigmaPSSR=" << sigmaPSSR<< G4endl;
466 466
467 //------------------------------------------ 467 //-----------------------------------------------------------------------------------------------------------------------
468 468
469 G4double pssDeltaK = (4./(systemMass*sigmaPS 469 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
470 // *** also called dzetaK*deltaK 470 // *** also called dzetaK*deltaK
471 // *** see Brandt, Phys Rev A23, p1727, f B2 471 // *** see Brandt, Phys Rev A23, p1727, f B2
472 472
473 if (verboseLevel>0) G4cout << " pssDeltaK 473 if (verboseLevel>0) G4cout << " pssDeltaK=" << pssDeltaK<< G4endl;
474 474
475 if (pssDeltaK>1) return 0.; 475 if (pssDeltaK>1) return 0.;
476 476
477 G4double energyLoss = std::pow(1-pssDeltaK,0 477 G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss
478 // *** also called zK 478 // *** also called zK
479 // *** see Brandt, Phys Rev A23, p1727, afte 479 // *** see Brandt, Phys Rev A23, p1727, after f B2
480 480
481 if (verboseLevel>0) G4cout << " energyLos 481 if (verboseLevel>0) G4cout << " energyLoss=" << energyLoss<< G4endl;
482 482
483 G4double energyLossFunction = (std::pow(2.,- 483 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));//energy loss function
484 // *** also called fs 484 // *** also called fs
485 // *** see Brandt, Phys Rev A23, p1718, f7 485 // *** see Brandt, Phys Rev A23, p1718, f7
486 486
487 if (verboseLevel>0) G4cout << " energyLos 487 if (verboseLevel>0) G4cout << " energyLossFunction=" << energyLossFunction<< G4endl;
488 488
489 //------------------------------------------ 489 //----------------------------------------------------------------------------------------------------------------------------------------------
490 490
491 G4double coulombDeflection = (4.*pi*zInciden 491 G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
492 // *** see Brandt, Phys Rev A23, p1727, f B3 492 // *** see Brandt, Phys Rev A23, p1727, f B3
493 493
494 if (verboseLevel>0) G4cout << " cParamete 494 if (verboseLevel>0) G4cout << " cParameter-short=" << coulombDeflection<< G4endl;
495 495
496 G4double cParameter = 2.*coulombDeflection/( 496 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
497 // *** see Brandt, Phys Rev A23, p1727, f B4 497 // *** see Brandt, Phys Rev A23, p1727, f B4
498 498
499 if (verboseLevel>0) G4cout << " cParamete 499 if (verboseLevel>0) G4cout << " cParameter-full=" << cParameter<< G4endl;
500 500
501 G4double coulombDeflectionFunction = 9.*ExpI 501 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter); //this function describes Coulomb-deflection effect
502 // *** see Brandt, Phys Rev A23, p1727 502 // *** see Brandt, Phys Rev A23, p1727
503 503
504 if (verboseLevel>0) G4cout << " ExpIntFun 504 if (verboseLevel>0) G4cout << " ExpIntFunction(10,cParameter) =" << ExpIntFunction(10,cParameter) << G4endl;
505 505
506 if (verboseLevel>0) G4cout << " coulombDe 506 if (verboseLevel>0) G4cout << " coulombDeflectionFunction =" << coulombDeflectionFunction << G4endl;
507 507
508 //------------------------------------------ 508 //--------------------------------------------------------------------------------------------------------------------------------------------------
509 509
510 G4double crossSection = 0; 510 G4double crossSection = 0;
511 511
512 crossSection = energyLossFunction* coulombDe 512 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR; //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
513 513 //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
514 514 //and the relativity(R) effects
515 515
516 //------------------------------------------ 516 //--------------------------------------------------------------------------------------------------------------------------------------------------
517 517
518 if (crossSection >= 0) { 518 if (crossSection >= 0) {
519 return crossSection * barn; 519 return crossSection * barn;
520 } 520 }
521 else {return 0;} 521 else {return 0;}
522 522
523 } 523 }
524 524
525 //....oooOO0OOooo........oooOO0OOooo........oo 525 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526 526
527 G4double G4ecpssrBaseKxsModel::FunctionFK(G4do 527 G4double G4ecpssrBaseKxsModel::FunctionFK(G4double k, G4double theta)
528 { 528 {
529 529
530 G4double sigma = 0.; 530 G4double sigma = 0.;
531 G4double valueT1 = 0; 531 G4double valueT1 = 0;
532 G4double valueT2 = 0; 532 G4double valueT2 = 0;
533 G4double valueE21 = 0; 533 G4double valueE21 = 0;
534 G4double valueE22 = 0; 534 G4double valueE22 = 0;
535 G4double valueE12 = 0; 535 G4double valueE12 = 0;
536 G4double valueE11 = 0; 536 G4double valueE11 = 0;
537 G4double xs11 = 0; 537 G4double xs11 = 0;
538 G4double xs12 = 0; 538 G4double xs12 = 0;
539 G4double xs21 = 0; 539 G4double xs21 = 0;
540 G4double xs22 = 0; 540 G4double xs22 = 0;
541 541
542 // PROTECTION TO ALLOW INTERPOLATION AT MINI 542 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values
543 // (in particular for FK computation at 8.66 543 // (in particular for FK computation at 8.66EXX for high velocity formula)
544 544
545 if ( 545 if (
546 theta==8.66e-3 || 546 theta==8.66e-3 ||
547 theta==8.66e-2 || 547 theta==8.66e-2 ||
548 theta==8.66e-1 || 548 theta==8.66e-1 ||
549 theta==8.66e+0 || 549 theta==8.66e+0 ||
550 theta==8.66e+1 550 theta==8.66e+1
551 ) theta=theta-1e-12; 551 ) theta=theta-1e-12;
552 552
553 if ( 553 if (
554 theta==1.e-3 || 554 theta==1.e-3 ||
555 theta==1.e-2 || 555 theta==1.e-2 ||
556 theta==1.e-1 || 556 theta==1.e-1 ||
557 theta==1.e+00 || 557 theta==1.e+00 ||
558 theta==1.e+01 558 theta==1.e+01
559 ) theta=theta+1e-12; 559 ) theta=theta+1e-12;
560 560
561 // END PROTECTION 561 // END PROTECTION
562 562
563 auto t2 = std::upper_bound(dummyVec.begin(), 563 auto t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
564 auto t1 = t2-1; 564 auto t1 = t2-1;
565 565
566 auto e12 = std::upper_bound(aVecMap[(*t1)].b 566 auto e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
567 auto e11 = e12-1; 567 auto e11 = e12-1;
568 568
569 auto e22 = std::upper_bound(aVecMap[(*t2)].b 569 auto e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
570 auto e21 = e22-1; 570 auto e21 = e22-1;
571 571
572 valueT1 =*t1; 572 valueT1 =*t1;
573 valueT2 =*t2; 573 valueT2 =*t2;
574 valueE21 =*e21; 574 valueE21 =*e21;
575 valueE22 =*e22; 575 valueE22 =*e22;
576 valueE12 =*e12; 576 valueE12 =*e12;
577 valueE11 =*e11; 577 valueE11 =*e11;
578 578
579 xs11 = FKData[valueT1][valueE11]; 579 xs11 = FKData[valueT1][valueE11];
580 xs12 = FKData[valueT1][valueE12]; 580 xs12 = FKData[valueT1][valueE12];
581 xs21 = FKData[valueT2][valueE21]; 581 xs21 = FKData[valueT2][valueE21];
582 xs22 = FKData[valueT2][valueE22]; 582 xs22 = FKData[valueT2][valueE22];
583 583
584 G4double xsProduct = xs11 * xs12 * xs21 * xs 584 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
585 585
586 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 586 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
587 587
588 if (xsProduct != 0.) 588 if (xsProduct != 0.)
589 { 589 {
590 sigma = QuadInterpolator( valueE11, value 590 sigma = QuadInterpolator( valueE11, valueE12,
591 valueE21, valueE22, 591 valueE21, valueE22,
592 xs11, xs12, 592 xs11, xs12,
593 xs21, xs22, 593 xs21, xs22,
594 valueT1, valueT2, 594 valueT1, valueT2,
595 k, theta ); 595 k, theta );
596 } 596 }
597 597
598 return sigma; 598 return sigma;
599 } 599 }
600 600
601 //....oooOO0OOooo........oooOO0OOooo........oo 601 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
602 602
603 G4double G4ecpssrBaseKxsModel::LinLogInterpola 603 G4double G4ecpssrBaseKxsModel::LinLogInterpolate(G4double e1,
604 G4double e2, 604 G4double e2,
605 G4double e, 605 G4double e,
606 G4double xs1, 606 G4double xs1,
607 G4double xs2) 607 G4double xs2)
608 { 608 {
609 G4double d1 = std::log(xs1); 609 G4double d1 = std::log(xs1);
610 G4double d2 = std::log(xs2); 610 G4double d2 = std::log(xs2);
611 G4double value = G4Exp(d1 + (d2 - d1)*(e - e 611 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
612 return value; 612 return value;
613 } 613 }
614 614
615 //....oooOO0OOooo........oooOO0OOooo........oo 615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
616 616
617 G4double G4ecpssrBaseKxsModel::LogLogInterpola 617 G4double G4ecpssrBaseKxsModel::LogLogInterpolate(G4double e1,
618 G4double e2, 618 G4double e2,
619 G4double e, 619 G4double e,
620 G4double xs1, 620 G4double xs1,
621 G4double xs2) 621 G4double xs2)
622 { 622 {
623 G4double a = (std::log10(xs2)-std::log10(xs1 623 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
624 G4double b = std::log10(xs2) - a*std::log10( 624 G4double b = std::log10(xs2) - a*std::log10(e2);
625 G4double sigma = a*std::log10(e) + b; 625 G4double sigma = a*std::log10(e) + b;
626 G4double value = (std::pow(10.,sigma)); 626 G4double value = (std::pow(10.,sigma));
627 return value; 627 return value;
628 } 628 }
629 629
630 //....oooOO0OOooo........oooOO0OOooo........oo 630 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
631 631
632 G4double G4ecpssrBaseKxsModel::QuadInterpolato 632 G4double G4ecpssrBaseKxsModel::QuadInterpolator(G4double e11, G4double e12,
633 G4double e21, G4double e22, 633 G4double e21, G4double e22,
634 G4double xs11, G4double xs1 634 G4double xs11, G4double xs12,
635 G4double xs21, G4double xs2 635 G4double xs21, G4double xs22,
636 G4double t1, G4double t2, 636 G4double t1, G4double t2,
637 G4double t, G4double e) 637 G4double t, G4double e)
638 { 638 {
639 // Log-Log 639 // Log-Log
640 G4double interpolatedvalue1 = LogLogInterpol 640 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
641 G4double interpolatedvalue2 = LogLogInterpol 641 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
642 G4double value = LogLogInterpolate(t1, t2, t 642 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
643 643
644 return value; 644 return value;
645 } 645 }
646 646