00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "G4NeutronHPInterpolator.hh"
00031
00032 G4double G4NeutronHPInterpolator::
00033 GetBinIntegral(const G4InterpolationScheme & aScheme,
00034 const G4double x1,const G4double x2,const G4double y1,const G4double y2)
00035 {
00036 G4double result = 0;
00037 if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
00038 {
00039 result = y1*(x2-x1);
00040 }
00041 else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
00042 {
00043 result = 0.5*(y2+y1)*(x2-x1);
00044 }
00045 else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
00046 {
00047 if(x1==0) result = y1;
00048 else if(x2==0) result = y2;
00049 else
00050 {
00051 G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
00052 G4double a = y1 - b*std::log(x1);
00053 result = (a-b)*(x2-x1) + b*(x2*std::log(x2)-x1*std::log(x1));
00054 }
00055 }
00056 else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
00057 {
00058 if(y1==0||y2==0) result =0;
00059 else
00060 {
00061 G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
00062 G4double a = std::log(y1) - b*x1;
00063
00064
00065
00066
00067 if(b!=0){
00068 result = (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
00069 }
00070 else{
00071 result=y2*(x2-x1);
00072 }
00073
00074 }
00075 }
00076 else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
00077 {
00078 if(x1==0) result = y1;
00079 else if(x2==0) result = y2;
00080 else if(y1==0||y2==0) result =0;
00081 else
00082 {
00083 G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
00084 G4double a = std::log(y1) - b*std::log(x1);;
00085 result = (std::exp(a)/(b+1))*(std::pow(x2,b+1)-std::pow(x1,b+1));
00086 }
00087 }
00088 else
00089 {
00090 throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
00091 }
00092 return result;
00093 }
00094 G4double G4NeutronHPInterpolator::
00095 GetWeightedBinIntegral(const G4InterpolationScheme & aScheme,
00096 const G4double x1,const G4double x2,const G4double y1,const G4double y2)
00097 {
00098 G4double result = 0;
00099 if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
00100 {
00101 result = 0.5*y1*(x2*x2-x1*x1);
00102 }
00103 else if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
00104 {
00105
00106
00107
00108
00109
00110 result = (y1*x2 - y2*x1)*(x2 + x1)/2. + (y2-y1)*(x2*x2 + x2*x1 + x1*x1)/3.;
00111 }
00112 else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
00113 {
00114 if(x1==0) result = y1;
00115 else if(x2==0) result = y2;
00116 else
00117 {
00118 G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
00119 G4double a = y1 - b*std::log(x1);
00120 result = ( x2*x2/2. * (a-b/2.+b*std::log(x2)) )
00121 -( x1*x1/2. * (a-b/2.+b*std::log(x1)) );
00122 }
00123 }
00124 else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
00125 {
00126 if(y1==0||y2==0) result = 0;
00127 else
00128 {
00129 G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
00130 G4double a = std::log(y1) - b*x1;
00131 result = std::exp(a)/(b*b)*( std::exp(b*x2)*(b*x2-1.) - std::exp(b*x1)*(b*x1-1.) );
00132 }
00133 }
00134 else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
00135 {
00136 if(x1==0) result = y1;
00137 else if(x2==0) result = y2;
00138 if(y1==0||y2==0) result = 0;
00139 else
00140 {
00141 G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
00142 G4double a = std::log(y1) - b*std::log(x1);;
00143 result = std::exp(a)/(b+2.)*( std::pow(x2, b+2.) - std::pow(x1, b+2) );
00144 }
00145 }
00146 else
00147 {
00148 throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
00149 }
00150 return result;
00151 }