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