85 std::istringstream& dataStream )
86 : Isotope_(WhichIsotope),
87 MetaState_(WhichMetaState),
89 YieldType_(WhichYieldType),
98 }
catch (std::exception e)
113 std::istringstream& dataStream)
114 : Isotope_(WhichIsotope),
115 MetaState_(WhichMetaState),
117 YieldType_(WhichYieldType),
118 Verbosity_(Verbosity)
126 }
catch (std::exception e)
177 }
catch (std::exception& e)
202 std::vector< G4ReactionProduct* >* Alphas =
new std::vector< G4ReactionProduct* >;
203 std::vector< G4ReactionProduct* >* Neutrons =
new std::vector< G4ReactionProduct* >;
204 std::vector< G4ReactionProduct* >* Gammas =
new std::vector< G4ReactionProduct* >;
250 G4cout <<
" -- First daughter product sampled" <<
G4endl;
267 if(
Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO)
272 G4cout <<
" -- Second daughter product sampled" <<
G4endl;
336 G4double NumeratorSqrt, NumeratorOther, Denominator;
341 if(Alphas->size() > 0)
366 MassRatio = 1 / MassRatio;
374 const G4double MeanAlphaAngle = 0.3644 * MassRatio * MassRatio * MassRatio
375 - 1.9766 * MassRatio * MassRatio
382 const G4double MeanAlphaAngleStdDev = 0.0523598776;
385 for(
unsigned int i = 0; i < Alphas->size(); ++i)
388 Theta = MeanAlphaAngle + PlusMinus;
392 }
else if(Theta >
pi)
394 Theta = (2 *
pi - Theta);
398 Direction.setRThetaPhi(1.0, Theta, Phi);
399 Magnitude = std::sqrt(2 * Alphas->at(i)->GetKineticEnergy() * Alphas->at(i)->GetDefinition()->GetPDGMass());
400 Alphas->at(i)->SetMomentum(Direction * Magnitude);
401 ResultantVector += Alphas->at(i)->GetMomentum();
406 if(Neutrons->size() != 0)
408 for(
unsigned int i = 0; i < Neutrons->size(); ++i)
413 Direction.setRThetaPhi(1.0, Theta, Phi);
414 Magnitude = std::sqrt(2 * Neutrons->at(i)->GetKineticEnergy() * Neutrons->at(i)->GetDefinition()->GetPDGMass());
415 Neutrons->at(i)->SetMomentum(Direction * Magnitude);
416 ResultantVector += Neutrons->at(i)->GetMomentum();
421 if(Gammas->size() != 0)
423 for(
unsigned int i = 0; i < Gammas->size(); ++i)
428 Direction.setRThetaPhi(1.0, Theta, Phi);
429 Magnitude = Gammas->at(i)->GetKineticEnergy() / CLHEP::c_light;
430 Gammas->at(i)->SetMomentum(Direction * Magnitude);
431 ResultantVector += Gammas->at(i)->GetMomentum();
440 G4double ResultantX, ResultantY, ResultantZ;
441 ResultantX = ResultantVector.getX();
442 ResultantY = ResultantVector.getY();
443 ResultantZ = ResultantVector.getZ();
447 LightFragment = FirstDaughter;
448 HeavyFragment = SecondDaughter;
451 LightFragment = SecondDaughter;
452 HeavyFragment = FirstDaughter;
457 LightFragmentDirection.setRThetaPhi(1.0, 0, 0);
492 NumeratorSqrt = std::sqrt(LightFragmentMass *
493 (LightFragmentMass * (2 * HeavyFragmentMass * FragmentsKE
494 - ResultantX * ResultantX
495 - ResultantY * ResultantY)
496 + HeavyFragmentMass * (2 * HeavyFragmentMass * FragmentsKE
497 - ResultantX * ResultantX
498 - ResultantY * ResultantY
499 - ResultantZ - ResultantZ)));
502 NumeratorOther = LightFragmentMass * ResultantZ;
505 Denominator = LightFragmentMass + HeavyFragmentMass;
508 const G4double LightFragmentMomentum = (NumeratorSqrt + NumeratorOther) / Denominator;
509 const G4double HeavyFragmentMomentum = std::sqrt(ResultantX * ResultantX
510 + ResultantY * ResultantY
511 + std::pow(LightFragmentMomentum + ResultantZ, 2));
514 LightFragment->
SetMomentum(LightFragmentDirection * LightFragmentMomentum);
515 HeavyFragmentDirection.setX(-ResultantX);
516 HeavyFragmentDirection.setY(-ResultantY);
517 HeavyFragmentDirection.setZ((-LightFragmentMomentum - ResultantZ));
519 HeavyFragmentDirection.setR(1.0);
520 HeavyFragment->
SetMomentum(HeavyFragmentDirection * HeavyFragmentMomentum);
527 G4cout <<
" -- Daugher product momenta finalized" <<
G4endl;
539 for(
unsigned int i = 0; i < Neutrons->size(); i++)
544 for(
unsigned int i = 0; i < Gammas->size(); i++)
549 for(
unsigned int i = 0; i < Alphas->size(); i++)
561 RotationAxis.setRThetaPhi(1.0, Theta, Phi);
564 ResultantVector.set(0.0, 0.0, 0.0);
565 for(
unsigned int i = 0; i < FissionProducts->size(); i++)
567 Direction = FissionProducts->at(i)->GetMomentumDirection();
568 Direction.rotateUz(RotationAxis);
569 FissionProducts->at(i)->SetMomentumDirection(Direction);
570 ResultantVector += FissionProducts->at(i)->GetMomentum();
575 G4double PossibleImbalance = ResultantVector.mag();
576 if(PossibleImbalance > 0.01)
578 std::ostringstream Temp;
579 Temp <<
"Momenta imbalance of ";
580 Temp << PossibleImbalance / (
MeV / CLHEP::c_light);
581 Temp <<
" MeV/c in the system";
582 G4Exception(
"G4FissionProductYieldDist::G4GetFission()",
585 "Results may not be valid");
589 delete FirstDaughter;
590 delete SecondDaughter;
596 return FissionProducts;
683 G4bool lowerExists =
false;
684 G4bool higherExists =
false;
699 }
else if(energyGroup == YieldEnergyGroups_ - 1)
718 G4Ions* FoundParticle = NULL;
719 if(isExact || YieldEnergyGroups_ == 1)
726 if(RandomParticle <=
Trees_[tree].ProbabilityRangeEnd[energyGroup])
737 while((RangeIsSmaller = (RandomParticle < Branch->ProbabilityRangeBottom[energyGroup]))
742 Branch = Branch->
Left;
744 Branch = Branch->
Right;
749 }
else if(lowerExists && higherExists)
761 return FoundParticle;
766 G4bool LowerEnergyGroupExists )
770 G4Ions* FoundParticle = NULL;
772 G4int NextNearestEnergy;
775 if(LowerEnergyGroupExists ==
true)
778 NextNearestEnergy = NearestEnergy - 1;
782 NextNearestEnergy = 1;
785 for(
G4int Tree = 0; Tree <
TreeCount_ && FoundParticle == NULL; Tree++)
794 return FoundParticle;
799 G4int LowerEnergyGroup )
803 G4Ions* FoundParticle = NULL;
804 G4int HigherEnergyGroup = LowerEnergyGroup + 1;
806 for(
G4int Tree = 0; Tree <
TreeCount_ && FoundParticle == NULL; Tree++)
815 return FoundParticle;
834 || EnergyGroup1 == EnergyGroup2
842 G4Ions* FoundParticle = NULL;
854 if(RandomParticle < RangeAtIncidentEnergy)
868 if(RandomParticle > RangeAtIncidentEnergy)
881 Particle = FoundParticle;
897 G4int NumberOfAlphasToProduce;
912 for(
int i = 0; i < NumberOfAlphasToProduce; i++)
931 G4int NeutronProduction;
936 for(
int i = 0; i < NeutronProduction; i++)
960 G4int Z = (Product -
A) / 1000;
1032 std::ostringstream DirectoryName;
1037 return DirectoryName.str();
1048 std::ostringstream FileName;
1051 if(Isotope < 100000)
1060 return FileName.str();
1071 return DynamicParticle;
1081 G4int A = Isotope % 1000;
1082 G4int Z = (Isotope -
A) / 1000;
1085 std::ostringstream IsotopeName;
1087 IsotopeName << Z <<
"_" <<
A;
1104 return IsotopeName.str();
1153 for(
G4int i = 0; i < ProductCount; i++)
1216 TotalAlphaEnergy = 0;
1220 for(
unsigned int i = 0; i < Alphas->size(); i++)
1226 Alphas->at(i)->SetKineticEnergy(AlphaEnergy);
1229 TotalAlphaEnergy += AlphaEnergy;
1233 MeanAlphaEnergy -= 0.1;
1277 Gammas->back()->SetTotalEnergy(SampleEnergy);
1294 Gammas->back()->SetTotalEnergy(SampleEnergy);
1318 TotalNeutronEnergy = 0;
1323 for(
unsigned int i = 0; i < Neutrons->size(); i++)
1327 Neutrons->at(i)->SetKineticEnergy(NeutronEnergy);
1330 TotalNeutronEnergy +=NeutronEnergy;
1359 XFactor = std::pow(10.0, -13.0);
1360 BFactor = std::pow(10.0, -4.0);
1362 + *(WhichNubar + 2) * BFactor;
1363 while(*WhichNubar != -1)
1368 + *(WhichNubar + 2) * BFactor;
1375 XFactor = std::pow((
G4double)10, -6);
1377 while(*WhichNubar != -1)
1399 NewBranch->Left = NULL;
1400 NewBranch->Right = NULL;
1453 while(BranchPosition > 1)
1455 if(BranchPosition & 1)
1458 WhichBranch = &((*WhichBranch)->Right);
1462 WhichBranch = &((*WhichBranch)->Left);
1465 BranchPosition >>= 1;
1468 *WhichBranch = NewBranch;
1480 G4int WhichTree = 0;
1511 delete Branch->
Left;
1513 delete Branch->
Right;
#define G4FFG_RECURSIVE_FUNCTIONLEAVE__
virtual ~G4FissionProductYieldDist(void)
Default deconstructor.
const G4FFGEnumerations::YieldType YieldType_
The type of yield to be used: INDEPENDET or CUMULATIVE.
G4ENDFYieldDataContainer * G4GetYield(G4int WhichYield)
Returns the data for the requested fission product.
G4FFGEnumerations::MetaState GetMetaState(void)
Get the meta state.
MetaState
ENDF format provides for 3 isomers - 1 ground state and 2 meta states.
G4double G4SampleGaussian(G4double Mean, G4double StdDev)
Returns a double value taken from a Gaussian distribution about Mean and with a standard deviation of...
void Copy(G4int Elements, T *To, T *From)
Copy values from one array to another.
G4double G4SampleWatt(G4int WhatIsotope, G4FFGEnumerations::FissionCause WhatCause, G4double WhatEnergy)
Samples the Watt fission spectrum for the selected isotope, using an algorithm adopted from Ref...
G4int IncidentEnergiesCount
Number of different incident energies are available.
void Divide(G4int Elements, T *To, T *Numerator, T *Denominator=NULL)
Divide an array by another.
G4double * ProbabilityRangeEnd
The value of the highest probability segment stored in this tree.
G4double G4SampleUniform(void)
Returns a double value evenly distributed in the range (0, 1].
CLHEP::Hep3Vector G4ThreeVector
void CheckAlphaSanity(void)
Checks to make sure that alpha overpopulation will not occur, which could result in an unsolvable zer...
G4Ions * LargestA_
Defines the largest Z particle in the field of trees.
G4double * GetYieldProbability(void)
Get the yield probability.
G4FissionProductYieldDist(G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, std::istringstream &dataStream)
Default constructor.
G4Ions * FindParticle(G4double RandomParticle)
Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0...
G4Ions * Particle
Pointer to the G4Ions definition that contains the data for the isotope/isomer for this ProbabilityBr...
G4FPYSamplingOps * RandomEngine_
Pointer to the CLHEP library random engine.
G4double * YieldEnergies_
Energy values of each energy.
void G4SetTernaryProbability(G4double TernaryProbability)
Sets the probability of ternary fission.
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4FFG_LOCATION__
G4FFG_LOCATION__ outputs the current location in the code.
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Alpha * Definition()
ProbabilityBranch * Right
Pointer to the branch that has a higher probability.
G4int Verbosity_
Verbosity level.
static const G4int SpontaneousNubar_[][3]
Evaluated nubar values for neutron induced fission.
G4DynamicParticleVector * G4GetFission(void)
Generates a fission event using default sampling and returns the pointer to that fission event...
G4String MakeDirectoryName(void)
Generates the directory location for the data file referenced by G4FissionProductYieldDist.
G4int GetPDGEncoding() const
G4ENDFTapeRead * ENDFData_
Name of the fission yield product data file that G4FissionProductYieldDist references.
void G4SetAlphaProduction(G4double WhatAlphaProduction)
Set the alpha production behavior for fission event generation.
void Initialize(std::istringstream &dataStream)
Initialize is a common function called by all constructors.
G4String MakeFileName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
Generates the appropriate file name for the isotope requested.
G4double RemainingEnergy_
Container for the energy remaining to be assigned in the fission generation.
YieldType
The two types of fission data available.
G4double * ProbabilityRangeBottom
The upper limit of the probability segment, indexed by incident energy.
#define G4FFG_DATA_FUNCTIONENTER__
virtual void SortProbability(G4ENDFYieldDataContainer *YieldData)
Sorts information for a potential new particle into the correct tree.
G4NeutronHPNames * ElementNames_
Pointer to G4NeutronHPNames Provides access to the list of element names included in Geant4...
void G4SetVerbosity(G4int WhatVerbosity)
Sets the verbosity levels.
G4double NubarWidth_
Width of the gaussian distribution that samples nubar for the isotope and incident neutron energy tha...
const G4String & GetParticleName() const
G4int RemainingA_
Counter for the number of nucleons available to the fission event.
G4Gamma * GammaDefinition_
Contains the g4ParticleDefinition pointer to a gamma particle.
G4double * ProbabilityRangeTop
The lower limit of the probability segment, indexed by incident energy.
FissionCause
Causes of fission.
G4int GetAtomicNumber() const
G4Ions * NeutronDefinition_
Contains the G4ParticleDefinition pointer to a neutron, cast as a G4Ion for compatibility.
Verbosity
These are the verbosity levels.
G4double TernaryProbability_
Sets the ternary fission probability.
static const G4int SpontaneousNubarWidth_[][2]
Recommended Gaussian widths for neutron induced fission.
virtual void GenerateNeutrons(std::vector< G4ReactionProduct * > *Neutrons)
Generate a linked chain of neutrons and return the pointer to the last neutron in the chain...
virtual void GenerateAlphas(std::vector< G4ReactionProduct * > *Alphas)
Generates a G4DynamicParticleVector with the fission alphas.
G4DynamicParticle * MakeG4DynamicParticle(G4ReactionProduct *)
Creates a G4DynamicParticle from an existing G4ReactionProduct.
const G4ParticleDefinition * GetDefinition() const
void SampleAlphaEnergies(std::vector< G4ReactionProduct * > *Alphas)
Sample the energy of the alpha particles.
G4GLOB_DLL std::ostream G4cout
virtual void MakeTrees(void)
Dynamically allocates and initializes the 'field' of 'trees' with the 'trunks'.
ProbabilityTree is the 'root' for each tree hierarchy, and contains pointer to the first ProbabilityB...
G4double Nubar_
Nubar for the isotope and incident neutron energy that G4FissionProductYieldDist references.
#define G4FFG_DATA_FUNCTIONLEAVE__
G4int YieldEnergyGroups_
Number of specific energy groups.
G4double AlphaProduction_
Controls whether alpha particles are emitted, and how many.
G4double * DataTotal_
A running total of all the probabilities.
G4FFGDefaultValues is a one-stop shop for storing the default values to variables that configure how ...
G4int BranchCount_
A run-time counter for the total number of branches stored.
static const G4int NeutronInducedNubar_[][3]
Evaluated nubar values for neutron induced fission.
void Multiply(G4int Elements, T *To, T *M1, T *M2=NULL)
Multiply two arrays together.
ProbabilityTree * Trees_
An array, or 'field', of the probability trees.
G4Ions * SmallestA_
Defines the smallest A particle in the field of trees.
std::vector< G4DynamicParticle * > G4DynamicParticleVector
const G4int Isotope_
Number in ZZZAAA format of the isotope that G4FissionProductYieldDist references. ...
static const G4String theString[100]
G4Ions * GetParticleDefinition(G4int Product, G4FFGEnumerations::MetaState MetaState)
Returns the G4Ions definition pointer to the isotope defined by Product and MetaState.
G4Ions * FindParticleBranchSearch(ProbabilityBranch *Branch, G4double RandomParticle, G4int EnergyGroup1, G4int EnergyGroup2)
Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0...
void G4SetEnergy(G4double WhatIncidentEnergy)
Sets the energy of the incident particle.
static const G4double MeanGammaEnergy
Default mean gamma energy for gamma sampling.
static const char ENDFFissionDataLocation[]
ENDF data tape location, reference against G4HPNEUTRONDATA.
void DeleteVectorOfPointers(std::vector< T > &Vector)
static const G4int Isotope
Default Isotope.
static const G4double A[nN]
G4ENDFYieldDataContainer is a simple data storage class that handles the memory management internally...
G4int GetAtomicMass() const
static G4IonTable * GetIonTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void Renormalize(ProbabilityBranch *Branch)
Renormalizes the data in a ProbabilityTree.
G4int GetProduct(void)
Get the product.
G4int G4GetNumberOfEnergyGroups(void)
Returns the number of energy yield groups that were extracted from the ENDF tape file.
ProbabilityBranch * Trunk
First branch, or 'trunk', in the tree.
G4Ions * SmallestZ_
Defines the smallest Z particle in the field of trees.
G4double GetPDGMass() const
G4ENDFTapeRead is a class designed to read in data from unformatted ENDF data tapes for MT = 454 or M...
void SampleGammaEnergies(std::vector< G4ReactionProduct * > *Gammas)
Samples the energy of the gamma rays.
G4Ions * FindParticleExtrapolation(G4double RandomParticle, G4bool LowerEnergyGroupExists)
Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0...
G4bool IsEnd
Variable to identify if this is the last tree in or not.
G4int BranchCount
Counter for the number of branches that this tree has.
void SetNubar(void)
Sets the nubar values for the isotope referenced by G4FissionProductYieldDistdefined from the data se...
G4double IncidentEnergy_
Kinetic energy, if any, of the incident particle in GeV.
G4Ions * FindParticleInterpolation(G4double RandomParticle, G4int LowerEnergyGroup)
Returns the G4Ions definitions pointer for the particle whose probability segment contains the (0...
virtual void ReadProbabilities(void)
Reads in the probability data from the data file.
const G4FFGEnumerations::FissionCause Cause_
The cause of fission: SPONTANEOUS or N_INDUCED.
static const G4int NeutronInducedNubarWidth_[][2]
Recommended Gaussian widths for neutron induced fission.
G4double * IncidentEnergies
The different energies available.
static G4Neutron * Definition()
ProbabilityBranch is a tree hierarchy storage method.
void Set(G4int Elements, T *To, T Value)
Set's all the values in an array to a constant.
void G4SetVerbosity(G4int WhatVerbosity)
Sets the verbosity levels.
G4ThreeVector GetMomentum() const
#define G4FFG_FUNCTIONLEAVE__
G4FPYSamplingOps performs all the uniform and Gaussian distribution sampling operations.
G4double * MaintainNormalizedData_
Variable for ensuring that the input data is normalized.
#define G4FFG_SPACING__
G4FFG_SPACING__ indents the debug messages according to the debugging depth.
void SampleNeutronEnergies(std::vector< G4ReactionProduct * > *Neutrons)
Sample the energy of the neutrons using the Watt fission spectrum.
#define G4FFG_RECURSIVE_FUNCTIONENTER__
G4int G4GetNumberOfFissionProducts(void)
Returns the number of fission products that were extracted from the ENDF tape file.
ProbabilityBranch * Left
Pointer to the branch that has a lower probability.
void BurnTree(ProbabilityBranch *Branch)
Recursively burns each branch in a probability tree.
G4IonTable * IonTable_
Pointer to G4IonTable All G4Ions are created using G4IonTable.
void G4SetVerbosity(G4int WhatVerbosity)
Sets the verbosity levels.
G4Ions * AlphaDefinition_
Contains the G4Ions pointer to an alpha particle.
G4ThreeVector G4ParticleMomentum
G4int TreeCount_
The number of trees in the field.
virtual G4Ions * GetFissionProduct(void)=0
Selects a fission product from the probability tree, limited by the number of nucleons available to t...
G4Ions * G4GetFissionProduct(void)
Selects a fission fragment at random from the probability tree and returns the G4Ions pointer...
#define G4FFG_FUNCTIONENTER__
G4String MakeIsotopeName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
Generates the unique name for an isotope/isomer defined by Isotope\ and MetaState in the following fo...
G4Ions * LargestZ_
Defines the largest Z particle in the field of trees.
void Add(G4int Elements, T *To, T *A1, T *A2=NULL)
Add two arrays together.
G4double * G4GetEnergyGroupValues(void)
Returns and array containing the values of each of the energy groups.
static G4Gamma * Definition()
G4int RemainingZ_
Counter for the number of protons available to the fission event.
G4int G4SampleIntegerGaussian(G4double Mean, G4double StdDev)
Returns an integer value taken from a Gaussian distribution.