| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 080612 SampleDiscreteTwoBody contribution f
31 // #3
32 //
33 // P. Arce, June-2014 Conversion neutron_hp to
34 //
35 #include "G4ParticleHPLegendreStore.hh"
36
37 #include "G4ParticleHPFastLegendre.hh"
38 #include "G4ParticleHPInterpolator.hh"
39 #include "G4ParticleHPVector.hh"
40 #include "Randomize.hh"
41
42 #include <iostream>
43
44 // 080612TK contribution from Benoit Pirard an
45 G4double G4ParticleHPLegendreStore::SampleDisc
46 {
47 G4double result = 0.;
48
49 G4int i0;
50 G4int low(0), high(0);
51 G4ParticleHPFastLegendre theLeg;
52 for (i0 = 0; i0 < nEnergy; i0++) {
53 high = i0;
54 if (theCoeff[i0].GetEnergy() > anEnergy) b
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;
66 max2 = 0;
67 G4int l, m_tmp;
68 for (i0 = 0; i0 < 601; i0++) {
69 costh = G4double(i0 - 300) / 300.;
70 try01 = 0.5;
71 for (m_tmp = 0; m_tmp < theCoeff[low].GetN
72 l = m_tmp + 1;
73 try01 += (2. * l + 1) / 2. * theCoeff[lo
74 }
75 if (try01 > max1) max1 = try01;
76 try02 = 0.5;
77 for (m_tmp = 0; m_tmp < theCoeff[high].Get
78 l = m_tmp + 1;
79 try02 += (2. * l + 1) / 2. * theCoeff[hi
80 }
81 if (try02 > max2) max2 = try02;
82 }
83 theNorm = theInt.Interpolate(theManager.GetS
84
85 G4double value, random;
86 G4double v1, v2;
87 G4int icounter = 0;
88 G4int icounter_max = 1024;
89 do {
90 icounter++;
91 if (icounter > icounter_max) {
92 G4cout << "Loop-counter exceeded the thr
93 << __FILE__ << "." << G4endl;
94 break;
95 }
96 v1 = 0.5;
97 v2 = 0.5;
98 result = 2. * G4UniformRand() - 1.;
99 for (m_tmp = 0; m_tmp < theCoeff[low].GetN
100 l = m_tmp + 1;
101 G4double legend = theLeg.Evaluate(l, res
102 v1 += (2. * l + 1) / 2. * theCoeff[low].
103 }
104 for (m_tmp = 0; m_tmp < theCoeff[high].Get
105 l = m_tmp + 1;
106 G4double legend = theLeg.Evaluate(l, res
107 v2 += (2. * l + 1) / 2. * theCoeff[high]
108 }
109 // v1 = std::max(0.,v1); // Workaround in
110 // v2 = std::max(0.,v2);
111 value = theInt.Interpolate(theManager.GetS
112 random = G4UniformRand();
113 if (0 >= theNorm) break; // Workaround fo
114 } while (random > value / theNorm); // Loop
115
116 return result;
117 }
118
119 G4double G4ParticleHPLegendreStore::SampleMax(
120 {
121 G4double result = 0.;
122
123 G4int i0;
124 G4int low(0), high(0);
125 G4ParticleHPFastLegendre theLeg;
126 for (i0 = 0; i0 < nEnergy; i0++) {
127 high = i0;
128 if (theCoeff[i0].GetEnergy() > anEnergy) b
129 }
130 low = std::max(0, high - 1);
131 G4ParticleHPInterpolator theInt;
132 G4double x, x1, x2;
133 x = anEnergy;
134 x1 = theCoeff[low].GetEnergy();
135 x2 = theCoeff[high].GetEnergy();
136 G4double theNorm = 0;
137 G4double try01 = 0, try02 = 0;
138 G4double max1, max2, costh;
139 max1 = 0;
140 max2 = 0;
141 G4int l;
142 for (i0 = 0; i0 < 601; i0++) {
143 costh = G4double(i0 - 300) / 300.;
144 try01 = 0;
145 for (l = 0; l < theCoeff[low].GetNumberOfP
146 try01 += (2. * l + 1) / 2. * theCoeff[lo
147 }
148 if (try01 > max1) max1 = try01;
149 try02 = 0;
150 for (l = 0; l < theCoeff[high].GetNumberOf
151 try02 += (2. * l + 1) / 2. * theCoeff[hi
152 }
153 if (try02 > max2) max2 = try02;
154 }
155 theNorm = theInt.Interpolate(theManager.GetS
156
157 G4double value, random;
158 G4double v1, v2;
159 G4int icounter = 0;
160 G4int icounter_max = 1024;
161 do {
162 icounter++;
163 if (icounter > icounter_max) {
164 G4cout << "Loop-counter exceeded the thr
165 << __FILE__ << "." << G4endl;
166 break;
167 }
168 v1 = 0;
169 v2 = 0;
170 result = 2. * G4UniformRand() - 1.;
171 for (l = 0; l < theCoeff[low].GetNumberOfP
172 G4double legend = theLeg.Evaluate(l, res
173 v1 += (2. * l + 1) / 2. * theCoeff[low].
174 }
175 for (l = 0; l < theCoeff[high].GetNumberOf
176 G4double legend = theLeg.Evaluate(l, res
177 v2 += (2. * l + 1) / 2. * theCoeff[high]
178 }
179 v1 = std::max(0., v1); // Workaround in c
180 v2 = std::max(0., v2);
181 value = theInt.Interpolate(theManager.GetS
182 random = G4UniformRand();
183 if (0 >= theNorm) break; // Workaround fo
184 } while (random > value / theNorm); // Loop
185 return result;
186 }
187
188 G4double G4ParticleHPLegendreStore::SampleElas
189 {
190 G4double result = 0.;
191
192 G4int i0;
193 G4int low(0), high(0);
194 G4ParticleHPFastLegendre theLeg;
195 for (i0 = 0; i0 < nEnergy; i0++) {
196 high = i0;
197 if (theCoeff[i0].GetEnergy() > anEnergy) b
198 }
199 low = std::max(0, high - 1);
200 G4ParticleHPInterpolator theInt;
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, tr
207 G4double try1, try2;
208 G4int l;
209 for (l = 0; l < theCoeff[low].GetNumberOfPol
210 try01 += (2. * l + 1) / 2. * theCoeff[low]
211 try11 += (2. * l + 1) / 2. * theCoeff[low]
212 }
213 for (l = 0; l < theCoeff[high].GetNumberOfPo
214 try02 += (2. * l + 1) / 2. * theCoeff[high
215 try12 += (2. * l + 1) / 2. * theCoeff[high
216 }
217 try1 = theInt.Interpolate(theManager.GetSche
218 try2 = theInt.Interpolate(theManager.GetSche
219 theNorm = std::max(try1, try2);
220
221 G4double value, random;
222 G4double v1, v2;
223 G4int icounter = 0;
224 G4int icounter_max = 1024;
225 do {
226 icounter++;
227 if (icounter > icounter_max) {
228 G4cout << "Loop-counter exceeded the thr
229 << __FILE__ << "." << G4endl;
230 break;
231 }
232 v1 = 0;
233 v2 = 0;
234 result = 2. * G4UniformRand() - 1.;
235 for (l = 0; l < theCoeff[low].GetNumberOfP
236 G4double legend = theLeg.Evaluate(l, res
237 v1 += (2. * l + 1) / 2. * theCoeff[low].
238 }
239 for (l = 0; l < theCoeff[high].GetNumberOf
240 G4double legend = theLeg.Evaluate(l, res
241 v2 += (2. * l + 1) / 2. * theCoeff[high]
242 }
243 value = theInt.Interpolate(theManager.GetS
244 random = G4UniformRand();
245 } while (random > value / theNorm); // Loop
246
247 return result;
248 }
249
250 G4double G4ParticleHPLegendreStore::Sample(G4d
251 {
252 G4int i0;
253 G4int low(0), high(0);
254 // G4cout << "G4ParticleHPLegendreStore::Sa
255 for (i0 = 0; i0 < nEnergy; i0++) {
256 // G4cout <<"theCoeff["<<i0<<"].GetEne
257 high = i0;
258 if (theCoeff[i0].GetEnergy() > energy) bre
259 }
260 low = std::max(0, high - 1);
261 // G4cout << "G4ParticleHPLegendreStore::Sa
262 G4ParticleHPVector theBuffer;
263 G4ParticleHPInterpolator theInt;
264 G4double x1, x2, y1, y2, y;
265 x1 = theCoeff[low].GetEnergy();
266 x2 = theCoeff[high].GetEnergy();
267 // G4cout << "the xes "<<x1<<" "<<x2<<G4end
268 G4double costh = 0;
269 for (i0 = 0; i0 < 601; i0++) {
270 costh = G4double(i0 - 300) / 300.;
271 y1 = Integrate(low, costh);
272 y2 = Integrate(high, costh);
273 y = theInt.Interpolate(theManager.GetSchem
274 theBuffer.SetData(i0, costh, y);
275 // G4cout << "Integration "<<low<<" "<
276 }
277 G4double rand = G4UniformRand();
278 G4int it;
279 for (i0 = 1; i0 < 601; i0++) {
280 it = i0;
281 if (rand < theBuffer.GetY(i0) / theBuffer.
282 // G4cout <<"sampling now "<<i0<<" "
283 // << theBuffer.GetY(i0)<<" "
284 // << theBuffer.GetY(600)<<" "
285 // << rand<<" "
286 // << theBuffer.GetY(i0)/theBuffer
287 }
288 if (it == 601) it = 600;
289 // G4cout << "G4ParticleHPLegendreStore::Sa
290 G4double norm = theBuffer.GetY(600);
291 if (norm == 0) return -DBL_MAX;
292 x1 = theBuffer.GetY(it) / norm;
293 x2 = theBuffer.GetY(it - 1) / norm;
294 y1 = theBuffer.GetX(it);
295 y2 = theBuffer.GetX(it - 1);
296 // G4cout << "G4ParticleHPLegendreStore::Sa
297 return theInt.Interpolate(theManager.GetSche
298 }
299
300 G4double
301 G4ParticleHPLegendreStore::Integrate(G4int k,
302 G4double
303 {
304 G4double result = 0.;
305 G4ParticleHPFastLegendre theLeg;
306 // G4cout <<"the COEFFS "<<k<<" ";
307 // G4cout <<theCoeff[k].GetNumberOfPoly()<<
308 for (G4int l = 0; l < theCoeff[k].GetNumberO
309 result += theCoeff[k].GetCoeff(l) * theLeg
310 // G4cout << theCoeff[k].GetCoeff(l)<<"
311 }
312 // G4cout <<G4endl;
313 return result;
314 }
315