#include <G4FissionProductYieldDist.hh>

Inheritance diagram for G4FissionProductYieldDist:

Public Member Functions
	G4FissionProductYieldDist (G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, std::istringstream &dataStream)

	G4FissionProductYieldDist (G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, G4int Verbosity, std::istringstream &dataStream)

G4DynamicParticleVector *	G4GetFission ()

G4Ions *	G4GetFissionProduct ()

void	G4SetAlphaProduction (G4double WhatAlphaProduction)

void	G4SetEnergy (G4double WhatIncidentEnergy)

void	G4SetTernaryProbability (G4double TernaryProbability)

void	G4SetVerbosity (G4int WhatVerbosity)

virtual	~G4FissionProductYieldDist ()

Protected Member Functions
void	CheckAlphaSanity ()

G4Ions *	FindParticle (G4double RandomParticle)

G4Ions *	FindParticleExtrapolation (G4double RandomParticle, G4bool LowerEnergyGroupExists)

G4Ions *	FindParticleInterpolation (G4double RandomParticle, G4int LowerEnergyGroup)

G4Ions *	FindParticleBranchSearch (ProbabilityBranch *Branch, G4double RandomParticle, G4int EnergyGroup1, G4int EnergyGroup2)

virtual void	GenerateAlphas (std::vector< G4ReactionProduct * > *Alphas)

virtual void	GenerateNeutrons (std::vector< G4ReactionProduct * > *Neutrons)

virtual G4Ions *	GetFissionProduct ()=0

G4Ions *	GetParticleDefinition (G4int Product, G4FFGEnumerations::MetaState MetaState)

G4String	MakeDirectoryName ()

G4String	MakeFileName (G4int Isotope, G4FFGEnumerations::MetaState MetaState)

G4DynamicParticle *	MakeG4DynamicParticle (G4ReactionProduct *)

G4String	MakeIsotopeName (G4int Isotope, G4FFGEnumerations::MetaState MetaState)

virtual void	MakeTrees ()

virtual void	ReadProbabilities ()

void	Renormalize (ProbabilityBranch *Branch)

void	SampleAlphaEnergies (std::vector< G4ReactionProduct * > *Alphas)

void	SampleGammaEnergies (std::vector< G4ReactionProduct * > *Gammas)

void	SampleNeutronEnergies (std::vector< G4ReactionProduct * > *Neutrons)

void	SetNubar ()

virtual void	SortProbability (G4ENDFYieldDataContainer *YieldData)

void	BurnTree (ProbabilityBranch *Branch)

Protected Attributes
const G4int	Isotope_

const G4FFGEnumerations::MetaState	MetaState_

const G4FFGEnumerations::FissionCause	Cause_

const G4FFGEnumerations::YieldType	YieldType_

G4ENDFTapeRead *	ENDFData_

G4Ions *	AlphaDefinition_

G4double	AlphaProduction_

G4double	TernaryProbability_

G4Gamma *	GammaDefinition_

G4double	IncidentEnergy_

G4double	MeanGammaEnergy_

G4Ions *	NeutronDefinition_

G4double	Nubar_

G4double	NubarWidth_

G4int	RemainingZ_

G4int	RemainingA_

G4double	RemainingEnergy_

G4int	Verbosity_

ProbabilityTree *	Trees_

G4Ions *	SmallestZ_

G4Ions *	SmallestA_

G4Ions *	LargestZ_

G4Ions *	LargestA_

G4int	YieldEnergyGroups_

G4double *	YieldEnergies_

G4double *	MaintainNormalizedData_

G4double *	DataTotal_

G4int	TreeCount_

G4int	BranchCount_

G4IonTable *	IonTable_

G4ParticleHPNames *	ElementNames_

G4FPYSamplingOps *	RandomEngine_

Detailed Description

G4FissionProductYieldDist is the base class for storing all the fission data and generating fission events.

Definition at line 53 of file G4FissionProductYieldDist.hh.

Constructor & Destructor Documentation

◆ G4FissionProductYieldDist() [1/2]

G4FissionProductYieldDist::G4FissionProductYieldDist	(	G4int	WhichIsotope,
		G4FFGEnumerations::MetaState	WhichMetaState,
		G4FFGEnumerations::FissionCause	WhichCause,
		G4FFGEnumerations::YieldType	WhichYieldType,
		std::istringstream &	dataStream )

Default constructor

Usage:
- WhichIsotope: Isotope number of the element in ZZZAAA form
- WhichMetaState: GROUND_STATE, META_1, or META_2
- WhichCause: SPONTANEOUS or N_INDUCED
- WhichYieldType: INDEPENDENT or CUMULATIVE
Notes:

Definition at line 85 of file G4FissionProductYieldDist.cc.

  : 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() [2/2]

G4FissionProductYieldDist::G4FissionProductYieldDist	(	G4int	WhichIsotope,
		G4FFGEnumerations::MetaState	WhichMetaState,
		G4FFGEnumerations::FissionCause	WhichCause,
		G4FFGEnumerations::YieldType	WhichYieldType,
		G4int	Verbosity,
		std::istringstream &	dataStream )

Overloaded constructor

Usage:
- WhichIsotope: Isotope number of the element in ZZZAAA form
- WhichMetaState: GROUND_STATE, META_1, or META_2
- WhichCause: SPONTANEOUS or N_INDUCED
- WhichYieldType: INDEPENDENT or CUMULATIVE
- Verbosity: Verbosity level
Notes:

Definition at line 110 of file G4FissionProductYieldDist.cc.

  : 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__
}

◆ ~G4FissionProductYieldDist()

G4FissionProductYieldDist::~G4FissionProductYieldDist ( )

virtual

Default deconstructor. It is a virtual function since G4FissionProductYieldDist is a parent class

Definition at line 1403 of file G4FissionProductYieldDist.cc.

{
  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__
}

Member Function Documentation

◆ BurnTree()

void G4FissionProductYieldDist::BurnTree ( ProbabilityBranch * Branch )

protected

Recursively burns each branch in a probability tree.

Definition at line 1428 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by BurnTree(), and ~G4FissionProductYieldDist().

◆ CheckAlphaSanity()

void G4FissionProductYieldDist::CheckAlphaSanity ( )

protected

Checks to make sure that alpha overpopulation will not occur, which could result in an unsolvable zero momentum in the LAB system.

Definition at line 637 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ FindParticle()

G4Ions * G4FissionProductYieldDist::FindParticle ( G4double RandomParticle )

protected

Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0, 1] random number RandomParticle

Definition at line 652 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by G4GetFissionProduct(), G4FPYBiasedLightFragmentDist::GetFissionProduct(), and G4FPYNormalFragmentDist::GetFissionProduct().

◆ FindParticleBranchSearch()

G4Ions * G4FissionProductYieldDist::FindParticleBranchSearch	(	ProbabilityBranch *	Branch,
		G4double	RandomParticle,
		G4int	EnergyGroup1,
		G4int	EnergyGroup2 )

protected

Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0, 1] random number RandomParticle by searching through a branch. Both the extrapolation and interpolation schemes currently use this function to identify the particle.

Definition at line 776 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by FindParticleBranchSearch(), FindParticleExtrapolation(), and FindParticleInterpolation().

◆ FindParticleExtrapolation()

G4Ions * G4FissionProductYieldDist::FindParticleExtrapolation	(	G4double	RandomParticle,
		G4bool	LowerEnergyGroupExists )

protected

Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0, 1] random number RandomParticle by extrapolating values using the current data set. This function exists so that that different models of extrapolation may be more easily implemented in the future.

Definition at line 731 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by FindParticle().

◆ FindParticleInterpolation()

G4Ions * G4FissionProductYieldDist::FindParticleInterpolation	(	G4double	RandomParticle,
		G4int	LowerEnergyGroup )

protected

Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0, 1] random number RandomParticle by interpolating values in the current data set. This function exists so that that different models of interpolation may be more easily implemented in the future.

Definition at line 759 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by FindParticle().

◆ G4GetFission()

G4DynamicParticleVector * G4FissionProductYieldDist::G4GetFission ( )

Generates a fission event using default sampling and returns the pointer to that fission event.

Usage: No arguments required
Notes:
- The fission products are loaded into FissionContainer in the following order:
  - First daughter product
  - Second daughter product
  - Alpha particles
  - Neutrons
  - Gamma rays
- The last particle will have a NULL NextFragment pointer

Definition at line 185 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by G4FissionFragmentGenerator::G4GenerateFission().

◆ G4GetFissionProduct()

G4Ions * G4FissionProductYieldDist::G4GetFissionProduct ( )

Selects a fission fragment at random from the probability tree and returns the G4Ions pointer.

Usage: No arguments required
Notes:

Definition at line 583 of file G4FissionProductYieldDist.cc.

{
  G4FFG_FUNCTIONENTER__
 
  G4Ions* Product = FindParticle(RandomEngine_->G4SampleUniform());
 
  G4FFG_FUNCTIONLEAVE__
  return Product;
}

Referenced by G4FissionFragmentGenerator::G4GenerateFissionProduct().

◆ G4SetAlphaProduction()

void G4FissionProductYieldDist::G4SetAlphaProduction ( G4double WhatAlphaProduction )

Set the alpha production behavior for fission event generation.

Usage:
- if AlphaProduction is negative then alpha particles are sampled randomly.
Notes:
- The maximum number of alpha particles that may be created is physically limited by the nucleons present in the parent nucleus. Setting the AlphaProduction too high will have unpredictable results on the sampling of the fission products.

Definition at line 593 of file G4FissionProductYieldDist.cc.

{
  G4FFG_FUNCTIONENTER__
 
  AlphaProduction_ = WhatAlphaProduction;
 
  G4FFG_FUNCTIONLEAVE__
}

Referenced by G4FissionFragmentGenerator::G4SetAlphaProduction(), and G4FissionFragmentGenerator::InitializeFissionProductYieldClass().

◆ G4SetEnergy()

void G4FissionProductYieldDist::G4SetEnergy ( G4double WhatIncidentEnergy )

Sets the energy of the incident particle

Usage:
- WhatIncidentEnergy: Kinetic energy, if any, of the incident neutron in GeV
Notes:

Definition at line 602 of file G4FissionProductYieldDist.cc.

{
  G4FFG_FUNCTIONENTER__
 
  if (Cause_ != G4FFGEnumerations::SPONTANEOUS) {
    IncidentEnergy_ = WhatIncidentEnergy;
  }
  else {
    IncidentEnergy_ = 0 * GeV;
  }
 
  G4FFG_FUNCTIONLEAVE__
}

Referenced by G4FissionFragmentGenerator::G4SetIncidentEnergy().

◆ G4SetTernaryProbability()

void G4FissionProductYieldDist::G4SetTernaryProbability ( G4double TernaryProbability )

Sets the probability of ternary fission

Usage:
- WhatTernaryProbability: Probability of generating a ternary fission event.
Notes:

Definition at line 616 of file G4FissionProductYieldDist.cc.

{
  G4FFG_FUNCTIONENTER__
 
  TernaryProbability_ = WhatTernaryProbability;
 
  G4FFG_FUNCTIONLEAVE__
}

Referenced by G4FissionFragmentGenerator::G4SetTernaryProbability(), and G4FissionFragmentGenerator::InitializeFissionProductYieldClass().

◆ G4SetVerbosity()

void G4FissionProductYieldDist::G4SetVerbosity ( G4int WhatVerbosity )

Sets the verbosity levels

Usage:
- WhichVerbosity: Combination of levels
Notes:
- SILENT: All verbose output is repressed
- UPDATES: Only high-level internal changes are reported
- DAUGHTER_INFO: Displays information about daughter product sampling
- NEUTRON_INFO: Displays information about neutron sampling
- GAMMA_INFO: Displays information about gamma sampling
- ALPHA_INFO: Displays information about alpha sampling
- MOMENTUM_INFO: Displays information about momentum balancing
- EXTRAPOLATION_INTERPOLATION_INFO: Displays information about any data extrapolation or interpolation that occurs
- DEBUG: Reports program flow as it steps through functions
- PRINT_ALL: Displays any and all output

Definition at line 625 of file G4FissionProductYieldDist.cc.

{
  G4FFG_FUNCTIONENTER__
 
  Verbosity_ = Verbosity;
 
  ENDFData_->G4SetVerbosity(Verbosity_);
  RandomEngine_->G4SetVerbosity(Verbosity_);
 
  G4FFG_FUNCTIONLEAVE__
}

Referenced by G4FissionFragmentGenerator::G4SetVerbosity().

◆ GenerateAlphas()

void G4FissionProductYieldDist::GenerateAlphas ( std::vector< G4ReactionProduct * > * Alphas )

protectedvirtual

Generates a G4DynamicParticleVector with the fission alphas

Definition at line 844 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ GenerateNeutrons()

void G4FissionProductYieldDist::GenerateNeutrons ( std::vector< G4ReactionProduct * > * Neutrons )

protectedvirtual

Generate a linked chain of neutrons and return the pointer to the last neutron in the chain.

Definition at line 877 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ GetFissionProduct()

virtual G4Ions * G4FissionProductYieldDist::GetFissionProduct ( )

protectedpure virtual

Selects a fission product from the probability tree, limited by the number of nucleons available to the system

Implemented in G4FPYBiasedLightFragmentDist, and G4FPYNormalFragmentDist.

Referenced by G4GetFission().

◆ GetParticleDefinition()

G4Ions * G4FissionProductYieldDist::GetParticleDefinition	(	G4int	Product,
		G4FFGEnumerations::MetaState	MetaState )

protected

Returns the G4Ions definition pointer to the isotope defined by Product and MetaState. Searches the ParticleTable for the particle defined by Product (ZZZAAA) and MetaState and returns the G4Ions pointer to that particle. If the particle does not exist then it is created in G4ParticleTable and the pointer to the new particle is returned.

Definition at line 898 of file G4FissionProductYieldDist.cc.

{
  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;
}

Referenced by G4GetFission(), and SortProbability().

◆ MakeDirectoryName()

G4String G4FissionProductYieldDist::MakeDirectoryName ( )

protected

Generates the directory location for the data file referenced by G4FissionProductYieldDist

Definition at line 975 of file G4FissionProductYieldDist.cc.

{
  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();
}

◆ MakeFileName()

G4String G4FissionProductYieldDist::MakeFileName	(	G4int	Isotope,
		G4FFGEnumerations::MetaState	MetaState )

protected

Generates the appropriate file name for the isotope requested

Definition at line 988 of file G4FissionProductYieldDist.cc.

{
  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();
}

◆ MakeG4DynamicParticle()

G4DynamicParticle * G4FissionProductYieldDist::MakeG4DynamicParticle ( G4ReactionProduct * ReactionProduct )

protected

Creates a G4DynamicParticle from an existing G4ReactionProduct

Definition at line 1009 of file G4FissionProductYieldDist.cc.

{
  G4FFG_DATA_FUNCTIONENTER__
 
  auto DynamicParticle =
    new G4DynamicParticle(ReactionProduct->GetDefinition(), ReactionProduct->GetMomentum());
 
  G4FFG_DATA_FUNCTIONLEAVE__
  return DynamicParticle;
}

Referenced by G4GetFission().

◆ MakeIsotopeName()

G4String G4FissionProductYieldDist::MakeIsotopeName	(	G4int	Isotope,
		G4FFGEnumerations::MetaState	MetaState )

protected

Generates the unique name for an isotope/isomer defined by Isotope\ and MetaState in the following format: ZZZ_AAAmX_NAME

Definition at line 1020 of file G4FissionProductYieldDist.cc.

{
  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();
}

Referenced by MakeFileName().

◆ MakeTrees()

void G4FissionProductYieldDist::MakeTrees ( )

protectedvirtual

Dynamically allocates and initializes the 'field' of 'trees' with the 'trunks'

Definition at line 1050 of file G4FissionProductYieldDist.cc.

{
  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__
}

◆ ReadProbabilities()

void G4FissionProductYieldDist::ReadProbabilities ( )

protectedvirtual

Reads in the probability data from the data file

Definition at line 1084 of file G4FissionProductYieldDist.cc.

{
  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__
}

◆ Renormalize()

void G4FissionProductYieldDist::Renormalize ( ProbabilityBranch * Branch )

protected

Renormalizes the data in a ProbabilityTree. Traverses the tree structure and renormalizes all the probability data into probability segments, ensuring that no segment overlaps the other.

Definition at line 1114 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by ReadProbabilities(), and Renormalize().

◆ SampleAlphaEnergies()

void G4FissionProductYieldDist::SampleAlphaEnergies ( std::vector< G4ReactionProduct * > * Alphas )

protected

Sample the energy of the alpha particles. The energy used by the alpha particles is subtracted from the available energy

Definition at line 1138 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ SampleGammaEnergies()

void G4FissionProductYieldDist::SampleGammaEnergies ( std::vector< G4ReactionProduct * > * Gammas )

protected

Samples the energy of the gamma rays

Definition at line 1174 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ SampleNeutronEnergies()

void G4FissionProductYieldDist::SampleNeutronEnergies ( std::vector< G4ReactionProduct * > * Neutrons )

protected

Sample the energy of the neutrons using the Watt fission spectrum. The kinetic energy consumed is returned.

Definition at line 1240 of file G4FissionProductYieldDist.cc.

{
  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__
}

Referenced by G4GetFission().

◆ SetNubar()

void G4FissionProductYieldDist::SetNubar ( )

protected

Sets the nubar values for the isotope referenced by G4FissionProductYieldDistdefined from the data sets defined in SpecialOps.hh

Definition at line 1281 of file G4FissionProductYieldDist.cc.

{
  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__
}

◆ SortProbability()

void G4FissionProductYieldDist::SortProbability ( G4ENDFYieldDataContainer * YieldData )

protectedvirtual

Sorts information for a potential new particle into the correct tree

Definition at line 1326 of file G4FissionProductYieldDist.cc.

{
  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__
}