87 std::istringstream& dataStream )
88 : Isotope_(WhichIsotope),
89 MetaState_(WhichMetaState),
91 YieldType_(WhichYieldType),
100 }
catch (std::exception&
e)
115 std::istringstream& dataStream)
116 : Isotope_(WhichIsotope),
117 MetaState_(WhichMetaState),
119 YieldType_(WhichYieldType),
120 Verbosity_(Verbosity)
128 }
catch (std::exception&
e)
179 }
catch (std::exception&
e)
204 std::vector< G4ReactionProduct* >* Alphas =
new std::vector< G4ReactionProduct* >;
205 std::vector< G4ReactionProduct* >* Neutrons =
new std::vector< G4ReactionProduct* >;
206 std::vector< G4ReactionProduct* >* Gammas =
new std::vector< G4ReactionProduct* >;
252 G4cout <<
" -- First daughter product sampled" <<
G4endl;
269 if(
Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO)
274 G4cout <<
" -- Second daughter product sampled" <<
G4endl;
317 G4int icounter_max=1024;
321 if ( icounter > icounter_max ) {
322 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
345 G4double NumeratorSqrt, NumeratorOther, Denominator;
350 if(Alphas->size() > 0)
375 MassRatio = 1 / MassRatio;
383 const G4double MeanAlphaAngle = 0.3644 * MassRatio * MassRatio * MassRatio
384 - 1.9766 * MassRatio * MassRatio
391 const G4double MeanAlphaAngleStdDev = 0.0523598776;
394 for(
unsigned int i = 0; i < Alphas->size(); ++i)
397 Theta = MeanAlphaAngle + PlusMinus;
401 }
else if(Theta >
pi)
403 Theta = (2 *
pi - Theta);
408 Magnitude = std::sqrt(2 * Alphas->at(i)->GetKineticEnergy() * Alphas->at(i)->GetDefinition()->GetPDGMass());
409 Alphas->at(i)->SetMomentum(Direction * Magnitude);
410 ResultantVector += Alphas->at(i)->GetMomentum();
415 if(Neutrons->size() != 0)
417 for(
unsigned int i = 0; i < Neutrons->size(); ++i)
423 Magnitude = std::sqrt(2 * Neutrons->at(i)->GetKineticEnergy() * Neutrons->at(i)->GetDefinition()->GetPDGMass());
424 Neutrons->at(i)->SetMomentum(Direction * Magnitude);
425 ResultantVector += Neutrons->at(i)->GetMomentum();
430 if(Gammas->size() != 0)
432 for(
unsigned int i = 0; i < Gammas->size(); ++i)
439 Gammas->at(i)->SetMomentum(Direction * Magnitude);
440 ResultantVector += Gammas->at(i)->GetMomentum();
449 G4double ResultantX, ResultantY, ResultantZ;
450 ResultantX = ResultantVector.
getX();
451 ResultantY = ResultantVector.
getY();
452 ResultantZ = ResultantVector.
getZ();
456 LightFragment = FirstDaughter;
457 HeavyFragment = SecondDaughter;
460 LightFragment = SecondDaughter;
461 HeavyFragment = FirstDaughter;
501 NumeratorSqrt = std::sqrt(LightFragmentMass *
502 (LightFragmentMass * (2 * HeavyFragmentMass * FragmentsKE
503 - ResultantX * ResultantX
504 - ResultantY * ResultantY)
505 + HeavyFragmentMass * (2 * HeavyFragmentMass * FragmentsKE
506 - ResultantX * ResultantX
507 - ResultantY * ResultantY
508 - ResultantZ - ResultantZ)));
511 NumeratorOther = LightFragmentMass * ResultantZ;
514 Denominator = LightFragmentMass + HeavyFragmentMass;
517 const G4double LightFragmentMomentum = (NumeratorSqrt + NumeratorOther) / Denominator;
518 const G4double HeavyFragmentMomentum = std::sqrt(ResultantX * ResultantX
519 + ResultantY * ResultantY
523 LightFragment->
SetMomentum(LightFragmentDirection * LightFragmentMomentum);
524 HeavyFragmentDirection.
setX(-ResultantX);
525 HeavyFragmentDirection.
setY(-ResultantY);
526 HeavyFragmentDirection.
setZ((-LightFragmentMomentum - ResultantZ));
528 HeavyFragmentDirection.
setR(1.0);
529 HeavyFragment->
SetMomentum(HeavyFragmentDirection * HeavyFragmentMomentum);
536 G4cout <<
" -- Daugher product momenta finalized" <<
G4endl;
548 for(
unsigned int i = 0; i < Neutrons->size(); i++)
553 for(
unsigned int i = 0; i < Gammas->size(); i++)
558 for(
unsigned int i = 0; i < Alphas->size(); i++)
573 ResultantVector.
set(0.0, 0.0, 0.0);
574 for(
unsigned int i = 0; i < FissionProducts->size(); i++)
576 Direction = FissionProducts->at(i)->GetMomentumDirection();
578 FissionProducts->at(i)->SetMomentumDirection(Direction);
579 ResultantVector += FissionProducts->at(i)->GetMomentum();
584 G4double PossibleImbalance = ResultantVector.
mag();
585 if(PossibleImbalance > 0.01)
587 std::ostringstream Temp;
588 Temp <<
"Momenta imbalance of ";
590 Temp <<
" MeV/c in the system";
591 G4Exception(
"G4FissionProductYieldDist::G4GetFission()",
594 "Results may not be valid");
598 delete FirstDaughter;
599 delete SecondDaughter;
605 return FissionProducts;
692 G4bool lowerExists =
false;
693 G4bool higherExists =
false;
708 }
else if(energyGroup == YieldEnergyGroups_ - 1)
727 G4Ions* FoundParticle = NULL;
728 if(isExact || YieldEnergyGroups_ == 1)
735 if(RandomParticle <=
Trees_[tree].ProbabilityRangeEnd[energyGroup])
746 while((RangeIsSmaller = (RandomParticle < Branch->ProbabilityRangeBottom[energyGroup]))
752 Branch = Branch->
Left;
754 Branch = Branch->
Right;
759 }
else if(lowerExists && higherExists)
771 return FoundParticle;
776 G4bool LowerEnergyGroupExists )
780 G4Ions* FoundParticle = NULL;
782 G4int NextNearestEnergy;
785 if(LowerEnergyGroupExists ==
true)
788 NextNearestEnergy = NearestEnergy - 1;
792 NextNearestEnergy = 1;
795 for(
G4int Tree = 0; Tree <
TreeCount_ && FoundParticle == NULL; Tree++)
804 return FoundParticle;
809 G4int LowerEnergyGroup )
813 G4Ions* FoundParticle = NULL;
814 G4int HigherEnergyGroup = LowerEnergyGroup + 1;
816 for(
G4int Tree = 0; Tree <
TreeCount_ && FoundParticle == NULL; Tree++)
825 return FoundParticle;
844 || EnergyGroup1 == EnergyGroup2
852 G4Ions* FoundParticle = NULL;
864 if(RandomParticle < RangeAtIncidentEnergy)
878 if(RandomParticle > RangeAtIncidentEnergy)
891 Particle = FoundParticle;
907 G4int NumberOfAlphasToProduce;
922 for(
int i = 0; i < NumberOfAlphasToProduce; i++)
941 G4int NeutronProduction;
946 for(
int i = 0; i < NeutronProduction; i++)
970 G4int Z = (Product -
A) / 1000;
1042 std::ostringstream DirectoryName;
1047 return DirectoryName.str();
1058 std::ostringstream FileName;
1061 if(Isotope < 100000)
1070 return FileName.str();
1081 return DynamicParticle;
1091 G4int A = Isotope % 1000;
1092 G4int Z = (Isotope -
A) / 1000;
1095 std::ostringstream IsotopeName;
1097 IsotopeName << Z <<
"_" <<
A;
1114 return IsotopeName.str();
1163 for(
G4int i = 0; i < ProductCount; i++)
1226 TotalAlphaEnergy = 0;
1230 for(
unsigned int i = 0; i < Alphas->size(); i++)
1236 Alphas->at(i)->SetKineticEnergy(AlphaEnergy);
1239 TotalAlphaEnergy += AlphaEnergy;
1243 MeanAlphaEnergy -= 0.1;
1268 G4int icounter_max=1024;
1272 if ( icounter > icounter_max ) {
1273 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
1294 Gammas->back()->SetTotalEnergy(SampleEnergy);
1311 Gammas->back()->SetTotalEnergy(SampleEnergy);
1334 G4int icounter_max=1024;
1338 if ( icounter > icounter_max ) {
1339 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
1342 TotalNeutronEnergy = 0;
1347 for(
unsigned int i = 0; i < Neutrons->size(); i++)
1351 Neutrons->at(i)->SetKineticEnergy(NeutronEnergy);
1354 TotalNeutronEnergy +=NeutronEnergy;
1386 + *(WhichNubar + 2) * BFactor;
1387 while(*WhichNubar != -1)
1392 + *(WhichNubar + 2) * BFactor;
1401 while(*WhichNubar != -1)
1423 NewBranch->Left = NULL;
1424 NewBranch->Right = NULL;
1477 while(BranchPosition > 1)
1479 if(BranchPosition & 1)
1482 WhichBranch = &((*WhichBranch)->Right);
1486 WhichBranch = &((*WhichBranch)->Left);
1489 BranchPosition >>= 1;
1492 *WhichBranch = NewBranch;
1504 G4int WhichTree = 0;
1535 delete Branch->
Left;
1537 delete Branch->
Right;
#define G4FFG_RECURSIVE_FUNCTIONLEAVE__
virtual void SortProbability(G4ENDFYieldDataContainer *YieldData)
void set(double x, double y, double z)
static const G4int NeutronInducedNubarWidth_[][2]
void Copy(G4int Elements, T *To, T *From)
static G4Gamma * Definition()
void CheckAlphaSanity(void)
G4double RemainingEnergy_
G4int GetAtomicNumber() const
#define G4FFG_RECURSIVE_FUNCTIONENTER__
G4double TernaryProbability_
void Initialize(std::istringstream &dataStream)
G4Ions * GetParticleDefinition(G4int Product, G4FFGEnumerations::MetaState MetaState)
static const G4int NeutronInducedNubar_[][3]
static G4Alpha * Definition()
static constexpr double MeV
G4ENDFYieldDataContainer * G4GetYield(G4int WhichYield)
static G4Neutron * Definition()
void SampleAlphaEnergies(std::vector< G4ReactionProduct * > *Alphas)
ProbabilityBranch * Right
static G4IonTable * GetIonTable()
static const G4int SpontaneousNubarWidth_[][2]
G4double * ProbabilityRangeTop
static constexpr double keV
G4Gamma * GammaDefinition_
#define G4FFG_FUNCTIONLEAVE__
#define G4FFG_DATA_FUNCTIONENTER__
void SetMomentum(const G4double x, const G4double y, const G4double z)
virtual void GenerateAlphas(std::vector< G4ReactionProduct * > *Alphas)
G4Ions * FindParticleBranchSearch(ProbabilityBranch *Branch, G4double RandomParticle, G4int EnergyGroup1, G4int EnergyGroup2)
const G4String & GetParticleName() const
std::vector< G4DynamicParticle * > G4DynamicParticleVector
void SampleGammaEnergies(std::vector< G4ReactionProduct * > *Gammas)
G4String MakeIsotopeName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
static const G4String theString[100]
void Add(G4int Elements, T *To, T *A1, T *A2=NULL)
Hep3Vector & rotateUz(const Hep3Vector &)
const G4ParticleDefinition * GetDefinition() const
G4double * IncidentEnergies
G4Ions * FindParticleInterpolation(G4double RandomParticle, G4int LowerEnergyGroup)
G4double GetPDGMass() const
virtual G4Ions * GetFissionProduct(void)=0
static const char ENDFFissionDataLocation[]
G4double powN(G4double x, G4int n) const
ProbabilityBranch * Trunk
G4int IncidentEnergiesCount
G4ENDFTapeRead * ENDFData_
G4Ions * FindParticleExtrapolation(G4double RandomParticle, G4bool LowerEnergyGroupExists)
virtual void MakeTrees(void)
static G4Pow * GetInstance()
const G4FFGEnumerations::YieldType YieldType_
G4double G4SampleWatt(G4int WhatIsotope, G4FFGEnumerations::FissionCause WhatCause, G4double WhatEnergy)
G4Ions * FindParticle(G4double RandomParticle)
G4double powA(G4double A, G4double y) const
G4String MakeDirectoryName(void)
void G4SetTernaryProbability(G4double TernaryProbability)
virtual void GenerateNeutrons(std::vector< G4ReactionProduct * > *Neutrons)
void G4SetVerbosity(G4int WhatVerbosity)
double A(double temperature)
G4double * YieldEnergies_
G4FPYSamplingOps * RandomEngine_
G4double * MaintainNormalizedData_
void SampleNeutronEnergies(std::vector< G4ReactionProduct * > *Neutrons)
G4Ions * AlphaDefinition_
void Divide(G4int Elements, T *To, T *Numerator, T *Denominator=NULL)
G4ThreeVector GetMomentum() const
virtual ~G4FissionProductYieldDist(void)
G4int G4GetNumberOfEnergyGroups(void)
const G4FFGEnumerations::FissionCause Cause_
void G4SetAlphaProduction(G4double WhatAlphaProduction)
G4double * GetYieldProbability(void)
G4int G4GetNumberOfFissionProducts(void)
void G4SetEnergy(G4double WhatIncidentEnergy)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void BurnTree(ProbabilityBranch *Branch)
static constexpr double c_light
void DeleteVectorOfPointers(std::vector< T > &Vector)
G4String MakeFileName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
G4double * G4GetEnergyGroupValues(void)
G4double AlphaProduction_
void setRThetaPhi(double r, double theta, double phi)
void Multiply(G4int Elements, T *To, T *M1, T *M2=NULL)
#define G4FFG_DATA_FUNCTIONLEAVE__
G4int GetAtomicMass() const
G4ParticleHPNames * ElementNames_
void Renormalize(ProbabilityBranch *Branch)
G4int GetPDGEncoding() const
G4GLOB_DLL std::ostream G4cout
G4Ions * NeutronDefinition_
G4DynamicParticle * MakeG4DynamicParticle(G4ReactionProduct *)
G4int G4SampleIntegerGaussian(G4double Mean, G4double StdDev)
void G4SetVerbosity(G4int WhatVerbosity)
static constexpr double pi
G4FFGEnumerations::MetaState GetMetaState(void)
static const G4double MeanGammaEnergy
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4double * ProbabilityRangeEnd
G4Ions * G4GetFissionProduct(void)
G4double G4SampleUniform(void)
void Set(G4int Elements, T *To, T Value)
G4DynamicParticleVector * G4GetFission(void)
G4double * ProbabilityRangeBottom
static constexpr double GeV
static const G4int SpontaneousNubar_[][3]
G4FissionProductYieldDist(G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, std::istringstream &dataStream)
static const G4int Isotope
#define G4FFG_FUNCTIONENTER__
static constexpr double pi
G4double G4SampleGaussian(G4double Mean, G4double StdDev)
void G4SetVerbosity(G4int WhatVerbosity)
virtual void ReadProbabilities(void)