ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FPYSamplingOps.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FPYSamplingOps.cc
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 }