193{
194
195 if (zTarget <=4) return 0.;
196
197
198
199
201
203
207
209
211
212 else
213 {
215
217
218 else
219 {
220 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
222 return 0;
223 }
224 }
225
227
229
230 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
231
233
234
235 G4double screenedzTarget = zTarget-zlshell;
236
237 const G4double rydbergMeV= 13.6056923e-6;
238
240
241
242 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
243
244
245 if (verboseLevel>0)
G4cout <<
" tetal1=" << tetal1<<
G4endl;
246
247 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
248
249
250
251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
252
253 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
254
255
256
258
259 if (verboseLevel>0)
G4cout <<
" velocityl1=" << velocityl1<<
G4endl;
260
261 const G4double l1AnalyticalApproximation= 1.5;
262 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
263
264
265
267
268 G4double electrIonizationEnergyl1=0.;
269
270
271
272
273 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*
pi*(std::log(1./(x1*x1))-1.);
274 else
275 {
276 if ( x1<=3.)
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));
278 else
279 {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
280 }
281
282 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3));
283
284
285 if (verboseLevel>0)
G4cout <<
" hFunctionl1=" << hFunctionl1<<
G4endl;
286
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.);
288
289
290 if (verboseLevel>0)
G4cout <<
" gFunctionl1=" << gFunctionl1<<
G4endl;
291
292 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1));
293
294
295
296
297 if (verboseLevel>0)
G4cout <<
"sigmaPSS_l1 =" << sigmaPSS_l1<<
G4endl;
298
300
301 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
302
303
304
305
306 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;
307
308
309
310
311
313
315
317
318
319
320 if ( velocityl1 <20. )
321 {
322
323 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
324
325
326
327
328
329 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
330
331 universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
332
333 if (verboseLevel>0)
G4cout <<
"at low velocity range, universalFunction_l1 =" << universalFunction_l1 <<
G4endl;
334
335 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;
336
337
338 if (verboseLevel>0)
G4cout <<
" at low velocity range, sigma PWBA L1 CS = " << sigmaPSSR_l1<<
G4endl;
339
340 }
341
342 else
343
344 {
345
346 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
347
348
349
350
351
352 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
353
354 universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
355
356 if (verboseLevel>0)
G4cout <<
"at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 <<
G4endl;
357
358 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;
359
360
361 if (verboseLevel>0)
G4cout <<
" sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<<
G4endl;
362 }
363
364 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
365
366
367
368 if (verboseLevel>0)
G4cout <<
" pssDeltal1=" << pssDeltal1<<
G4endl;
369
370 if (pssDeltal1>1) return 0.;
371
372 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
373
374
375
376 if (verboseLevel>0)
G4cout <<
" energyLossl1=" << energyLossl1<<
G4endl;
377
379 (8.*
pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
380
381
382
383 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
384
385
387
388 if (verboseLevel>0)
G4cout <<
" coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 <<
G4endl;
389
390 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
391
392
393
394
395 if (verboseLevel>0)
G4cout <<
" crossSection_L1 =" << crossSection_L1 <<
G4endl;
396
397 if (crossSection_L1 >= 0)
398 {
399 return crossSection_L1 * barn;
400 }
401
402 else {return 0;}
403}
G4DLLIMPORT std::ostream G4cout
G4double BindingEnergy() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4AtomicTransitionManager * Instance()
static G4NistManager * Instance()
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4double GetPDGCharge() const
static G4Proton * Proton()
G4double ExpIntFunction(G4int n, G4double x)
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)