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