Geant4  10.03.p03
 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 #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 = 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 }
G4Exp.hh
C2
const double C2
Definition: MCGIDI_KalbachMann.cc:15
C1
const double C1
Definition: MCGIDI_KalbachMann.cc:15
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4NistManager.hh
d2
static const G4double d2
Definition: G4ElectroNuclearCrossSection.cc:89
G4ecpssrBaseKxsModel.hh
GT
Definition: Evaluator.cc:66
eps
static const G4double eps
Definition: G4ElectroNuclearCrossSection.cc:99
test.x
tuple x
Definition: test.py:50
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:80
G4int
int G4int
Definition: G4Types.hh:78
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
C3
#define C3
test.b
tuple b
Definition: test.py:12
G4cout
G4GLOB_DLL std::ostream G4cout
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4Proton.hh
eplus
static constexpr double eplus
Definition: G4SIunits.hh:199
G4CrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4CrossSectionDataSet.cc:482
G4Alpha.hh
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4ecpssrBaseKxsModel::CalculateCrossSection
G4double CalculateCrossSection(G4int, G4double, G4double)
Definition: G4ecpssrBaseKxsModel.cc:197
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4ecpssrBaseKxsModel::ExpIntFunction
G4double ExpIntFunction(G4int n, G4double x)
Definition: G4ecpssrBaseKxsModel.cc:125
d1
static const G4double d1
Definition: G4ElectroNuclearCrossSection.cc:88
print
void print(G4double elem)
Definition: G4ecpssrBaseKxsModel.cc:108
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4AtomicTransitionManager.hh
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
plottest35.t1
tuple t1
Definition: plottest35.py:33
FatalException
Definition: G4ExceptionSeverity.hh:60
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:359
G4endl
#define G4endl
Definition: G4ios.hh:61
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4SemiLogInterpolation.hh
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
python.hepunit.amu_c2
float amu_c2
Definition: hepunit.py:277
test.c
tuple c
Definition: test.py:13
barn
static constexpr double barn
Definition: G4SIunits.hh:105
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:45
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:95
G4CrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4CrossSectionDataSet.cc:84