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