| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 Verbo
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 gene
72 // RandomEngine_ = CLHEP::HepRandom::getTheE
73 RandomEngine_ = G4Random::getTheEngine();
74
75 // Initialize the data members
76 ShiftedGaussianValues_ = new G4ShiftedGaussi
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(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;
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(G4
110 G4
111 {
112 G4FFG_SAMPLING_FUNCTIONENTER__
113
114 if (Range == G4FFGEnumerations::ALL) {
115 // Call the overloaded function
116 G4double Sample = G4SampleGaussian(Mean, S
117
118 G4FFG_SAMPLING_FUNCTIONLEAVE__
119 return Sample;
120 }
121
122 // Determine if the parameters have changed
123 G4bool ParametersChanged = (Mean_ != Mean ||
124 if (static_cast<int>(ParametersChanged) == T
125 if (Mean <= 0) {
126 std::ostringstream Temp;
127 Temp << "Mean value of " << Mean << " ou
128 G4Exception("G4FPYGaussianOps::G4SampleI
129 "A value of '0' will be used
130
131 G4FFG_SAMPLING_FUNCTIONLEAVE__
132 return 0;
133 }
134
135 // Set the new values if the parameters ha
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
148
149 G4FFG_SAMPLING_FUNCTIONLEAVE__
150 return Rand;
151 }
152
153 G4int G4FPYSamplingOps::G4SampleIntegerGaussia
154 {
155 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
167 // Return the sample integer value
168 auto Sample = (G4int)(std::floor(SampleGauss
169
170 G4FFG_SAMPLING_FUNCTIONLEAVE__
171 return Sample;
172 }
173
174 G4int G4FPYSamplingOps::G4SampleIntegerGaussia
175
176 {
177 G4FFG_SAMPLING_FUNCTIONENTER__
178
179 if (Range == G4FFGEnumerations::ALL) {
180 // Call the overloaded function
181 G4int Sample = G4SampleIntegerGaussian(Mea
182
183 G4FFG_SAMPLING_FUNCTIONLEAVE__
184 return Sample;
185 }
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
195 // return 0;
196 }
197
198 if (Mean / StdDev < 2) {
199 // Temp << "Non-ideal conditions:\tMean:"
200 // << StdDev;
201 // G4Exception("G4FPYGaussianOps::G4Sample
202 // Temp.str().c_str(),
203 // JustWarning,
204 // "Results not guaranteed: Co
205 }
206
207 // Determine if the parameters have changed
208 G4bool ParametersChanged = (Mean_ != Mean ||
209 if (static_cast<int>(ParametersChanged) == T
210 // Set the new values if the parameters ha
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 n
220 // accepted
221 do {
222 RandInt = (G4int)floor(SampleGaussian());
223 } while (RandInt < 0); // Loop checking, 11
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(G4d
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() * Differ
248
249 G4FFG_SAMPLING_FUNCTIONLEAVE__
250 return Sample;
251 }
252
253 G4double G4FPYSamplingOps::G4SampleWatt(G4int
254 G4FFGE
255 G4doub
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 /
263 || WattConstants_->Energy != WhatEnergy)
264 {
265 // If the parameters are different or have
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 - WattCo
281 > WattConstants_->B * WattConstants_-
282 {
283 icounter++;
284 if (icounter > icounter_max) {
285 G4cout << "Loop-counter exceeded the thr
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 Ve
298 {
299 G4FFG_SAMPLING_FUNCTIONENTER__
300
301 Verbosity_ = Verbosity;
302
303 ShiftedGaussianValues_->G4SetVerbosity(Verbo
304
305 G4FFG_SAMPLING_FUNCTIONLEAVE__
306 }
307
308 G4bool G4FPYSamplingOps::CheckAndSetParameters
309 {
310 G4FFG_SAMPLING_FUNCTIONENTER__
311
312 G4double ShiftedMean = ShiftedGaussianValues
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
331 G4int IsotopeIndex = 0;
332
333 if (WattConstants_->Cause == G4FFGEnumeratio
334 // Determine if the isotope requested exis
335 for (G4int i = 0; SpontaneousWattIsotopesI
336 if (SpontaneousWattIsotopesIndex[i] == W
337 IsotopeIndex = i;
338
339 break;
340 }
341 }
342
343 // Get A and B
344 A = SpontaneousWattConstants[IsotopeIndex]
345 WattConstants_->B = SpontaneousWattConstan
346 }
347 else if (WattConstants_->Cause == G4FFGEnume
348 // Determine if the isotope requested exis
349 for (G4int i = 0; NeutronInducedWattIsotop
350 if (NeutronInducedWattIsotopesIndex[i] =
351 IsotopeIndex = i;
352 break;
353 }
354 }
355
356 // Determine the Watt fission spectrum con
357 // the incident neutron
358 if (WattConstants_->Energy == G4FFGDefault
359 A = NeutronInducedWattConstants[IsotopeI
360 WattConstants_->B = NeutronInducedWattCo
361 }
362 else if (WattConstants_->Energy > 14.0 * C
363 G4Exception("G4FPYSamplingOps::G4SampleW
364 "Incident neutron energy abo
365 "Using Watt fission constant
366
367 A = NeutronInducedWattConstants[IsotopeI
368 WattConstants_->B = NeutronInducedWattCo
369 }
370 else {
371 G4int EnergyIndex = 0;
372 G4double EnergyDifference = 0;
373 G4double RangeDifference, ConstantDiffer
374
375 for (G4int i = 1; IncidentEnergyBins[i]
376 if (WattConstants_->Energy <= Incident
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
393 RangeDifference = IncidentEnergyBins[Ene
394
395 // Interpolate the value for 'a'
396 ConstantDifference = NeutronInducedWattC
397 - NeutronInducedWat
398 A = (EnergyDifference / RangeDifference)
399 + NeutronInducedWattConstants[Isotop
400
401 // Interpolate the value for 'b'
402 ConstantDifference = NeutronInducedWattC
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 }
421 G4Exception("G4FPYSamplingOps::G4SampleWat
422 "Fission events will not be sa
423 }
424
425 // Calculate the constants
426 K = 1 + (WattConstants_->B / (8.0 * A));
427 WattConstants_->L = (K + G4Pow::GetInstance(
428 WattConstants_->M = A * WattConstants_->L -
429
430 G4FFG_SAMPLING_FUNCTIONLEAVE__
431 }
432
433 G4double G4FPYSamplingOps::SampleGaussian()
434 {
435 G4FFG_SAMPLING_FUNCTIONENTER__
436
437 if (static_cast<int>(NextGaussianIsStoredInM
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% 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
546 // Container = (erf(HighErf) + erf(LowErf))/2.
547 #if defined WIN32 - VC
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
562 #if defined WIN32 - VC
563 // TK modified to avoid problem caused
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 }
608
609 G4FPYSamplingOps::~G4FPYSamplingOps()
610 {
611 G4FFG_FUNCTIONENTER__
612
613 delete ShiftedGaussianValues_;
614 delete WattConstants_;
615
616 G4FFG_FUNCTIONLEAVE__
617 }
618