Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPLegendreStore Class Reference

#include <G4ParticleHPLegendreStore.hh>

Public Member Functions

 G4ParticleHPLegendreStore (G4int n)
 
 ~G4ParticleHPLegendreStore ()
 
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, G4ParticleHPLegendreTable *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 37 of file G4ParticleHPLegendreStore.hh.

Constructor & Destructor Documentation

G4ParticleHPLegendreStore::G4ParticleHPLegendreStore ( G4int  n)
inline

Definition at line 41 of file G4ParticleHPLegendreStore.hh.

42  {
43  theCoeff = new G4ParticleHPLegendreTable[n];
44  nEnergy = n;
45  }
const G4int n
G4ParticleHPLegendreStore::~G4ParticleHPLegendreStore ( )
inline

Definition at line 47 of file G4ParticleHPLegendreStore.hh.

48  {
49  delete [] theCoeff;
50  }

Member Function Documentation

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

Definition at line 68 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPLegendreStore::GetEnergy ( G4int  i)
inline

Definition at line 69 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4ParticleHPLegendreStore::GetNumberOfPoly ( G4int  i)
inline

Definition at line 71 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

G4double G4ParticleHPLegendreStore::GetTemperature ( G4int  i)
inline

Definition at line 70 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

void G4ParticleHPLegendreStore::Init ( G4int  i,
G4double  e,
G4int  n 
)
inline

Definition at line 52 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPLegendreStore::InitInterpolation ( std::istream &  aDataFile)
inline

Definition at line 79 of file G4ParticleHPLegendreStore.hh.

80  {
81  theManager.Init(aDataFile);
82  }
void Init(G4int aScheme, G4int aRange)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 323 of file G4ParticleHPLegendreStore.cc.

324 {
325  G4double result=0;
327 // G4cout <<"the COEFFS "<<k<<" ";
328 // G4cout <<theCoeff[k].GetNumberOfPoly()<<" ";
329  for(G4int l=0; l<theCoeff[k].GetNumberOfPoly() ; l++)
330  {
331  result += theCoeff[k].GetCoeff(l)*theLeg.Integrate(l, costh);
332 // G4cout << theCoeff[k].GetCoeff(l)<<" ";
333  }
334 // G4cout <<G4endl;
335  return result;
336 }
G4double G4ParticleHPJENDLHEData::G4double result
G4double Integrate(G4int l, G4double costh)
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPLegendreStore::Sample ( G4double  energy)

Definition at line 270 of file G4ParticleHPLegendreStore.cc.

271 {
272  G4int i0;
273  G4int low(0), high(0);
274 // G4cout << "G4ParticleHPLegendreStore::Sample "<<energy<<" "<<energy<<" "<<nEnergy<<G4endl;
275  for (i0=0; i0<nEnergy; i0++)
276  {
277 // G4cout <<"theCoeff["<<i0<<"].GetEnergy() = "<<theCoeff[i0].GetEnergy()<<G4endl;
278  high = i0;
279  if(theCoeff[i0].GetEnergy()>energy) break;
280  }
281  low = std::max(0, high-1);
282 // G4cout << "G4ParticleHPLegendreStore::Sample high, low: "<<high<<", "<<low<<G4endl;
283  G4ParticleHPVector theBuffer;
285  G4double x1, x2, y1, y2, y;
286  x1 = theCoeff[low].GetEnergy();
287  x2 = theCoeff[high].GetEnergy();
288 // G4cout << "the xes "<<x1<<" "<<x2<<G4endl;
289  G4double costh=0;
290  for(i0=0; i0<601; i0++)
291  {
292  costh = G4double(i0-300)/300.;
293  y1 = Integrate(low, costh);
294  y2 = Integrate(high, costh);
295  y = theInt.Interpolate(theManager.GetScheme(high), energy, x1, x2, y1, y2);
296  theBuffer.SetData(i0, costh, y);
297 // G4cout << "Integration "<<low<<" "<<costh<<" "<<y1<<" "<<y2<<" "<<y<<G4endl;
298  }
299  G4double rand = G4UniformRand();
300  G4int it;
301  for (i0=1; i0<601; i0++)
302  {
303  it = i0;
304  if(rand < theBuffer.GetY(i0)/theBuffer.GetY(600)) break;
305 // G4cout <<"sampling now "<<i0<<" "
306 // << theBuffer.GetY(i0)<<" "
307 // << theBuffer.GetY(600)<<" "
308 // << rand<<" "
309 // << theBuffer.GetY(i0)/theBuffer.GetY(600)<<G4endl;;
310  }
311  if(it==601) it=600;
312 // G4cout << "G4ParticleHPLegendreStore::Sample it "<<rand<<" "<<it<<G4endl;
313  G4double norm = theBuffer.GetY(600);
314  if(norm==0) return -DBL_MAX;
315  x1 = theBuffer.GetY(it)/norm;
316  x2 = theBuffer.GetY(it-1)/norm;
317  y1 = theBuffer.GetX(it);
318  y2 = theBuffer.GetX(it-1);
319 // G4cout << "G4ParticleHPLegendreStore::Sample x y "<<x1<<" "<<y1<<" "<<x2<<" "<<y2<<G4endl;
320  return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2);
321 }
G4double Integrate(G4int k, G4double costh)
void SetData(G4int i, G4double x, G4double y)
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
G4InterpolationScheme GetScheme(G4int index) const
G4double GetX(G4int i) const
G4double GetY(G4double x)
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

G4double G4ParticleHPLegendreStore::SampleDiscreteTwoBody ( G4double  anEnergy)

Definition at line 44 of file G4ParticleHPLegendreStore.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPLegendreStore::SampleElastic ( G4double  anEnergy)

Definition at line 202 of file G4ParticleHPLegendreStore.cc.

203 {
205 
206  G4int i0;
207  G4int low(0), high(0);
209  for (i0=0; i0<nEnergy; i0++)
210  {
211  high = i0;
212  if(theCoeff[i0].GetEnergy()>anEnergy) break;
213  }
214  low = std::max(0, high-1);
216  G4double x, x1, x2;
217  x = anEnergy;
218  x1 = theCoeff[low].GetEnergy();
219  x2 = theCoeff[high].GetEnergy();
220  G4double theNorm = 0;
221  G4double try01=0, try02=0, try11=0, try12=0;
222  G4double try1, try2;
223  G4int l;
224  for(l=0; l<theCoeff[low].GetNumberOfPoly(); l++)
225  {
226  try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, -1.);
227  try11 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, +1.);
228  }
229  for(l=0; l<theCoeff[high].GetNumberOfPoly(); l++)
230  {
231  try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, -1.);
232  try12 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, +1.);
233  }
234  try1 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try01, try02);
235  try2 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try11, try12);
236  theNorm = std::max(try1, try2);
237 
238  G4double value, random;
239  G4double v1, v2;
240  G4int icounter=0;
241  G4int icounter_max=1024;
242  do
243  {
244  icounter++;
245  if ( icounter > icounter_max ) {
246  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
247  break;
248  }
249  v1 = 0;
250  v2 = 0;
251  result = 2.*G4UniformRand()-1.;
252  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
253  {
254  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
255  v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
256  }
257  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
258  {
259  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
260  v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
261  }
262  value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
263  random = G4UniformRand();
264  }
265  while(random>value/theNorm); // Loop checking, 11.05.2015, T. Koi
266 
267  return result;
268 }
G4double G4ParticleHPJENDLHEData::G4double result
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double GetCoeff(G4int i, G4int l)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4InterpolationScheme GetScheme(G4int index) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double Evaluate(G4int l, G4double costh)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPLegendreStore::SampleMax ( G4double  energy)

Definition at line 126 of file G4ParticleHPLegendreStore.cc.

127 {
129 
130  G4int i0;
131  G4int low(0), high(0);
133  for (i0=0; i0<nEnergy; i0++)
134  {
135  high = i0;
136  if(theCoeff[i0].GetEnergy()>anEnergy) break;
137  }
138  low = std::max(0, high-1);
140  G4double x, x1, x2;
141  x = anEnergy;
142  x1 = theCoeff[low].GetEnergy();
143  x2 = theCoeff[high].GetEnergy();
144  G4double theNorm = 0;
145  G4double try01=0, try02=0;
146  G4double max1, max2, costh;
147  max1 = 0; max2 = 0;
148  G4int l;
149  for(i0=0; i0<601; i0++)
150  {
151  costh = G4double(i0-300)/300.;
152  try01 = 0;
153  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
154  {
155  try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, costh);
156  }
157  if(try01>max1) max1=try01;
158  try02 = 0;
159  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
160  {
161  try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, costh);
162  }
163  if(try02>max2) max2=try02;
164  }
165  theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2);
166 
167  G4double value, random;
168  G4double v1, v2;
169  G4int icounter=0;
170  G4int icounter_max=1024;
171  do
172  {
173  icounter++;
174  if ( icounter > icounter_max ) {
175  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
176  break;
177  }
178  v1 = 0;
179  v2 = 0;
180  result = 2.*G4UniformRand()-1.;
181  for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++)
182  {
183  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
184  v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend;
185  }
186  for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++)
187  {
188  G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN
189  v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend;
190  }
191  v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical.
192  v2 = std::max(0.,v2);
193  value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2);
194  random = G4UniformRand();
195  if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000
196  }
197  while(random>value/theNorm); // Loop checking, 11.05.2015, T. Koi
198  return result;
199 }
G4double G4ParticleHPJENDLHEData::G4double result
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double GetCoeff(G4int i, G4int l)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4InterpolationScheme GetScheme(G4int index) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double Evaluate(G4int l, G4double costh)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPLegendreStore::SetCoeff ( G4int  i,
G4int  l,
G4double  coeff 
)
inline

Definition at line 59 of file G4ParticleHPLegendreStore.hh.

59 {theCoeff[i].SetCoeff(l, coeff); }
void SetCoeff(G4int l, G4double coeff)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPLegendreStore::SetCoeff ( G4int  i,
G4ParticleHPLegendreTable theTable 
)
inline

Definition at line 60 of file G4ParticleHPLegendreStore.hh.

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

Definition at line 57 of file G4ParticleHPLegendreStore.hh.

57 { theCoeff[i].SetEnergy(energy); }
G4double energy(const ThreeVector &p, const G4double m)

Here is the call graph for this function:

void G4ParticleHPLegendreStore::SetManager ( G4InterpolationManager aManager)
inline

Definition at line 84 of file G4ParticleHPLegendreStore.hh.

85  {
86  theManager = aManager;
87  }

Here is the caller graph for this function:

void G4ParticleHPLegendreStore::SetNPoints ( G4int  n)
inline

Definition at line 56 of file G4ParticleHPLegendreStore.hh.

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

Definition at line 58 of file G4ParticleHPLegendreStore.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:


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