Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FPYSamplingOps.cc
Go to the documentation of this file.
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 <iostream>
34 
35 #include <CLHEP/Random/Stat.h>
36 #include "Randomize.hh"
37 #include "globals.hh"
38 #include "G4Log.hh"
39 #include "G4Pow.hh"
40 
41 #include "G4FFGDebuggingMacros.hh"
42 #include "G4FFGDefaultValues.hh"
43 #include "G4FFGEnumerations.hh"
44 #include "G4FPYSamplingOps.hh"
45 #include "G4ShiftedGaussian.hh"
47 
50 {
51  // Set the default verbosity
53 
54  // Initialize the class
55  Initialize();
56 }
57 
60 {
61  // Set the default verbosity
63 
64  // Initialize the class
65  Initialize();
66 }
67 
69 Initialize( void )
70 {
72 
73  // Get the pointer to the random number generator
74  //RandomEngine_ = CLHEP::HepRandom::getTheEngine();
75  RandomEngine_ = G4Random::getTheEngine();
76 
77  // Initialize the data members
79  Mean_ = 0;
80  StdDev_ = 0;
82  GaussianOne_ = 0;
83  GaussianTwo_ = 0;
84  Tolerance_ = 0.000001;
87 
89 }
90 
93  G4double StdDev )
94 {
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
103 
104  Mean_ = Mean;
105  StdDev_ = StdDev;
106  }
107 
108  G4double Sample = SampleGaussian();
109 
111  return Sample;
112 }
113 
116  G4double StdDev,
118 {
120 
121  if(Range == G4FFGEnumerations::ALL)
122  {
123  // Call the overloaded function
124  G4double Sample = G4SampleGaussian(Mean, StdDev);
125 
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.");
143 
145  return 0;
146  }
147 
148  // Set the new values if the parameters have changed and then perform
149  // the shift
150  Mean_ = Mean;
151  StdDev_ = StdDev;
152 
154  }
155 
156  // Sample the Gaussian distribution
157  G4double Rand;
158  do
159  {
160  Rand = SampleGaussian();
161  } while(Rand < 0); // Loop checking, 11.05.2015, T. Koi
162 
164  return Rand;
165 }
166 
169  G4double StdDev )
170 {
172 
173  // Determine if the parameters have changed
174  G4bool ParametersChanged = (Mean_ != Mean || StdDev_ != StdDev);
175  if(ParametersChanged == TRUE)
176  {
177  // Set the new values if the parameters have changed
179 
180  Mean_ = Mean;
181  StdDev_ = StdDev;
182  }
183 
184  // Return the sample integer value
185  G4int Sample = (G4int)(std::floor(SampleGaussian()));
186 
188  return Sample;
189 }
190 
193  G4double StdDev,
195 {
197 
198  if(Range == G4FFGEnumerations::ALL)
199  {
200  // Call the overloaded function
201  G4int Sample = G4SampleIntegerGaussian(Mean, StdDev);
202 
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.");
216 
217  // return 0;
218  }
219 
220  if(Mean / StdDev < 2)
221  {
222  //Temp << "Non-ideal conditions:\tMean:" << Mean << "\tStdDev: "
223  // << StdDev;
224  //G4Exception("G4FPYGaussianOps::G4SampleIntegerGaussian()",
225  // Temp.str().c_str(),
226  // JustWarning,
227  // "Results not guaranteed: Consider tightening the standard deviation");
228  }
229 
230  // Determine if the parameters have changed
231  G4bool ParametersChanged = (Mean_ != Mean ||
232  StdDev_ != StdDev);
233  if(ParametersChanged == TRUE)
234  {
235  // Set the new values if the parameters have changed and then perform
236  // the shift
237  Mean_ = Mean;
238  StdDev_ = StdDev;
239 
241  }
242 
243  G4int RandInt;
244  // Sample the Gaussian distribution - only non-negative values are
245  // accepted
246  do
247  {
248  RandInt = (G4int)floor(SampleGaussian());
249  } while (RandInt < 0); // Loop checking, 11.05.2015, T. Koi
250 
252  return RandInt;
253  }
254 }
255 
258 {
260 
261  G4double Sample = RandomEngine_->flat();
262 
264  return Sample;
265 }
266 
269  G4double Upper )
270 {
272 
273  // Calculate the difference
274  G4double Difference = Upper - Lower;
275 
276  // Scale appropriately and return the value
277  G4double Sample = G4SampleUniform() * Difference + Lower;
278 
280  return Sample;
281 }
282 
284 G4SampleWatt( G4int WhatIsotope,
286  G4double WhatEnergy )
287 {
289 
290  // Determine if the parameters are different
291  //TK modified 131108
292  //if(WattConstants_->Product != WhatIsotope
293  if(WattConstants_->Product != WhatIsotope/10
294  || WattConstants_->Cause != WhatCause
295  || WattConstants_->Energy!= WhatEnergy )
296  {
297  // If the parameters are different or have not yet been defined then we
298  // need to re-evaluate the constants
299  //TK modified 131108
300  //WattConstants_->Product = WhatIsotope;
301  WattConstants_->Product = WhatIsotope/10;
302  WattConstants_->Cause = WhatCause;
303  WattConstants_->Energy = WhatEnergy;
304 
306  }
307 
310  G4int icounter=0;
311  G4int icounter_max=1024;
312  while(G4Pow::GetInstance()->powN(Y - WattConstants_->M*(X + 1), 2)
313  > WattConstants_->B * WattConstants_->L * X) // Loop checking, 11.05.2015, T. Koi
314  {
315  icounter++;
316  if ( icounter > icounter_max ) {
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());
322  }
323 
325  return WattConstants_->L * X;
326 }
327 
330 {
332 
334 
336 
338 }
339 
342 {
344 
346  if(ShiftedMean == 0)
347  {
349  return FALSE;
350  }
351 
352  Mean_ = ShiftedMean;
353 
355  return TRUE;
356 }
357 
360 {
362 
363  G4double A, K;
364  A = K = 0;
365  // Use the default values if IsotopeIndex is not reset
366  G4int IsotopeIndex = 0;
367 
369  {
370  // Determine if the isotope requested exists in the lookup tables
371  for(G4int i = 0; SpontaneousWattIsotopesIndex[i] != -1; i++)
372  {
375  {
376  IsotopeIndex = i;
377 
378  break;
379  }
380  }
381 
382  // Get A and B
383  A = SpontaneousWattConstants[IsotopeIndex][0];
384  WattConstants_->B = SpontaneousWattConstants[IsotopeIndex][1];
386  {
387  // Determine if the isotope requested exists in the lookup tables
388  for(G4int i = 0; NeutronInducedWattIsotopesIndex[i] != -1; i++)
389  {
391  {
392  IsotopeIndex = i;
393  break;
394  }
395  }
396 
397  // Determine the Watt fission spectrum constants based on the energy of
398  // the incident neutron
400  {
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  {
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 ";
466  {
467  Temp += "proton induced fission.";
469  {
470  Temp += "gamma induced fission.";
471  } else
472  {
473  Temp += "!Warning! unknown cause.";
474  }
475  G4Exception("G4FPYSamplingOps::G4SampleWatt()",
476  Temp,
478  "Fission events will not be sampled in this run.");
479  }
480 
481  // Calculate the constants
482  K = 1 + (WattConstants_->B / (8.0 * A));
483  WattConstants_->L = (K + G4Pow::GetInstance()->powA((K * K - 1), 0.5)) / A;
484  WattConstants_->M = A * WattConstants_->L - 1;
485 
487 }
488 
491 {
493 
495  {
497 
499  return GaussianTwo_;
500  }
501 
502  // Define the variables to be used
503  G4double Radius;
504  G4double MappingFactor;
505 
506  // Sample from the unit circle (21.4% rejection probability)
507  do
508  {
509  GaussianOne_ = 2.0 * G4SampleUniform() - 1.0;
510  GaussianTwo_ = 2.0 * G4SampleUniform() - 1.0;
512  } while (Radius > 1.0); // Loop checking, 11.05.2015, T. Koi
513 
514  // Translate the values to Gaussian space
515  MappingFactor = std::sqrt(-2.0*G4Log(Radius)/Radius) * StdDev_;
516  GaussianOne_ = Mean_ + GaussianOne_*MappingFactor;
517  GaussianTwo_ = Mean_ + GaussianTwo_*MappingFactor;
518 
519  // Set the flag that a value is now stored in memory
521 
523  return GaussianOne_;
524 }
525 
528 {
530 
531  // Set the flag that any second value stored is no longer valid
533 
534  // Check if the requested parameters have already been calculated
535  if(CheckAndSetParameters() == TRUE)
536  {
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  {
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 
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 
676 }
677 
679 {
681 
682  delete ShiftedGaussianValues_;
683  delete WattConstants_;
684 
686 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double G4SampleGaussian(G4double Mean, G4double StdDev)
static const G4double NeutronInducedWattConstants[][3][2]
G4double SampleGaussian(void)
double Y(double density)
void EvaluateWattConstants(void)
G4double G4SampleWatt(G4int WhatIsotope, G4FFGEnumerations::FissionCause WhatCause, G4double WhatEnergy)
#define G4FFG_SAMPLING_FUNCTIONENTER__
G4double G4SampleUniform(void)
static const G4double SpontaneousWattConstants[][2]
G4double G4FindShiftedMean(G4double RequestedMean, G4double RequestedStdDev)
G4FFGEnumerations::FissionCause Cause
virtual double flat()=0
G4ShiftedGaussian * ShiftedGaussianValues_
int G4int
Definition: G4Types.hh:78
WattSpectrumConstants * WattConstants_
void G4SetVerbosity(G4int WhatVerbosity)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static const G4int SpontaneousWattIsotopesIndex[]
bool G4bool
Definition: G4Types.hh:79
void G4InsertShiftedMean(G4double ShiftedMean, G4double RequestedMean, G4double RequestedStdDev)
void ShiftParameters(G4FFGEnumerations::GaussianReturnType Type)
#define FALSE
Definition: globals.hh:52
static constexpr double MeV
#define TRUE
Definition: globals.hh:55
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4int NeutronInducedWattIsotopesIndex[]
G4double G4Log(G4double x)
Definition: G4Log.hh:230
CLHEP::HepRandomEngine * RandomEngine_
#define G4FFG_FUNCTIONLEAVE__
#define G4endl
Definition: G4ios.hh:61
static double erf(double x)
#define G4FFG_SAMPLING_FUNCTIONLEAVE__
double G4double
Definition: G4Types.hh:76
void G4SetVerbosity(G4int WhatVerbosity)
static const G4double IncidentEnergyBins[]
static const G4double ThermalNeutronEnergy
G4bool CheckAndSetParameters(void)
#define G4FFG_FUNCTIONENTER__
static const G4FFGEnumerations::Verbosity Verbosity
G4bool NextGaussianIsStoredInMemory_
G4int G4SampleIntegerGaussian(G4double Mean, G4double StdDev)