G4FissionProductYieldDist Class Reference

#include <G4FissionProductYieldDist.hh>

Inheritance diagram for G4FissionProductYieldDist:

Collaboration 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::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 83 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::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 ( 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__
}

Here is the call graph for this function:

Member Function Documentation

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

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 192 of file G4FissionProductYieldDist.cc.

{
G4FFG_FUNCTIONENTER__

    // 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_ % 1000;
        //const G4int ParentZ = (Isotope_ - ParentA) / 1000;
        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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.

Here is the caller graph for this function:

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);
//
//        // retreive  /control/verbose
//        G4UImanager::GetUIpointer()->SetVerboseLevel(tempVerboseLevel);
//    }

G4FFG_DATA_FUNCTIONLEAVE__
    return Temp;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 << getenv("G4NEUTRONHPDATA") << G4FFGDefaultValues::ENDFFissionDataLocation;

    // Return the directory structure
G4FFG_FUNCTIONLEAVE__
    return DirectoryName.str();
}

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

Here is the call graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

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

Here is the call graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

G4Ions* G4FissionProductYieldDist::AlphaDefinition_

protected

Contains the G4Ions pointer to an alpha particle

Definition at line 187 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::AlphaProduction_

protected

Controls whether alpha particles are emitted, and how many

Definition at line 189 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::BranchCount_

protected

A run-time counter for the total number of branches stored

Definition at line 243 of file G4FissionProductYieldDist.hh.

const G4FFGEnumerations::FissionCause G4FissionProductYieldDist::Cause_

protected

The cause of fission: SPONTANEOUS or N_INDUCED.

Definition at line 175 of file G4FissionProductYieldDist.hh.

G4double* G4FissionProductYieldDist::DataTotal_

protected

A running total of all the probabilities

Definition at line 239 of file G4FissionProductYieldDist.hh.

G4ParticleHPNames* G4FissionProductYieldDist::ElementNames_

protected

Pointer to G4NeutronHPNames
Provides access to the list of element names included in Geant4

Definition at line 255 of file G4FissionProductYieldDist.hh.

G4ENDFTapeRead* G4FissionProductYieldDist::ENDFData_

protected

Name of the fission yield product data file that G4FissionProductYieldDist references

Definition at line 183 of file G4FissionProductYieldDist.hh.

G4Gamma* G4FissionProductYieldDist::GammaDefinition_

protected

Contains the g4ParticleDefinition pointer to a gamma particle

Definition at line 193 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::IncidentEnergy_

protected

Kinetic energy, if any, of the incident particle in GeV.

Definition at line 195 of file G4FissionProductYieldDist.hh.

G4IonTable* G4FissionProductYieldDist::IonTable_

protected

Pointer to G4IonTable
All G4Ions are created using G4IonTable

Definition at line 251 of file G4FissionProductYieldDist.hh.

const G4int G4FissionProductYieldDist::Isotope_

protected

Number in ZZZAAA format of the isotope that G4FissionProductYieldDist references

Definition at line 168 of file G4FissionProductYieldDist.hh.

G4Ions* G4FissionProductYieldDist::LargestA_

protected

Defines the largest Z particle in the field of trees

Definition at line 231 of file G4FissionProductYieldDist.hh.

G4Ions* G4FissionProductYieldDist::LargestZ_

protected

Defines the largest Z particle in the field of trees.

Definition at line 229 of file G4FissionProductYieldDist.hh.

G4double* G4FissionProductYieldDist::MaintainNormalizedData_

protected

Variable for ensuring that the input data is normalized

Definition at line 237 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::MeanGammaEnergy_

protected

Sets the mean gamma energy, in MeV, produced by the fission of the isotope described by Isotope_

Definition at line 199 of file G4FissionProductYieldDist.hh.

const G4FFGEnumerations::MetaState G4FissionProductYieldDist::MetaState_

protected

MetaState information of the isotope that G4FissionProductYieldDist references
Possible values are GROUND_STATE, META_1, or META_2

Definition at line 173 of file G4FissionProductYieldDist.hh.

G4Ions* G4FissionProductYieldDist::NeutronDefinition_

protected

Contains the G4ParticleDefinition pointer to a neutron, cast as a G4Ion for compatibility

Definition at line 202 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::Nubar_

protected

Nubar for the isotope and incident neutron energy that G4FissionProductYieldDist references.

Definition at line 206 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::NubarWidth_

protected

Width of the gaussian distribution that samples nubar for the isotope and incident neutron energy that G4FissionProductYieldDist references.

Definition at line 211 of file G4FissionProductYieldDist.hh.

G4FPYSamplingOps* G4FissionProductYieldDist::RandomEngine_

protected

Pointer to the CLHEP library random engine

Definition at line 257 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::RemainingA_

protected

Counter for the number of nucleons available to the fission event

Definition at line 215 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::RemainingEnergy_

protected

Container for the energy remaining to be assigned in the fission generation

Definition at line 217 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::RemainingZ_

protected

Counter for the number of protons available to the fission event

Definition at line 213 of file G4FissionProductYieldDist.hh.

G4Ions* G4FissionProductYieldDist::SmallestA_

protected

Defines the smallest A particle in the field of trees

Definition at line 227 of file G4FissionProductYieldDist.hh.

G4Ions* G4FissionProductYieldDist::SmallestZ_

protected

Defines the smallest Z particle in the field of trees

Definition at line 225 of file G4FissionProductYieldDist.hh.

G4double G4FissionProductYieldDist::TernaryProbability_

protected

Sets the ternary fission probability. Valid ranges are [0, 1]

Definition at line 191 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::TreeCount_

protected

The number of trees in the field

Definition at line 241 of file G4FissionProductYieldDist.hh.

ProbabilityTree* G4FissionProductYieldDist::Trees_

protected

An array, or 'field', of the probability trees

Definition at line 223 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::Verbosity_

protected

Verbosity level

Definition at line 219 of file G4FissionProductYieldDist.hh.

G4double* G4FissionProductYieldDist::YieldEnergies_

protected

Energy values of each energy

Definition at line 235 of file G4FissionProductYieldDist.hh.

G4int G4FissionProductYieldDist::YieldEnergyGroups_

protected

Number of specific energy groups

Definition at line 233 of file G4FissionProductYieldDist.hh.

const G4FFGEnumerations::YieldType G4FissionProductYieldDist::YieldType_

protected

The type of yield to be used: INDEPENDET or CUMULATIVE

Definition at line 177 of file G4FissionProductYieldDist.hh.

---

The documentation for this class was generated from the following files:

• geant4.10.03.p01/source/processes/hadronic/models/particle_hp/include/G4FissionProductYieldDist.hh
• geant4.10.03.p01/source/processes/hadronic/models/particle_hp/src/G4FissionProductYieldDist.cc