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 #include "G4NeutronHPPartial.hh"
00033 #include "G4NeutronHPInterpolator.hh"
00034 #include "Randomize.hh"
00035
00036 G4NeutronHPVector * G4NeutronHPPartial::GetY(G4double e1)
00037 {
00038 G4NeutronHPVector * aBuffer = new G4NeutronHPVector();
00039 G4int i;
00040 if(nData==1)
00041 {
00042 for(i=0; i<data[0].GetVectorLength(); i++)
00043 {
00044 aBuffer->SetInterpolationManager(data[0].GetInterpolationManager());
00045 aBuffer->SetData(i , data[0].GetX(i), data[0].GetY(i));
00046 }
00047 return aBuffer;
00048 }
00049 for (i=0; i<nData; i++)
00050 {
00051 if(X[i]>e1) break;
00052 }
00053 if(i==nData) i--;
00054 if(0==i) i=1;
00055 G4double x1,x2,y1,y2,y;
00056 G4int i1=0, ib=0;
00057 G4double E1 = X[i-1];
00058 G4double E2 = X[i];
00059 for(G4int ii=0; ii<data[i].GetVectorLength(); ii++)
00060 {
00061 x1 = data[i-1].GetX(std::min(i1, data[i-1].GetVectorLength()-1));
00062 x2 = data[i].GetX(ii);
00063 if(x1<x2&&i1<data[i-1].GetVectorLength())
00064 {
00065 y1 = data[i-1].GetY(i1);
00066 y2 = data[i].GetY(x1);
00067 if(E2-E1!=0)
00068 {
00069 y = theInt.Interpolate(theManager.GetScheme(i), e1, E1, E2, y1, y2);
00070 }
00071 else
00072 {
00073 y = 0.5*(y1+y2);
00074 }
00075 aBuffer->SetData(ib, x1, y);
00076 aBuffer->SetScheme(ib++, data[i-1].GetScheme(i1));
00077 i1++;
00078 if(x2-x1>0.001*x1)
00079 {
00080 ii--;
00081 }
00082 }
00083 else
00084 {
00085 y1 = data[i-1].GetY(x2);
00086 y2 = data[i].GetY(ii);
00087 if(E2-E1!=0)
00088 {
00089 y = theInt.Interpolate(theManager.GetScheme(i), e1, E1, E2, y1, y2);
00090 }
00091 else
00092 {
00093 y = 0.5*(y1+y2);
00094 }
00095 aBuffer->SetData(ib, x2, y);
00096 aBuffer->SetScheme(ib++, data[i].GetScheme(ii));
00097 if(x1-x2<0.001*x2) i1++;
00098 }
00099 }
00100 return aBuffer;
00101 }
00102
00103 G4double G4NeutronHPPartial::Sample(G4double x)
00104 {
00105 G4int i;
00106 for (i=0; i<nData; i++)
00107 {
00108 if(x<X[i]) break;
00109 }
00110 G4NeutronHPVector theBuff;
00111 if(i==0)
00112 {
00113 theBuff.SetInterpolationManager(data[0].GetInterpolationManager());
00114 for(G4int ii=0;ii<GetNEntries(0);ii++)
00115 {
00116 theBuff.SetX(ii, GetX(0,ii));
00117 theBuff.SetY(ii, GetY(0,ii));
00118 }
00119 }
00120
00121 else if ( i == nData )
00122 {
00123 for(i=0;i<GetNEntries(nData-1);i++)
00124 {
00125 theBuff.SetX(i, GetX(nData-1,i));
00126 theBuff.SetY(i, GetY(nData-1,i));
00127 theBuff.SetInterpolationManager(data[nData-1].GetInterpolationManager());
00128 }
00129 }
00130 else
00131 {
00132 G4int low = i-1;
00133 G4int high = low+1;
00134 G4double x1,x2,y1,y2;
00135 G4int i1=0, i2=0, ii=0;
00136 x1 = X[low];
00137 x2 = X[high];
00138 while(i1<GetNEntries(low)||i2<GetNEntries(high))
00139 {
00140 if( (GetX(low,i1)<GetX(high,i2) && i1<GetNEntries(low))
00141 ||(i2==GetNEntries(high)) )
00142 {
00143 theBuff.SetX(ii, GetX(low,i1));
00144 y1 = GetY(low,i1);
00145 y2 = GetY(high, GetX(low,i1));
00146 theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high),
00147 x, x1, x2, y1, y2));
00148 theBuff.SetScheme(ii, data[low].GetScheme(i1));
00149 if(std::abs(GetX(low,i1)-GetX(high,i2))<0.001) i2++;
00150 i1++;
00151 ii++;
00152 }
00153 else
00154 {
00155 theBuff.SetX(ii, GetX(high,i2));
00156
00157
00158
00159
00160 y2 = GetY(high,i2);
00161 y1 = GetY(low, GetX(high,i2));
00162
00163 theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high),
00164 x, x1, x2, y1, y2));
00165 theBuff.SetScheme(ii, data[high].GetScheme(i2));
00166 if(std::abs(GetX(low,i1)-GetX(high,i2))<0.001) i1++;
00167 i2++;
00168 ii++;
00169 }
00170 }
00171 }
00172
00173 return theBuff.Sample();
00174 }