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
00031
00032 #ifndef G4NeutronHPInterpolator_h
00033 #define G4NeutronHPInterpolator_h 1
00034
00035 #include "globals.hh"
00036 #include "G4InterpolationScheme.hh"
00037 #include "Randomize.hh"
00038 #include "G4ios.hh"
00039 #include "G4HadronicException.hh"
00040
00041
00042 class G4NeutronHPInterpolator
00043 {
00044 public:
00045
00046 G4NeutronHPInterpolator(){}
00047 ~G4NeutronHPInterpolator()
00048 {
00049
00050 }
00051
00052 inline G4double Lin(G4double x,G4double x1,G4double x2,G4double y1,G4double y2)
00053 {
00054 G4double slope=0, off=0;
00055 if(x2-x1==0) return (y2+y1)/2.;
00056 slope = (y2-y1)/(x2-x1);
00057 off = y2-x2*slope;
00058 G4double y = x*slope+off;
00059 return y;
00060 }
00061
00062 inline G4double Interpolate(G4InterpolationScheme aScheme,
00063 G4double x, G4double x1, G4double x2,
00064 G4double y1, G4double y2) const;
00065
00066 G4double
00067 GetBinIntegral(const G4InterpolationScheme & aScheme,
00068 const G4double x1,const G4double x2,const G4double y1,const G4double y2);
00069
00070 G4double
00071 GetWeightedBinIntegral(const G4InterpolationScheme & aScheme,
00072 const G4double x1,const G4double x2,const G4double y1,const G4double y2);
00073
00074 private:
00075
00076 inline G4double Histogram(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00077 inline G4double LinearLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00078 inline G4double LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00079 inline G4double LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00080 inline G4double LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00081 inline G4double Random(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const;
00082
00083 };
00084
00085 inline G4double G4NeutronHPInterpolator::
00086 Interpolate(G4InterpolationScheme aScheme,
00087 G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
00088 {
00089 G4double result(0);
00090 G4int theScheme = aScheme;
00091 theScheme = theScheme%CSTART_;
00092 switch(theScheme)
00093 {
00094 case 1:
00095
00096
00097 result = LinearLinear(x, x1, x2, y1, y2);
00098 break;
00099 case 2:
00100 result = LinearLinear(x, x1, x2, y1, y2);
00101 break;
00102 case 3:
00103 result = LinearLogarithmic(x, x1, x2, y1, y2);
00104 break;
00105 case 4:
00106 result = LogarithmicLinear(x, x1, x2, y1, y2);
00107 break;
00108 case 5:
00109 result = LogarithmicLogarithmic(x, x1, x2, y1, y2);
00110 break;
00111 case 6:
00112 result = Random(x, x1, x2, y1, y2);
00113 break;
00114 default:
00115 G4cout << "theScheme = "<<theScheme<<G4endl;
00116 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPInterpolator::Carthesian Invalid InterpolationScheme");
00117 break;
00118 }
00119 return result;
00120 }
00121
00122 inline G4double G4NeutronHPInterpolator::
00123 Histogram(G4double , G4double , G4double , G4double y1, G4double ) const
00124 {
00125 G4double result;
00126 result = y1;
00127 return result;
00128 }
00129
00130 inline G4double G4NeutronHPInterpolator::
00131 LinearLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
00132 {
00133 G4double slope=0, off=0;
00134 if(x2-x1==0) return (y2+y1)/2.;
00135 slope = (y2-y1)/(x2-x1);
00136 off = y2-x2*slope;
00137 G4double y = x*slope+off;
00138 return y;
00139 }
00140
00141 inline G4double G4NeutronHPInterpolator::
00142 LinearLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
00143 {
00144 G4double result;
00145 if(x==0) result = y1+y2/2.;
00146 else if(x1==0) result = y1;
00147 else if(x2==0) result = y2;
00148 else result = LinearLinear(std::log(x), std::log(x1), std::log(x2), y1, y2);
00149 return result;
00150 }
00151
00152 inline G4double G4NeutronHPInterpolator::
00153 LogarithmicLinear(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
00154 {
00155 G4double result;
00156 if(y1==0||y2==0) result = 0;
00157 else
00158 {
00159 result = LinearLinear(x, x1, x2, std::log(y1), std::log(y2));
00160 result = std::exp(result);
00161 }
00162 return result;
00163 }
00164
00165 inline G4double G4NeutronHPInterpolator::
00166 LogarithmicLogarithmic(G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
00167 {
00168 G4double result;
00169 if(x==0) result = y1+y2/2.;
00170 else if(x1==0) result = y1;
00171 else if(x2==0) result = y2;
00172 if(y1==0||y2==0) result = 0;
00173 else
00174 {
00175 result = LinearLinear(std::log(x), std::log(x1), std::log(x2), std::log(y1), std::log(y2));
00176 result = std::exp(result);
00177 }
00178 return result;
00179 }
00180
00181 inline G4double G4NeutronHPInterpolator::
00182 Random(G4double , G4double , G4double , G4double y1, G4double y2) const
00183 {
00184 G4double result;
00185 result = y1+G4UniformRand()*(y2-y1);
00186 return result;
00187 }
00188
00189 #endif