G4NeutronHPLegendreStore Class Reference

#include <G4NeutronHPLegendreStore.hh>


Public Member Functions

 G4NeutronHPLegendreStore (G4int n)
 ~G4NeutronHPLegendreStore ()
void Init (G4int i, G4double e, G4int n)
void SetNPoints (G4int n)
void SetEnergy (G4int i, G4double energy)
void SetTemperature (G4int i, G4double temp)
void SetCoeff (G4int i, G4int l, G4double coeff)
void SetCoeff (G4int i, G4NeutronHPLegendreTable *theTable)
G4double GetCoeff (G4int i, G4int l)
G4double GetEnergy (G4int i)
G4double GetTemperature (G4int i)
G4int GetNumberOfPoly (G4int i)
G4double SampleDiscreteTwoBody (G4double anEnergy)
G4double SampleElastic (G4double anEnergy)
G4double Sample (G4double energy)
G4double SampleMax (G4double energy)
G4double Integrate (G4int k, G4double costh)
void InitInterpolation (std::ifstream &aDataFile)
void SetManager (G4InterpolationManager &aManager)


Detailed Description

Definition at line 37 of file G4NeutronHPLegendreStore.hh.


Constructor & Destructor Documentation

G4NeutronHPLegendreStore::G4NeutronHPLegendreStore ( G4int  n  )  [inline]

Definition at line 41 of file G4NeutronHPLegendreStore.hh.

00042   {
00043     theCoeff = new G4NeutronHPLegendreTable[n];
00044     nEnergy = n;
00045   }

G4NeutronHPLegendreStore::~G4NeutronHPLegendreStore (  )  [inline]

Definition at line 47 of file G4NeutronHPLegendreStore.hh.

00048   {
00049     delete [] theCoeff;
00050   }


Member Function Documentation

G4double G4NeutronHPLegendreStore::GetCoeff ( G4int  i,
G4int  l 
) [inline]

Definition at line 68 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetCoeff().

Referenced by Integrate(), SampleDiscreteTwoBody(), SampleElastic(), and SampleMax().

00068 {return theCoeff[i].GetCoeff(l);}

G4double G4NeutronHPLegendreStore::GetEnergy ( G4int  i  )  [inline]

Definition at line 69 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetEnergy().

Referenced by Sample(), SampleDiscreteTwoBody(), SampleElastic(), and SampleMax().

00069 {return theCoeff[i].GetEnergy();}

G4int G4NeutronHPLegendreStore::GetNumberOfPoly ( G4int  i  )  [inline]

Definition at line 71 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetNumberOfPoly().

00071 {return theCoeff[i].GetNumberOfPoly();}

G4double G4NeutronHPLegendreStore::GetTemperature ( G4int  i  )  [inline]

Definition at line 70 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetTemperature().

00070 {return theCoeff[i].GetTemperature();}

void G4NeutronHPLegendreStore::Init ( G4int  i,
G4double  e,
G4int  n 
) [inline]

Definition at line 52 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::Init().

Referenced by G4NeutronHPElasticFS::Init(), G4NeutronHPAngular::Init(), and G4NeutronHPContAngularPar::Sample().

00053   {
00054     theCoeff[i].Init(e, n);
00055   }

void G4NeutronHPLegendreStore::InitInterpolation ( std::ifstream &  aDataFile  )  [inline]

Definition at line 79 of file G4NeutronHPLegendreStore.hh.

References G4InterpolationManager::Init().

Referenced by G4NeutronHPElasticFS::Init(), and G4NeutronHPAngular::Init().

00080   {
00081     theManager.Init(aDataFile);
00082   }

G4double G4NeutronHPLegendreStore::Integrate ( G4int  k,
G4double  costh 
)

Definition at line 301 of file G4NeutronHPLegendreStore.cc.

References GetCoeff(), G4NeutronHPLegendreTable::GetNumberOfPoly(), and G4NeutronHPFastLegendre::Integrate().

Referenced by Sample().

00302 {
00303   G4double result=0;
00304   G4NeutronHPFastLegendre theLeg;
00305 //  G4cout <<"the COEFFS "<<k<<" ";
00306 //  G4cout <<theCoeff[k].GetNumberOfPoly()<<" ";
00307   for(G4int l=0; l<theCoeff[k].GetNumberOfPoly() ; l++)
00308   {
00309     result += theCoeff[k].GetCoeff(l)*theLeg.Integrate(l, costh);
00310 //    G4cout << theCoeff[k].GetCoeff(l)<<" ";
00311   } 
00312 //  G4cout <<G4endl;
00313   return result;
00314 }

G4double G4NeutronHPLegendreStore::Sample ( G4double  energy  ) 

Definition at line 248 of file G4NeutronHPLegendreStore.cc.

References DBL_MAX, G4UniformRand, G4NeutronHPLegendreTable::GetEnergy(), GetEnergy(), G4InterpolationManager::GetScheme(), G4NeutronHPVector::GetX(), G4NeutronHPVector::GetY(), Integrate(), G4NeutronHPInterpolator::Interpolate(), and G4NeutronHPVector::SetData().

00249 {
00250   G4int i0;
00251   G4int low(0), high(0);
00252 //  G4cout << "G4NeutronHPLegendreStore::Sample "<<energy<<" "<<energy<<" "<<nEnergy<<G4endl;
00253   for (i0=0; i0<nEnergy; i0++)
00254   {
00255 //     G4cout <<"theCoeff["<<i0<<"].GetEnergy() = "<<theCoeff[i0].GetEnergy()<<G4endl;
00256     high = i0;
00257     if(theCoeff[i0].GetEnergy()>energy) break;
00258   }
00259   low = std::max(0, high-1);
00260 //  G4cout << "G4NeutronHPLegendreStore::Sample high, low: "<<high<<", "<<low<<G4endl;
00261   G4NeutronHPVector theBuffer;
00262   G4NeutronHPInterpolator theInt;
00263   G4double x1, x2, y1, y2, y;
00264   x1 = theCoeff[low].GetEnergy();
00265   x2 = theCoeff[high].GetEnergy();
00266 //  G4cout << "the xes "<<x1<<" "<<x2<<G4endl;
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 //     G4cout << "Integration "<<low<<" "<<costh<<" "<<y1<<" "<<y2<<" "<<y<<G4endl;
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 //    G4cout <<"sampling now "<<i0<<" "
00284 //         << theBuffer.GetY(i0)<<" "
00285 //         << theBuffer.GetY(600)<<" "
00286 //         << rand<<" "
00287 //         << theBuffer.GetY(i0)/theBuffer.GetY(600)<<G4endl;;
00288   }
00289   if(it==601) it=600;
00290 //  G4cout << "G4NeutronHPLegendreStore::Sample it "<<rand<<" "<<it<<G4endl;
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 //  G4cout << "G4NeutronHPLegendreStore::Sample x y "<<x1<<" "<<y1<<" "<<x2<<" "<<y2<<G4endl;
00298   return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2);
00299 }

G4double G4NeutronHPLegendreStore::SampleDiscreteTwoBody ( G4double  anEnergy  ) 

Definition at line 42 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), G4NeutronHPLegendreTable::GetEnergy(), GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), and G4NeutronHPInterpolator::Interpolate().

Referenced by G4NeutronHPDiscreteTwoBody::Sample().

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); // @@@ done to avoid optimization error on SUN
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); // @@@ done to avoid optimization error on SUN
00102         v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*legend;
00103     } 
00104     // v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
00105     // v2 = std::max(0.,v2); 
00106     value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
00107     random = G4UniformRand();
00108     if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
00109   }
00110   while(random>value/theNorm);
00111 
00112   return result;
00113 }

G4double G4NeutronHPLegendreStore::SampleElastic ( G4double  anEnergy  ) 

Definition at line 187 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), G4NeutronHPLegendreTable::GetCoeff(), G4NeutronHPLegendreTable::GetEnergy(), GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), and G4NeutronHPInterpolator::Interpolate().

Referenced by G4NeutronHPElasticFS::ApplyYourself().

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); // @@@ done to avoid optimization error on SUN
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); // @@@ done to avoid optimization error on SUN
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 }

G4double G4NeutronHPLegendreStore::SampleMax ( G4double  energy  ) 

Definition at line 117 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), G4NeutronHPLegendreTable::GetCoeff(), G4NeutronHPLegendreTable::GetEnergy(), GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), and G4NeutronHPInterpolator::Interpolate().

Referenced by G4NeutronHPPhotonDist::GetPhotons(), G4NeutronHPContAngularPar::Sample(), and G4NeutronHPAngular::SampleAndUpdate().

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); // @@@ done to avoid optimization error on SUN
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); // @@@ done to avoid optimization error on SUN
00173       v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
00174     } 
00175     v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
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; // Workaround for negative cross-section values. @@@@ 31 May 2000
00180   }
00181   while(random>value/theNorm);
00182   
00183   return result;
00184 }

void G4NeutronHPLegendreStore::SetCoeff ( G4int  i,
G4NeutronHPLegendreTable theTable 
) [inline]

Definition at line 60 of file G4NeutronHPLegendreStore.hh.

00061   {
00062     if(i>nEnergy) throw G4HadronicException(__FILE__, __LINE__, "LegendreTableIndex out of range");
00063     theCoeff[i] = *theTable;
00064 // not here -- see G4NeutronHPPhotonDist.cc line 275
00065 //    delete theTable;
00066   }

void G4NeutronHPLegendreStore::SetCoeff ( G4int  i,
G4int  l,
G4double  coeff 
) [inline]

Definition at line 59 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::SetCoeff().

Referenced by G4NeutronHPPhotonDist::GetPhotons(), G4NeutronHPElasticFS::Init(), G4NeutronHPAngular::Init(), G4NeutronHPDiscreteTwoBody::Sample(), and G4NeutronHPContAngularPar::Sample().

00059 {theCoeff[i].SetCoeff(l, coeff); }

void G4NeutronHPLegendreStore::SetEnergy ( G4int  i,
G4double  energy 
) [inline]

Definition at line 57 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::SetEnergy().

00057 { theCoeff[i].SetEnergy(energy); }

void G4NeutronHPLegendreStore::SetManager ( G4InterpolationManager aManager  )  [inline]

Definition at line 84 of file G4NeutronHPLegendreStore.hh.

Referenced by G4NeutronHPDiscreteTwoBody::Sample(), and G4NeutronHPContAngularPar::Sample().

00085   {
00086     theManager = aManager;
00087   }

void G4NeutronHPLegendreStore::SetNPoints ( G4int  n  )  [inline]

Definition at line 56 of file G4NeutronHPLegendreStore.hh.

00056 { nEnergy = n; }

void G4NeutronHPLegendreStore::SetTemperature ( G4int  i,
G4double  temp 
) [inline]

Definition at line 58 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::SetTemperature().

Referenced by G4NeutronHPElasticFS::Init(), and G4NeutronHPAngular::Init().

00058 { theCoeff[i].SetTemperature(temp); }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:39 2013 for Geant4 by  doxygen 1.4.7