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