Geant4_10
G4NeutronHPLegendreStore.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 080612 SampleDiscreteTwoBody contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
31 //
33 #include "G4NeutronHPVector.hh"
36 #include "Randomize.hh"
37 #include <iostream>
38 
39 
40 
41 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
43 {
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 }
114 
115 
116 
118 {
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 }
185 
186 
188 {
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 }
247 
248 G4double G4NeutronHPLegendreStore::Sample (G4double energy) // still in interpolation; do not use
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 }
300 
301 G4double G4NeutronHPLegendreStore::Integrate(G4int k, G4double costh) // still in interpolation; not used anymore
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 }
Double_t y2[nxs]
Definition: Style.C:21
G4double GetY(G4double x)
TLegend * legend
Definition: egs.C:37
Double_t y1[nxs]
Definition: Style.C:20
G4double SampleMax(G4double energy)
Double_t x2[nxs]
Definition: Style.C:19
G4double Integrate(G4int k, G4double costh)
Float_t norm
G4double G4NeutronHPJENDLHEData::G4double result
G4double Evaluate(G4int l, G4double costh)
G4double GetX(G4int i) const
tuple x
Definition: test.py:50
void SetData(G4int i, G4double x, G4double y)
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
Double_t y
Definition: plot.C:279
#define G4UniformRand()
Definition: Randomize.hh:87
double energy
Definition: plottest35.C:25
G4double SampleElastic(G4double anEnergy)
G4InterpolationScheme GetScheme(G4int index) const
Double_t x1[nxs]
Definition: Style.C:18
G4double SampleDiscreteTwoBody(G4double anEnergy)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Sample(G4double energy)
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double Integrate(G4int l, G4double costh)
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83