G4FissionProductYieldDist.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
/*
 * File:   G4FissionProductYieldDist.cc
 * Author: B. Wendt ([email protected])
 *
 * Created on May 11, 2011, 12:04 PM
 */
 
/* * * * * * * * * * * * * * * *   References   * * * * * * * * * * * * * * * *
 *                                                                            *
 *  1.  "Systematics of fission fragment total kinetic energy release",       *
 *      V. E. Viola, K. Kwiatkowski, and M. Walker, Physical Review C, 31.4,  *
 *      April 1985                                                            *
 *  2.  "Reactor Handbook", United States Atomic Energy Commission,           *
 *      III.A:Physics, 1962                                                   *
 *  3.  "Properties of the Alpha Particles Emitted in the Spontaneous Fission *
 *      of Cf252", Z. Fraenkel and S. G. Thompson, Physical Review Letters,   *
 *      13.14, October 1964                                                   *
 *                                                                            *
 * * * * * * * * * * * * * * * *   References   * * * * * * * * * * * * * * * */
 
// #include <ios>
// #include <iostream>
 
#include "G4Alpha.hh"
#include "G4Gamma.hh"
#include "G4Ions.hh"
#include "G4Neutron.hh"
// #include "G4NeutronHPNames.hh"
#include "G4ArrayOps.hh"
#include "G4DynamicParticle.hh"
#include "G4DynamicParticleVector.hh"
#include "G4ENDFTapeRead.hh"
#include "G4ENDFYieldDataContainer.hh"
#include "G4Exp.hh"
#include "G4FFGDebuggingMacros.hh"
#include "G4FFGDefaultValues.hh"
#include "G4FFGEnumerations.hh"
#include "G4FPYNubarValues.hh"
#include "G4FPYSamplingOps.hh"
#include "G4FPYTreeStructures.hh"
#include "G4FissionProductYieldDist.hh"
#include "G4HadFinalState.hh"
#include "G4NucleiProperties.hh"
#include "G4ParticleTable.hh"
#include "G4Pow.hh"
#include "G4ReactionProduct.hh"
#include "G4SystemOfUnits.hh"
#include "G4ThreeVector.hh"
#include "G4UImanager.hh"
#include "Randomize.hh"
#include "globals.hh"
 
using CLHEP::pi;
 
#ifdef G4MULTITHREADED
#  include "G4AutoLock.hh"
G4Mutex G4FissionProductYieldDist::fissprodMutex = G4MUTEX_INITIALIZER;
#endif
 
G4FissionProductYieldDist::G4FissionProductYieldDist(G4int WhichIsotope,
                                                     G4FFGEnumerations::MetaState WhichMetaState,
                                                     G4FFGEnumerations::FissionCause WhichCause,
                                                     G4FFGEnumerations::YieldType WhichYieldType,
                                                     std::istringstream& dataStream)
  : Isotope_(WhichIsotope),
    MetaState_(WhichMetaState),
    Cause_(WhichCause),
    YieldType_(WhichYieldType),
    Verbosity_(G4FFGDefaultValues::Verbosity)
{
  G4FFG_FUNCTIONENTER__
 
  try {
    // Initialize the class
    Initialize(dataStream);
  }
  catch (std::exception& e) {
    G4FFG_FUNCTIONLEAVE__
    throw e;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
G4FissionProductYieldDist::G4FissionProductYieldDist(G4int WhichIsotope,
                                                     G4FFGEnumerations::MetaState WhichMetaState,
                                                     G4FFGEnumerations::FissionCause WhichCause,
                                                     G4FFGEnumerations::YieldType WhichYieldType,
                                                     G4int Verbosity,
                                                     std::istringstream& dataStream)
  : Isotope_(WhichIsotope),
    MetaState_(WhichMetaState),
    Cause_(WhichCause),
    YieldType_(WhichYieldType),
    Verbosity_(Verbosity)
{
  G4FFG_FUNCTIONENTER__
 
  try {
    // Initialize the class
    Initialize(dataStream);
  }
  catch (std::exception& e) {
    G4FFG_FUNCTIONLEAVE__
    throw e;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::Initialize(std::istringstream& dataStream)
{
  G4FFG_FUNCTIONENTER__
 
  IncidentEnergy_ = 0.0;
  TernaryProbability_ = 0;
  AlphaProduction_ = 0;
  SetNubar();
 
  // Set miscellaneous variables
  AlphaDefinition_ = static_cast<G4Ions*>(G4Alpha::Definition());
  NeutronDefinition_ = static_cast<G4Ions*>(G4Neutron::Definition());
  GammaDefinition_ = G4Gamma::Definition();
  SmallestZ_ = SmallestA_ = LargestZ_ = LargestA_ = nullptr;
 
  // Construct G4NeutronHPNames: provides access to the element names
  ElementNames_ = new G4ParticleHPNames;
  // Get the pointer to G4ParticleTable: stores all G4Ions
  IonTable_ = G4IonTable::GetIonTable();
  // Construct the pointer to the random engine
  // TODO Make G4FPSamplingOps a singleton so that only one instance is used across all classes
  RandomEngine_ = new G4FPYSamplingOps;
 
  try {
    // Read in and sort the probability data
    ENDFData_ = new G4ENDFTapeRead(dataStream, YieldType_, Cause_, Verbosity_);
    //        ENDFData_ = new G4ENDFTapeRead(MakeDirectoryName(),
    //                                       MakeFileName(Isotope_, MetaState_),
    //                                       YieldType_,
    //                                       Cause_);
    YieldEnergyGroups_ = ENDFData_->G4GetNumberOfEnergyGroups();
    DataTotal_ = new G4double[YieldEnergyGroups_];
    MaintainNormalizedData_ = new G4double[YieldEnergyGroups_];
    YieldEnergies_ = new G4double[YieldEnergyGroups_];
    G4ArrayOps::Copy(YieldEnergyGroups_, YieldEnergies_, ENDFData_->G4GetEnergyGroupValues());
    MakeTrees();
    ReadProbabilities();
  }
  catch (std::exception& e) {
    delete ElementNames_;
    delete RandomEngine_;
 
    G4FFG_FUNCTIONLEAVE__
    throw e;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
G4DynamicParticleVector* G4FissionProductYieldDist::G4GetFission()
{
  G4FFG_FUNCTIONENTER__
 
#ifdef G4MULTITHREADED
  G4AutoLock lk(&G4FissionProductYieldDist::fissprodMutex);
#endif
 
  // Check to see if the user has set the alpha production to a somewhat
  // reasonable level
  CheckAlphaSanity();
 
  // Generate the new G4DynamicParticle pointers to identify key locations in
  // the G4DynamicParticle chain that will be passed to the G4FissionEvent
  G4ReactionProduct* FirstDaughter = nullptr;
  G4ReactionProduct* SecondDaughter = nullptr;
  auto Alphas = new std::vector<G4ReactionProduct*>;
  auto Neutrons = new std::vector<G4ReactionProduct*>;
  auto Gammas = new std::vector<G4ReactionProduct*>;
 
  // Generate all the nucleonic fission products
  // How many nucleons do we have to work with?
  // TK modified 131108
  const G4int ParentA = (Isotope_ / 10) % 1000;
  const G4int ParentZ = ((Isotope_ / 10) - ParentA) / 1000;
  RemainingA_ = ParentA;
  RemainingZ_ = ParentZ;
 
  // Don't forget the extra nucleons depending on the fission cause
  switch (Cause_) {
    case G4FFGEnumerations::NEUTRON_INDUCED:
      ++RemainingA_;
      break;
 
    case G4FFGEnumerations::PROTON_INDUCED:
      ++RemainingZ_;
      break;
 
    case G4FFGEnumerations::GAMMA_INDUCED:
    case G4FFGEnumerations::SPONTANEOUS:
    default:
      // Nothing to do here
      break;
  }
 
  // Ternary fission can be set by the user. Thus, it is necessary to
  // sample the alpha particle first and the first daughter product
  // second. See the discussion in
  // G4FissionProductYieldDist::G4GetFissionProduct() for more information
  // as to why the fission events are sampled this way.
  GenerateAlphas(Alphas);
 
  // Generate the first daughter product
  FirstDaughter = new G4ReactionProduct(GetFissionProduct());
  RemainingA_ -= FirstDaughter->GetDefinition()->GetAtomicMass();
  RemainingZ_ -= FirstDaughter->GetDefinition()->GetAtomicNumber();
  if ((Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO) != 0) {
    G4FFG_SPACING__
    G4FFG_LOCATION__
 
    G4cout << " -- First daughter product sampled" << G4endl;
    G4FFG_SPACING__
    G4cout << "  Name:       " << FirstDaughter->GetDefinition()->GetParticleName() << G4endl;
    G4FFG_SPACING__
    G4cout << "  Z:          " << FirstDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
    G4FFG_SPACING__
    G4cout << "  A:          " << FirstDaughter->GetDefinition()->GetAtomicMass() << G4endl;
    G4FFG_SPACING__
    G4cout << "  Meta State: " << (FirstDaughter->GetDefinition()->GetPDGEncoding() % 10) << G4endl;
  }
 
  GenerateNeutrons(Neutrons);
 
  // Now that all the nucleonic particles have been generated, we can
  // calculate the composition of the second daughter product.
  G4int NewIsotope = RemainingZ_ * 1000 + RemainingA_;
  SecondDaughter =
    new G4ReactionProduct(GetParticleDefinition(NewIsotope, G4FFGEnumerations::GROUND_STATE));
  if ((Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO) != 0) {
    G4FFG_SPACING__
    G4FFG_LOCATION__
 
    G4cout << " -- Second daughter product sampled" << G4endl;
    G4FFG_SPACING__
    G4cout << "  Name:       " << SecondDaughter->GetDefinition()->GetParticleName() << G4endl;
    G4FFG_SPACING__
    G4cout << "  Z:          " << SecondDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
    G4FFG_SPACING__
    G4cout << "  A:          " << SecondDaughter->GetDefinition()->GetAtomicMass() << G4endl;
    G4FFG_SPACING__
    G4cout << "  Meta State: " << (SecondDaughter->GetDefinition()->GetPDGEncoding() % 10)
           << G4endl;
  }
 
  // Calculate how much kinetic energy will be available
  // 195 to 205 MeV are available in a fission reaction, but about 20 MeV
  // are from delayed sources. We are concerned only with prompt sources,
  // so sample a Gaussian distribution about 20 MeV and subtract the
  // result from the total available energy. Also, the energy of fission
  // neutrinos is neglected. Fission neutrinos would add ~11 MeV
  // additional energy to the fission. (Ref 2)
  // Finally, add in the kinetic energy contribution of the fission
  // inducing particle, if any.
  const G4double TotalKE = RandomEngine_->G4SampleUniform(195.0, 205.0) * MeV
                           - RandomEngine_->G4SampleGaussian(20.0, 3.0) * MeV + IncidentEnergy_;
  RemainingEnergy_ = TotalKE;
 
  // Calculate the energies of the alpha particles and neutrons
  // Once again, since the alpha particles are user defined, we must
  // sample their kinetic energy first. SampleAlphaEnergies() returns the
  // amount of energy consumed by the alpha particles, so remove the total
  // energy alloted to the alpha particles from the available energy
  SampleAlphaEnergies(Alphas);
 
  // Second, the neutrons are sampled from the Watt fission spectrum.
  SampleNeutronEnergies(Neutrons);
 
  // Calculate the total energy available to the daughter products
  // A Gaussian distribution about the average calculated energy with
  // a standard deviation of 1.5 MeV (Ref. 2) is used. Since the energy
  // distribution is dependant on the alpha particle generation and the
  // Watt fission sampling for neutrons, we only have the left-over energy
  // to work with for the fission daughter products.
  G4double FragmentsKE = 0.;
  G4int icounter = 0;
  G4int icounter_max = 1024;
  do {
    icounter++;
    if (icounter > icounter_max) {
      G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
             << __FILE__ << "." << G4endl;
      break;
    }
    FragmentsKE = RandomEngine_->G4SampleGaussian(RemainingEnergy_, 1.5 * MeV);
  } while (FragmentsKE > RemainingEnergy_);  // Loop checking, 11.05.2015, T. Koi
 
  // Make sure that we don't produce any sub-gamma photons
  if ((RemainingEnergy_ - FragmentsKE) / (100 * keV) < 1.0) {
    FragmentsKE = RemainingEnergy_;
  }
 
  // This energy has now been allotted to the fission fragments.
  // Subtract FragmentsKE from the total available fission energy.
  RemainingEnergy_ -= FragmentsKE;
 
  // Sample the energies of the gamma rays
  // Assign the remainder, if any, of the energy to the gamma rays
  SampleGammaEnergies(Gammas);
 
  // Prepare to balance the momenta of the system
  // We will need these for sampling the angles and balancing the momenta
  // of all the fission products
  G4double NumeratorSqrt, NumeratorOther, Denominator;
  G4double Theta, Phi, Magnitude;
  G4ThreeVector Direction;
  G4ParticleMomentum ResultantVector(0, 0, 0);
 
  if (!Alphas->empty()) {
    // Sample the angles of the alpha particles and neutrons, then calculate
    // the total moment contribution to the system
    // The average angle of the alpha particles with respect to the
    // light fragment is dependent on the ratio of the kinetic energies.
    // This equation was determined by performing a fit on data from
    // Ref. 3 using the website:
    // http://soft.arquimedex.com/linear_regression.php
    //
    // RatioOfKE    Angle (rad)     Angle (degrees)
    // 1.05         1.257           72
    // 1.155        1.361           78
    // 1.28         1.414           81
    // 1.5          1.518           87
    // 1.75         1.606           92
    // 1.9          1.623           93
    // 2.2          1.728           99
    // This equation generates the angle in radians. If the RatioOfKE is
    // greater than 2.25 the angle defaults to 1.3963 rad (100 degrees)
    G4double MassRatio =
      FirstDaughter->GetDefinition()->GetPDGMass() / SecondDaughter->GetDefinition()->GetPDGMass();
 
    // Invert the mass ratio if the first daughter product is the lighter fragment
    if (MassRatio < 1) {
      MassRatio = 1 / MassRatio;
    }
 
    // The empirical equation is valid for mass ratios up to 2.75
    if (MassRatio > 2.75) {
      MassRatio = 2.75;
    }
    const G4double MeanAlphaAngle = 0.3644 * MassRatio * MassRatio * MassRatio
                                    - 1.9766 * MassRatio * MassRatio + 3.8207 * MassRatio - 0.9917;
 
    // Sample the directions of the alpha particles with respect to the
    // light fragment. For the moment we will assume that the light
    // fragment is traveling along the z-axis in the positive direction.
    const G4double MeanAlphaAngleStdDev = 0.0523598776;
    G4double PlusMinus;
 
    for (auto& Alpha : *Alphas) {
      PlusMinus = std::acos(RandomEngine_->G4SampleGaussian(0, MeanAlphaAngleStdDev)) - (pi / 2);
      Theta = MeanAlphaAngle + PlusMinus;
      if (Theta < 0) {
        Theta = 0.0 - Theta;
      }
      else if (Theta > pi) {
        Theta = (2 * pi - Theta);
      }
      Phi = RandomEngine_->G4SampleUniform(-pi, pi);
 
      Direction.setRThetaPhi(1.0, Theta, Phi);
      Magnitude = std::sqrt(2 * Alpha->GetKineticEnergy() * Alpha->GetDefinition()->GetPDGMass());
      Alpha->SetMomentum(Direction * Magnitude);
      ResultantVector += Alpha->GetMomentum();
    }
  }
 
  // Sample the directions of the neutrons.
  if (!Neutrons->empty()) {
    for (auto& Neutron : *Neutrons) {
      Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
      Phi = RandomEngine_->G4SampleUniform(-pi, pi);
 
      Direction.setRThetaPhi(1.0, Theta, Phi);
      Magnitude =
        std::sqrt(2 * Neutron->GetKineticEnergy() * Neutron->GetDefinition()->GetPDGMass());
      Neutron->SetMomentum(Direction * Magnitude);
      ResultantVector += Neutron->GetMomentum();
    }
  }
 
  // Sample the directions of the gamma rays
  if (!Gammas->empty()) {
    for (auto& Gamma : *Gammas) {
      Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
      Phi = RandomEngine_->G4SampleUniform(-pi, pi);
 
      Direction.setRThetaPhi(1.0, Theta, Phi);
      Magnitude = Gamma->GetKineticEnergy() / CLHEP::c_light;
      Gamma->SetMomentum(Direction * Magnitude);
      ResultantVector += Gamma->GetMomentum();
    }
  }
 
  // Calculate the momenta of the two daughter products
  G4ReactionProduct* LightFragment;
  G4ReactionProduct* HeavyFragment;
  G4ThreeVector LightFragmentDirection;
  G4ThreeVector HeavyFragmentDirection;
  G4double ResultantX, ResultantY, ResultantZ;
  ResultantX = ResultantVector.getX();
  ResultantY = ResultantVector.getY();
  ResultantZ = ResultantVector.getZ();
 
  if (FirstDaughter->GetDefinition()->GetPDGMass() < SecondDaughter->GetDefinition()->GetPDGMass())
  {
    LightFragment = FirstDaughter;
    HeavyFragment = SecondDaughter;
  }
  else {
    LightFragment = SecondDaughter;
    HeavyFragment = FirstDaughter;
  }
  const G4double LightFragmentMass = LightFragment->GetDefinition()->GetPDGMass();
  const G4double HeavyFragmentMass = HeavyFragment->GetDefinition()->GetPDGMass();
 
  LightFragmentDirection.setRThetaPhi(1.0, 0, 0);
 
  // Fit the momenta of the daughter products to the resultant vector of
  // the remaining fission products. This will be done in the Cartesian
  // coordinate system, not spherical. This is done using the following
  // table listing the system momenta and the corresponding equations:
  //              X               Y               Z
  //
  //      A       0               0               P
  //
  //      B       -R_x            -R_y            -P - R_z
  //
  //      R       R_x             R_y             R_z
  //
  // v = sqrt(2*m*k)  ->  k = v^2/(2*m)
  // tk = k_A + k_B
  // k_L = P^2/(2*m_L)
  // k_H = ((-R_x)^2 + (-R_y)^2 + (-P - R_z)^2)/(2*m_H)
  // where:
  // P: momentum of the light daughter product
  // R: the remaining fission products' resultant vector
  // v: momentum
  // m: mass
  // k: kinetic energy
  // tk: total kinetic energy available to the daughter products
  //
  // Below is the solved form for P, with the solution generated using
  // the WolframAlpha website:
  // http://www.wolframalpha.com/input/?i=
  // solve+((-x)^2+%2B+(-y)^2+%2B+(-P-z)^2)%2F(2*B)+%2B+L^2%2F(2*A)+%3D+k+
  // for+P
  //
  //
  // nsqrt = sqrt(m_L*(m_L*(2*m_H*tk - R_x^2 - R_y^2) +
  //                   m_H*(2*m_H*tk - R_x^2 - R_y^2 - R_z^2))
  NumeratorSqrt = std::sqrt(
    LightFragmentMass
    * (LightFragmentMass
         * (2 * HeavyFragmentMass * FragmentsKE - ResultantX * ResultantX - ResultantY * ResultantY)
       + HeavyFragmentMass
           * (2 * HeavyFragmentMass * FragmentsKE - ResultantX * ResultantX
              - ResultantY * ResultantY - ResultantZ - ResultantZ)));
 
  // nother = m_L*R_z
  NumeratorOther = LightFragmentMass * ResultantZ;
 
  // denom = m_L + m_H
  Denominator = LightFragmentMass + HeavyFragmentMass;
 
  // P = (nsqrt + nother) / denom
  const G4double LightFragmentMomentum = (NumeratorSqrt + NumeratorOther) / Denominator;
  const G4double HeavyFragmentMomentum =
    std::sqrt(ResultantX * ResultantX + ResultantY * ResultantY
              + G4Pow::GetInstance()->powN(LightFragmentMomentum + ResultantZ, 2));
 
  // Finally! We now have everything we need for the daughter products
  LightFragment->SetMomentum(LightFragmentDirection * LightFragmentMomentum);
  HeavyFragmentDirection.setX(-ResultantX);
  HeavyFragmentDirection.setY(-ResultantY);
  HeavyFragmentDirection.setZ((-LightFragmentMomentum - ResultantZ));
  // Don't forget to normalize the vector to the unit sphere
  HeavyFragmentDirection.setR(1.0);
  HeavyFragment->SetMomentum(HeavyFragmentDirection * HeavyFragmentMomentum);
 
  if ((Verbosity_ & (G4FFGEnumerations::DAUGHTER_INFO | G4FFGEnumerations::MOMENTUM_INFO)) != 0) {
    G4FFG_SPACING__
    G4FFG_LOCATION__
 
    G4cout << " -- Daugher product momenta finalized" << G4endl;
    G4FFG_SPACING__
  }
 
  // Load all the particles into a contiguous set
  // TK modifed 131108
  // G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector(2 + Alphas->size() +
  // Neutrons->size() + Gammas->size());
  auto FissionProducts = new G4DynamicParticleVector();
  // Load the fission fragments
  FissionProducts->push_back(MakeG4DynamicParticle(LightFragment));
  FissionProducts->push_back(MakeG4DynamicParticle(HeavyFragment));
  // Load the neutrons
  for (auto& Neutron : *Neutrons) {
    FissionProducts->push_back(MakeG4DynamicParticle(Neutron));
  }
  // Load the gammas
  for (auto& Gamma : *Gammas) {
    FissionProducts->push_back(MakeG4DynamicParticle(Gamma));
  }
  // Load the alphas
  for (auto& Alpha : *Alphas) {
    FissionProducts->push_back(MakeG4DynamicParticle(Alpha));
  }
 
  // Rotate the system to a random location so that the light fission fragment
  // is not always traveling along the positive z-axis
  // Sample Theta and Phi.
  G4ThreeVector RotationAxis;
 
  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
  Phi = RandomEngine_->G4SampleUniform(-pi, pi);
  RotationAxis.setRThetaPhi(1.0, Theta, Phi);
 
  // We will also check the net momenta
  ResultantVector.set(0.0, 0.0, 0.0);
  for (auto& FissionProduct : *FissionProducts) {
    Direction = FissionProduct->GetMomentumDirection();
    Direction.rotateUz(RotationAxis);
    FissionProduct->SetMomentumDirection(Direction);
    ResultantVector += FissionProduct->GetMomentum();
  }
 
  // Warn if the sum momenta of the system is not within a reasonable
  // tolerance
  G4double PossibleImbalance = ResultantVector.mag();
  if (PossibleImbalance > 0.01) {
    std::ostringstream Temp;
    Temp << "Momenta imbalance of ";
    Temp << PossibleImbalance / (MeV / CLHEP::c_light);
    Temp << " MeV/c in the system";
    G4Exception("G4FissionProductYieldDist::G4GetFission()", Temp.str().c_str(), JustWarning,
                "Results may not be valid");
  }
 
  // Clean up
  delete FirstDaughter;
  delete SecondDaughter;
  G4ArrayOps::DeleteVectorOfPointers(*Neutrons);
  G4ArrayOps::DeleteVectorOfPointers(*Gammas);
  G4ArrayOps::DeleteVectorOfPointers(*Alphas);
 
  G4FFG_FUNCTIONLEAVE__
  return FissionProducts;
}
 
G4Ions* G4FissionProductYieldDist::G4GetFissionProduct()
{
  G4FFG_FUNCTIONENTER__
 
  G4Ions* Product = FindParticle(RandomEngine_->G4SampleUniform());
 
  G4FFG_FUNCTIONLEAVE__
  return Product;
}
 
void G4FissionProductYieldDist::G4SetAlphaProduction(G4double WhatAlphaProduction)
{
  G4FFG_FUNCTIONENTER__
 
  AlphaProduction_ = WhatAlphaProduction;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::G4SetEnergy(G4double WhatIncidentEnergy)
{
  G4FFG_FUNCTIONENTER__
 
  if (Cause_ != G4FFGEnumerations::SPONTANEOUS) {
    IncidentEnergy_ = WhatIncidentEnergy;
  }
  else {
    IncidentEnergy_ = 0 * GeV;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::G4SetTernaryProbability(G4double WhatTernaryProbability)
{
  G4FFG_FUNCTIONENTER__
 
  TernaryProbability_ = WhatTernaryProbability;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::G4SetVerbosity(G4int Verbosity)
{
  G4FFG_FUNCTIONENTER__
 
  Verbosity_ = Verbosity;
 
  ENDFData_->G4SetVerbosity(Verbosity_);
  RandomEngine_->G4SetVerbosity(Verbosity_);
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::CheckAlphaSanity()
{
  G4FFG_FUNCTIONENTER__
 
  // This provides comfortable breathing room at 16 MeV per alpha
  if (AlphaProduction_ > 10) {
    AlphaProduction_ = 10;
  }
  else if (AlphaProduction_ < -7) {
    AlphaProduction_ = -7;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
G4Ions* G4FissionProductYieldDist::FindParticle(G4double RandomParticle)
{
  G4FFG_FUNCTIONENTER__
 
  // Determine which energy group is currently in use
  G4bool isExact = false;
  G4bool lowerExists = false;
  G4bool higherExists = false;
  G4int energyGroup;
  for (energyGroup = 0; energyGroup < YieldEnergyGroups_; energyGroup++) {
    if (IncidentEnergy_ == YieldEnergies_[energyGroup]) {
      isExact = true;
      break;
    }
 
    if (energyGroup == 0 && IncidentEnergy_ < YieldEnergies_[energyGroup]) {
      // Break if the energy is less than the lowest energy
      higherExists = true;
      break;
    }
    if (energyGroup == YieldEnergyGroups_ - 1) {
      // The energy is greater than any values in the yield data.
      lowerExists = true;
      break;
    }
    // Break if the energy is less than the lowest energy
    if (IncidentEnergy_ > YieldEnergies_[energyGroup]) {
      energyGroup--;
      lowerExists = true;
      higherExists = true;
      break;
    }
  }
 
  // Determine which particle it is
  G4Ions* FoundParticle = nullptr;
  if (isExact || YieldEnergyGroups_ == 1) {
    // Determine which tree contains the random value
    G4int tree;
    for (tree = 0; tree < TreeCount_; tree++) {
      // Break if a tree is identified as containing the random particle
      if (RandomParticle <= Trees_[tree].ProbabilityRangeEnd[energyGroup]) {
        break;
      }
    }
    ProbabilityBranch* Branch = Trees_[tree].Trunk;
 
    // Iteratively traverse the tree until the particle addressed by the random
    // variable is found
    G4bool RangeIsSmaller;
    G4bool RangeIsGreater;
    while ((RangeIsSmaller = (RandomParticle < Branch->ProbabilityRangeBottom[energyGroup]))
           || (RangeIsGreater = (RandomParticle > Branch->ProbabilityRangeTop[energyGroup])))
    // Loop checking, 11.05.2015, T. Koi
    {
      if (RangeIsSmaller) {
        Branch = Branch->Left;
      }
      else {
        Branch = Branch->Right;
      }
    }
 
    FoundParticle = Branch->Particle;
  }
  else if (lowerExists && higherExists) {
    // We need to do some interpolation
    FoundParticle = FindParticleInterpolation(RandomParticle, energyGroup);
  }
  else {
    // We need to do some extrapolation
    FoundParticle = FindParticleExtrapolation(RandomParticle, lowerExists);
  }
 
  // Return the particle
  G4FFG_FUNCTIONLEAVE__
  return FoundParticle;
}
 
G4Ions* G4FissionProductYieldDist::FindParticleExtrapolation(G4double RandomParticle,
                                                             G4bool LowerEnergyGroupExists)
{
  G4FFG_FUNCTIONENTER__
 
  G4Ions* FoundParticle = nullptr;
  G4int NearestEnergy;
  G4int NextNearestEnergy;
 
  // Check to see if we are extrapolating above or below the data set
  if (LowerEnergyGroupExists) {
    NearestEnergy = YieldEnergyGroups_ - 1;
    NextNearestEnergy = NearestEnergy - 1;
  }
  else {
    NearestEnergy = 0;
    NextNearestEnergy = 1;
  }
 
  for (G4int Tree = 0; Tree < TreeCount_ && FoundParticle == nullptr; Tree++) {
    FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk, RandomParticle, NearestEnergy,
                                             NextNearestEnergy);
  }
 
  G4FFG_FUNCTIONLEAVE__
  return FoundParticle;
}
 
G4Ions* G4FissionProductYieldDist::FindParticleInterpolation(G4double RandomParticle,
                                                             G4int LowerEnergyGroup)
{
  G4FFG_FUNCTIONENTER__
 
  G4Ions* FoundParticle = nullptr;
  G4int HigherEnergyGroup = LowerEnergyGroup + 1;
 
  for (G4int Tree = 0; Tree < TreeCount_ && FoundParticle == nullptr; Tree++) {
    FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk, RandomParticle, LowerEnergyGroup,
                                             HigherEnergyGroup);
  }
 
  G4FFG_FUNCTIONLEAVE__
  return FoundParticle;
}
 
G4Ions* G4FissionProductYieldDist::FindParticleBranchSearch(ProbabilityBranch* Branch,
                                                            G4double RandomParticle,
                                                            G4int EnergyGroup1, G4int EnergyGroup2)
{
  G4FFG_RECURSIVE_FUNCTIONENTER__
 
  G4Ions* Particle;
 
  // Verify that the branch exists
  if (Branch == nullptr) {
    Particle = nullptr;
  }
  else if (EnergyGroup1 >= Branch->IncidentEnergiesCount
           || EnergyGroup2 >= Branch->IncidentEnergiesCount || EnergyGroup1 == EnergyGroup2
           || Branch->IncidentEnergies[EnergyGroup1] == Branch->IncidentEnergies[EnergyGroup2])
  {
    // Set NULL if any invalid conditions exist
    Particle = nullptr;
  }
  else {
    // Everything check out - proceed
    G4Ions* FoundParticle = nullptr;
    G4double Intercept;
    G4double Slope;
    G4double RangeAtIncidentEnergy;
    G4double Denominator =
      Branch->IncidentEnergies[EnergyGroup1] - Branch->IncidentEnergies[EnergyGroup2];
 
    // Calculate the lower probability bounds
    Slope =
      (Branch->ProbabilityRangeBottom[EnergyGroup1] - Branch->ProbabilityRangeBottom[EnergyGroup2])
      / Denominator;
    Intercept =
      Branch->ProbabilityRangeBottom[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
    RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
 
    // Go right if the particle is below the probability bounds
    if (RandomParticle < RangeAtIncidentEnergy) {
      FoundParticle =
        FindParticleBranchSearch(Branch->Left, RandomParticle, EnergyGroup1, EnergyGroup2);
    }
    else {
      // Calculate the upper probability bounds
      Slope =
        (Branch->ProbabilityRangeTop[EnergyGroup1] - Branch->ProbabilityRangeTop[EnergyGroup2])
        / Denominator;
      Intercept =
        Branch->ProbabilityRangeTop[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
      RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
 
      // Go left if the particle is above the probability bounds
      if (RandomParticle > RangeAtIncidentEnergy) {
        FoundParticle =
          FindParticleBranchSearch(Branch->Right, RandomParticle, EnergyGroup1, EnergyGroup2);
      }
      else {
        // If the particle is bounded then we found it!
        FoundParticle = Branch->Particle;
      }
    }
 
    Particle = FoundParticle;
  }
 
  G4FFG_RECURSIVE_FUNCTIONLEAVE__
  return Particle;
}
 
void G4FissionProductYieldDist::GenerateAlphas(std::vector<G4ReactionProduct*>* Alphas)
{
  G4FFG_FUNCTIONENTER__
 
  // Throw the dice to determine if ternary fission occurs
  G4bool MakeAlphas = RandomEngine_->G4SampleUniform() <= TernaryProbability_;
  if (MakeAlphas) {
    G4int NumberOfAlphasToProduce;
 
    // Determine how many alpha particles to produce for the ternary fission
    if (AlphaProduction_ < 0) {
      NumberOfAlphasToProduce = RandomEngine_->G4SampleIntegerGaussian(AlphaProduction_ * -1, 1,
                                                                       G4FFGEnumerations::POSITIVE);
    }
    else {
      NumberOfAlphasToProduce = (G4int)AlphaProduction_;
    }
 
    // TK modifed 131108
    // Alphas->resize(NumberOfAlphasToProduce);
    for (int i = 0; i < NumberOfAlphasToProduce; i++) {
      // Set the G4Ions as an alpha particle
      Alphas->push_back(new G4ReactionProduct(AlphaDefinition_));
 
      // Remove 4 nucleons (2 protons and 2 neutrons) for each alpha added
      RemainingZ_ -= 2;
      RemainingA_ -= 4;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::GenerateNeutrons(std::vector<G4ReactionProduct*>* Neutrons)
{
  G4FFG_FUNCTIONENTER__
 
  G4int NeutronProduction;
  NeutronProduction =
    RandomEngine_->G4SampleIntegerGaussian(Nubar_, NubarWidth_, G4FFGEnumerations::POSITIVE);
 
  // TK modifed 131108
  // Neutrons->resize(NeutronProduction);
  for (int i = 0; i < NeutronProduction; i++) {
    // Define the fragment as a neutron
    Neutrons->push_back(new G4ReactionProduct(NeutronDefinition_));
 
    // Remove 1 nucleon for each neutron added
    RemainingA_--;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
G4Ions* G4FissionProductYieldDist::GetParticleDefinition(G4int Product,
                                                         // TK modified 131108
                                                         // G4FFGEnumerations::MetaState MetaState )
                                                         G4FFGEnumerations::MetaState /*MetaState*/)
{
  G4FFG_DATA_FUNCTIONENTER__
 
  G4Ions* Temp;
 
  // Break Product down into its A and Z components
  G4int A = Product % 1000;  // Extract A
  G4int Z = (Product - A) / 1000;  // Extract Z
 
  // Check to see if the particle is registered using the PDG code
  // TODO Add metastable state when supported by G4IonTable::GetIon()
  Temp = static_cast<G4Ions*>(IonTable_->GetIon(Z, A));
 
  // Removed in favor of the G4IonTable::GetIon() method
  //    // Register the particle if it does not exist
  //    if(Temp == NULL)
  //    {
  //        // Define the particle properties
  //        G4String Name = MakeIsotopeName(Product, MetaState);
  //        // Calculate the rest mass using a function already in Geant4
  //        G4double Mass = G4NucleiProperties::
  //                        GetNuclearMass((double)A, (double)Z );
  //        G4double Charge = Z*eplus;
  //        G4int BaryonNum = A;
  //        G4bool Stable = TRUE;
  //
  //        // I am unsure about the following properties:
  //        //     2*Spin, Parity, C-conjugation, 2*Isospin, 2*Isospin3, G-parity.
  //        // Perhaps is would be a good idea to have a physicist familiar with
  //        // Geant4 nomenclature to review and correct these properties.
  //        Temp = new G4Ions (
  //            // Name           Mass           Width          Charge
  //               Name,          Mass,          0.0,           Charge,
  //
  //            // 2*Spin         Parity         C-conjugation  2*Isospin
  //               0,             1,             0,             0,
  //
  //            // 2*Isospin3     G-parity       Type           Lepton number
  //               0,             0,             "nucleus",     0,
  //
  //            // Baryon number  PDG encoding   Stable         Lifetime
  //               BaryonNum,     PDGCode,       Stable,        -1,
  //
  //            // Decay table    Shortlived     SubType        Anti_encoding
  //               NULL,          FALSE,        "generic",     0,
  //
  //            // Excitation
  //               0.0);
  //        Temp->SetPDGMagneticMoment(0.0);
  //
  //        // Declare that there is no anti-particle
  //        Temp->SetAntiPDGEncoding(0);
  //
  //        // Define the processes to use in transporting the particles
  //        std::ostringstream osAdd;
  //        osAdd << "/run/particle/addProcManager " << Name;
  //        G4String cmdAdd = osAdd.str();
  //
  //        // set /control/verbose 0
  //        G4int tempVerboseLevel = G4UImanager::GetUIpointer()->GetVerboseLevel();
  //        G4UImanager::GetUIpointer()->SetVerboseLevel(0);
  //
  //        // issue /run/particle/addProcManage
  //        G4UImanager::GetUIpointer()->ApplyCommand(cmdAdd);
  //
  //        // retrieve  /control/verbose
  //        G4UImanager::GetUIpointer()->SetVerboseLevel(tempVerboseLevel);
  //    }
 
  G4FFG_DATA_FUNCTIONLEAVE__
  return Temp;
}
 
G4String G4FissionProductYieldDist::MakeDirectoryName()
{
  G4FFG_FUNCTIONENTER__
 
  // Generate the file location starting in the Geant4 data directory
  std::ostringstream DirectoryName;
  DirectoryName << G4FindDataDir("G4NEUTRONHPDATA") << G4FFGDefaultValues::ENDFFissionDataLocation;
 
  // Return the directory structure
  G4FFG_FUNCTIONLEAVE__
  return DirectoryName.str();
}
 
G4String G4FissionProductYieldDist::MakeFileName(G4int Isotope,
                                                 G4FFGEnumerations::MetaState MetaState)
{
  G4FFG_FUNCTIONENTER__
 
  // Create the unique identifying name for the particle
  std::ostringstream FileName;
 
  // Determine if a leading 0 is needed (ZZZAAA or 0ZZAAA)
  if (Isotope < 100000) {
    FileName << "0";
  }
 
  // Add the name of the element and the extension
  FileName << MakeIsotopeName(Isotope, MetaState) << ".fpy";
 
  G4FFG_FUNCTIONLEAVE__
  return FileName.str();
}
 
G4DynamicParticle*
G4FissionProductYieldDist::MakeG4DynamicParticle(G4ReactionProduct* ReactionProduct)
{
  G4FFG_DATA_FUNCTIONENTER__
 
  auto DynamicParticle =
    new G4DynamicParticle(ReactionProduct->GetDefinition(), ReactionProduct->GetMomentum());
 
  G4FFG_DATA_FUNCTIONLEAVE__
  return DynamicParticle;
}
 
G4String G4FissionProductYieldDist::MakeIsotopeName(G4int Isotope,
                                                    G4FFGEnumerations::MetaState MetaState)
{
  G4FFG_DATA_FUNCTIONENTER__
 
  // Break Product down into its A and Z components
  G4int A = Isotope % 1000;
  G4int Z = (Isotope - A) / 1000;
 
  // Create the unique identifying name for the particle
  std::ostringstream IsotopeName;
 
  IsotopeName << Z << "_" << A;
 
  // If it is metastable then append "m" to the name
  if (MetaState != G4FFGEnumerations::GROUND_STATE) {
    IsotopeName << "m";
 
    // If it is a second isomeric state then append "2" to the name
    if (MetaState == G4FFGEnumerations::META_2) {
      IsotopeName << "2";
    }
  }
  // Add the name of the element and the extension
  IsotopeName << "_" << ElementNames_->GetName(Z - 1);
 
  G4FFG_DATA_FUNCTIONLEAVE__
  return IsotopeName.str();
}
 
void G4FissionProductYieldDist::MakeTrees()
{
  G4FFG_FUNCTIONENTER__
 
  // Allocate the space
  // We will make each tree a binary search
  // The maximum number of iterations required to find a single fission product
  // based on it's probability is defined by the following:
  //      x = number of fission products
  //      Trees       = T(x)  = ceil( ln(x) )
  //      Rows/Tree   = R(x)  = ceil(( sqrt( (8 * x / T(x)) + 1) - 1) / 2)
  //      Maximum     = M(x)  = T(x) + R(x)
  //      Results: x    =>    M(x)
  //               10         5
  //               100        10
  //               1000       25
  //               10000      54
  //               100000     140
  TreeCount_ = (G4int)ceil((G4double)log((G4double)ENDFData_->G4GetNumberOfFissionProducts()));
  Trees_ = new ProbabilityTree[TreeCount_];
 
  // Initialize the range of each node
  for (G4int i = 0; i < TreeCount_; i++) {
    Trees_[i].ProbabilityRangeEnd = new G4double[YieldEnergyGroups_];
    Trees_[i].Trunk = nullptr;
    Trees_[i].BranchCount = 0;
    Trees_[i].IsEnd = FALSE;
  }
  // Mark the last tree as the ending tree
  Trees_[TreeCount_ - 1].IsEnd = TRUE;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::ReadProbabilities()
{
  G4FFG_DATA_FUNCTIONENTER__
 
  G4int ProductCount = ENDFData_->G4GetNumberOfFissionProducts();
  BranchCount_ = 0;
  G4ArrayOps::Set(YieldEnergyGroups_, DataTotal_, 0.0);
 
  // Loop through all the products
  for (G4int i = 0; i < ProductCount; i++) {
    // Acquire the data and sort it
    SortProbability(ENDFData_->G4GetYield(i));
  }
 
  // Generate the true normalization factor, since round-off errors may result
  // in non-singular normalization of the data files. Also, reset DataTotal_
  // since it is used by Renormalize() to set the probability segments.
  G4ArrayOps::Divide(YieldEnergyGroups_, MaintainNormalizedData_, 1.0, DataTotal_);
  G4ArrayOps::Set(YieldEnergyGroups_, DataTotal_, 0.0);
 
  // Go through all the trees one at a time
  for (G4int i = 0; i < TreeCount_; i++) {
    Renormalize(Trees_[i].Trunk);
    // Set the max range of the tree to DataTotal
    G4ArrayOps::Copy(YieldEnergyGroups_, Trees_[i].ProbabilityRangeEnd, DataTotal_);
  }
 
  G4FFG_DATA_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::Renormalize(ProbabilityBranch* Branch)
{
  G4FFG_RECURSIVE_FUNCTIONENTER__
 
  // Check to see if Branch exists. Branch will be a null pointer if it
  // doesn't exist
  if (Branch != nullptr) {
    // Call the lower branch to set the probability segment first, since it
    // supposed to have a lower probability segment that this node
    Renormalize(Branch->Left);
 
    // Set this node as the next sequential probability segment
    G4ArrayOps::Copy(YieldEnergyGroups_, Branch->ProbabilityRangeBottom, DataTotal_);
    G4ArrayOps::Multiply(YieldEnergyGroups_, Branch->ProbabilityRangeTop, MaintainNormalizedData_);
    G4ArrayOps::Add(YieldEnergyGroups_, Branch->ProbabilityRangeTop, DataTotal_);
    G4ArrayOps::Copy(YieldEnergyGroups_, DataTotal_, Branch->ProbabilityRangeTop);
 
    // Now call the upper branch to set those probability segments
    Renormalize(Branch->Right);
  }
 
  G4FFG_RECURSIVE_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::SampleAlphaEnergies(std::vector<G4ReactionProduct*>* Alphas)
{
  G4FFG_FUNCTIONENTER__
 
  // The condition of sampling more energy from the fission products than is
  // alloted is statistically unfavorable, but it could still happen. The
  // do-while loop prevents such an occurrence from happening
  G4double MeanAlphaEnergy = 16.0;
  G4double TotalAlphaEnergy;
 
  do {
    G4double AlphaEnergy;
    TotalAlphaEnergy = 0;
 
    // Walk through the alpha particles one at a time and sample each's
    // energy
    for (auto& Alpha : *Alphas) {
      AlphaEnergy =
        RandomEngine_->G4SampleGaussian(MeanAlphaEnergy, 2.35, G4FFGEnumerations::POSITIVE) * MeV;
      // Assign the energy to the alpha particle
      Alpha->SetKineticEnergy(AlphaEnergy);
 
      // Add up the total amount of kinetic energy consumed.
      TotalAlphaEnergy += AlphaEnergy;
    }
 
    // If true, decrement the mean alpha energy by 0.1 and try again.
    MeanAlphaEnergy -= 0.1;
  } while (TotalAlphaEnergy >= RemainingEnergy_);  // Loop checking, 11.05.2015, T. Koi
 
  // Subtract the total amount of energy that was assigned.
  RemainingEnergy_ -= TotalAlphaEnergy;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::SampleGammaEnergies(std::vector<G4ReactionProduct*>* Gammas)
{
  G4FFG_FUNCTIONENTER__
 
  // Make sure that there is energy to assign to the gamma rays
  if (RemainingEnergy_ != 0) {
    G4double SampleEnergy;
 
    // Sample from RemainingEnergy until it is all gone. Also,
    // RemainingEnergy should not be smaller than
    // G4FFGDefaultValues::MeanGammaEnergy. This will prevent the
    // sampling of a fractional portion of the Gaussian distribution
    // in an attempt to find a new gamma ray energy.
    G4int icounter = 0;
    G4int icounter_max = 1024;
    while (RemainingEnergy_
           >= G4FFGDefaultValues::MeanGammaEnergy)  // Loop checking, 11.05.2015, T. Koi
    {
      icounter++;
      if (icounter > icounter_max) {
        G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
               << __FILE__ << "." << G4endl;
        break;
      }
      SampleEnergy = RandomEngine_->G4SampleGaussian(G4FFGDefaultValues::MeanGammaEnergy, 1.0 * MeV,
                                                     G4FFGEnumerations::POSITIVE);
      // Make sure that we didn't sample more energy than was available
      if (SampleEnergy <= RemainingEnergy_) {
        // If this energy assignment would leave less energy than the
        // 'intrinsic' minimal energy of a gamma ray then just assign
        // all of the remaining energy
        if (RemainingEnergy_ - SampleEnergy < 100 * keV) {
          SampleEnergy = RemainingEnergy_;
        }
 
        // Create the new particle
        Gammas->push_back(new G4ReactionProduct());
 
        // Set the properties
        Gammas->back()->SetDefinition(GammaDefinition_);
        Gammas->back()->SetTotalEnergy(SampleEnergy);
 
        // Calculate how much is left
        RemainingEnergy_ -= SampleEnergy;
      }
    }
 
    // If there is anything left over, the energy must be above 100 keV but
    // less than G4FFGDefaultValues::MeanGammaEnergy. Arbitrarily assign
    // RemainingEnergy to a new particle
    if (RemainingEnergy_ > 0) {
      SampleEnergy = RemainingEnergy_;
      Gammas->push_back(new G4ReactionProduct());
 
      // Set the properties
      Gammas->back()->SetDefinition(GammaDefinition_);
      Gammas->back()->SetTotalEnergy(SampleEnergy);
 
      // Calculate how much is left
      RemainingEnergy_ -= SampleEnergy;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::SampleNeutronEnergies(std::vector<G4ReactionProduct*>* Neutrons)
{
  G4FFG_FUNCTIONENTER__
 
  // The condition of sampling more energy from the fission products than is
  // alloted is statistically unfavorable, but it could still happen. The
  // do-while loop prevents such an occurrence from happening
  G4double TotalNeutronEnergy = 0.;
  G4double NeutronEnergy = 0.;
 
  // Make sure that we don't sample more energy than is available
  G4int icounter = 0;
  G4int icounter_max = 1024;
  do {
    icounter++;
    if (icounter > icounter_max) {
      G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
             << __FILE__ << "." << G4endl;
      break;
    }
    TotalNeutronEnergy = 0;
 
    // Walk through the neutrons one at a time and sample the energies.
    // The gamma rays have not yet been sampled, so the last neutron will
    // have a NULL value for NextFragment
    for (auto& Neutron : *Neutrons) {
      // Assign the energy to the neutron
      NeutronEnergy = RandomEngine_->G4SampleWatt(Isotope_, Cause_, IncidentEnergy_);
      Neutron->SetKineticEnergy(NeutronEnergy);
 
      // Add up the total amount of kinetic energy consumed.
      TotalNeutronEnergy += NeutronEnergy;
    }
  } while (TotalNeutronEnergy > RemainingEnergy_);  // Loop checking, 11.05.2015, T. Koi
 
  // Subtract the total amount of energy that was assigned.
  RemainingEnergy_ -= TotalNeutronEnergy;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::SetNubar()
{
  G4FFG_FUNCTIONENTER__
 
  G4int* WhichNubar;
  G4int* NubarWidth;
  G4double XFactor, BFactor;
 
  if (Cause_ == G4FFGEnumerations::SPONTANEOUS) {
    WhichNubar = const_cast<G4int*>(&SpontaneousNubar_[0][0]);
    NubarWidth = const_cast<G4int*>(&SpontaneousNubarWidth_[0][0]);
  }
  else {
    WhichNubar = const_cast<G4int*>(&NeutronInducedNubar_[0][0]);
    NubarWidth = const_cast<G4int*>(&NeutronInducedNubarWidth_[0][0]);
  }
 
  XFactor = G4Pow::GetInstance()->powA(10.0, -13.0);
  BFactor = G4Pow::GetInstance()->powA(10.0, -4.0);
  Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor + *(WhichNubar + 2) * BFactor;
  while (*WhichNubar != -1)  // Loop checking, 11.05.2015, T. Koi
  {
    if (*WhichNubar == Isotope_) {
      Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor + *(WhichNubar + 2) * BFactor;
 
      break;
    }
    WhichNubar += 3;
  }
 
  XFactor = G4Pow::GetInstance()->powN((G4double)10, -6);
  NubarWidth_ = *(NubarWidth + 1) * XFactor;
  while (*WhichNubar != -1)  // Loop checking, 11.05.2015, T. Koi
  {
    if (*WhichNubar == Isotope_) {
      NubarWidth_ = *(NubarWidth + 1) * XFactor;
 
      break;
    }
    WhichNubar += 2;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::SortProbability(G4ENDFYieldDataContainer* YieldData)
{
  G4FFG_DATA_FUNCTIONENTER__
 
  // Initialize the new branch
  auto NewBranch = new ProbabilityBranch;
  NewBranch->IncidentEnergiesCount = YieldEnergyGroups_;
  NewBranch->Left = nullptr;
  NewBranch->Right = nullptr;
  NewBranch->Particle = GetParticleDefinition(YieldData->GetProduct(), YieldData->GetMetaState());
  NewBranch->IncidentEnergies = new G4double[YieldEnergyGroups_];
  NewBranch->ProbabilityRangeTop = new G4double[YieldEnergyGroups_];
  NewBranch->ProbabilityRangeBottom = new G4double[YieldEnergyGroups_];
  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->ProbabilityRangeTop,
                   YieldData->GetYieldProbability());
  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->IncidentEnergies, YieldEnergies_);
  G4ArrayOps::Add(YieldEnergyGroups_, DataTotal_, YieldData->GetYieldProbability());
 
  // Check to see if the this is the smallest/largest particle. First, check
  // to see if this is the first particle in the system
  if (SmallestZ_ == nullptr) {
    SmallestZ_ = SmallestA_ = LargestZ_ = LargestA_ = NewBranch->Particle;
  }
  else {
    G4bool IsSmallerZ = NewBranch->Particle->GetAtomicNumber() < SmallestZ_->GetAtomicNumber();
    G4bool IsSmallerA = NewBranch->Particle->GetAtomicMass() < SmallestA_->GetAtomicMass();
    G4bool IsLargerZ = NewBranch->Particle->GetAtomicNumber() > LargestZ_->GetAtomicNumber();
    G4bool IsLargerA = NewBranch->Particle->GetAtomicMass() > LargestA_->GetAtomicMass();
 
    if (IsSmallerZ) {
      SmallestZ_ = NewBranch->Particle;
    }
 
    if (IsLargerZ) {
      LargestA_ = NewBranch->Particle;
    }
 
    if (IsSmallerA) {
      SmallestA_ = NewBranch->Particle;
    }
 
    if (IsLargerA) {
      LargestA_ = NewBranch->Particle;
    }
  }
 
  // Place the new branch
  // Determine which tree the new branch goes into
  auto WhichTree = (G4int)floor((G4double)(BranchCount_ % TreeCount_));
  ProbabilityBranch** WhichBranch = &(Trees_[WhichTree].Trunk);
  Trees_[WhichTree].BranchCount++;
 
  // Search for the position
  // Determine where the branch goes
  G4int BranchPosition = (G4int)floor((G4double)(BranchCount_ / TreeCount_)) + 1;
 
  // Run through the tree until the end branch is reached
  while (BranchPosition > 1)  // Loop checking, 11.05.2015, T. Koi
  {
    if ((BranchPosition & 1) != 0) {
      // If the 1's bit is on then move to the next 'right' branch
      WhichBranch = &((*WhichBranch)->Right);
    }
    else {
      // If the 1's bit is off then move to the next 'down' branch
      WhichBranch = &((*WhichBranch)->Left);
    }
 
    BranchPosition >>= 1;
  }
 
  *WhichBranch = NewBranch;
  BranchCount_++;
 
  G4FFG_DATA_FUNCTIONLEAVE__
}
 
G4FissionProductYieldDist::~G4FissionProductYieldDist()
{
  G4FFG_FUNCTIONENTER__
 
  // Burn each tree, one by one
  G4int WhichTree = 0;
  while (static_cast<int>(Trees_[WhichTree].IsEnd) != TRUE)  // Loop checking, 11.05.2015, T. Koi
  {
    BurnTree(Trees_[WhichTree].Trunk);
    delete Trees_[WhichTree].Trunk;
    delete[] Trees_[WhichTree].ProbabilityRangeEnd;
    WhichTree++;
  }
 
  // Delete each dynamically allocated variable
  delete ENDFData_;
  delete[] Trees_;
  delete[] DataTotal_;
  delete[] MaintainNormalizedData_;
  delete ElementNames_;
  delete RandomEngine_;
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionProductYieldDist::BurnTree(ProbabilityBranch* Branch)
{
  G4FFG_RECURSIVE_FUNCTIONENTER__
 
  // Check to see it Branch exists. Branch will be a null pointer if it
  // doesn't exist
  if (Branch != nullptr) {
    // Burn down before you burn up
    BurnTree(Branch->Left);
    delete Branch->Left;
    BurnTree(Branch->Right);
    delete Branch->Right;
 
    delete[] Branch->IncidentEnergies;
    delete[] Branch->ProbabilityRangeTop;
    delete[] Branch->ProbabilityRangeBottom;
  }
 
  G4FFG_RECURSIVE_FUNCTIONLEAVE__
}