85 std::istringstream& dataStream )
86 : Isotope_(WhichIsotope),
87 MetaState_(WhichMetaState),
89 YieldType_(WhichYieldType),
97 Initialize(dataStream);
98 }
catch (std::exception e)
113 std::istringstream& dataStream)
114 : Isotope_(WhichIsotope),
115 MetaState_(WhichMetaState),
117 YieldType_(WhichYieldType),
118 Verbosity_(Verbosity)
125 Initialize(dataStream);
126 }
catch (std::exception e)
135 void G4FissionProductYieldDist::
136 Initialize( std::istringstream& dataStream )
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);
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)
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)
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;
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++)
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();
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;
627 IncidentEnergy_ = WhatIncidentEnergy;
630 IncidentEnergy_ = 0 *
GeV;
683 G4bool lowerExists =
false;
684 G4bool higherExists =
false;
694 if(energyGroup == 0 && IncidentEnergy_ <
YieldEnergies_[energyGroup])
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;
851 RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
854 if(RandomParticle < RangeAtIncidentEnergy)
865 RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
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++)
959 G4int A = Product % 1000;
960 G4int Z = (Product - A) / 1000;
1032 std::ostringstream DirectoryName;
1033 DirectoryName << getenv(
"G4NEUTRONHPDATA") << G4FFGDefaultValues::ENDFFissionDataLocation;
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;
1351 WhichNubar =
const_cast<G4int*
>(&SpontaneousNubar_[0][0]);
1352 NubarWidth =
const_cast<G4int*
>(&SpontaneousNubarWidth_[0][0]);
1355 WhichNubar =
const_cast<G4int*
>(&NeutronInducedNubar_[0][0]);
1356 NubarWidth =
const_cast<G4int*
>(&NeutronInducedNubarWidth_[0][0]);
1359 XFactor = std::pow(10.0, -13.0);
1360 BFactor = std::pow(10.0, -4.0);
1361 Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor
1362 + *(WhichNubar + 2) * BFactor;
1363 while(*WhichNubar != -1)
1367 Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor
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__
void set(double x, double y, double z)
virtual ~G4FissionProductYieldDist(void)
const G4FFGEnumerations::YieldType YieldType_
G4ENDFYieldDataContainer * G4GetYield(G4int WhichYield)
G4FFGEnumerations::MetaState GetMetaState(void)
G4double G4SampleGaussian(G4double Mean, G4double StdDev)
void Copy(G4int Elements, T *To, T *From)
G4double G4SampleWatt(G4int WhatIsotope, G4FFGEnumerations::FissionCause WhatCause, G4double WhatEnergy)
G4int IncidentEnergiesCount
void Divide(G4int Elements, T *To, T *Numerator, T *Denominator=NULL)
G4double * ProbabilityRangeEnd
G4double G4SampleUniform(void)
void CheckAlphaSanity(void)
G4double * GetYieldProbability(void)
G4FissionProductYieldDist(G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, std::istringstream &dataStream)
G4Ions * FindParticle(G4double RandomParticle)
G4FPYSamplingOps * RandomEngine_
G4double * YieldEnergies_
void G4SetTernaryProbability(G4double TernaryProbability)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Alpha * Definition()
ProbabilityBranch * Right
G4DynamicParticleVector * G4GetFission(void)
G4String MakeDirectoryName(void)
G4int GetPDGEncoding() const
G4ENDFTapeRead * ENDFData_
void G4SetAlphaProduction(G4double WhatAlphaProduction)
G4String MakeFileName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
G4double RemainingEnergy_
G4double * ProbabilityRangeBottom
#define G4FFG_DATA_FUNCTIONENTER__
virtual void SortProbability(G4ENDFYieldDataContainer *YieldData)
G4NeutronHPNames * ElementNames_
void G4SetVerbosity(G4int WhatVerbosity)
const G4String & GetParticleName() const
G4double * ProbabilityRangeTop
G4int GetAtomicNumber() const
G4Ions * NeutronDefinition_
G4double TernaryProbability_
G4ParticleDefinition * GetDefinition() const
virtual void GenerateNeutrons(std::vector< G4ReactionProduct * > *Neutrons)
virtual void GenerateAlphas(std::vector< G4ReactionProduct * > *Alphas)
G4DynamicParticle * MakeG4DynamicParticle(G4ReactionProduct *)
void SampleAlphaEnergies(std::vector< G4ReactionProduct * > *Alphas)
G4GLOB_DLL std::ostream G4cout
virtual void MakeTrees(void)
#define G4FFG_DATA_FUNCTIONLEAVE__
G4double AlphaProduction_
Hep3Vector & rotateUz(const Hep3Vector &)
void Multiply(G4int Elements, T *To, T *M1, T *M2=NULL)
std::vector< G4DynamicParticle * > G4DynamicParticleVector
static const G4String theString[100]
G4Ions * GetParticleDefinition(G4int Product, G4FFGEnumerations::MetaState MetaState)
G4Ions * FindParticleBranchSearch(ProbabilityBranch *Branch, G4double RandomParticle, G4int EnergyGroup1, G4int EnergyGroup2)
void G4SetEnergy(G4double WhatIncidentEnergy)
G4Ions * GammaDefinition_
void DeleteVectorOfPointers(std::vector< T > &Vector)
void setRThetaPhi(double r, double theta, double phi)
G4int GetAtomicMass() const
static G4IonTable * GetIonTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void Renormalize(ProbabilityBranch *Branch)
G4int G4GetNumberOfEnergyGroups(void)
ProbabilityBranch * Trunk
G4double GetPDGMass() const
void SampleGammaEnergies(std::vector< G4ReactionProduct * > *Gammas)
G4Ions * FindParticleExtrapolation(G4double RandomParticle, G4bool LowerEnergyGroupExists)
G4Ions * FindParticleInterpolation(G4double RandomParticle, G4int LowerEnergyGroup)
virtual void ReadProbabilities(void)
const G4FFGEnumerations::FissionCause Cause_
G4double * IncidentEnergies
static G4Neutron * Definition()
void Set(G4int Elements, T *To, T Value)
void G4SetVerbosity(G4int WhatVerbosity)
G4ThreeVector GetMomentum() const
#define G4FFG_FUNCTIONLEAVE__
G4double * MaintainNormalizedData_
void SampleNeutronEnergies(std::vector< G4ReactionProduct * > *Neutrons)
#define G4FFG_RECURSIVE_FUNCTIONENTER__
G4int G4GetNumberOfFissionProducts(void)
void BurnTree(ProbabilityBranch *Branch)
void G4SetVerbosity(G4int WhatVerbosity)
G4Ions * AlphaDefinition_
virtual G4Ions * GetFissionProduct(void)=0
G4Ions * G4GetFissionProduct(void)
#define G4FFG_FUNCTIONENTER__
G4String MakeIsotopeName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
void Add(G4int Elements, T *To, T *A1, T *A2=NULL)
G4double * G4GetEnergyGroupValues(void)
static G4Gamma * Definition()
G4int G4SampleIntegerGaussian(G4double Mean, G4double StdDev)