Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FissionProductYieldDist.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /*
27  * File: G4FissionProductYieldDist.cc
28  * Author: B. Wendt (wendbryc@isu.edu)
29  *
30  * Created on May 11, 2011, 12:04 PM
31  */
32 
33 /* * * * * * * * * * * * * * * * References * * * * * * * * * * * * * * * *
34  * *
35  * 1. "Systematics of fission fragment total kinetic energy release", *
36  * V. E. Viola, K. Kwiatkowski, and M. Walker, Physical Review C, 31.4, *
37  * April 1985 *
38  * 2. "Reactor Handbook", United States Atomic Energy Commission, *
39  * III.A:Physics, 1962 *
40  * 3. "Properties of the Alpha Particles Emitted in the Spontaneous Fission *
41  * of Cf252", Z. Fraenkel and S. G. Thompson, Physical Review Letters, *
42  * 13.14, October 1964 *
43  * *
44  * * * * * * * * * * * * * * * * References * * * * * * * * * * * * * * * */
45 
46 
47 //#include <ios>
48 //#include <iostream>
49 
50 #include "G4Alpha.hh"
51 #include "G4Gamma.hh"
52 #include "G4Ions.hh"
53 #include "G4Neutron.hh"
54 //#include "G4NeutronHPNames.hh"
55 #include "G4NucleiProperties.hh"
56 #include "G4ParticleTable.hh"
57 #include "G4ThreeVector.hh"
58 #include "Randomize.hh"
59 #include "G4UImanager.hh"
60 #include "globals.hh"
61 #include "G4Exp.hh"
62 #include "G4SystemOfUnits.hh"
63 #include "G4HadFinalState.hh"
64 #include "G4DynamicParticle.hh"
66 #include "G4ReactionProduct.hh"
67 
68 #include "G4ArrayOps.hh"
69 #include "G4ENDFTapeRead.hh"
71 #include "G4FFGDebuggingMacros.hh"
72 #include "G4FFGDefaultValues.hh"
73 #include "G4FFGEnumerations.hh"
74 #include "G4FPYNubarValues.hh"
75 #include "G4FPYSamplingOps.hh"
76 #include "G4FPYTreeStructures.hh"
78 #include "G4Pow.hh"
79 
80 using CLHEP::pi;
81 
84  G4FFGEnumerations::MetaState WhichMetaState,
86  G4FFGEnumerations::YieldType WhichYieldType,
87  std::istringstream& dataStream )
88 : Isotope_(WhichIsotope),
89  MetaState_(WhichMetaState),
90  Cause_(WhichCause),
91  YieldType_(WhichYieldType),
92  Verbosity_(G4FFGDefaultValues::Verbosity)
93 {
95 
96  try
97  {
98  // Initialize the class
99  Initialize(dataStream);
100  } catch (std::exception e)
101  {
103  throw e;
104  }
105 
107 }
108 
111  G4FFGEnumerations::MetaState WhichMetaState,
113  G4FFGEnumerations::YieldType WhichYieldType,
115  std::istringstream& dataStream)
116 : Isotope_(WhichIsotope),
117  MetaState_(WhichMetaState),
118  Cause_(WhichCause),
119  YieldType_(WhichYieldType),
120  Verbosity_(Verbosity)
121 {
123 
124  try
125  {
126  // Initialize the class
127  Initialize(dataStream);
128  } catch (std::exception e)
129  {
131  throw e;
132  }
133 
135 }
136 
137 void G4FissionProductYieldDist::
138 Initialize( std::istringstream& dataStream )
139 {
141 
142  IncidentEnergy_ = 0.0;
144  AlphaProduction_ = 0;
145  SetNubar();
146 
147  // Set miscellaneous variables
148  AlphaDefinition_ = reinterpret_cast<G4Ions*>(G4Alpha::Definition());
149  NeutronDefinition_ = reinterpret_cast<G4Ions*>(G4Neutron::Definition());
152 
153  // Construct G4NeutronHPNames: provides access to the element names
155  // Get the pointer to G4ParticleTable: stores all G4Ions
157  // Construct the pointer to the random engine
158  // TODO Make G4FPSamplingOps a singleton so that only one instance is used across all classes
160 
161  try
162  {
163  // Read in and sort the probability data
164  ENDFData_ = new G4ENDFTapeRead(dataStream,
165  YieldType_,
166  Cause_,
167  Verbosity_);
168 // ENDFData_ = new G4ENDFTapeRead(MakeDirectoryName(),
169 // MakeFileName(Isotope_, MetaState_),
170 // YieldType_,
171 // Cause_);
177  MakeTrees();
179  } catch (std::exception& e)
180  {
181  delete ElementNames_;
182  delete RandomEngine_;
183 
185  throw e;
186  }
187 
189 }
190 
193 {
195 
196  // Check to see if the user has set the alpha production to a somewhat
197  // reasonable level
199 
200  // Generate the new G4DynamicParticle pointers to identify key locations in
201  // the G4DynamicParticle chain that will be passed to the G4FissionEvent
202  G4ReactionProduct* FirstDaughter = NULL;
203  G4ReactionProduct* SecondDaughter = NULL;
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* >;
207 
208  // Generate all the nucleonic fission products
209  // How many nucleons do we have to work with?
210  //TK modified 131108
211  //const G4int ParentA = Isotope_ % 1000;
212  //const G4int ParentZ = (Isotope_ - ParentA) / 1000;
213  const G4int ParentA = (Isotope_/10) % 1000;
214  const G4int ParentZ = ((Isotope_/10) - ParentA) / 1000;
215  RemainingA_ = ParentA;
216  RemainingZ_ = ParentZ;
217 
218  // Don't forget the extra nucleons depending on the fission cause
219  switch(Cause_)
220  {
222  ++RemainingA_;
223  break;
224 
226  ++RemainingZ_;
227  break;
228 
231  default:
232  // Nothing to do here
233  break;
234  }
235 
236  // Ternary fission can be set by the user. Thus, it is necessary to
237  // sample the alpha particle first and the first daughter product
238  // second. See the discussion in
239  // G4FissionProductYieldDist::G4GetFissionProduct() for more information
240  // as to why the fission events are sampled this way.
241  GenerateAlphas(Alphas);
242 
243  // Generate the first daughter product
244  FirstDaughter = new G4ReactionProduct(GetFissionProduct());
245  RemainingA_ -= FirstDaughter->GetDefinition()->GetAtomicMass();
246  RemainingZ_ -= FirstDaughter->GetDefinition()->GetAtomicNumber();
248  {
251 
252  G4cout << " -- First daughter product sampled" << G4endl;
254  G4cout << " Name: " << FirstDaughter->GetDefinition()->GetParticleName() << G4endl;
256  G4cout << " Z: " << FirstDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
258  G4cout << " A: " << FirstDaughter->GetDefinition()->GetAtomicMass() << G4endl;
260  G4cout << " Meta State: " << (FirstDaughter->GetDefinition()->GetPDGEncoding() % 10) << G4endl;
261  }
262 
263  GenerateNeutrons(Neutrons);
264 
265  // Now that all the nucleonic particles have been generated, we can
266  // calculate the composition of the second daughter product.
267  G4int NewIsotope = RemainingZ_ * 1000 + RemainingA_;
269  if(Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO)
270  {
273 
274  G4cout << " -- Second daughter product sampled" << G4endl;
276  G4cout << " Name: " << SecondDaughter->GetDefinition()->GetParticleName() << G4endl;
278  G4cout << " Z: " << SecondDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
280  G4cout << " A: " << SecondDaughter->GetDefinition()->GetAtomicMass() << G4endl;
282  G4cout << " Meta State: " << (SecondDaughter->GetDefinition()->GetPDGEncoding() % 10) << G4endl;
283  }
284 
285  // Calculate how much kinetic energy will be available
286  // 195 to 205 MeV are available in a fission reaction, but about 20 MeV
287  // are from delayed sources. We are concerned only with prompt sources,
288  // so sample a Gaussian distribution about 20 MeV and subtract the
289  // result from the total available energy. Also, the energy of fission
290  // neutrinos is neglected. Fission neutrinos would add ~11 MeV
291  // additional energy to the fission. (Ref 2)
292  // Finally, add in the kinetic energy contribution of the fission
293  // inducing particle, if any.
294  const G4double TotalKE = RandomEngine_->G4SampleUniform(195.0, 205.0) * MeV
295  - RandomEngine_->G4SampleGaussian(20.0, 3.0) * MeV
296  + IncidentEnergy_;
297  RemainingEnergy_ = TotalKE;
298 
299  // Calculate the energies of the alpha particles and neutrons
300  // Once again, since the alpha particles are user defined, we must
301  // sample their kinetic energy first. SampleAlphaEnergies() returns the
302  // amount of energy consumed by the alpha particles, so remove the total
303  // energy alloted to the alpha particles from the available energy
304  SampleAlphaEnergies(Alphas);
305 
306  // Second, the neutrons are sampled from the Watt fission spectrum.
307  SampleNeutronEnergies(Neutrons);
308 
309  // Calculate the total energy available to the daughter products
310  // A Gaussian distribution about the average calculated energy with
311  // a standard deviation of 1.5 MeV (Ref. 2) is used. Since the energy
312  // distribution is dependant on the alpha particle generation and the
313  // Watt fission sampling for neutrons, we only have the left-over energy
314  // to work with for the fission daughter products.
315  G4double FragmentsKE;
316  G4int icounter=0;
317  G4int icounter_max=1024;
318  do
319  {
320  icounter++;
321  if ( icounter > icounter_max ) {
322  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
323  break;
324  }
325  FragmentsKE = RandomEngine_->G4SampleGaussian(RemainingEnergy_, 1.5 *MeV);
326  } while(FragmentsKE > RemainingEnergy_); // Loop checking, 11.05.2015, T. Koi
327 
328  // Make sure that we don't produce any sub-gamma photons
329  if((RemainingEnergy_ - FragmentsKE) / (100 * keV) < 1.0)
330  {
331  FragmentsKE = RemainingEnergy_;
332  }
333 
334  // This energy has now been allotted to the fission fragments.
335  // Subtract FragmentsKE from the total available fission energy.
336  RemainingEnergy_ -= FragmentsKE;
337 
338  // Sample the energies of the gamma rays
339  // Assign the remainder, if any, of the energy to the gamma rays
340  SampleGammaEnergies(Gammas);
341 
342  // Prepare to balance the momenta of the system
343  // We will need these for sampling the angles and balancing the momenta
344  // of all the fission products
345  G4double NumeratorSqrt, NumeratorOther, Denominator;
346  G4double Theta, Phi, Magnitude;
347  G4ThreeVector Direction;
348  G4ParticleMomentum ResultantVector(0, 0, 0);
349 
350  if(Alphas->size() > 0)
351  {
352  // Sample the angles of the alpha particles and neutrons, then calculate
353  // the total moment contribution to the system
354  // The average angle of the alpha particles with respect to the
355  // light fragment is dependent on the ratio of the kinetic energies.
356  // This equation was determined by performing a fit on data from
357  // Ref. 3 using the website:
358  // http://soft.arquimedex.com/linear_regression.php
359  //
360  // RatioOfKE Angle (rad) Angle (degrees)
361  // 1.05 1.257 72
362  // 1.155 1.361 78
363  // 1.28 1.414 81
364  // 1.5 1.518 87
365  // 1.75 1.606 92
366  // 1.9 1.623 93
367  // 2.2 1.728 99
368  // This equation generates the angle in radians. If the RatioOfKE is
369  // greater than 2.25 the angle defaults to 1.3963 rad (100 degrees)
370  G4double MassRatio = FirstDaughter->GetDefinition()->GetPDGMass() / SecondDaughter->GetDefinition()->GetPDGMass();
371 
372  // Invert the mass ratio if the first daughter product is the lighter fragment
373  if(MassRatio < 1)
374  {
375  MassRatio = 1 / MassRatio;
376  }
377 
378  // The empirical equation is valid for mass ratios up to 2.75
379  if(MassRatio > 2.75)
380  {
381  MassRatio = 2.75;
382  }
383  const G4double MeanAlphaAngle = 0.3644 * MassRatio * MassRatio * MassRatio
384  - 1.9766 * MassRatio * MassRatio
385  + 3.8207 * MassRatio
386  - 0.9917;
387 
388  // Sample the directions of the alpha particles with respect to the
389  // light fragment. For the moment we will assume that the light
390  // fragment is traveling along the z-axis in the positive direction.
391  const G4double MeanAlphaAngleStdDev = 0.0523598776;
392  G4double PlusMinus;
393 
394  for(unsigned int i = 0; i < Alphas->size(); ++i)
395  {
396  PlusMinus = std::acos(RandomEngine_->G4SampleGaussian(0, MeanAlphaAngleStdDev)) - (pi / 2);
397  Theta = MeanAlphaAngle + PlusMinus;
398  if(Theta < 0)
399  {
400  Theta = 0.0 - Theta;
401  } else if(Theta > pi)
402  {
403  Theta = (2 * pi - Theta);
404  }
406 
407  Direction.setRThetaPhi(1.0, Theta, Phi);
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();
411  }
412  }
413 
414  // Sample the directions of the neutrons.
415  if(Neutrons->size() != 0)
416  {
417  for(unsigned int i = 0; i < Neutrons->size(); ++i)
418  {
419  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
421 
422  Direction.setRThetaPhi(1.0, Theta, Phi);
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();
426  }
427  }
428 
429  // Sample the directions of the gamma rays
430  if(Gammas->size() != 0)
431  {
432  for(unsigned int i = 0; i < Gammas->size(); ++i)
433  {
434  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
436 
437  Direction.setRThetaPhi(1.0, Theta, Phi);
438  Magnitude = Gammas->at(i)->GetKineticEnergy() / CLHEP::c_light;
439  Gammas->at(i)->SetMomentum(Direction * Magnitude);
440  ResultantVector += Gammas->at(i)->GetMomentum();
441  }
442  }
443 
444  // Calculate the momenta of the two daughter products
445  G4ReactionProduct* LightFragment;
446  G4ReactionProduct* HeavyFragment;
447  G4ThreeVector LightFragmentDirection;
448  G4ThreeVector HeavyFragmentDirection;
449  G4double ResultantX, ResultantY, ResultantZ;
450  ResultantX = ResultantVector.getX();
451  ResultantY = ResultantVector.getY();
452  ResultantZ = ResultantVector.getZ();
453 
454  if(FirstDaughter->GetDefinition()->GetPDGMass() < SecondDaughter->GetDefinition()->GetPDGMass())
455  {
456  LightFragment = FirstDaughter;
457  HeavyFragment = SecondDaughter;
458  } else
459  {
460  LightFragment = SecondDaughter;
461  HeavyFragment = FirstDaughter;
462  }
463  const G4double LightFragmentMass = LightFragment->GetDefinition()->GetPDGMass();
464  const G4double HeavyFragmentMass = HeavyFragment->GetDefinition()->GetPDGMass();
465 
466  LightFragmentDirection.setRThetaPhi(1.0, 0, 0);
467 
468  // Fit the momenta of the daughter products to the resultant vector of
469  // the remaining fission products. This will be done in the Cartesian
470  // coordinate system, not spherical. This is done using the following
471  // table listing the system momenta and the corresponding equations:
472  // X Y Z
473  //
474  // A 0 0 P
475  //
476  // B -R_x -R_y -P - R_z
477  //
478  // R R_x R_y R_z
479  //
480  // v = sqrt(2*m*k) -> k = v^2/(2*m)
481  // tk = k_A + k_B
482  // k_L = P^2/(2*m_L)
483  // k_H = ((-R_x)^2 + (-R_y)^2 + (-P - R_z)^2)/(2*m_H)
484  // where:
485  // P: momentum of the light daughter product
486  // R: the remaining fission products' resultant vector
487  // v: momentum
488  // m: mass
489  // k: kinetic energy
490  // tk: total kinetic energy available to the daughter products
491  //
492  // Below is the solved form for P, with the solution generated using
493  // the WolframAlpha website:
494  // http://www.wolframalpha.com/input/?i=
495  // solve+((-x)^2+%2B+(-y)^2+%2B+(-P-z)^2)%2F(2*B)+%2B+L^2%2F(2*A)+%3D+k+
496  // for+P
497  //
498  //
499  // nsqrt = sqrt(m_L*(m_L*(2*m_H*tk - R_x^2 - R_y^2) +
500  // m_H*(2*m_H*tk - R_x^2 - R_y^2 - R_z^2))
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)));
509 
510  // nother = m_L*R_z
511  NumeratorOther = LightFragmentMass * ResultantZ;
512 
513  // denom = m_L + m_H
514  Denominator = LightFragmentMass + HeavyFragmentMass;
515 
516  // P = (nsqrt + nother) / denom
517  const G4double LightFragmentMomentum = (NumeratorSqrt + NumeratorOther) / Denominator;
518  const G4double HeavyFragmentMomentum = std::sqrt(ResultantX * ResultantX
519  + ResultantY * ResultantY
520  + G4Pow::GetInstance()->powN(LightFragmentMomentum + ResultantZ, 2));
521 
522  // Finally! We now have everything we need for the daughter products
523  LightFragment->SetMomentum(LightFragmentDirection * LightFragmentMomentum);
524  HeavyFragmentDirection.setX(-ResultantX);
525  HeavyFragmentDirection.setY(-ResultantY);
526  HeavyFragmentDirection.setZ((-LightFragmentMomentum - ResultantZ));
527  // Don't forget to normalize the vector to the unit sphere
528  HeavyFragmentDirection.setR(1.0);
529  HeavyFragment->SetMomentum(HeavyFragmentDirection * HeavyFragmentMomentum);
530 
531  if(Verbosity_ & (G4FFGEnumerations::DAUGHTER_INFO | G4FFGEnumerations::MOMENTUM_INFO))
532  {
535 
536  G4cout << " -- Daugher product momenta finalized" << G4endl;
538  }
539 
540  // Load all the particles into a contiguous set
541  //TK modifed 131108
542  //G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector(2 + Alphas->size() + Neutrons->size() + Gammas->size());
543  G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector();
544  // Load the fission fragments
545  FissionProducts->push_back(MakeG4DynamicParticle(LightFragment));
546  FissionProducts->push_back(MakeG4DynamicParticle(HeavyFragment));
547  // Load the neutrons
548  for(unsigned int i = 0; i < Neutrons->size(); i++)
549  {
550  FissionProducts->push_back(MakeG4DynamicParticle(Neutrons->at(i)));
551  }
552  // Load the gammas
553  for(unsigned int i = 0; i < Gammas->size(); i++)
554  {
555  FissionProducts->push_back(MakeG4DynamicParticle(Gammas->at(i)));
556  }
557  // Load the alphas
558  for(unsigned int i = 0; i < Alphas->size(); i++)
559  {
560  FissionProducts->push_back(MakeG4DynamicParticle(Alphas->at(i)));
561  }
562 
563  // Rotate the system to a random location so that the light fission fragment
564  // is not always traveling along the positive z-axis
565  // Sample Theta and Phi.
566  G4ThreeVector RotationAxis;
567 
568  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
570  RotationAxis.setRThetaPhi(1.0, Theta, Phi);
571 
572  // We will also check the net momenta
573  ResultantVector.set(0.0, 0.0, 0.0);
574  for(unsigned int i = 0; i < FissionProducts->size(); i++)
575  {
576  Direction = FissionProducts->at(i)->GetMomentumDirection();
577  Direction.rotateUz(RotationAxis);
578  FissionProducts->at(i)->SetMomentumDirection(Direction);
579  ResultantVector += FissionProducts->at(i)->GetMomentum();
580  }
581 
582  // Warn if the sum momenta of the system is not within a reasonable
583  // tolerance
584  G4double PossibleImbalance = ResultantVector.mag();
585  if(PossibleImbalance > 0.01)
586  {
587  std::ostringstream Temp;
588  Temp << "Momenta imbalance of ";
589  Temp << PossibleImbalance / (MeV / CLHEP::c_light);
590  Temp << " MeV/c in the system";
591  G4Exception("G4FissionProductYieldDist::G4GetFission()",
592  Temp.str().c_str(),
593  JustWarning,
594  "Results may not be valid");
595  }
596 
597  // Clean up
598  delete FirstDaughter;
599  delete SecondDaughter;
603 
605  return FissionProducts;
606 }
607 
610 {
612 
614 
616  return Product;
617 }
618 
620 G4SetAlphaProduction(G4double WhatAlphaProduction)
621 {
623 
624  AlphaProduction_ = WhatAlphaProduction;
625 
627 }
628 
630 G4SetEnergy( G4double WhatIncidentEnergy )
631 {
633 
635  {
636  IncidentEnergy_ = WhatIncidentEnergy;
637  } else
638  {
639  IncidentEnergy_ = 0 * GeV;
640  }
641 
643 }
644 
646 G4SetTernaryProbability( G4double WhatTernaryProbability )
647 {
649 
650  TernaryProbability_ = WhatTernaryProbability;
651 
653 }
654 
657 {
659 
661 
662  ENDFData_->G4SetVerbosity(Verbosity_);
663  RandomEngine_->G4SetVerbosity(Verbosity_);
664 
666 }
667 
670 {
672 
673  // This provides comfortable breathing room at 16 MeV per alpha
674  if(AlphaProduction_ > 10)
675  {
676  AlphaProduction_ = 10;
677  } else if(AlphaProduction_ < -7)
678  {
679  AlphaProduction_ = -7;
680  }
681 
683 }
684 
686 FindParticle( G4double RandomParticle )
687 {
689 
690  // Determine which energy group is currently in use
691  G4bool isExact = false;
692  G4bool lowerExists = false;
693  G4bool higherExists = false;
694  G4int energyGroup;
695  for(energyGroup = 0; energyGroup < YieldEnergyGroups_; energyGroup++)
696  {
697  if(IncidentEnergy_ == YieldEnergies_[energyGroup])
698  {
699  isExact = true;
700  break;
701  }
702 
703  if(energyGroup == 0 && IncidentEnergy_ < YieldEnergies_[energyGroup])
704  {
705  // Break if the energy is less than the lowest energy
706  higherExists = true;
707  break;
708  } else if(energyGroup == YieldEnergyGroups_ - 1)
709  {
710  // The energy is greater than any values in the yield data.
711  lowerExists = true;
712  break;
713  } else
714  {
715  // Break if the energy is less than the lowest energy
716  if(IncidentEnergy_ > YieldEnergies_[energyGroup])
717  {
718  energyGroup--;
719  lowerExists = true;
720  higherExists = true;
721  break;
722  }
723  }
724  }
725 
726  // Determine which particle it is
727  G4Ions* FoundParticle = NULL;
728  if(isExact || YieldEnergyGroups_ == 1)
729  {
730  // Determine which tree contains the random value
731  G4int tree;
732  for(tree = 0; tree < TreeCount_; tree++)
733  {
734  // Break if a tree is identified as containing the random particle
735  if(RandomParticle <= Trees_[tree].ProbabilityRangeEnd[energyGroup])
736  {
737  break;
738  }
739  }
740  ProbabilityBranch* Branch = Trees_[tree].Trunk;
741 
742  // Iteratively traverse the tree until the particle addressed by the random
743  // variable is found
744  G4bool RangeIsSmaller;
745  G4bool RangeIsGreater;
746  while((RangeIsSmaller = (RandomParticle < Branch->ProbabilityRangeBottom[energyGroup]))
747  || (RangeIsGreater = (RandomParticle > Branch->ProbabilityRangeTop[energyGroup])))
748  // Loop checking, 11.05.2015, T. Koi
749  {
750  if(RangeIsSmaller)
751  {
752  Branch = Branch->Left;
753  } else {
754  Branch = Branch->Right;
755  }
756  }
757 
758  FoundParticle = Branch->Particle;
759  } else if(lowerExists && higherExists)
760  {
761  // We need to do some interpolation
762  FoundParticle = FindParticleInterpolation(RandomParticle, energyGroup);
763  } else
764  {
765  // We need to do some extrapolation
766  FoundParticle = FindParticleExtrapolation(RandomParticle, lowerExists);
767  }
768 
769  // Return the particle
771  return FoundParticle;
772 }
773 
776  G4bool LowerEnergyGroupExists )
777 {
779 
780  G4Ions* FoundParticle = NULL;
781  G4int NearestEnergy;
782  G4int NextNearestEnergy;
783 
784  // Check to see if we are extrapolating above or below the data set
785  if(LowerEnergyGroupExists == true)
786  {
787  NearestEnergy = YieldEnergyGroups_ - 1;
788  NextNearestEnergy = NearestEnergy - 1;
789  } else
790  {
791  NearestEnergy = 0;
792  NextNearestEnergy = 1;
793  }
794 
795  for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; Tree++)
796  {
797  FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk,
798  RandomParticle,
799  NearestEnergy,
800  NextNearestEnergy);
801  }
802 
804  return FoundParticle;
805 }
806 
809  G4int LowerEnergyGroup )
810 {
812 
813  G4Ions* FoundParticle = NULL;
814  G4int HigherEnergyGroup = LowerEnergyGroup + 1;
815 
816  for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; Tree++)
817  {
818  FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk,
819  RandomParticle,
820  LowerEnergyGroup,
821  HigherEnergyGroup);
822  }
823 
825  return FoundParticle;
826 }
827 
830  G4double RandomParticle,
831  G4int EnergyGroup1,
832  G4int EnergyGroup2 )
833 {
835 
836  G4Ions* Particle;
837 
838  // Verify that the branch exists
839  if(Branch == NULL)
840  {
841  Particle = NULL;
842  } else if(EnergyGroup1 >= Branch->IncidentEnergiesCount
843  || EnergyGroup2 >= Branch->IncidentEnergiesCount
844  || EnergyGroup1 == EnergyGroup2
845  || Branch->IncidentEnergies[EnergyGroup1] == Branch->IncidentEnergies[EnergyGroup2])
846  {
847  // Set NULL if any invalid conditions exist
848  Particle = NULL;
849  } else
850  {
851  // Everything check out - proceed
852  G4Ions* FoundParticle = NULL;
853  G4double Intercept;
854  G4double Slope;
855  G4double RangeAtIncidentEnergy;
856  G4double Denominator = Branch->IncidentEnergies[EnergyGroup1] - Branch->IncidentEnergies[EnergyGroup2];
857 
858  // Calculate the lower probability bounds
859  Slope = (Branch->ProbabilityRangeBottom[EnergyGroup1] - Branch->ProbabilityRangeBottom[EnergyGroup2]) / Denominator;
860  Intercept = Branch->ProbabilityRangeBottom[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
861  RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
862 
863  // Go right if the particle is below the probability bounds
864  if(RandomParticle < RangeAtIncidentEnergy)
865  {
866  FoundParticle = FindParticleBranchSearch(Branch->Left,
867  RandomParticle,
868  EnergyGroup1,
869  EnergyGroup2);
870  } else
871  {
872  // Calculate the upper probability bounds
873  Slope = (Branch->ProbabilityRangeTop[EnergyGroup1] - Branch->ProbabilityRangeTop[EnergyGroup2]) / Denominator;
874  Intercept = Branch->ProbabilityRangeTop[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
875  RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
876 
877  // Go left if the particle is above the probability bounds
878  if(RandomParticle > RangeAtIncidentEnergy)
879  {
880  FoundParticle = FindParticleBranchSearch(Branch->Right,
881  RandomParticle,
882  EnergyGroup1,
883  EnergyGroup2);
884  } else
885  {
886  // If the particle is bounded then we found it!
887  FoundParticle = Branch->Particle;
888  }
889  }
890 
891  Particle = FoundParticle;
892  }
893 
895  return Particle;
896 }
897 
899 GenerateAlphas( std::vector< G4ReactionProduct* >* Alphas )
900 {
902 
903  // Throw the dice to determine if ternary fission occurs
905  if(MakeAlphas)
906  {
907  G4int NumberOfAlphasToProduce;
908 
909  // Determine how many alpha particles to produce for the ternary fission
910  if(AlphaProduction_ < 0)
911  {
912  NumberOfAlphasToProduce = RandomEngine_->G4SampleIntegerGaussian(AlphaProduction_ * -1,
913  1,
915  } else
916  {
917  NumberOfAlphasToProduce = (G4int)AlphaProduction_;
918  }
919 
920  //TK modifed 131108
921  //Alphas->resize(NumberOfAlphasToProduce);
922  for(int i = 0; i < NumberOfAlphasToProduce; i++)
923  {
924  // Set the G4Ions as an alpha particle
925  Alphas->push_back(new G4ReactionProduct(AlphaDefinition_));
926 
927  // Remove 4 nucleons (2 protons and 2 neutrons) for each alpha added
928  RemainingZ_ -= 2;
929  RemainingA_ -= 4;
930  }
931  }
932 
934 }
935 
937 GenerateNeutrons( std::vector< G4ReactionProduct* >* Neutrons )
938 {
940 
941  G4int NeutronProduction;
943 
944  //TK modifed 131108
945  //Neutrons->resize(NeutronProduction);
946  for(int i = 0; i < NeutronProduction; i++)
947  {
948  // Define the fragment as a neutron
949  Neutrons->push_back(new G4ReactionProduct(NeutronDefinition_));
950 
951  // Remove 1 nucleon for each neutron added
952  RemainingA_--;
953  }
954 
956 }
957 
960  //TK modified 131108
961  //G4FFGEnumerations::MetaState MetaState )
962  G4FFGEnumerations::MetaState /*MetaState*/ )
963 {
965 
966  G4Ions* Temp;
967 
968  // Break Product down into its A and Z components
969  G4int A = Product % 1000; // Extract A
970  G4int Z = (Product - A) / 1000; // Extract Z
971 
972  // Check to see if the particle is registered using the PDG code
973  // TODO Add metastable state when supported by G4IonTable::GetIon()
974  Temp = reinterpret_cast<G4Ions*>(IonTable_->GetIon(Z, A));
975 
976  // Removed in favor of the G4IonTable::GetIon() method
977 // // Register the particle if it does not exist
978 // if(Temp == NULL)
979 // {
980 // // Define the particle properties
981 // G4String Name = MakeIsotopeName(Product, MetaState);
982 // // Calculate the rest mass using a function already in Geant4
983 // G4double Mass = G4NucleiProperties::
984 // GetNuclearMass((double)A, (double)Z );
985 // G4double Charge = Z*eplus;
986 // G4int BaryonNum = A;
987 // G4bool Stable = TRUE;
988 //
989 // // I am unsure about the following properties:
990 // // 2*Spin, Parity, C-conjugation, 2*Isospin, 2*Isospin3, G-parity.
991 // // Perhaps is would be a good idea to have a physicist familiar with
992 // // Geant4 nomenclature to review and correct these properties.
993 // Temp = new G4Ions (
994 // // Name Mass Width Charge
995 // Name, Mass, 0.0, Charge,
996 //
997 // // 2*Spin Parity C-conjugation 2*Isospin
998 // 0, 1, 0, 0,
999 //
1000 // // 2*Isospin3 G-parity Type Lepton number
1001 // 0, 0, "nucleus", 0,
1002 //
1003 // // Baryon number PDG encoding Stable Lifetime
1004 // BaryonNum, PDGCode, Stable, -1,
1005 //
1006 // // Decay table Shortlived SubType Anti_encoding
1007 // NULL, FALSE, "generic", 0,
1008 //
1009 // // Excitation
1010 // 0.0);
1011 // Temp->SetPDGMagneticMoment(0.0);
1012 //
1013 // // Declare that there is no anti-particle
1014 // Temp->SetAntiPDGEncoding(0);
1015 //
1016 // // Define the processes to use in transporting the particles
1017 // std::ostringstream osAdd;
1018 // osAdd << "/run/particle/addProcManager " << Name;
1019 // G4String cmdAdd = osAdd.str();
1020 //
1021 // // set /control/verbose 0
1022 // G4int tempVerboseLevel = G4UImanager::GetUIpointer()->GetVerboseLevel();
1023 // G4UImanager::GetUIpointer()->SetVerboseLevel(0);
1024 //
1025 // // issue /run/particle/addProcManage
1026 // G4UImanager::GetUIpointer()->ApplyCommand(cmdAdd);
1027 //
1028 // // retreive /control/verbose
1029 // G4UImanager::GetUIpointer()->SetVerboseLevel(tempVerboseLevel);
1030 // }
1031 
1033  return Temp;
1034 }
1035 
1038 {
1040 
1041  // Generate the file location starting in the Geant4 data directory
1042  std::ostringstream DirectoryName;
1043  DirectoryName << getenv("G4NEUTRONHPDATA") << G4FFGDefaultValues::ENDFFissionDataLocation;
1044 
1045  // Return the directory structure
1047  return DirectoryName.str();
1048 }
1049 
1053 {
1055 
1056 
1057  // Create the unique identifying name for the particle
1058  std::ostringstream FileName;
1059 
1060  // Determine if a leading 0 is needed (ZZZAAA or 0ZZAAA)
1061  if(Isotope < 100000)
1062  {
1063  FileName << "0";
1064  }
1065 
1066  // Add the name of the element and the extension
1067  FileName << MakeIsotopeName(Isotope, MetaState) << ".fpy";
1068 
1070  return FileName.str();
1071 }
1072 
1075 {
1077 
1078  G4DynamicParticle* DynamicParticle = new G4DynamicParticle(ReactionProduct->GetDefinition(), ReactionProduct->GetMomentum());
1079 
1081  return DynamicParticle;
1082 }
1083 
1087 {
1089 
1090  // Break Product down into its A and Z components
1091  G4int A = Isotope % 1000;
1092  G4int Z = (Isotope - A) / 1000;
1093 
1094  // Create the unique identifying name for the particle
1095  std::ostringstream IsotopeName;
1096 
1097  IsotopeName << Z << "_" << A;
1098 
1099  // If it is metastable then append "m" to the name
1100  if(MetaState != G4FFGEnumerations::GROUND_STATE)
1101  {
1102  IsotopeName << "m";
1103 
1104  // If it is a second isomeric state then append "2" to the name
1105  if(MetaState == G4FFGEnumerations::META_2)
1106  {
1107  IsotopeName << "2";
1108  }
1109  }
1110  // Add the name of the element and the extension
1111  IsotopeName << "_" << ElementNames_->theString[Z - 1];
1112 
1114  return IsotopeName.str();
1115 }
1116 
1118 MakeTrees( void )
1119 {
1121 
1122  // Allocate the space
1123  // We will make each tree a binary search
1124  // The maximum number of iterations required to find a single fission product
1125  // based on it's probability is defined by the following:
1126  // x = number of fission products
1127  // Trees = T(x) = ceil( ln(x) )
1128  // Rows/Tree = R(x) = ceil(( sqrt( (8 * x / T(x)) + 1) - 1) / 2)
1129  // Maximum = M(x) = T(x) + R(x)
1130  // Results: x => M(x)
1131  // 10 5
1132  // 100 10
1133  // 1000 25
1134  // 10000 54
1135  // 100000 140
1138 
1139  // Initialize the range of each node
1140  for(G4int i = 0; i < TreeCount_; i++)
1141  {
1143  Trees_[i].Trunk = NULL;
1144  Trees_[i].BranchCount = 0;
1145  Trees_[i].IsEnd = FALSE;
1146  }
1147  // Mark the last tree as the ending tree
1148  Trees_[TreeCount_ - 1].IsEnd = TRUE;
1149 
1151 }
1152 
1155 {
1157 
1158  G4int ProductCount = ENDFData_->G4GetNumberOfFissionProducts();
1159  BranchCount_ = 0;
1161 
1162  // Loop through all the products
1163  for(G4int i = 0; i < ProductCount; i++)
1164  {
1165  // Acquire the data and sort it
1167  }
1168 
1169  // Generate the true normalization factor, since round-off errors may result
1170  // in non-singular normalization of the data files. Also, reset DataTotal_
1171  // since it is used by Renormalize() to set the probability segments.
1174 
1175  // Go through all the trees one at a time
1176  for(G4int i = 0; i < TreeCount_; i++)
1177  {
1178  Renormalize(Trees_[i].Trunk);
1179  // Set the max range of the tree to DataTotal
1180  G4ArrayOps::Copy(YieldEnergyGroups_, Trees_[i].ProbabilityRangeEnd, DataTotal_);
1181  }
1182 
1184 }
1185 
1188 {
1190 
1191  // Check to see if Branch exists. Branch will be a null pointer if it
1192  // doesn't exist
1193  if(Branch != NULL)
1194  {
1195  // Call the lower branch to set the probability segment first, since it
1196  // supposed to have a lower probability segment that this node
1197  Renormalize(Branch->Left);
1198 
1199  // Set this node as the next sequential probability segment
1204 
1205  // Now call the upper branch to set those probability segments
1206  Renormalize(Branch->Right);
1207  }
1208 
1210 }
1211 
1213 SampleAlphaEnergies( std::vector< G4ReactionProduct* >* Alphas )
1214 {
1216 
1217  // The condition of sampling more energy from the fission products than is
1218  // alloted is statistically unfavorable, but it could still happen. The
1219  // do-while loop prevents such an occurrence from happening
1220  G4double MeanAlphaEnergy = 16.0;
1221  G4double TotalAlphaEnergy;
1222 
1223  do
1224  {
1225  G4double AlphaEnergy;
1226  TotalAlphaEnergy = 0;
1227 
1228  // Walk through the alpha particles one at a time and sample each's
1229  // energy
1230  for(unsigned int i = 0; i < Alphas->size(); i++)
1231  {
1232  AlphaEnergy = RandomEngine_->G4SampleGaussian(MeanAlphaEnergy,
1233  2.35,
1235  // Assign the energy to the alpha particle
1236  Alphas->at(i)->SetKineticEnergy(AlphaEnergy);
1237 
1238  // Add up the total amount of kinetic energy consumed.
1239  TotalAlphaEnergy += AlphaEnergy;
1240  }
1241 
1242  // If true, decrement the mean alpha energy by 0.1 and try again.
1243  MeanAlphaEnergy -= 0.1;
1244  } while(TotalAlphaEnergy >= RemainingEnergy_); // Loop checking, 11.05.2015, T. Koi
1245 
1246  // Subtract the total amount of energy that was assigned.
1247  RemainingEnergy_ -= TotalAlphaEnergy;
1248 
1250 }
1251 
1253 SampleGammaEnergies( std::vector< G4ReactionProduct* >* Gammas )
1254 {
1256 
1257  // Make sure that there is energy to assign to the gamma rays
1258  if(RemainingEnergy_ != 0)
1259  {
1260  G4double SampleEnergy;
1261 
1262  // Sample from RemainingEnergy until it is all gone. Also,
1263  // RemainingEnergy should not be smaller than
1264  // G4FFGDefaultValues::MeanGammaEnergy. This will prevent the
1265  // sampling of a fractional portion of the Gaussian distribution
1266  // in an attempt to find a new gamma ray energy.
1267  G4int icounter=0;
1268  G4int icounter_max=1024;
1269  while(RemainingEnergy_ >= G4FFGDefaultValues::MeanGammaEnergy ) // Loop checking, 11.05.2015, T. Koi
1270  {
1271  icounter++;
1272  if ( icounter > icounter_max ) {
1273  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
1274  break;
1275  }
1276  SampleEnergy = RandomEngine_->
1278  // Make sure that we didn't sample more energy than was available
1279  if(SampleEnergy <= RemainingEnergy_)
1280  {
1281  // If this energy assignment would leave less energy than the
1282  // 'intrinsic' minimal energy of a gamma ray then just assign
1283  // all of the remaining energy
1284  if(RemainingEnergy_ - SampleEnergy < 100 * keV)
1285  {
1286  SampleEnergy = RemainingEnergy_;
1287  }
1288 
1289  // Create the new particle
1290  Gammas->push_back(new G4ReactionProduct());
1291 
1292  // Set the properties
1293  Gammas->back()->SetDefinition(GammaDefinition_);
1294  Gammas->back()->SetTotalEnergy(SampleEnergy);
1295 
1296  // Calculate how much is left
1297  RemainingEnergy_ -= SampleEnergy;
1298  }
1299  }
1300 
1301  // If there is anything left over, the energy must be above 100 keV but
1302  // less than G4FFGDefaultValues::MeanGammaEnergy. Arbitrarily assign
1303  // RemainingEnergy to a new particle
1304  if(RemainingEnergy_ > 0)
1305  {
1306  SampleEnergy = RemainingEnergy_;
1307  Gammas->push_back(new G4ReactionProduct());
1308 
1309  // Set the properties
1310  Gammas->back()->SetDefinition(GammaDefinition_);
1311  Gammas->back()->SetTotalEnergy(SampleEnergy);
1312 
1313  // Calculate how much is left
1314  RemainingEnergy_ -= SampleEnergy;
1315  }
1316  }
1317 
1319 }
1320 
1322 SampleNeutronEnergies( std::vector< G4ReactionProduct* >* Neutrons )
1323 {
1325 
1326  // The condition of sampling more energy from the fission products than is
1327  // alloted is statistically unfavorable, but it could still happen. The
1328  // do-while loop prevents such an occurrence from happening
1329  G4double TotalNeutronEnergy;
1330  G4double NeutronEnergy;
1331 
1332  // Make sure that we don't sample more energy than is available
1333  G4int icounter=0;
1334  G4int icounter_max=1024;
1335  do
1336  {
1337  icounter++;
1338  if ( icounter > icounter_max ) {
1339  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
1340  break;
1341  }
1342  TotalNeutronEnergy = 0;
1343 
1344  // Walk through the neutrons one at a time and sample the energies.
1345  // The gamma rays have not yet been sampled, so the last neutron will
1346  // have a NULL value for NextFragment
1347  for(unsigned int i = 0; i < Neutrons->size(); i++)
1348  {
1349  // Assign the energy to the neutron
1350  NeutronEnergy = RandomEngine_->G4SampleWatt(Isotope_, Cause_, IncidentEnergy_);
1351  Neutrons->at(i)->SetKineticEnergy(NeutronEnergy);
1352 
1353  // Add up the total amount of kinetic energy consumed.
1354  TotalNeutronEnergy +=NeutronEnergy;
1355  }
1356  } while (TotalNeutronEnergy > RemainingEnergy_); // Loop checking, 11.05.2015, T. Koi
1357 
1358  // Subtract the total amount of energy that was assigned.
1359  RemainingEnergy_ -= TotalNeutronEnergy;
1360 
1362 }
1363 
1365 SetNubar( void )
1366 {
1368 
1369  G4int* WhichNubar;
1370  G4int* NubarWidth;
1371  G4double XFactor, BFactor;
1372 
1374  {
1375  WhichNubar = const_cast<G4int*>(&SpontaneousNubar_[0][0]);
1376  NubarWidth = const_cast<G4int*>(&SpontaneousNubarWidth_[0][0]);
1377  } else
1378  {
1379  WhichNubar = const_cast<G4int*>(&NeutronInducedNubar_[0][0]);
1380  NubarWidth = const_cast<G4int*>(&NeutronInducedNubarWidth_[0][0]);
1381  }
1382 
1383  XFactor = G4Pow::GetInstance()->powA(10.0, -13.0);
1384  BFactor = G4Pow::GetInstance()->powA(10.0, -4.0);
1385  Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor
1386  + *(WhichNubar + 2) * BFactor;
1387  while(*WhichNubar != -1) // Loop checking, 11.05.2015, T. Koi
1388  {
1389  if(*WhichNubar == Isotope_)
1390  {
1391  Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor
1392  + *(WhichNubar + 2) * BFactor;
1393 
1394  break;
1395  }
1396  WhichNubar += 3;
1397  }
1398 
1399  XFactor = G4Pow::GetInstance()->powN((G4double)10, -6);
1400  NubarWidth_ = *(NubarWidth + 1) * XFactor;
1401  while(*WhichNubar != -1) // Loop checking, 11.05.2015, T. Koi
1402  {
1403  if(*WhichNubar == Isotope_)
1404  {
1405  NubarWidth_ = *(NubarWidth + 1) * XFactor;
1406 
1407  break;
1408  }
1409  WhichNubar += 2;
1410  }
1411 
1413 }
1414 
1417 {
1419 
1420  // Initialize the new branch
1421  ProbabilityBranch* NewBranch = new ProbabilityBranch;
1423  NewBranch->Left = NULL;
1424  NewBranch->Right = NULL;
1425  NewBranch->Particle = GetParticleDefinition(YieldData->GetProduct(), YieldData->GetMetaState());
1426  NewBranch->IncidentEnergies = new G4double[YieldEnergyGroups_];
1427  NewBranch->ProbabilityRangeTop = new G4double[YieldEnergyGroups_];
1428  NewBranch->ProbabilityRangeBottom = new G4double[YieldEnergyGroups_];
1429  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->ProbabilityRangeTop, YieldData->GetYieldProbability());
1430  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->IncidentEnergies, YieldEnergies_);
1432 
1433  // Check to see if the this is the smallest/largest particle. First, check
1434  // to see if this is the first particle in the system
1435  if(SmallestZ_ == NULL)
1436  {
1437  SmallestZ_ = SmallestA_ = LargestZ_ = LargestA_ = NewBranch->Particle;
1438  } else
1439  {
1440  G4bool IsSmallerZ = NewBranch->Particle->GetAtomicNumber() < SmallestZ_->GetAtomicNumber();
1441  G4bool IsSmallerA = NewBranch->Particle->GetAtomicMass() < SmallestA_->GetAtomicMass();
1442  G4bool IsLargerZ = NewBranch->Particle->GetAtomicNumber() > LargestZ_->GetAtomicNumber();
1443  G4bool IsLargerA = NewBranch->Particle->GetAtomicMass() > LargestA_->GetAtomicMass();
1444 
1445  if(IsSmallerZ)
1446  {
1447  SmallestZ_ = NewBranch->Particle;
1448  }
1449 
1450  if(IsLargerZ)
1451  {
1452  LargestA_ = NewBranch->Particle;
1453  }
1454 
1455  if(IsSmallerA)
1456  {
1457  SmallestA_ = NewBranch->Particle;
1458  }
1459 
1460  if(IsLargerA)
1461  {
1462  LargestA_ = NewBranch->Particle;
1463  }
1464  }
1465 
1466  // Place the new branch
1467  // Determine which tree the new branch goes into
1468  G4int WhichTree = (G4int)floor((G4double)(BranchCount_ % TreeCount_));
1469  ProbabilityBranch** WhichBranch = &(Trees_[WhichTree].Trunk);
1470  Trees_[WhichTree].BranchCount++;
1471 
1472  // Search for the position
1473  // Determine where the branch goes
1474  G4int BranchPosition = (G4int)floor((G4double)(BranchCount_ / TreeCount_)) + 1;
1475 
1476  // Run through the tree until the end branch is reached
1477  while(BranchPosition > 1) // Loop checking, 11.05.2015, T. Koi
1478  {
1479  if(BranchPosition & 1)
1480  {
1481  // If the 1's bit is on then move to the next 'right' branch
1482  WhichBranch = &((*WhichBranch)->Right);
1483  } else
1484  {
1485  // If the 1's bit is off then move to the next 'down' branch
1486  WhichBranch = &((*WhichBranch)->Left);
1487  }
1488 
1489  BranchPosition >>= 1;
1490  }
1491 
1492  *WhichBranch = NewBranch;
1493  BranchCount_++;
1494 
1496 }
1497 
1500 {
1502 
1503  // Burn each tree, one by one
1504  G4int WhichTree = 0;
1505  while(Trees_[WhichTree].IsEnd != TRUE) // Loop checking, 11.05.2015, T. Koi
1506  {
1507  BurnTree(Trees_[WhichTree].Trunk);
1508  delete Trees_[WhichTree].Trunk;
1509  delete[] Trees_[WhichTree].ProbabilityRangeEnd;
1510  WhichTree++;
1511  }
1512 
1513  // Delete each dynamically allocated variable
1514  delete ENDFData_;
1515  delete[] Trees_;
1516  delete[] DataTotal_;
1517  delete[] MaintainNormalizedData_;
1518  delete ElementNames_;
1519  delete RandomEngine_;
1520 
1522 }
1523 
1526 {
1528 
1529  // Check to see it Branch exists. Branch will be a null pointer if it
1530  // doesn't exist
1531  if(Branch)
1532  {
1533  // Burn down before you burn up
1534  BurnTree(Branch->Left);
1535  delete Branch->Left;
1536  BurnTree(Branch->Right);
1537  delete Branch->Right;
1538 
1539  delete[] Branch->IncidentEnergies;
1540  delete[] Branch->ProbabilityRangeTop;
1541  delete[] Branch->ProbabilityRangeBottom;
1542  }
1543 
1545 }
1546 
#define G4FFG_RECURSIVE_FUNCTIONLEAVE__
void set(double x, double y, double z)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
const G4FFGEnumerations::YieldType YieldType_
G4ENDFYieldDataContainer * G4GetYield(G4int WhichYield)
G4FFGEnumerations::MetaState GetMetaState(void)
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double G4SampleGaussian(G4double Mean, G4double StdDev)
void Copy(G4int Elements, T *To, T *From)
Definition: G4ArrayOps.hh:63
G4double G4SampleWatt(G4int WhatIsotope, G4FFGEnumerations::FissionCause WhatCause, G4double WhatEnergy)
void Divide(G4int Elements, T *To, T *Numerator, T *Denominator=NULL)
Definition: G4ArrayOps.hh:178
G4double * ProbabilityRangeEnd
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
G4double G4SampleUniform(void)
G4FissionProductYieldDist(G4int WhichIsotope, G4FFGEnumerations::MetaState WhichMetaState, G4FFGEnumerations::FissionCause WhichCause, G4FFGEnumerations::YieldType WhichYieldType, std::istringstream &dataStream)
G4Ions * FindParticle(G4double RandomParticle)
void G4SetTernaryProbability(G4double TernaryProbability)
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4FFG_LOCATION__
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
static G4Alpha * Definition()
Definition: G4Alpha.cc:49
ProbabilityBranch * Right
void setR(double s)
static const G4int SpontaneousNubar_[][3]
G4DynamicParticleVector * G4GetFission(void)
double getY() const
void G4SetAlphaProduction(G4double WhatAlphaProduction)
G4String MakeFileName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
static const G4String theString[100]
G4double * ProbabilityRangeBottom
#define G4FFG_DATA_FUNCTIONENTER__
virtual void SortProbability(G4ENDFYieldDataContainer *YieldData)
int G4int
Definition: G4Types.hh:78
void setY(double)
void G4SetVerbosity(G4int WhatVerbosity)
const G4String & GetParticleName() const
void setZ(double)
G4double * ProbabilityRangeTop
void setX(double)
G4int GetAtomicNumber() const
static const G4int SpontaneousNubarWidth_[][2]
double getX() const
virtual void GenerateNeutrons(std::vector< G4ReactionProduct * > *Neutrons)
virtual void GenerateAlphas(std::vector< G4ReactionProduct * > *Alphas)
G4DynamicParticle * MakeG4DynamicParticle(G4ReactionProduct *)
const G4ParticleDefinition * GetDefinition() const
void SampleAlphaEnergies(std::vector< G4ReactionProduct * > *Alphas)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
tuple tree
Definition: gammaraytel.py:4
Definition: G4Ions.hh:51
#define G4FFG_DATA_FUNCTIONLEAVE__
bool G4bool
Definition: G4Types.hh:79
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
#define FALSE
Definition: globals.hh:52
static const G4int NeutronInducedNubar_[][3]
void Multiply(G4int Elements, T *To, T *M1, T *M2=NULL)
Definition: G4ArrayOps.hh:138
std::vector< G4DynamicParticle * > G4DynamicParticleVector
#define TRUE
Definition: globals.hh:55
G4Ions * GetParticleDefinition(G4int Product, G4FFGEnumerations::MetaState MetaState)
G4Ions * FindParticleBranchSearch(ProbabilityBranch *Branch, G4double RandomParticle, G4int EnergyGroup1, G4int EnergyGroup2)
void G4SetEnergy(G4double WhatIncidentEnergy)
static const G4double MeanGammaEnergy
static const char ENDFFissionDataLocation[]
void DeleteVectorOfPointers(std::vector< T > &Vector)
Definition: G4ArrayOps.hh:216
static const G4int Isotope
void setRThetaPhi(double r, double theta, double phi)
G4int GetAtomicMass() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
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)
double getZ() const
const G4FFGEnumerations::FissionCause Cause_
static constexpr double c_light
static const G4int NeutronInducedNubarWidth_[][2]
static G4Neutron * Definition()
Definition: G4Neutron.cc:54
static constexpr double GeV
Definition: G4SIunits.hh:217
void Set(G4int Elements, T *To, T Value)
Definition: G4ArrayOps.hh:51
void G4SetVerbosity(G4int WhatVerbosity)
G4ThreeVector GetMomentum() const
#define G4FFG_FUNCTIONLEAVE__
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
#define G4FFG_SPACING__
void SampleNeutronEnergies(std::vector< G4ReactionProduct * > *Neutrons)
static constexpr double pi
Definition: G4SIunits.hh:75
#define G4FFG_RECURSIVE_FUNCTIONENTER__
G4int G4GetNumberOfFissionProducts(void)
ProbabilityBranch * Left
void BurnTree(ProbabilityBranch *Branch)
double G4double
Definition: G4Types.hh:76
void G4SetVerbosity(G4int WhatVerbosity)
virtual G4Ions * GetFissionProduct(void)=0
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
#define G4FFG_FUNCTIONENTER__
G4String MakeIsotopeName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
static constexpr double pi
Definition: SystemOfUnits.h:54
void Add(G4int Elements, T *To, T *A1, T *A2=NULL)
Definition: G4ArrayOps.hh:77
G4double * G4GetEnergyGroupValues(void)
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
G4int G4SampleIntegerGaussian(G4double Mean, G4double StdDev)