ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPLegendreStore.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
34 #include "G4ParticleHPLegendreStore.hh"
35 #include "G4ParticleHPVector.hh"
36 #include "G4ParticleHPInterpolator.hh"
37 #include "G4ParticleHPFastLegendre.hh"
38 #include "Randomize.hh"
39 #include <iostream>
40 
41 
42 
43 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
44 G4double G4ParticleHPLegendreStore::SampleDiscreteTwoBody (G4double anEnergy)
45 {
46  G4double result;
47 
48  G4int i0;
49  G4int low(0), high(0);
50  G4ParticleHPFastLegendre theLeg;
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);
57  G4ParticleHPInterpolator theInt;
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 
126 G4double G4ParticleHPLegendreStore::SampleMax (G4double anEnergy)
127 {
128  G4double result;
129 
130  G4int i0;
131  G4int low(0), high(0);
132  G4ParticleHPFastLegendre theLeg;
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);
139  G4ParticleHPInterpolator theInt;
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 
202 G4double G4ParticleHPLegendreStore::SampleElastic (G4double anEnergy)
203 {
204  G4double result;
205 
206  G4int i0;
207  G4int low(0), high(0);
208  G4ParticleHPFastLegendre theLeg;
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);
215  G4ParticleHPInterpolator theInt;
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;
284  G4ParticleHPInterpolator theInt;
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;
326  G4ParticleHPFastLegendre theLeg;
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 }