Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
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::istream &aDataFile)
 
void SetManager (G4InterpolationManager &aManager)
 

Detailed Description

Definition at line 36 of file G4NeutronHPLegendreStore.hh.

Constructor & Destructor Documentation

G4NeutronHPLegendreStore::G4NeutronHPLegendreStore ( G4int  n)
inline

Definition at line 40 of file G4NeutronHPLegendreStore.hh.

References n.

41  {
42  theCoeff = new G4NeutronHPLegendreTable[n];
43  nEnergy = n;
44  }
const G4int n
G4NeutronHPLegendreStore::~G4NeutronHPLegendreStore ( )
inline

Definition at line 46 of file G4NeutronHPLegendreStore.hh.

47  {
48  delete [] theCoeff;
49  }

Member Function Documentation

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

Definition at line 67 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetCoeff().

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

67 {return theCoeff[i].GetCoeff(l);}
G4double G4NeutronHPLegendreStore::GetEnergy ( G4int  i)
inline

Definition at line 68 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetEnergy().

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

68 {return theCoeff[i].GetEnergy();}
G4int G4NeutronHPLegendreStore::GetNumberOfPoly ( G4int  i)
inline

Definition at line 70 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetNumberOfPoly().

70 {return theCoeff[i].GetNumberOfPoly();}
G4double G4NeutronHPLegendreStore::GetTemperature ( G4int  i)
inline

Definition at line 69 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::GetTemperature().

69 {return theCoeff[i].GetTemperature();}
void G4NeutronHPLegendreStore::Init ( G4int  i,
G4double  e,
G4int  n 
)
inline

Definition at line 51 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::Init().

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

52  {
53  theCoeff[i].Init(e, n);
54  }
void Init(std::istream &aDataFile)
const G4int n
void G4NeutronHPLegendreStore::InitInterpolation ( std::istream &  aDataFile)
inline

Definition at line 78 of file G4NeutronHPLegendreStore.hh.

References G4InterpolationManager::Init().

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

79  {
80  theManager.Init(aDataFile);
81  }
void Init(G4int aScheme, G4int aRange)
G4double G4NeutronHPLegendreStore::Integrate ( G4int  k,
G4double  costh 
)

Definition at line 301 of file G4NeutronHPLegendreStore.cc.

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

Referenced by Sample().

302 {
303  G4double result=0;
305 // G4cout <<"the COEFFS "<<k<<" ";
306 // G4cout <<theCoeff[k].GetNumberOfPoly()<<" ";
307  for(G4int l=0; l<theCoeff[k].GetNumberOfPoly() ; l++)
308  {
309  result += theCoeff[k].GetCoeff(l)*theLeg.Integrate(l, costh);
310 // G4cout << theCoeff[k].GetCoeff(l)<<" ";
311  }
312 // G4cout <<G4endl;
313  return result;
314 }
int G4int
Definition: G4Types.hh:78
G4double Integrate(G4int l, G4double costh)
double G4double
Definition: G4Types.hh:76
G4double G4NeutronHPLegendreStore::Sample ( G4double  energy)

Definition at line 248 of file G4NeutronHPLegendreStore.cc.

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

249 {
250  G4int i0;
251  G4int low(0), high(0);
252 // G4cout << "G4NeutronHPLegendreStore::Sample "<<energy<<" "<<energy<<" "<<nEnergy<<G4endl;
253  for (i0=0; i0<nEnergy; i0++)
254  {
255 // G4cout <<"theCoeff["<<i0<<"].GetEnergy() = "<<theCoeff[i0].GetEnergy()<<G4endl;
256  high = i0;
257  if(theCoeff[i0].GetEnergy()>energy) break;
258  }
259  low = std::max(0, high-1);
260 // G4cout << "G4NeutronHPLegendreStore::Sample high, low: "<<high<<", "<<low<<G4endl;
261  G4NeutronHPVector theBuffer;
263  G4double x1, x2, y1, y2, y;
264  x1 = theCoeff[low].GetEnergy();
265  x2 = theCoeff[high].GetEnergy();
266 // G4cout << "the xes "<<x1<<" "<<x2<<G4endl;
267  G4double costh=0;
268  for(i0=0; i0<601; i0++)
269  {
270  costh = G4double(i0-300)/300.;
271  y1 = Integrate(low, costh);
272  y2 = Integrate(high, costh);
273  y = theInt.Interpolate(theManager.GetScheme(high), energy, x1, x2, y1, y2);
274  theBuffer.SetData(i0, costh, y);
275 // G4cout << "Integration "<<low<<" "<<costh<<" "<<y1<<" "<<y2<<" "<<y<<G4endl;
276  }
277  G4double rand = G4UniformRand();
278  G4int it;
279  for (i0=1; i0<601; i0++)
280  {
281  it = i0;
282  if(rand < theBuffer.GetY(i0)/theBuffer.GetY(600)) break;
283 // G4cout <<"sampling now "<<i0<<" "
284 // << theBuffer.GetY(i0)<<" "
285 // << theBuffer.GetY(600)<<" "
286 // << rand<<" "
287 // << theBuffer.GetY(i0)/theBuffer.GetY(600)<<G4endl;;
288  }
289  if(it==601) it=600;
290 // G4cout << "G4NeutronHPLegendreStore::Sample it "<<rand<<" "<<it<<G4endl;
291  G4double norm = theBuffer.GetY(600);
292  if(norm==0) return -DBL_MAX;
293  x1 = theBuffer.GetY(it)/norm;
294  x2 = theBuffer.GetY(it-1)/norm;
295  y1 = theBuffer.GetX(it);
296  y2 = theBuffer.GetX(it-1);
297 // G4cout << "G4NeutronHPLegendreStore::Sample x y "<<x1<<" "<<y1<<" "<<x2<<" "<<y2<<G4endl;
298  return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2);
299 }
G4double GetY(G4double x)
G4double Integrate(G4int k, G4double costh)
G4double GetX(G4int i) const
void SetData(G4int i, G4double x, G4double y)
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
#define G4UniformRand()
Definition: Randomize.hh:87
G4InterpolationScheme GetScheme(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4double G4NeutronHPLegendreStore::SampleDiscreteTwoBody ( G4double  anEnergy)

Definition at line 42 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), GetEnergy(), G4NeutronHPLegendreTable::GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), G4NeutronHPInterpolator::Interpolate(), G4INCL::Math::max(), and test::x.

Referenced by G4NeutronHPDiscreteTwoBody::Sample().

43 {
44  G4double result;
45 
46  G4int i0;
47  G4int low(0), high(0);
49  for (i0=0; i0<nEnergy; i0++)
50  {
51  high = i0;
52  if(theCoeff[i0].GetEnergy()>anEnergy) break;
53  }
54  low = std::max(0, high-1);
56  G4double x, x1, x2;
57  x = anEnergy;
58  x1 = theCoeff[low].GetEnergy();
59  x2 = theCoeff[high].GetEnergy();
60  G4double theNorm = 0;
61  G4double try01=0, try02=0;
62  G4double max1, max2, costh;
63  max1 = 0; max2 = 0;
64  G4int l,m_tmp;
65  for(i0=0; i0<601; i0++)
66  {
67  costh = G4double(i0-300)/300.;
68  try01 = 0.5;
69  for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++)
70  {
71  l=m_tmp+1;
72  try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh);
73  }
74  if(try01>max1) max1=try01;
75  try02 = 0.5;
76  for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++)
77  {
78  l=m_tmp+1;
79  try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh);
80  }
81  if(try02>max2) max2=try02;
82  }
83  theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
84 
85  G4double value, random;
86  G4double v1, v2;
87  do
88  {
89  v1 = 0.5;
90  v2 = 0.5;
91  result = 2.*G4UniformRand()-1.;
92  for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++)
93  {
94  l=m_tmp+1;
95  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
96  v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*legend;
97  }
98  for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++)
99  {
100  l=m_tmp+1;
101  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
102  v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*legend;
103  }
104  // v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
105  // v2 = std::max(0.,v2);
106  value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
107  random = G4UniformRand();
108  if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
109  }
110  while(random>value/theNorm);
111 
112  return result;
113 }
G4double Evaluate(G4int l, G4double costh)
G4double GetCoeff(G4int i, G4int l)
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4InterpolationScheme GetScheme(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const XML_Char int const XML_Char * value
double G4double
Definition: G4Types.hh:76
G4double G4NeutronHPLegendreStore::SampleElastic ( G4double  anEnergy)

Definition at line 187 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), GetEnergy(), G4NeutronHPLegendreTable::GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), G4NeutronHPInterpolator::Interpolate(), G4INCL::Math::max(), and test::x.

Referenced by G4NeutronHPElasticFS::ApplyYourself().

188 {
189  G4double result;
190 
191  G4int i0;
192  G4int low(0), high(0);
194  for (i0=0; i0<nEnergy; i0++)
195  {
196  high = i0;
197  if(theCoeff[i0].GetEnergy()>anEnergy) break;
198  }
199  low = std::max(0, high-1);
201  G4double x, x1, x2;
202  x = anEnergy;
203  x1 = theCoeff[low].GetEnergy();
204  x2 = theCoeff[high].GetEnergy();
205  G4double theNorm = 0;
206  G4double try01=0, try02=0, try11=0, try12=0;
207  G4double try1, try2;
208  G4int l;
209  for(l=0; l<theCoeff[low].GetNumberOfPoly(); l++)
210  {
211  try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, -1.);
212  try11 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, +1.);
213  }
214  for(l=0; l<theCoeff[high].GetNumberOfPoly(); l++)
215  {
216  try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, -1.);
217  try12 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, +1.);
218  }
219  try1 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try01, try02);
220  try2 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try11, try12);
221  theNorm = std::max(try1, try2);
222 
223  G4double value, random;
224  G4double v1, v2;
225  do
226  {
227  v1 = 0;
228  v2 = 0;
229  result = 2.*G4UniformRand()-1.;
230  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
231  {
232  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
233  v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
234  }
235  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
236  {
237  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
238  v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
239  }
240  value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
241  random = G4UniformRand();
242  }
243  while(random>value/theNorm);
244 
245  return result;
246 }
G4double Evaluate(G4int l, G4double costh)
G4double GetCoeff(G4int i, G4int l)
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4InterpolationScheme GetScheme(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const XML_Char int const XML_Char * value
double G4double
Definition: G4Types.hh:76
G4double G4NeutronHPLegendreStore::SampleMax ( G4double  energy)

Definition at line 117 of file G4NeutronHPLegendreStore.cc.

References G4NeutronHPFastLegendre::Evaluate(), G4UniformRand, GetCoeff(), GetEnergy(), G4NeutronHPLegendreTable::GetEnergy(), G4NeutronHPLegendreTable::GetNumberOfPoly(), G4InterpolationManager::GetScheme(), G4NeutronHPInterpolator::Interpolate(), G4INCL::Math::max(), and test::x.

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

118 {
119  G4double result;
120 
121  G4int i0;
122  G4int low(0), high(0);
124  for (i0=0; i0<nEnergy; i0++)
125  {
126  high = i0;
127  if(theCoeff[i0].GetEnergy()>anEnergy) break;
128  }
129  low = std::max(0, high-1);
131  G4double x, x1, x2;
132  x = anEnergy;
133  x1 = theCoeff[low].GetEnergy();
134  x2 = theCoeff[high].GetEnergy();
135  G4double theNorm = 0;
136  G4double try01=0, try02=0;
137  G4double max1, max2, costh;
138  max1 = 0; max2 = 0;
139  G4int l;
140  for(i0=0; i0<601; i0++)
141  {
142  costh = G4double(i0-300)/300.;
143  try01 = 0;
144  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
145  {
146  try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, costh);
147  }
148  if(try01>max1) max1=try01;
149  try02 = 0;
150  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
151  {
152  try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, costh);
153  }
154  if(try02>max2) max2=try02;
155  }
156  theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
157 
158  G4double value, random;
159  G4double v1, v2;
160  do
161  {
162  v1 = 0;
163  v2 = 0;
164  result = 2.*G4UniformRand()-1.;
165  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
166  {
167  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
168  v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
169  }
170  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
171  {
172  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
173  v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
174  }
175  v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
176  v2 = std::max(0.,v2);
177  value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
178  random = G4UniformRand();
179  if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
180  }
181  while(random>value/theNorm);
182 
183  return result;
184 }
G4double Evaluate(G4int l, G4double costh)
G4double GetCoeff(G4int i, G4int l)
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4InterpolationScheme GetScheme(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const XML_Char int const XML_Char * value
double G4double
Definition: G4Types.hh:76
void G4NeutronHPLegendreStore::SetCoeff ( G4int  i,
G4int  l,
G4double  coeff 
)
inline
void G4NeutronHPLegendreStore::SetCoeff ( G4int  i,
G4NeutronHPLegendreTable theTable 
)
inline

Definition at line 59 of file G4NeutronHPLegendreStore.hh.

60  {
61  if(i>nEnergy) throw G4HadronicException(__FILE__, __LINE__, "LegendreTableIndex out of range");
62  theCoeff[i] = *theTable;
63 // not here -- see G4NeutronHPPhotonDist.cc line 275
64 // delete theTable;
65  }
void G4NeutronHPLegendreStore::SetEnergy ( G4int  i,
G4double  energy 
)
inline

Definition at line 56 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::SetEnergy().

56 { theCoeff[i].SetEnergy(energy); }
double precision function energy(A, Z)
Definition: dpm25nuc6.f:4106
void G4NeutronHPLegendreStore::SetManager ( G4InterpolationManager aManager)
inline

Definition at line 83 of file G4NeutronHPLegendreStore.hh.

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

84  {
85  theManager = aManager;
86  }
void G4NeutronHPLegendreStore::SetNPoints ( G4int  n)
inline

Definition at line 55 of file G4NeutronHPLegendreStore.hh.

References n.

55 { nEnergy = n; }
const G4int n
void G4NeutronHPLegendreStore::SetTemperature ( G4int  i,
G4double  temp 
)
inline

Definition at line 57 of file G4NeutronHPLegendreStore.hh.

References G4NeutronHPLegendreTable::SetTemperature().

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

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

The documentation for this class was generated from the following files: