| 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 /* 26 /*
27 * File: G4FPYSamplingOps.cc 27 * File: G4FPYSamplingOps.cc
28 * Author: B. Wendt (wendbryc@isu.edu) 28 * Author: B. Wendt (wendbryc@isu.edu)
29 * 29 *
30 * Created on June 30, 2011, 11:10 AM 30 * Created on June 30, 2011, 11:10 AM
31 */ 31 */
32 32
33 #include "G4FPYSamplingOps.hh" 33 #include "G4FPYSamplingOps.hh"
34 34
35 #include "G4FFGDebuggingMacros.hh" 35 #include "G4FFGDebuggingMacros.hh"
36 #include "G4FFGDefaultValues.hh" 36 #include "G4FFGDefaultValues.hh"
37 #include "G4FFGEnumerations.hh" 37 #include "G4FFGEnumerations.hh"
38 #include "G4Log.hh" 38 #include "G4Log.hh"
39 #include "G4Pow.hh" 39 #include "G4Pow.hh"
40 #include "G4ShiftedGaussian.hh" 40 #include "G4ShiftedGaussian.hh"
41 #include "G4WattFissionSpectrumValues.hh" 41 #include "G4WattFissionSpectrumValues.hh"
42 #include "Randomize.hh" 42 #include "Randomize.hh"
43 #include "globals.hh" 43 #include "globals.hh"
44 44
45 #include <CLHEP/Random/Stat.h> 45 #include <CLHEP/Random/Stat.h>
46 46
47 #include <iostream> 47 #include <iostream>
48 48
49 G4FPYSamplingOps::G4FPYSamplingOps() 49 G4FPYSamplingOps::G4FPYSamplingOps()
50 { 50 {
51 // Set the default verbosity 51 // Set the default verbosity
52 Verbosity_ = G4FFGDefaultValues::Verbosity; 52 Verbosity_ = G4FFGDefaultValues::Verbosity;
53 53
54 // Initialize the class 54 // Initialize the class
55 Initialize(); 55 Initialize();
56 } 56 }
57 57
58 G4FPYSamplingOps::G4FPYSamplingOps(G4int Verbo 58 G4FPYSamplingOps::G4FPYSamplingOps(G4int Verbosity)
59 { 59 {
60 // Set the default verbosity 60 // Set the default verbosity
61 Verbosity_ = Verbosity; 61 Verbosity_ = Verbosity;
62 62
63 // Initialize the class 63 // Initialize the class
64 Initialize(); 64 Initialize();
65 } 65 }
66 66
67 void G4FPYSamplingOps::Initialize() 67 void G4FPYSamplingOps::Initialize()
68 { 68 {
69 G4FFG_FUNCTIONENTER__ 69 G4FFG_FUNCTIONENTER__
70 70
71 // Get the pointer to the random number gene 71 // Get the pointer to the random number generator
72 // RandomEngine_ = CLHEP::HepRandom::getTheE 72 // RandomEngine_ = CLHEP::HepRandom::getTheEngine();
73 RandomEngine_ = G4Random::getTheEngine(); 73 RandomEngine_ = G4Random::getTheEngine();
74 74
75 // Initialize the data members 75 // Initialize the data members
76 ShiftedGaussianValues_ = new G4ShiftedGaussi 76 ShiftedGaussianValues_ = new G4ShiftedGaussian;
77 Mean_ = 0; 77 Mean_ = 0;
78 StdDev_ = 0; 78 StdDev_ = 0;
79 NextGaussianIsStoredInMemory_ = FALSE; 79 NextGaussianIsStoredInMemory_ = FALSE;
80 GaussianOne_ = 0; 80 GaussianOne_ = 0;
81 GaussianTwo_ = 0; 81 GaussianTwo_ = 0;
82 Tolerance_ = 0.000001; 82 Tolerance_ = 0.000001;
83 WattConstants_ = new WattSpectrumConstants; 83 WattConstants_ = new WattSpectrumConstants;
84 WattConstants_->Product = 0; 84 WattConstants_->Product = 0;
85 85
86 G4FFG_FUNCTIONLEAVE__ 86 G4FFG_FUNCTIONLEAVE__
87 } 87 }
88 88
89 G4double G4FPYSamplingOps::G4SampleGaussian(G4 89 G4double G4FPYSamplingOps::G4SampleGaussian(G4double Mean, G4double StdDev)
90 { 90 {
91 G4FFG_SAMPLING_FUNCTIONENTER__ 91 G4FFG_SAMPLING_FUNCTIONENTER__
92 92
93 // Determine if the parameters have changed 93 // Determine if the parameters have changed
94 G4bool ParametersChanged = (Mean_ != Mean || 94 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
95 if (static_cast<int>(ParametersChanged) == T 95 if (static_cast<int>(ParametersChanged) == TRUE) {
96 // Set the new values if the parameters ha 96 // Set the new values if the parameters have changed
97 NextGaussianIsStoredInMemory_ = FALSE; 97 NextGaussianIsStoredInMemory_ = FALSE;
98 98
99 Mean_ = Mean; 99 Mean_ = Mean;
100 StdDev_ = StdDev; 100 StdDev_ = StdDev;
101 } 101 }
102 102
103 G4double Sample = SampleGaussian(); 103 G4double Sample = SampleGaussian();
104 104
105 G4FFG_SAMPLING_FUNCTIONLEAVE__ 105 G4FFG_SAMPLING_FUNCTIONLEAVE__
106 return Sample; 106 return Sample;
107 } 107 }
108 108
109 G4double G4FPYSamplingOps::G4SampleGaussian(G4 109 G4double G4FPYSamplingOps::G4SampleGaussian(G4double Mean, G4double StdDev,
110 G4 110 G4FFGEnumerations::GaussianRange Range)
111 { 111 {
112 G4FFG_SAMPLING_FUNCTIONENTER__ 112 G4FFG_SAMPLING_FUNCTIONENTER__
113 113
114 if (Range == G4FFGEnumerations::ALL) { 114 if (Range == G4FFGEnumerations::ALL) {
115 // Call the overloaded function 115 // Call the overloaded function
116 G4double Sample = G4SampleGaussian(Mean, S 116 G4double Sample = G4SampleGaussian(Mean, StdDev);
117 117
118 G4FFG_SAMPLING_FUNCTIONLEAVE__ 118 G4FFG_SAMPLING_FUNCTIONLEAVE__
119 return Sample; 119 return Sample;
120 } 120 }
121 121
122 // Determine if the parameters have changed 122 // Determine if the parameters have changed
123 G4bool ParametersChanged = (Mean_ != Mean || 123 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
124 if (static_cast<int>(ParametersChanged) == T 124 if (static_cast<int>(ParametersChanged) == TRUE) {
125 if (Mean <= 0) { 125 if (Mean <= 0) {
126 std::ostringstream Temp; 126 std::ostringstream Temp;
127 Temp << "Mean value of " << Mean << " ou 127 Temp << "Mean value of " << Mean << " out of range";
128 G4Exception("G4FPYGaussianOps::G4SampleI 128 G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()", Temp.str().c_str(), JustWarning,
129 "A value of '0' will be used 129 "A value of '0' will be used instead.");
130 130
131 G4FFG_SAMPLING_FUNCTIONLEAVE__ 131 G4FFG_SAMPLING_FUNCTIONLEAVE__
132 return 0; 132 return 0;
133 } 133 }
134 134
135 // Set the new values if the parameters ha 135 // Set the new values if the parameters have changed and then perform
136 // the shift 136 // the shift
137 Mean_ = Mean; 137 Mean_ = Mean;
138 StdDev_ = StdDev; 138 StdDev_ = StdDev;
139 139
140 ShiftParameters(G4FFGEnumerations::DOUBLE) 140 ShiftParameters(G4FFGEnumerations::DOUBLE);
141 } 141 }
142 142
143 // Sample the Gaussian distribution 143 // Sample the Gaussian distribution
144 G4double Rand; 144 G4double Rand;
145 do { 145 do {
146 Rand = SampleGaussian(); 146 Rand = SampleGaussian();
147 } while (Rand < 0); // Loop checking, 11.05 147 } while (Rand < 0); // Loop checking, 11.05.2015, T. Koi
148 148
149 G4FFG_SAMPLING_FUNCTIONLEAVE__ 149 G4FFG_SAMPLING_FUNCTIONLEAVE__
150 return Rand; 150 return Rand;
151 } 151 }
152 152
153 G4int G4FPYSamplingOps::G4SampleIntegerGaussia 153 G4int G4FPYSamplingOps::G4SampleIntegerGaussian(G4double Mean, G4double StdDev)
154 { 154 {
155 G4FFG_SAMPLING_FUNCTIONENTER__ 155 G4FFG_SAMPLING_FUNCTIONENTER__
156 156
157 // Determine if the parameters have changed 157 // Determine if the parameters have changed
158 G4bool ParametersChanged = (Mean_ != Mean || 158 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
159 if (static_cast<int>(ParametersChanged) == T 159 if (static_cast<int>(ParametersChanged) == TRUE) {
160 // Set the new values if the parameters ha 160 // Set the new values if the parameters have changed
161 NextGaussianIsStoredInMemory_ = FALSE; 161 NextGaussianIsStoredInMemory_ = FALSE;
162 162
163 Mean_ = Mean; 163 Mean_ = Mean;
164 StdDev_ = StdDev; 164 StdDev_ = StdDev;
165 } 165 }
166 166
167 // Return the sample integer value 167 // Return the sample integer value
168 auto Sample = (G4int)(std::floor(SampleGauss 168 auto Sample = (G4int)(std::floor(SampleGaussian()));
169 169
170 G4FFG_SAMPLING_FUNCTIONLEAVE__ 170 G4FFG_SAMPLING_FUNCTIONLEAVE__
171 return Sample; 171 return Sample;
172 } 172 }
173 173
174 G4int G4FPYSamplingOps::G4SampleIntegerGaussia 174 G4int G4FPYSamplingOps::G4SampleIntegerGaussian(G4double Mean, G4double StdDev,
175 175 G4FFGEnumerations::GaussianRange Range)
176 { 176 {
177 G4FFG_SAMPLING_FUNCTIONENTER__ 177 G4FFG_SAMPLING_FUNCTIONENTER__
178 178
179 if (Range == G4FFGEnumerations::ALL) { 179 if (Range == G4FFGEnumerations::ALL) {
180 // Call the overloaded function 180 // Call the overloaded function
181 G4int Sample = G4SampleIntegerGaussian(Mea 181 G4int Sample = G4SampleIntegerGaussian(Mean, StdDev);
182 182
183 G4FFG_SAMPLING_FUNCTIONLEAVE__ 183 G4FFG_SAMPLING_FUNCTIONLEAVE__
184 return Sample; 184 return Sample;
185 } 185 }
186 // ParameterShift() locks up if the mean is 186 // ParameterShift() locks up if the mean is less than 1.
187 std::ostringstream Temp; 187 std::ostringstream Temp;
188 if (Mean < 1) { 188 if (Mean < 1) {
189 // Temp << "Mean value of " << Mean << 189 // Temp << "Mean value of " << Mean << " out of range";
190 // G4Exception("G4FPYGaussianOps::G4Sam 190 // G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()",
191 // Temp.str().c_str(), 191 // Temp.str().c_str(),
192 // JustWarning, 192 // JustWarning,
193 // "A value of '0' will be 193 // "A value of '0' will be used instead.");
194 194
195 // return 0; 195 // return 0;
196 } 196 }
197 197
198 if (Mean / StdDev < 2) { 198 if (Mean / StdDev < 2) {
199 // Temp << "Non-ideal conditions:\tMean:" 199 // Temp << "Non-ideal conditions:\tMean:" << Mean << "\tStdDev: "
200 // << StdDev; 200 // << StdDev;
201 // G4Exception("G4FPYGaussianOps::G4Sample 201 // G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()",
202 // Temp.str().c_str(), 202 // Temp.str().c_str(),
203 // JustWarning, 203 // JustWarning,
204 // "Results not guaranteed: Co 204 // "Results not guaranteed: Consider tightening the standard deviation");
205 } 205 }
206 206
207 // Determine if the parameters have changed 207 // Determine if the parameters have changed
208 G4bool ParametersChanged = (Mean_ != Mean || 208 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
209 if (static_cast<int>(ParametersChanged) == T 209 if (static_cast<int>(ParametersChanged) == TRUE) {
210 // Set the new values if the parameters ha 210 // Set the new values if the parameters have changed and then perform
211 // the shift 211 // the shift
212 Mean_ = Mean; 212 Mean_ = Mean;
213 StdDev_ = StdDev; 213 StdDev_ = StdDev;
214 214
215 ShiftParameters(G4FFGEnumerations::INT); 215 ShiftParameters(G4FFGEnumerations::INT);
216 } 216 }
217 217
218 G4int RandInt; 218 G4int RandInt;
219 // Sample the Gaussian distribution - only n 219 // Sample the Gaussian distribution - only non-negative values are
220 // accepted 220 // accepted
221 do { 221 do {
222 RandInt = (G4int)floor(SampleGaussian()); 222 RandInt = (G4int)floor(SampleGaussian());
223 } while (RandInt < 0); // Loop checking, 11 223 } while (RandInt < 0); // Loop checking, 11.05.2015, T. Koi
224 224
225 G4FFG_SAMPLING_FUNCTIONLEAVE__ 225 G4FFG_SAMPLING_FUNCTIONLEAVE__
226 return RandInt; 226 return RandInt;
227 } 227 }
228 228
229 G4double G4FPYSamplingOps::G4SampleUniform() 229 G4double G4FPYSamplingOps::G4SampleUniform()
230 { 230 {
231 G4FFG_SAMPLING_FUNCTIONENTER__ 231 G4FFG_SAMPLING_FUNCTIONENTER__
232 232
233 G4double Sample = RandomEngine_->flat(); 233 G4double Sample = RandomEngine_->flat();
234 234
235 G4FFG_SAMPLING_FUNCTIONLEAVE__ 235 G4FFG_SAMPLING_FUNCTIONLEAVE__
236 return Sample; 236 return Sample;
237 } 237 }
238 238
239 G4double G4FPYSamplingOps::G4SampleUniform(G4d 239 G4double G4FPYSamplingOps::G4SampleUniform(G4double Lower, G4double Upper)
240 { 240 {
241 G4FFG_SAMPLING_FUNCTIONENTER__ 241 G4FFG_SAMPLING_FUNCTIONENTER__
242 242
243 // Calculate the difference 243 // Calculate the difference
244 G4double Difference = Upper - Lower; 244 G4double Difference = Upper - Lower;
245 245
246 // Scale appropriately and return the value 246 // Scale appropriately and return the value
247 G4double Sample = G4SampleUniform() * Differ 247 G4double Sample = G4SampleUniform() * Difference + Lower;
248 248
249 G4FFG_SAMPLING_FUNCTIONLEAVE__ 249 G4FFG_SAMPLING_FUNCTIONLEAVE__
250 return Sample; 250 return Sample;
251 } 251 }
252 252
253 G4double G4FPYSamplingOps::G4SampleWatt(G4int 253 G4double G4FPYSamplingOps::G4SampleWatt(G4int WhatIsotope,
254 G4FFGE 254 G4FFGEnumerations::FissionCause WhatCause,
255 G4doub 255 G4double WhatEnergy)
256 { 256 {
257 G4FFG_SAMPLING_FUNCTIONENTER__ 257 G4FFG_SAMPLING_FUNCTIONENTER__
258 258
259 // Determine if the parameters are different 259 // Determine if the parameters are different
260 // TK modified 131108 260 // TK modified 131108
261 // if(WattConstants_->Product != WhatIsotope 261 // if(WattConstants_->Product != WhatIsotope
262 if (WattConstants_->Product != WhatIsotope / 262 if (WattConstants_->Product != WhatIsotope / 10 || WattConstants_->Cause != WhatCause
263 || WattConstants_->Energy != WhatEnergy) 263 || WattConstants_->Energy != WhatEnergy)
264 { 264 {
265 // If the parameters are different or have 265 // If the parameters are different or have not yet been defined then we
266 // need to re-evaluate the constants 266 // need to re-evaluate the constants
267 // TK modified 131108 267 // TK modified 131108
268 // WattConstants_->Product = WhatIsotope; 268 // WattConstants_->Product = WhatIsotope;
269 WattConstants_->Product = WhatIsotope / 10 269 WattConstants_->Product = WhatIsotope / 10;
270 WattConstants_->Cause = WhatCause; 270 WattConstants_->Cause = WhatCause;
271 WattConstants_->Energy = WhatEnergy; 271 WattConstants_->Energy = WhatEnergy;
272 272
273 EvaluateWattConstants(); 273 EvaluateWattConstants();
274 } 274 }
275 275
276 G4double X = -G4Log(G4SampleUniform()); 276 G4double X = -G4Log(G4SampleUniform());
277 G4double Y = -G4Log(G4SampleUniform()); 277 G4double Y = -G4Log(G4SampleUniform());
278 G4int icounter = 0; 278 G4int icounter = 0;
279 G4int icounter_max = 1024; 279 G4int icounter_max = 1024;
280 while (G4Pow::GetInstance()->powN(Y - WattCo 280 while (G4Pow::GetInstance()->powN(Y - WattConstants_->M * (X + 1), 2)
281 > WattConstants_->B * WattConstants_- 281 > WattConstants_->B * WattConstants_->L * X) // Loop checking, 11.05.2015, T. Koi
282 { 282 {
283 icounter++; 283 icounter++;
284 if (icounter > icounter_max) { 284 if (icounter > icounter_max) {
285 G4cout << "Loop-counter exceeded the thr 285 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
286 << __FILE__ << "." << G4endl; 286 << __FILE__ << "." << G4endl;
287 break; 287 break;
288 } 288 }
289 X = -G4Log(G4SampleUniform()); 289 X = -G4Log(G4SampleUniform());
290 Y = -G4Log(G4SampleUniform()); 290 Y = -G4Log(G4SampleUniform());
291 } 291 }
292 292
293 G4FFG_SAMPLING_FUNCTIONLEAVE__ 293 G4FFG_SAMPLING_FUNCTIONLEAVE__
294 return WattConstants_->L * X; 294 return WattConstants_->L * X;
295 } 295 }
296 296
297 void G4FPYSamplingOps::G4SetVerbosity(G4int Ve 297 void G4FPYSamplingOps::G4SetVerbosity(G4int Verbosity)
298 { 298 {
299 G4FFG_SAMPLING_FUNCTIONENTER__ 299 G4FFG_SAMPLING_FUNCTIONENTER__
300 300
301 Verbosity_ = Verbosity; 301 Verbosity_ = Verbosity;
302 302
303 ShiftedGaussianValues_->G4SetVerbosity(Verbo 303 ShiftedGaussianValues_->G4SetVerbosity(Verbosity_);
304 304
305 G4FFG_SAMPLING_FUNCTIONLEAVE__ 305 G4FFG_SAMPLING_FUNCTIONLEAVE__
306 } 306 }
307 307
308 G4bool G4FPYSamplingOps::CheckAndSetParameters 308 G4bool G4FPYSamplingOps::CheckAndSetParameters()
309 { 309 {
310 G4FFG_SAMPLING_FUNCTIONENTER__ 310 G4FFG_SAMPLING_FUNCTIONENTER__
311 311
312 G4double ShiftedMean = ShiftedGaussianValues 312 G4double ShiftedMean = ShiftedGaussianValues_->G4FindShiftedMean(Mean_, StdDev_);
313 if (ShiftedMean == 0) { 313 if (ShiftedMean == 0) {
314 G4FFG_SAMPLING_FUNCTIONLEAVE__ 314 G4FFG_SAMPLING_FUNCTIONLEAVE__
315 return FALSE; 315 return FALSE;
316 } 316 }
317 317
318 Mean_ = ShiftedMean; 318 Mean_ = ShiftedMean;
319 319
320 G4FFG_SAMPLING_FUNCTIONLEAVE__ 320 G4FFG_SAMPLING_FUNCTIONLEAVE__
321 return TRUE; 321 return TRUE;
322 } 322 }
323 323
324 void G4FPYSamplingOps::EvaluateWattConstants() 324 void G4FPYSamplingOps::EvaluateWattConstants()
325 { 325 {
326 G4FFG_SAMPLING_FUNCTIONENTER__ 326 G4FFG_SAMPLING_FUNCTIONENTER__
327 327
328 G4double A, K; 328 G4double A, K;
329 A = K = 0; 329 A = K = 0;
330 // Use the default values if IsotopeIndex is 330 // Use the default values if IsotopeIndex is not reset
331 G4int IsotopeIndex = 0; 331 G4int IsotopeIndex = 0;
332 332
333 if (WattConstants_->Cause == G4FFGEnumeratio 333 if (WattConstants_->Cause == G4FFGEnumerations::SPONTANEOUS) {
334 // Determine if the isotope requested exis 334 // Determine if the isotope requested exists in the lookup tables
335 for (G4int i = 0; SpontaneousWattIsotopesI 335 for (G4int i = 0; SpontaneousWattIsotopesIndex[i] != -1; i++) {
336 if (SpontaneousWattIsotopesIndex[i] == W 336 if (SpontaneousWattIsotopesIndex[i] == WattConstants_->Product) {
337 IsotopeIndex = i; 337 IsotopeIndex = i;
338 338
339 break; 339 break;
340 } 340 }
341 } 341 }
342 342
343 // Get A and B 343 // Get A and B
344 A = SpontaneousWattConstants[IsotopeIndex] 344 A = SpontaneousWattConstants[IsotopeIndex][0];
345 WattConstants_->B = SpontaneousWattConstan 345 WattConstants_->B = SpontaneousWattConstants[IsotopeIndex][1];
346 } 346 }
347 else if (WattConstants_->Cause == G4FFGEnume 347 else if (WattConstants_->Cause == G4FFGEnumerations::NEUTRON_INDUCED) {
348 // Determine if the isotope requested exis 348 // Determine if the isotope requested exists in the lookup tables
349 for (G4int i = 0; NeutronInducedWattIsotop 349 for (G4int i = 0; NeutronInducedWattIsotopesIndex[i] != -1; i++) {
350 if (NeutronInducedWattIsotopesIndex[i] = 350 if (NeutronInducedWattIsotopesIndex[i] == WattConstants_->Product) {
351 IsotopeIndex = i; 351 IsotopeIndex = i;
352 break; 352 break;
353 } 353 }
354 } 354 }
355 355
356 // Determine the Watt fission spectrum con 356 // Determine the Watt fission spectrum constants based on the energy of
357 // the incident neutron 357 // the incident neutron
358 if (WattConstants_->Energy == G4FFGDefault 358 if (WattConstants_->Energy == G4FFGDefaultValues::ThermalNeutronEnergy) {
359 A = NeutronInducedWattConstants[IsotopeI 359 A = NeutronInducedWattConstants[IsotopeIndex][0][0];
360 WattConstants_->B = NeutronInducedWattCo 360 WattConstants_->B = NeutronInducedWattConstants[IsotopeIndex][0][1];
361 } 361 }
362 else if (WattConstants_->Energy > 14.0 * C 362 else if (WattConstants_->Energy > 14.0 * CLHEP::MeV) {
363 G4Exception("G4FPYSamplingOps::G4SampleW 363 G4Exception("G4FPYSamplingOps::G4SampleWatt()",
364 "Incident neutron energy abo 364 "Incident neutron energy above 14 MeV requested.", JustWarning,
365 "Using Watt fission constant 365 "Using Watt fission constants for 14 Mev.");
366 366
367 A = NeutronInducedWattConstants[IsotopeI 367 A = NeutronInducedWattConstants[IsotopeIndex][2][0];
368 WattConstants_->B = NeutronInducedWattCo 368 WattConstants_->B = NeutronInducedWattConstants[IsotopeIndex][2][1];
369 } 369 }
370 else { 370 else {
371 G4int EnergyIndex = 0; 371 G4int EnergyIndex = 0;
372 G4double EnergyDifference = 0; 372 G4double EnergyDifference = 0;
373 G4double RangeDifference, ConstantDiffer 373 G4double RangeDifference, ConstantDifference;
374 374
375 for (G4int i = 1; IncidentEnergyBins[i] 375 for (G4int i = 1; IncidentEnergyBins[i] != -1; i++) {
376 if (WattConstants_->Energy <= Incident 376 if (WattConstants_->Energy <= IncidentEnergyBins[i]) {
377 EnergyIndex = i; 377 EnergyIndex = i;
378 EnergyDifference = IncidentEnergyBin 378 EnergyDifference = IncidentEnergyBins[EnergyIndex] - WattConstants_->Energy;
379 if (EnergyDifference != 0) { 379 if (EnergyDifference != 0) {
380 std::ostringstream Temp; 380 std::ostringstream Temp;
381 Temp << "Incident neutron energy o 381 Temp << "Incident neutron energy of ";
382 Temp << WattConstants_->Energy << 382 Temp << WattConstants_->Energy << " MeV is not ";
383 Temp << "explicitly listed in the 383 Temp << "explicitly listed in the data tables";
384 // G4Except 384 // G4Exception("G4FPYSamplingOps::G4SampleWatt()",
385 // 385 // Temp.str().c_str(),
386 // 386 // JustWarning,
387 // 387 // "Using linear interpolation.");
388 } 388 }
389 break; 389 break;
390 } 390 }
391 } 391 }
392 392
393 RangeDifference = IncidentEnergyBins[Ene 393 RangeDifference = IncidentEnergyBins[EnergyIndex] - IncidentEnergyBins[EnergyIndex - 1];
394 394
395 // Interpolate the value for 'a' 395 // Interpolate the value for 'a'
396 ConstantDifference = NeutronInducedWattC 396 ConstantDifference = NeutronInducedWattConstants[IsotopeIndex][EnergyIndex][0]
397 - NeutronInducedWat 397 - NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][0];
398 A = (EnergyDifference / RangeDifference) 398 A = (EnergyDifference / RangeDifference) * ConstantDifference
399 + NeutronInducedWattConstants[Isotop 399 + NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][0];
400 400
401 // Interpolate the value for 'b' 401 // Interpolate the value for 'b'
402 ConstantDifference = NeutronInducedWattC 402 ConstantDifference = NeutronInducedWattConstants[IsotopeIndex][EnergyIndex][1]
403 - NeutronInducedWat 403 - NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][1];
404 WattConstants_->B = (EnergyDifference / 404 WattConstants_->B = (EnergyDifference / RangeDifference) * ConstantDifference
405 + NeutronInducedWatt 405 + NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][1];
406 } 406 }
407 } 407 }
408 else { 408 else {
409 // Produce an error since an unsupported f 409 // Produce an error since an unsupported fission type was requested and
410 // abort the current run 410 // abort the current run
411 G4String Temp = "Watt fission spectra data 411 G4String Temp = "Watt fission spectra data not available for ";
412 if (WattConstants_->Cause == G4FFGEnumerat 412 if (WattConstants_->Cause == G4FFGEnumerations::PROTON_INDUCED) {
413 Temp += "proton induced fission."; 413 Temp += "proton induced fission.";
414 } 414 }
415 else if (WattConstants_->Cause == G4FFGEnu 415 else if (WattConstants_->Cause == G4FFGEnumerations::GAMMA_INDUCED) {
416 Temp += "gamma induced fission."; 416 Temp += "gamma induced fission.";
417 } 417 }
418 else { 418 else {
419 Temp += "!Warning! unknown cause."; 419 Temp += "!Warning! unknown cause.";
420 } 420 }
421 G4Exception("G4FPYSamplingOps::G4SampleWat 421 G4Exception("G4FPYSamplingOps::G4SampleWatt()", Temp, RunMustBeAborted,
422 "Fission events will not be sa 422 "Fission events will not be sampled in this run.");
423 } 423 }
424 424
425 // Calculate the constants 425 // Calculate the constants
426 K = 1 + (WattConstants_->B / (8.0 * A)); 426 K = 1 + (WattConstants_->B / (8.0 * A));
427 WattConstants_->L = (K + G4Pow::GetInstance( 427 WattConstants_->L = (K + G4Pow::GetInstance()->powA((K * K - 1), 0.5)) / A;
428 WattConstants_->M = A * WattConstants_->L - 428 WattConstants_->M = A * WattConstants_->L - 1;
429 429
430 G4FFG_SAMPLING_FUNCTIONLEAVE__ 430 G4FFG_SAMPLING_FUNCTIONLEAVE__
431 } 431 }
432 432
433 G4double G4FPYSamplingOps::SampleGaussian() 433 G4double G4FPYSamplingOps::SampleGaussian()
434 { 434 {
435 G4FFG_SAMPLING_FUNCTIONENTER__ 435 G4FFG_SAMPLING_FUNCTIONENTER__
436 436
437 if (static_cast<int>(NextGaussianIsStoredInM 437 if (static_cast<int>(NextGaussianIsStoredInMemory_) == TRUE) {
438 NextGaussianIsStoredInMemory_ = FALSE; 438 NextGaussianIsStoredInMemory_ = FALSE;
439 439
440 G4FFG_SAMPLING_FUNCTIONLEAVE__ 440 G4FFG_SAMPLING_FUNCTIONLEAVE__
441 return GaussianTwo_; 441 return GaussianTwo_;
442 } 442 }
443 443
444 // Define the variables to be used 444 // Define the variables to be used
445 G4double Radius; 445 G4double Radius;
446 G4double MappingFactor; 446 G4double MappingFactor;
447 447
448 // Sample from the unit circle (21.4% reject 448 // Sample from the unit circle (21.4% rejection probability)
449 do { 449 do {
450 GaussianOne_ = 2.0 * G4SampleUniform() - 1 450 GaussianOne_ = 2.0 * G4SampleUniform() - 1.0;
451 GaussianTwo_ = 2.0 * G4SampleUniform() - 1 451 GaussianTwo_ = 2.0 * G4SampleUniform() - 1.0;
452 Radius = GaussianOne_ * GaussianOne_ + Gau 452 Radius = GaussianOne_ * GaussianOne_ + GaussianTwo_ * GaussianTwo_;
453 } while (Radius > 1.0); // Loop checking, 1 453 } while (Radius > 1.0); // Loop checking, 11.05.2015, T. Koi
454 454
455 // Translate the values to Gaussian space 455 // Translate the values to Gaussian space
456 MappingFactor = std::sqrt(-2.0 * G4Log(Radiu 456 MappingFactor = std::sqrt(-2.0 * G4Log(Radius) / Radius) * StdDev_;
457 GaussianOne_ = Mean_ + GaussianOne_ * Mappin 457 GaussianOne_ = Mean_ + GaussianOne_ * MappingFactor;
458 GaussianTwo_ = Mean_ + GaussianTwo_ * Mappin 458 GaussianTwo_ = Mean_ + GaussianTwo_ * MappingFactor;
459 459
460 // Set the flag that a value is now stored i 460 // Set the flag that a value is now stored in memory
461 NextGaussianIsStoredInMemory_ = TRUE; 461 NextGaussianIsStoredInMemory_ = TRUE;
462 462
463 G4FFG_SAMPLING_FUNCTIONLEAVE__ 463 G4FFG_SAMPLING_FUNCTIONLEAVE__
464 return GaussianOne_; 464 return GaussianOne_;
465 } 465 }
466 466
467 void G4FPYSamplingOps::ShiftParameters(G4FFGEn 467 void G4FPYSamplingOps::ShiftParameters(G4FFGEnumerations::GaussianReturnType Type)
468 { 468 {
469 G4FFG_SAMPLING_FUNCTIONENTER__ 469 G4FFG_SAMPLING_FUNCTIONENTER__
470 470
471 // Set the flag that any second value stored 471 // Set the flag that any second value stored is no longer valid
472 NextGaussianIsStoredInMemory_ = FALSE; 472 NextGaussianIsStoredInMemory_ = FALSE;
473 473
474 // Check if the requested parameters have al 474 // Check if the requested parameters have already been calculated
475 if (static_cast<int>(CheckAndSetParameters() 475 if (static_cast<int>(CheckAndSetParameters()) == TRUE) {
476 G4FFG_SAMPLING_FUNCTIONLEAVE__ 476 G4FFG_SAMPLING_FUNCTIONLEAVE__
477 return; 477 return;
478 } 478 }
479 479
480 // If the requested type is INT, then perfor 480 // If the requested type is INT, then perform an iterative refinement of the
481 // mean that is going to be sampled 481 // mean that is going to be sampled
482 if (Type == G4FFGEnumerations::INT) { 482 if (Type == G4FFGEnumerations::INT) {
483 // Return if the mean is greater than 7 st 483 // Return if the mean is greater than 7 standard deviations away from 0
484 // since there is essentially 0 probabilit 484 // since there is essentially 0 probability that a sampled number will
485 // be less than 0 485 // be less than 0
486 if (Mean_ > 7 * StdDev_) { 486 if (Mean_ > 7 * StdDev_) {
487 G4FFG_SAMPLING_FUNCTIONLEAVE__ 487 G4FFG_SAMPLING_FUNCTIONLEAVE__
488 return; 488 return;
489 } 489 }
490 // Variables that contain the area and wei 490 // Variables that contain the area and weighted area information for
491 // calculating the statistical mean of the 491 // calculating the statistical mean of the Gaussian distribution when
492 // converted to a step function 492 // converted to a step function
493 G4double ErfContainer, AdjustedErfContaine 493 G4double ErfContainer, AdjustedErfContainer, Container;
494 494
495 // Variables that store lower and upper bo 495 // Variables that store lower and upper bounds information
496 G4double LowErf, HighErf; 496 G4double LowErf, HighErf;
497 497
498 // Values for determining the convergence 498 // Values for determining the convergence of the solution
499 G4double AdjMean = Mean_; 499 G4double AdjMean = Mean_;
500 G4double Delta = 1.0; 500 G4double Delta = 1.0;
501 G4bool HalfDelta = false; 501 G4bool HalfDelta = false;
502 G4bool ToleranceCheck = false; 502 G4bool ToleranceCheck = false;
503 503
504 // Normalization constant for use with erf 504 // Normalization constant for use with erf()
505 const G4double Normalization = StdDev_ * s 505 const G4double Normalization = StdDev_ * std::sqrt(2.0);
506 506
507 // Determine the upper limit of the estima 507 // Determine the upper limit of the estimates
508 const auto UpperLimit = (G4int)std::ceil(M 508 const auto UpperLimit = (G4int)std::ceil(Mean_ + 9 * StdDev_);
509 509
510 // Calculate the integral of the Gaussian 510 // Calculate the integral of the Gaussian distribution
511 511
512 G4int icounter = 0; 512 G4int icounter = 0;
513 G4int icounter_max = 1024; 513 G4int icounter_max = 1024;
514 do { 514 do {
515 icounter++; 515 icounter++;
516 if (icounter > icounter_max) { 516 if (icounter > icounter_max) {
517 G4cout << "Loop-counter exceeded the t 517 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
518 << __FILE__ << "." << G4endl; 518 << __FILE__ << "." << G4endl;
519 break; 519 break;
520 } 520 }
521 // Set the variables 521 // Set the variables
522 ErfContainer = 0; 522 ErfContainer = 0;
523 AdjustedErfContainer = 0; 523 AdjustedErfContainer = 0;
524 524
525 // Calculate the area and weighted area 525 // Calculate the area and weighted area
526 for (G4int i = 0; i <= UpperLimit; i++) 526 for (G4int i = 0; i <= UpperLimit; i++) {
527 // Determine the lower and upper bound 527 // Determine the lower and upper bounds
528 LowErf = ((AdjMean - i) / Normalizatio 528 LowErf = ((AdjMean - i) / Normalization);
529 HighErf = ((AdjMean - (i + 1.0)) / Nor 529 HighErf = ((AdjMean - (i + 1.0)) / Normalization);
530 530
531 // Correct the bounds for how they lie 531 // Correct the bounds for how they lie on the x-axis with
532 // respect to the mean 532 // respect to the mean
533 if (LowErf <= 0) { 533 if (LowErf <= 0) {
534 LowErf *= -1; 534 LowErf *= -1;
535 HighErf *= -1; 535 HighErf *= -1;
536 // Container = (erf(HighErf) - erf(LowErf))/2. 536 // Container = (erf(HighErf) - erf(LowErf))/2.0;
537 #if defined WIN32 - VC 537 #if defined WIN32 - VC
538 Container = (CLHEP::HepStat::erf(Hig 538 Container = (CLHEP::HepStat::erf(HighErf) - CLHEP::HepStat::erf(LowErf)) / 2.0;
539 #else 539 #else
540 Container = (erf(HighErf) - erf(LowE 540 Container = (erf(HighErf) - erf(LowErf)) / 2.0;
541 #endif 541 #endif
542 } 542 }
543 else if (HighErf < 0) { 543 else if (HighErf < 0) {
544 HighErf *= -1; 544 HighErf *= -1;
545 545
546 // Container = (erf(HighErf) + erf(LowErf))/2. 546 // Container = (erf(HighErf) + erf(LowErf))/2.0;
547 #if defined WIN32 - VC 547 #if defined WIN32 - VC
548 Container = (CLHEP::HepStat::erf(Hig 548 Container = (CLHEP::HepStat::erf(HighErf) + CLHEP::HepStat::erf(LowErf)) / 2.0;
549 #else 549 #else
550 Container = (erf(HighErf) + erf(LowE 550 Container = (erf(HighErf) + erf(LowErf)) / 2.0;
551 #endif 551 #endif
552 } 552 }
553 else { 553 else {
554 // Container = (erf(LowErf) - erf(HighErf))/2. 554 // Container = (erf(LowErf) - erf(HighErf))/2.0;
555 #if defined WIN32 - VC 555 #if defined WIN32 - VC
556 Container = (CLHEP::HepStat::erf(Low 556 Container = (CLHEP::HepStat::erf(LowErf) - CLHEP::HepStat::erf(HighErf)) / 2.0;
557 #else 557 #else
558 Container = (erf(LowErf) - erf(HighE 558 Container = (erf(LowErf) - erf(HighErf)) / 2.0;
559 #endif 559 #endif
560 } 560 }
561 561
562 #if defined WIN32 - VC 562 #if defined WIN32 - VC
563 // TK modified to avoid problem caused 563 // TK modified to avoid problem caused by QNAN
564 if (Container != Container) Container 564 if (Container != Container) Container = 0;
565 #endif 565 #endif
566 566
567 // Add up the weighted area 567 // Add up the weighted area
568 ErfContainer += Container; 568 ErfContainer += Container;
569 AdjustedErfContainer += Container * i; 569 AdjustedErfContainer += Container * i;
570 } 570 }
571 571
572 // Calculate the statistical mean 572 // Calculate the statistical mean
573 Container = AdjustedErfContainer / ErfCo 573 Container = AdjustedErfContainer / ErfContainer;
574 574
575 // Is it close enough to what we want? 575 // Is it close enough to what we want?
576 ToleranceCheck = (std::fabs(Mean_ - Cont 576 ToleranceCheck = (std::fabs(Mean_ - Container) < Tolerance_);
577 if (static_cast<int>(ToleranceCheck) == 577 if (static_cast<int>(ToleranceCheck) == TRUE) {
578 break; 578 break;
579 } 579 }
580 580
581 // Determine the step size 581 // Determine the step size
582 if (static_cast<int>(HalfDelta) == TRUE) 582 if (static_cast<int>(HalfDelta) == TRUE) {
583 Delta /= 2; 583 Delta /= 2;
584 } 584 }
585 585
586 // Step in the appropriate direction 586 // Step in the appropriate direction
587 if (Container > Mean_) { 587 if (Container > Mean_) {
588 AdjMean -= Delta; 588 AdjMean -= Delta;
589 } 589 }
590 else { 590 else {
591 HalfDelta = TRUE; 591 HalfDelta = TRUE;
592 AdjMean += Delta; 592 AdjMean += Delta;
593 } 593 }
594 } while (TRUE); // Loop checking, 11.05.2 594 } while (TRUE); // Loop checking, 11.05.2015, T. Koi
595 595
596 ShiftedGaussianValues_->G4InsertShiftedMea 596 ShiftedGaussianValues_->G4InsertShiftedMean(AdjMean, Mean_, StdDev_);
597 Mean_ = AdjMean; 597 Mean_ = AdjMean;
598 } 598 }
599 else if (Mean_ / 7 < StdDev_) { 599 else if (Mean_ / 7 < StdDev_) {
600 // If the requested type is double, then j 600 // If the requested type is double, then just re-define the standard
601 // deviation appropriately - chances are a 601 // deviation appropriately - chances are approximately 2.56E-12 that
602 // the value will be negative using this s 602 // the value will be negative using this sampling scheme
603 StdDev_ = Mean_ / 7; 603 StdDev_ = Mean_ / 7;
604 } 604 }
605 605
606 G4FFG_SAMPLING_FUNCTIONLEAVE__ 606 G4FFG_SAMPLING_FUNCTIONLEAVE__
607 } 607 }
608 608
609 G4FPYSamplingOps::~G4FPYSamplingOps() 609 G4FPYSamplingOps::~G4FPYSamplingOps()
610 { 610 {
611 G4FFG_FUNCTIONENTER__ 611 G4FFG_FUNCTIONENTER__
612 612
613 delete ShiftedGaussianValues_; 613 delete ShiftedGaussianValues_;
614 delete WattConstants_; 614 delete WattConstants_;
615 615
616 G4FFG_FUNCTIONLEAVE__ 616 G4FFG_FUNCTIONLEAVE__
617 } 617 }
618 618