Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
45 G4ecpssrBaseKxsModel::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 
107 void print (G4double elem)
108 {
109  G4cout << elem << " ";
110 }
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 
113 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel()
114 {
115 
116  delete tableC1;
117  delete tableC2;
118  delete tableC3;
119 
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 
124 G4double 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  }
190 return ans;
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194 
195 
196 G4double 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 
541 G4double 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 
645 G4double 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 
659 G4double 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 
674 G4double 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 
Generated on Sat May 25 2013 14:33:33 for Geant4 by   doxygen 1.8.4