#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 (void)

G4Ions *	G4GetFissionProduct (void)

void	G4SetAlphaProduction (G4double WhatAlphaProduction)

void	G4SetEnergy (G4double WhatIncidentEnergy)

void	G4SetTernaryProbability (G4double TernaryProbability)

void	G4SetVerbosity (G4int WhatVerbosity)

virtual	~G4FissionProductYieldDist (void)

Protected Member Functions
void	CheckAlphaSanity (void)

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 (void)=0

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

G4String	MakeDirectoryName (void)

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

G4DynamicParticle *	MakeG4DynamicParticle (G4ReactionProduct *)

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

virtual void	MakeTrees (void)

virtual void	ReadProbabilities (void)

void	Renormalize (ProbabilityBranch *Branch)

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

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

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

void	SetNubar (void)

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 54 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 87 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 114 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 ( void )

virtual

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

Definition at line 1499 of file G4FissionProductYieldDist.cc.

{
G4FFG_FUNCTIONENTER__
 
    // Burn each tree, one by one
    G4int WhichTree = 0;
    while(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 1525 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)
    {
        // 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 ( void )

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 669 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 686 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;
        } else if(energyGroup == YieldEnergyGroups_ - 1)
        {
            // The energy is greater than any values in the yield data.
            lowerExists = true;
            break;
        } else
        {
            // 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 = NULL;
    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 829 of file G4FissionProductYieldDist.cc.

{
G4FFG_RECURSIVE_FUNCTIONENTER__
 
    G4Ions* Particle;
 
    // Verify that the branch exists
    if(Branch == NULL)
    {
        Particle = NULL;
    } else if(EnergyGroup1 >= Branch->IncidentEnergiesCount
              || EnergyGroup2 >= Branch->IncidentEnergiesCount
              || EnergyGroup1 == EnergyGroup2
              || Branch->IncidentEnergies[EnergyGroup1] == Branch->IncidentEnergies[EnergyGroup2])
    {
        // Set NULL if any invalid conditions exist
        Particle = NULL;
    } else
    {
        // Everything check out - proceed
        G4Ions* FoundParticle = NULL;
        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 775 of file G4FissionProductYieldDist.cc.

{
G4FFG_FUNCTIONENTER__
 
    G4Ions* FoundParticle = NULL;
    G4int NearestEnergy;
    G4int NextNearestEnergy;
 
    // Check to see if we are extrapolating above or below the data set    
    if(LowerEnergyGroupExists == true)
    {
        NearestEnergy = YieldEnergyGroups_ - 1;
        NextNearestEnergy = NearestEnergy - 1;
    } else
    {
        NearestEnergy = 0;
        NextNearestEnergy = 1;
    }
 
    for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; 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 808 of file G4FissionProductYieldDist.cc.

{
G4FFG_FUNCTIONENTER__
 
    G4Ions* FoundParticle = NULL;
    G4int HigherEnergyGroup = LowerEnergyGroup + 1;
 
    for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; Tree++)
    {
        FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk,
                                                 RandomParticle,
                                                 LowerEnergyGroup,
                                                 HigherEnergyGroup);
    }
 
G4FFG_FUNCTIONLEAVE__
    return FoundParticle;
}

Referenced by FindParticle().

◆ G4GetFission()

G4DynamicParticleVector * G4FissionProductYieldDist::G4GetFission ( void )

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 193 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 = NULL;
  G4ReactionProduct* SecondDaughter = NULL;
  std::vector< G4ReactionProduct* >* Alphas = new std::vector< G4ReactionProduct* >;
  std::vector< G4ReactionProduct* >* Neutrons = new std::vector< G4ReactionProduct* >;
  std::vector< G4ReactionProduct* >* 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) {
    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) {
    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;
        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->size() > 0)
        {
        // 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(unsigned int i = 0; i < Alphas->size(); ++i)
            {
                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 * Alphas->at(i)->GetKineticEnergy() * Alphas->at(i)->GetDefinition()->GetPDGMass());
                Alphas->at(i)->SetMomentum(Direction * Magnitude);
                ResultantVector += Alphas->at(i)->GetMomentum();
            }
        }
 
        // Sample the directions of the neutrons.
        if(Neutrons->size() != 0)
        {
            for(unsigned int i = 0; i < Neutrons->size(); ++i)
            {
                Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
                Phi = RandomEngine_->G4SampleUniform(-pi, pi);
 
                Direction.setRThetaPhi(1.0, Theta, Phi);
                Magnitude = std::sqrt(2 * Neutrons->at(i)->GetKineticEnergy() * Neutrons->at(i)->GetDefinition()->GetPDGMass());
                Neutrons->at(i)->SetMomentum(Direction * Magnitude);
                ResultantVector += Neutrons->at(i)->GetMomentum();
            }
        }
 
        // Sample the directions of the gamma rays
        if(Gammas->size() != 0)
        {
            for(unsigned int i = 0; i < Gammas->size(); ++i)
            {
                Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
                Phi = RandomEngine_->G4SampleUniform(-pi, pi);
 
                Direction.setRThetaPhi(1.0, Theta, Phi);
                Magnitude = Gammas->at(i)->GetKineticEnergy() / CLHEP::c_light;
                Gammas->at(i)->SetMomentum(Direction * Magnitude);
                ResultantVector += Gammas->at(i)->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))
        {
            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());
        G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector();
        // Load the fission fragments
        FissionProducts->push_back(MakeG4DynamicParticle(LightFragment));
        FissionProducts->push_back(MakeG4DynamicParticle(HeavyFragment));
        // Load the neutrons
        for(unsigned int i = 0; i < Neutrons->size(); i++)
        {
            FissionProducts->push_back(MakeG4DynamicParticle(Neutrons->at(i)));
        }
        // Load the gammas
        for(unsigned int i = 0; i < Gammas->size(); i++)
        {
            FissionProducts->push_back(MakeG4DynamicParticle(Gammas->at(i)));
        }
        // Load the alphas
        for(unsigned int i = 0; i < Alphas->size(); i++)
        {
            FissionProducts->push_back(MakeG4DynamicParticle(Alphas->at(i)));
        }
 
    // 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(unsigned int i = 0; i < FissionProducts->size(); i++)
        {
            Direction = FissionProducts->at(i)->GetMomentumDirection();
            Direction.rotateUz(RotationAxis);
            FissionProducts->at(i)->SetMomentumDirection(Direction);
            ResultantVector += FissionProducts->at(i)->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 ( void )

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

Usage: No arguments required
Notes:

Definition at line 609 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 620 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 630 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 646 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 656 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 899 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 937 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 ( void )

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 959 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 = reinterpret_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 ( void )

protected

Generates the directory location for the data file referenced by G4FissionProductYieldDist

Definition at line 1037 of file G4FissionProductYieldDist.cc.

{
G4FFG_FUNCTIONENTER__
 
    // Generate the file location starting in the Geant4 data directory
    std::ostringstream DirectoryName;
    DirectoryName << std::getenv("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 1051 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 1074 of file G4FissionProductYieldDist.cc.

{
G4FFG_DATA_FUNCTIONENTER__
 
    G4DynamicParticle* 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 1085 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_->theString[Z - 1];
 
G4FFG_DATA_FUNCTIONLEAVE__
    return IsotopeName.str();
}

Referenced by MakeFileName().

◆ MakeTrees()

void G4FissionProductYieldDist::MakeTrees ( void )

protectedvirtual

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

Definition at line 1118 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 = NULL;
        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 ( void )

protectedvirtual

Reads in the probability data from the data file

Definition at line 1154 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 1187 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 != NULL)
    {
        // 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 1213 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(unsigned int i = 0; i < Alphas->size(); i++)
        {
            AlphaEnergy = RandomEngine_->G4SampleGaussian(MeanAlphaEnergy,
                                                          2.35,
                                                          G4FFGEnumerations::POSITIVE) * MeV;
            // Assign the energy to the alpha particle
            Alphas->at(i)->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 1253 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 1322 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;
    G4double NeutronEnergy;
 
    // 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(unsigned int i = 0; i < Neutrons->size(); i++)
        {
            // Assign the energy to the neutron
            NeutronEnergy = RandomEngine_->G4SampleWatt(Isotope_, Cause_, IncidentEnergy_);
            Neutrons->at(i)->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 ( void )

protected

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

Definition at line 1365 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 1416 of file G4FissionProductYieldDist.cc.

{
G4FFG_DATA_FUNCTIONENTER__
 
    // Initialize the new branch
    ProbabilityBranch* NewBranch = new ProbabilityBranch;
    NewBranch->IncidentEnergiesCount = YieldEnergyGroups_;
    NewBranch->Left = NULL;
    NewBranch->Right = NULL;
    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_ == NULL)
    {
        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
    G4int 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)
        {
            // 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__
}