| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
27 * File: G4FPYSamplingOps.cc
28 * Author: B. Wendt (wendbryc@isu.edu)
29 *
30 * Created on June 30, 2011, 11:10 AM
31 */
32
33 #include "G4FPYSamplingOps.hh"
34
35 #include "G4FFGDebuggingMacros.hh"
36 #include "G4FFGDefaultValues.hh"
37 #include "G4FFGEnumerations.hh"
38 #include "G4Log.hh"
39 #include "G4Pow.hh"
40 #include "G4ShiftedGaussian.hh"
41 #include "G4WattFissionSpectrumValues.hh"
42 #include "Randomize.hh"
43 #include "globals.hh"
44
45 #include <CLHEP/Random/Stat.h>
46
47 #include <iostream>
48
49 G4FPYSamplingOps::G4FPYSamplingOps()
50 {
51 // Set the default verbosity
52 Verbosity_ = G4FFGDefaultValues::Verbosity;
53
54 // Initialize the class
55 Initialize();
56 }
57
58 G4FPYSamplingOps::G4FPYSamplingOps(G4int Verbosity)
59 {
60 // Set the default verbosity
61 Verbosity_ = Verbosity;
62
63 // Initialize the class
64 Initialize();
65 }
66
67 void G4FPYSamplingOps::Initialize()
68 {
69 G4FFG_FUNCTIONENTER__
70
71 // Get the pointer to the random number generator
72 // RandomEngine_ = CLHEP::HepRandom::getTheEngine();
73 RandomEngine_ = G4Random::getTheEngine();
74
75 // Initialize the data members
76 ShiftedGaussianValues_ = new G4ShiftedGaussian;
77 Mean_ = 0;
78 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(G4double Mean, G4double StdDev)
90 {
91 G4FFG_SAMPLING_FUNCTIONENTER__
92
93 // Determine if the parameters have changed
94 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
95 if (static_cast<int>(ParametersChanged) == TRUE) {
96 // Set the new values if the parameters have changed
97 NextGaussianIsStoredInMemory_ = FALSE;
98
99 Mean_ = Mean;
100 StdDev_ = StdDev;
101 }
102
103 G4double Sample = SampleGaussian();
104
105 G4FFG_SAMPLING_FUNCTIONLEAVE__
106 return Sample;
107 }
108
109 G4double G4FPYSamplingOps::G4SampleGaussian(G4double Mean, G4double StdDev,
110 G4FFGEnumerations::GaussianRange Range)
111 {
112 G4FFG_SAMPLING_FUNCTIONENTER__
113
114 if (Range == G4FFGEnumerations::ALL) {
115 // Call the overloaded function
116 G4double Sample = G4SampleGaussian(Mean, StdDev);
117
118 G4FFG_SAMPLING_FUNCTIONLEAVE__
119 return Sample;
120 }
121
122 // Determine if the parameters have changed
123 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
124 if (static_cast<int>(ParametersChanged) == TRUE) {
125 if (Mean <= 0) {
126 std::ostringstream Temp;
127 Temp << "Mean value of " << Mean << " out of range";
128 G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()", Temp.str().c_str(), JustWarning,
129 "A value of '0' will be used instead.");
130
131 G4FFG_SAMPLING_FUNCTIONLEAVE__
132 return 0;
133 }
134
135 // Set the new values if the parameters have changed and then perform
136 // the shift
137 Mean_ = Mean;
138 StdDev_ = StdDev;
139
140 ShiftParameters(G4FFGEnumerations::DOUBLE);
141 }
142
143 // Sample the Gaussian distribution
144 G4double Rand;
145 do {
146 Rand = SampleGaussian();
147 } while (Rand < 0); // Loop checking, 11.05.2015, T. Koi
148
149 G4FFG_SAMPLING_FUNCTIONLEAVE__
150 return Rand;
151 }
152
153 G4int G4FPYSamplingOps::G4SampleIntegerGaussian(G4double Mean, G4double StdDev)
154 {
155 G4FFG_SAMPLING_FUNCTIONENTER__
156
157 // Determine if the parameters have changed
158 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
159 if (static_cast<int>(ParametersChanged) == TRUE) {
160 // Set the new values if the parameters have changed
161 NextGaussianIsStoredInMemory_ = FALSE;
162
163 Mean_ = Mean;
164 StdDev_ = StdDev;
165 }
166
167 // Return the sample integer value
168 auto Sample = (G4int)(std::floor(SampleGaussian()));
169
170 G4FFG_SAMPLING_FUNCTIONLEAVE__
171 return Sample;
172 }
173
174 G4int G4FPYSamplingOps::G4SampleIntegerGaussian(G4double Mean, G4double StdDev,
175 G4FFGEnumerations::GaussianRange Range)
176 {
177 G4FFG_SAMPLING_FUNCTIONENTER__
178
179 if (Range == G4FFGEnumerations::ALL) {
180 // Call the overloaded function
181 G4int Sample = G4SampleIntegerGaussian(Mean, StdDev);
182
183 G4FFG_SAMPLING_FUNCTIONLEAVE__
184 return Sample;
185 }
186 // ParameterShift() locks up if the mean is less than 1.
187 std::ostringstream Temp;
188 if (Mean < 1) {
189 // Temp << "Mean value of " << Mean << " out of range";
190 // G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()",
191 // Temp.str().c_str(),
192 // JustWarning,
193 // "A value of '0' will be used instead.");
194
195 // return 0;
196 }
197
198 if (Mean / StdDev < 2) {
199 // Temp << "Non-ideal conditions:\tMean:" << Mean << "\tStdDev: "
200 // << StdDev;
201 // G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()",
202 // Temp.str().c_str(),
203 // JustWarning,
204 // "Results not guaranteed: Consider tightening the standard deviation");
205 }
206
207 // Determine if the parameters have changed
208 G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
209 if (static_cast<int>(ParametersChanged) == TRUE) {
210 // Set the new values if the parameters have changed and then perform
211 // the shift
212 Mean_ = Mean;
213 StdDev_ = StdDev;
214
215 ShiftParameters(G4FFGEnumerations::INT);
216 }
217
218 G4int RandInt;
219 // Sample the Gaussian distribution - only non-negative values are
220 // accepted
221 do {
222 RandInt = (G4int)floor(SampleGaussian());
223 } while (RandInt < 0); // Loop checking, 11.05.2015, T. Koi
224
225 G4FFG_SAMPLING_FUNCTIONLEAVE__
226 return RandInt;
227 }
228
229 G4double G4FPYSamplingOps::G4SampleUniform()
230 {
231 G4FFG_SAMPLING_FUNCTIONENTER__
232
233 G4double Sample = RandomEngine_->flat();
234
235 G4FFG_SAMPLING_FUNCTIONLEAVE__
236 return Sample;
237 }
238
239 G4double G4FPYSamplingOps::G4SampleUniform(G4double Lower, G4double Upper)
240 {
241 G4FFG_SAMPLING_FUNCTIONENTER__
242
243 // Calculate the difference
244 G4double Difference = Upper - Lower;
245
246 // Scale appropriately and return the value
247 G4double Sample = G4SampleUniform() * Difference + Lower;
248
249 G4FFG_SAMPLING_FUNCTIONLEAVE__
250 return Sample;
251 }
252
253 G4double G4FPYSamplingOps::G4SampleWatt(G4int WhatIsotope,
254 G4FFGEnumerations::FissionCause WhatCause,
255 G4double WhatEnergy)
256 {
257 G4FFG_SAMPLING_FUNCTIONENTER__
258
259 // Determine if the parameters are different
260 // TK modified 131108
261 // if(WattConstants_->Product != WhatIsotope
262 if (WattConstants_->Product != WhatIsotope / 10 || WattConstants_->Cause != WhatCause
263 || WattConstants_->Energy != WhatEnergy)
264 {
265 // If the parameters are different or have not yet been defined then we
266 // need to re-evaluate the constants
267 // TK modified 131108
268 // WattConstants_->Product = WhatIsotope;
269 WattConstants_->Product = WhatIsotope / 10;
270 WattConstants_->Cause = WhatCause;
271 WattConstants_->Energy = WhatEnergy;
272
273 EvaluateWattConstants();
274 }
275
276 G4double X = -G4Log(G4SampleUniform());
277 G4double Y = -G4Log(G4SampleUniform());
278 G4int icounter = 0;
279 G4int icounter_max = 1024;
280 while (G4Pow::GetInstance()->powN(Y - WattConstants_->M * (X + 1), 2)
281 > WattConstants_->B * WattConstants_->L * X) // Loop checking, 11.05.2015, T. Koi
282 {
283 icounter++;
284 if (icounter > icounter_max) {
285 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
286 << __FILE__ << "." << G4endl;
287 break;
288 }
289 X = -G4Log(G4SampleUniform());
290 Y = -G4Log(G4SampleUniform());
291 }
292
293 G4FFG_SAMPLING_FUNCTIONLEAVE__
294 return WattConstants_->L * X;
295 }
296
297 void G4FPYSamplingOps::G4SetVerbosity(G4int Verbosity)
298 {
299 G4FFG_SAMPLING_FUNCTIONENTER__
300
301 Verbosity_ = Verbosity;
302
303 ShiftedGaussianValues_->G4SetVerbosity(Verbosity_);
304
305 G4FFG_SAMPLING_FUNCTIONLEAVE__
306 }
307
308 G4bool G4FPYSamplingOps::CheckAndSetParameters()
309 {
310 G4FFG_SAMPLING_FUNCTIONENTER__
311
312 G4double ShiftedMean = ShiftedGaussianValues_->G4FindShiftedMean(Mean_, StdDev_);
313 if (ShiftedMean == 0) {
314 G4FFG_SAMPLING_FUNCTIONLEAVE__
315 return FALSE;
316 }
317
318 Mean_ = ShiftedMean;
319
320 G4FFG_SAMPLING_FUNCTIONLEAVE__
321 return TRUE;
322 }
323
324 void G4FPYSamplingOps::EvaluateWattConstants()
325 {
326 G4FFG_SAMPLING_FUNCTIONENTER__
327
328 G4double A, K;
329 A = K = 0;
330 // Use the default values if IsotopeIndex is not reset
331 G4int IsotopeIndex = 0;
332
333 if (WattConstants_->Cause == G4FFGEnumerations::SPONTANEOUS) {
334 // Determine if the isotope requested exists in the lookup tables
335 for (G4int i = 0; SpontaneousWattIsotopesIndex[i] != -1; i++) {
336 if (SpontaneousWattIsotopesIndex[i] == WattConstants_->Product) {
337 IsotopeIndex = i;
338
339 break;
340 }
341 }
342
343 // Get A and B
344 A = SpontaneousWattConstants[IsotopeIndex][0];
345 WattConstants_->B = SpontaneousWattConstants[IsotopeIndex][1];
346 }
347 else if (WattConstants_->Cause == G4FFGEnumerations::NEUTRON_INDUCED) {
348 // Determine if the isotope requested exists in the lookup tables
349 for (G4int i = 0; NeutronInducedWattIsotopesIndex[i] != -1; i++) {
350 if (NeutronInducedWattIsotopesIndex[i] == WattConstants_->Product) {
351 IsotopeIndex = i;
352 break;
353 }
354 }
355
356 // Determine the Watt fission spectrum constants based on the energy of
357 // the incident neutron
358 if (WattConstants_->Energy == G4FFGDefaultValues::ThermalNeutronEnergy) {
359 A = NeutronInducedWattConstants[IsotopeIndex][0][0];
360 WattConstants_->B = NeutronInducedWattConstants[IsotopeIndex][0][1];
361 }
362 else if (WattConstants_->Energy > 14.0 * CLHEP::MeV) {
363 G4Exception("G4FPYSamplingOps::G4SampleWatt()",
364 "Incident neutron energy above 14 MeV requested.", JustWarning,
365 "Using Watt fission constants for 14 Mev.");
366
367 A = NeutronInducedWattConstants[IsotopeIndex][2][0];
368 WattConstants_->B = NeutronInducedWattConstants[IsotopeIndex][2][1];
369 }
370 else {
371 G4int EnergyIndex = 0;
372 G4double EnergyDifference = 0;
373 G4double RangeDifference, ConstantDifference;
374
375 for (G4int i = 1; IncidentEnergyBins[i] != -1; i++) {
376 if (WattConstants_->Energy <= IncidentEnergyBins[i]) {
377 EnergyIndex = i;
378 EnergyDifference = IncidentEnergyBins[EnergyIndex] - WattConstants_->Energy;
379 if (EnergyDifference != 0) {
380 std::ostringstream Temp;
381 Temp << "Incident neutron energy of ";
382 Temp << WattConstants_->Energy << " MeV is not ";
383 Temp << "explicitly listed in the data tables";
384 // G4Exception("G4FPYSamplingOps::G4SampleWatt()",
385 // Temp.str().c_str(),
386 // JustWarning,
387 // "Using linear interpolation.");
388 }
389 break;
390 }
391 }
392
393 RangeDifference = IncidentEnergyBins[EnergyIndex] - IncidentEnergyBins[EnergyIndex - 1];
394
395 // Interpolate the value for 'a'
396 ConstantDifference = NeutronInducedWattConstants[IsotopeIndex][EnergyIndex][0]
397 - NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][0];
398 A = (EnergyDifference / RangeDifference) * ConstantDifference
399 + NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][0];
400
401 // Interpolate the value for 'b'
402 ConstantDifference = NeutronInducedWattConstants[IsotopeIndex][EnergyIndex][1]
403 - NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][1];
404 WattConstants_->B = (EnergyDifference / RangeDifference) * ConstantDifference
405 + NeutronInducedWattConstants[IsotopeIndex][EnergyIndex - 1][1];
406 }
407 }
408 else {
409 // Produce an error since an unsupported fission type was requested and
410 // abort the current run
411 G4String Temp = "Watt fission spectra data not available for ";
412 if (WattConstants_->Cause == G4FFGEnumerations::PROTON_INDUCED) {
413 Temp += "proton induced fission.";
414 }
415 else if (WattConstants_->Cause == G4FFGEnumerations::GAMMA_INDUCED) {
416 Temp += "gamma induced fission.";
417 }
418 else {
419 Temp += "!Warning! unknown cause.";
420 }
421 G4Exception("G4FPYSamplingOps::G4SampleWatt()", Temp, RunMustBeAborted,
422 "Fission events will not be sampled in this run.");
423 }
424
425 // Calculate the constants
426 K = 1 + (WattConstants_->B / (8.0 * A));
427 WattConstants_->L = (K + G4Pow::GetInstance()->powA((K * K - 1), 0.5)) / A;
428 WattConstants_->M = A * WattConstants_->L - 1;
429
430 G4FFG_SAMPLING_FUNCTIONLEAVE__
431 }
432
433 G4double G4FPYSamplingOps::SampleGaussian()
434 {
435 G4FFG_SAMPLING_FUNCTIONENTER__
436
437 if (static_cast<int>(NextGaussianIsStoredInMemory_) == TRUE) {
438 NextGaussianIsStoredInMemory_ = FALSE;
439
440 G4FFG_SAMPLING_FUNCTIONLEAVE__
441 return GaussianTwo_;
442 }
443
444 // Define the variables to be used
445 G4double Radius;
446 G4double MappingFactor;
447
448 // Sample from the unit circle (21.4% rejection probability)
449 do {
450 GaussianOne_ = 2.0 * G4SampleUniform() - 1.0;
451 GaussianTwo_ = 2.0 * G4SampleUniform() - 1.0;
452 Radius = GaussianOne_ * GaussianOne_ + GaussianTwo_ * GaussianTwo_;
453 } while (Radius > 1.0); // Loop checking, 11.05.2015, T. Koi
454
455 // Translate the values to Gaussian space
456 MappingFactor = std::sqrt(-2.0 * G4Log(Radius) / Radius) * StdDev_;
457 GaussianOne_ = Mean_ + GaussianOne_ * MappingFactor;
458 GaussianTwo_ = Mean_ + GaussianTwo_ * MappingFactor;
459
460 // Set the flag that a value is now stored in memory
461 NextGaussianIsStoredInMemory_ = TRUE;
462
463 G4FFG_SAMPLING_FUNCTIONLEAVE__
464 return GaussianOne_;
465 }
466
467 void G4FPYSamplingOps::ShiftParameters(G4FFGEnumerations::GaussianReturnType Type)
468 {
469 G4FFG_SAMPLING_FUNCTIONENTER__
470
471 // Set the flag that any second value stored is no longer valid
472 NextGaussianIsStoredInMemory_ = FALSE;
473
474 // Check if the requested parameters have already been calculated
475 if (static_cast<int>(CheckAndSetParameters()) == TRUE) {
476 G4FFG_SAMPLING_FUNCTIONLEAVE__
477 return;
478 }
479
480 // If the requested type is INT, then perform an iterative refinement of the
481 // mean that is going to be sampled
482 if (Type == G4FFGEnumerations::INT) {
483 // Return if the mean is greater than 7 standard deviations away from 0
484 // since there is essentially 0 probability that a sampled number will
485 // be less than 0
486 if (Mean_ > 7 * StdDev_) {
487 G4FFG_SAMPLING_FUNCTIONLEAVE__
488 return;
489 }
490 // Variables that contain the area and weighted area information for
491 // calculating the statistical mean of the Gaussian distribution when
492 // converted to a step function
493 G4double ErfContainer, AdjustedErfContainer, Container;
494
495 // Variables that store lower and upper bounds information
496 G4double LowErf, HighErf;
497
498 // Values for determining the convergence of the solution
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_ * std::sqrt(2.0);
506
507 // Determine the upper limit of the estimates
508 const auto UpperLimit = (G4int)std::ceil(Mean_ + 9 * StdDev_);
509
510 // Calculate the integral of the Gaussian distribution
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 threshold value at " << __LINE__ << "th line of "
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 bounds
528 LowErf = ((AdjMean - i) / Normalization);
529 HighErf = ((AdjMean - (i + 1.0)) / Normalization);
530
531 // Correct the bounds for how they lie on the x-axis with
532 // respect to the mean
533 if (LowErf <= 0) {
534 LowErf *= -1;
535 HighErf *= -1;
536 // Container = (erf(HighErf) - erf(LowErf))/2.0;
537 #if defined WIN32 - VC
538 Container = (CLHEP::HepStat::erf(HighErf) - CLHEP::HepStat::erf(LowErf)) / 2.0;
539 #else
540 Container = (erf(HighErf) - erf(LowErf)) / 2.0;
541 #endif
542 }
543 else if (HighErf < 0) {
544 HighErf *= -1;
545
546 // Container = (erf(HighErf) + erf(LowErf))/2.0;
547 #if defined WIN32 - VC
548 Container = (CLHEP::HepStat::erf(HighErf) + CLHEP::HepStat::erf(LowErf)) / 2.0;
549 #else
550 Container = (erf(HighErf) + erf(LowErf)) / 2.0;
551 #endif
552 }
553 else {
554 // Container = (erf(LowErf) - erf(HighErf))/2.0;
555 #if defined WIN32 - VC
556 Container = (CLHEP::HepStat::erf(LowErf) - CLHEP::HepStat::erf(HighErf)) / 2.0;
557 #else
558 Container = (erf(LowErf) - erf(HighErf)) / 2.0;
559 #endif
560 }
561
562 #if defined WIN32 - VC
563 // TK modified to avoid problem caused by QNAN
564 if (Container != Container) Container = 0;
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 / ErfContainer;
574
575 // Is it close enough to what we want?
576 ToleranceCheck = (std::fabs(Mean_ - Container) < Tolerance_);
577 if (static_cast<int>(ToleranceCheck) == TRUE) {
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.2015, T. Koi
595
596 ShiftedGaussianValues_->G4InsertShiftedMean(AdjMean, Mean_, StdDev_);
597 Mean_ = AdjMean;
598 }
599 else if (Mean_ / 7 < StdDev_) {
600 // If the requested type is double, then just re-define the standard
601 // deviation appropriately - chances are approximately 2.56E-12 that
602 // the value will be negative using this sampling scheme
603 StdDev_ = Mean_ / 7;
604 }
605
606 G4FFG_SAMPLING_FUNCTIONLEAVE__
607 }
608
609 G4FPYSamplingOps::~G4FPYSamplingOps()
610 {
611 G4FFG_FUNCTIONENTER__
612
613 delete ShiftedGaussianValues_;
614 delete WattConstants_;
615
616 G4FFG_FUNCTIONLEAVE__
617 }
618