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