Geant4  10.02.p01
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"
38 #include "G4NistManager.hh"
39 #include "G4Proton.hh"
40 #include "G4Alpha.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 
46 {
47  verboseLevel=0;
48 
49  // Storing C coefficients for high velocity formula
50 
51  G4String fileC1("pixe/uf/c1");
53 
54  G4String fileC2("pixe/uf/c2");
56 
57  G4String fileC3("pixe/uf/c3");
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 
107 void print (G4double elem)
108 {
109  G4cout << elem << " ";
110 }
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 
114 {
115 
116  delete tableC1;
117  delete tableC2;
118  delete tableC3;
119 
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
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  }
190 return ans;
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194 
195 
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 
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 
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 
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 
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 
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 
const double C2
G4double FunctionFK(G4double k, G4double theta)
static const G4double d1
const double C1
const G4double fact
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static const G4double e2
static const G4double eps
G4double a
Definition: TRTMaterials.hh:39
G4double BindingEnergy() const
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
std::vector< double > dummyVec
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
#define C3
G4CrossSectionDataSet * tableC3
G4CrossSectionDataSet * tableC1
G4GLOB_DLL std::ostream G4cout
virtual G4double FindValue(G4double e, G4int componentId=0) const
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double CalculateCrossSection(G4int, G4double, G4double)
const G4int n
G4double ExpIntFunction(G4int n, G4double x)
void print(G4double elem)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Definition: G4SIunits.hh:212
G4double GetPDGMass() const
static const double pi
Definition: G4SIunits.hh:74
const G4double x[NPOINTSGL]
G4double GetAtomicMassAmu(const G4String &symb) const
G4CrossSectionDataSet * tableC2
#define G4endl
Definition: G4ios.hh:61
Definition: Evaluator.cc:66
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static const double barn
Definition: G4SIunits.hh:104
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:196
G4double GetPDGCharge() const
static const G4double d2
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
virtual G4bool LoadData(const G4String &argFileName)