Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPLegendreStore.cc
이 파일의 문서화 페이지로 가기
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 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4ParticleHPVector.hh"
38 #include "Randomize.hh"
39 #include <iostream>
40 
41 
42 
43 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
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 }
123 
124 
125 
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 }
200 
201 
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 }
269 
270 G4double G4ParticleHPLegendreStore::Sample (G4double energy) // still in interpolation; do not use
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 }
322 
323 G4double G4ParticleHPLegendreStore::Integrate(G4int k, G4double costh) // still in interpolation; not used anymore
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 }
Float_t x
Definition: compare.C:6
G4double Integrate(G4int k, G4double costh)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Integrate(G4int l, G4double costh)
G4ParticleHPLegendreTable * theCoeff
Float_t y1[n_points_granero]
Definition: compare.C:5
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4double GetX(G4int i) const
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
G4double Sample(G4double energy)
G4InterpolationScheme GetScheme(G4int index) const
Float_t y2[n_points_geant4]
Definition: compare.C:26
double G4double
Definition: G4Types.hh:76
G4double GetY(G4double x)
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double SampleDiscreteTwoBody(G4double anEnergy)
#define G4UniformRand()
Definition: Randomize.hh:53
double energy
Definition: plottest35.C:25
void SetData(G4int i, G4double x, G4double y)
TLegend * legend
Definition: egs.C:40
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
G4double GetCoeff(G4int i, G4int l)
Float_t norm
G4double SampleMax(G4double energy)
G4GLOB_DLL std::ostream G4cout
G4double SampleElastic(G4double anEnergy)
Float_t x2[n_points_geant4]
Definition: compare.C:26
#define DBL_MAX
Definition: templates.hh:83
G4double Evaluate(G4int l, G4double costh)