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 "G4NeutronHPLegendreStore.hh"
00033 #include "G4NeutronHPVector.hh"
00034 #include "G4NeutronHPInterpolator.hh"
00035 #include "G4NeutronHPFastLegendre.hh"
00036 #include "Randomize.hh"
00037 #include <iostream>
00038
00039
00040
00041
00042 G4double G4NeutronHPLegendreStore::SampleDiscreteTwoBody (G4double anEnergy)
00043 {
00044 G4double result;
00045
00046 G4int i0;
00047 G4int low(0), high(0);
00048 G4NeutronHPFastLegendre theLeg;
00049 for (i0=0; i0<nEnergy; i0++)
00050 {
00051 high = i0;
00052 if(theCoeff[i0].GetEnergy()>anEnergy) break;
00053 }
00054 low = std::max(0, high-1);
00055 G4NeutronHPInterpolator theInt;
00056 G4double x, x1, x2;
00057 x = anEnergy;
00058 x1 = theCoeff[low].GetEnergy();
00059 x2 = theCoeff[high].GetEnergy();
00060 G4double theNorm = 0;
00061 G4double try01=0, try02=0;
00062 G4double max1, max2, costh;
00063 max1 = 0; max2 = 0;
00064 G4int l,m_tmp;
00065 for(i0=0; i0<601; i0++)
00066 {
00067 costh = G4double(i0-300)/300.;
00068 try01 = 0.5;
00069 for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++)
00070 {
00071 l=m_tmp+1;
00072 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh);
00073 }
00074 if(try01>max1) max1=try01;
00075 try02 = 0.5;
00076 for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++)
00077 {
00078 l=m_tmp+1;
00079 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh);
00080 }
00081 if(try02>max2) max2=try02;
00082 }
00083 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
00084
00085 G4double value, random;
00086 G4double v1, v2;
00087 do
00088 {
00089 v1 = 0.5;
00090 v2 = 0.5;
00091 result = 2.*G4UniformRand()-1.;
00092 for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++)
00093 {
00094 l=m_tmp+1;
00095 G4double legend = theLeg.Evaluate(l, result);
00096 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*legend;
00097 }
00098 for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++)
00099 {
00100 l=m_tmp+1;
00101 G4double legend = theLeg.Evaluate(l, result);
00102 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*legend;
00103 }
00104
00105
00106 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
00107 random = G4UniformRand();
00108 if(0>=theNorm) break;
00109 }
00110 while(random>value/theNorm);
00111
00112 return result;
00113 }
00114
00115
00116
00117 G4double G4NeutronHPLegendreStore::SampleMax (G4double anEnergy)
00118 {
00119 G4double result;
00120
00121 G4int i0;
00122 G4int low(0), high(0);
00123 G4NeutronHPFastLegendre theLeg;
00124 for (i0=0; i0<nEnergy; i0++)
00125 {
00126 high = i0;
00127 if(theCoeff[i0].GetEnergy()>anEnergy) break;
00128 }
00129 low = std::max(0, high-1);
00130 G4NeutronHPInterpolator theInt;
00131 G4double x, x1, x2;
00132 x = anEnergy;
00133 x1 = theCoeff[low].GetEnergy();
00134 x2 = theCoeff[high].GetEnergy();
00135 G4double theNorm = 0;
00136 G4double try01=0, try02=0;
00137 G4double max1, max2, costh;
00138 max1 = 0; max2 = 0;
00139 G4int l;
00140 for(i0=0; i0<601; i0++)
00141 {
00142 costh = G4double(i0-300)/300.;
00143 try01 = 0;
00144 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
00145 {
00146 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, costh);
00147 }
00148 if(try01>max1) max1=try01;
00149 try02 = 0;
00150 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
00151 {
00152 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, costh);
00153 }
00154 if(try02>max2) max2=try02;
00155 }
00156 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
00157
00158 G4double value, random;
00159 G4double v1, v2;
00160 do
00161 {
00162 v1 = 0;
00163 v2 = 0;
00164 result = 2.*G4UniformRand()-1.;
00165 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
00166 {
00167 G4double legend = theLeg.Evaluate(l, result);
00168 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
00169 }
00170 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
00171 {
00172 G4double legend = theLeg.Evaluate(l, result);
00173 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
00174 }
00175 v1 = std::max(0.,v1);
00176 v2 = std::max(0.,v2);
00177 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
00178 random = G4UniformRand();
00179 if(0>=theNorm) break;
00180 }
00181 while(random>value/theNorm);
00182
00183 return result;
00184 }
00185
00186
00187 G4double G4NeutronHPLegendreStore::SampleElastic (G4double anEnergy)
00188 {
00189 G4double result;
00190
00191 G4int i0;
00192 G4int low(0), high(0);
00193 G4NeutronHPFastLegendre theLeg;
00194 for (i0=0; i0<nEnergy; i0++)
00195 {
00196 high = i0;
00197 if(theCoeff[i0].GetEnergy()>anEnergy) break;
00198 }
00199 low = std::max(0, high-1);
00200 G4NeutronHPInterpolator theInt;
00201 G4double x, x1, x2;
00202 x = anEnergy;
00203 x1 = theCoeff[low].GetEnergy();
00204 x2 = theCoeff[high].GetEnergy();
00205 G4double theNorm = 0;
00206 G4double try01=0, try02=0, try11=0, try12=0;
00207 G4double try1, try2;
00208 G4int l;
00209 for(l=0; l<theCoeff[low].GetNumberOfPoly(); l++)
00210 {
00211 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, -1.);
00212 try11 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, +1.);
00213 }
00214 for(l=0; l<theCoeff[high].GetNumberOfPoly(); l++)
00215 {
00216 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, -1.);
00217 try12 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, +1.);
00218 }
00219 try1 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try01, try02);
00220 try2 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try11, try12);
00221 theNorm = std::max(try1, try2);
00222
00223 G4double value, random;
00224 G4double v1, v2;
00225 do
00226 {
00227 v1 = 0;
00228 v2 = 0;
00229 result = 2.*G4UniformRand()-1.;
00230 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
00231 {
00232 G4double legend = theLeg.Evaluate(l, result);
00233 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
00234 }
00235 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
00236 {
00237 G4double legend = theLeg.Evaluate(l, result);
00238 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
00239 }
00240 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
00241 random = G4UniformRand();
00242 }
00243 while(random>value/theNorm);
00244
00245 return result;
00246 }
00247
00248 G4double G4NeutronHPLegendreStore::Sample (G4double energy)
00249 {
00250 G4int i0;
00251 G4int low(0), high(0);
00252
00253 for (i0=0; i0<nEnergy; i0++)
00254 {
00255
00256 high = i0;
00257 if(theCoeff[i0].GetEnergy()>energy) break;
00258 }
00259 low = std::max(0, high-1);
00260
00261 G4NeutronHPVector theBuffer;
00262 G4NeutronHPInterpolator theInt;
00263 G4double x1, x2, y1, y2, y;
00264 x1 = theCoeff[low].GetEnergy();
00265 x2 = theCoeff[high].GetEnergy();
00266
00267 G4double costh=0;
00268 for(i0=0; i0<601; i0++)
00269 {
00270 costh = G4double(i0-300)/300.;
00271 y1 = Integrate(low, costh);
00272 y2 = Integrate(high, costh);
00273 y = theInt.Interpolate(theManager.GetScheme(high), energy, x1, x2, y1, y2);
00274 theBuffer.SetData(i0, costh, y);
00275
00276 }
00277 G4double rand = G4UniformRand();
00278 G4int it;
00279 for (i0=1; i0<601; i0++)
00280 {
00281 it = i0;
00282 if(rand < theBuffer.GetY(i0)/theBuffer.GetY(600)) break;
00283
00284
00285
00286
00287
00288 }
00289 if(it==601) it=600;
00290
00291 G4double norm = theBuffer.GetY(600);
00292 if(norm==0) return -DBL_MAX;
00293 x1 = theBuffer.GetY(it)/norm;
00294 x2 = theBuffer.GetY(it-1)/norm;
00295 y1 = theBuffer.GetX(it);
00296 y2 = theBuffer.GetX(it-1);
00297
00298 return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2);
00299 }
00300
00301 G4double G4NeutronHPLegendreStore::Integrate(G4int k, G4double costh)
00302 {
00303 G4double result=0;
00304 G4NeutronHPFastLegendre theLeg;
00305
00306
00307 for(G4int l=0; l<theCoeff[k].GetNumberOfPoly() ; l++)
00308 {
00309 result += theCoeff[k].GetCoeff(l)*theLeg.Integrate(l, costh);
00310
00311 }
00312
00313 return result;
00314 }