Geant4.10
 All Data Structures 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 "G4SystemOfUnits.hh"
62 #include "G4HadFinalState.hh"
63 #include "G4DynamicParticle.hh"
65 #include "G4ReactionProduct.hh"
66 
67 #include "G4ArrayOps.hh"
68 #include "G4ENDFTapeRead.hh"
70 #include "G4FFGDebuggingMacros.hh"
71 #include "G4FFGDefaultValues.hh"
72 #include "G4FFGEnumerations.hh"
73 #include "G4FPYNubarValues.hh"
74 #include "G4FPYSamplingOps.hh"
75 #include "G4FPYTreeStructures.hh"
77 
78 using CLHEP::pi;
79 
82  G4FFGEnumerations::MetaState WhichMetaState,
84  G4FFGEnumerations::YieldType WhichYieldType,
85  std::istringstream& dataStream )
86 : Isotope_(WhichIsotope),
87  MetaState_(WhichMetaState),
88  Cause_(WhichCause),
89  YieldType_(WhichYieldType),
90  Verbosity_(G4FFGDefaultValues::Verbosity)
91 {
93 
94  try
95  {
96  // Initialize the class
97  Initialize(dataStream);
98  } catch (std::exception e)
99  {
101  throw e;
102  }
103 
105 }
106 
109  G4FFGEnumerations::MetaState WhichMetaState,
111  G4FFGEnumerations::YieldType WhichYieldType,
113  std::istringstream& dataStream)
114 : Isotope_(WhichIsotope),
115  MetaState_(WhichMetaState),
116  Cause_(WhichCause),
117  YieldType_(WhichYieldType),
118  Verbosity_(Verbosity)
119 {
121 
122  try
123  {
124  // Initialize the class
125  Initialize(dataStream);
126  } catch (std::exception e)
127  {
129  throw e;
130  }
131 
133 }
134 
135 void G4FissionProductYieldDist::
136 Initialize( std::istringstream& dataStream )
137 {
139 
140  IncidentEnergy_ = 0.0;
142  AlphaProduction_ = 0;
143  SetNubar();
144 
145  // Set miscellaneous variables
146  AlphaDefinition_ = reinterpret_cast<G4Ions*>(G4Alpha::Definition());
147  NeutronDefinition_ = reinterpret_cast<G4Ions*>(G4Neutron::Definition());
148  GammaDefinition_ = reinterpret_cast<G4Ions*>(G4Gamma::Definition());
150 
151  // Construct G4NeutronHPNames: provides access to the element names
153  // Get the pointer to G4ParticleTable: stores all G4Ions
155  // Construct the pointer to the random engine
156  // TODO Make G4FPSamplingOps a singleton so that only one instance is used across all classes
158 
159  try
160  {
161  // Read in and sort the probability data
162  ENDFData_ = new G4ENDFTapeRead(dataStream,
163  YieldType_,
164  Cause_,
165  Verbosity_);
166 // ENDFData_ = new G4ENDFTapeRead(MakeDirectoryName(),
167 // MakeFileName(Isotope_, MetaState_),
168 // YieldType_,
169 // Cause_);
175  MakeTrees();
177  } catch (std::exception& e)
178  {
179  delete ElementNames_;
180  delete RandomEngine_;
181 
183  throw e;
184  }
185 
187 }
188 
191 {
193 
194  // Check to see if the user has set the alpha production to a somewhat
195  // reasonable level
197 
198  // Generate the new G4DynamicParticle pointers to identify key locations in
199  // the G4DynamicParticle chain that will be passed to the G4FissionEvent
200  G4ReactionProduct* FirstDaughter = NULL;
201  G4ReactionProduct* SecondDaughter = NULL;
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* >;
205 
206  // Generate all the nucleonic fission products
207  // How many nucleons do we have to work with?
208  //TK modified 131108
209  //const G4int ParentA = Isotope_ % 1000;
210  //const G4int ParentZ = (Isotope_ - ParentA) / 1000;
211  const G4int ParentA = (Isotope_/10) % 1000;
212  const G4int ParentZ = ((Isotope_/10) - ParentA) / 1000;
213  RemainingA_ = ParentA;
214  RemainingZ_ = ParentZ;
215 
216  // Don't forget the extra nucleons depending on the fission cause
217  switch(Cause_)
218  {
220  ++RemainingA_;
221  break;
222 
224  ++RemainingZ_;
225  break;
226 
229  default:
230  // Nothing to do here
231  break;
232  }
233 
234  // Ternary fission can be set by the user. Thus, it is necessary to
235  // sample the alpha particle first and the first daughter product
236  // second. See the discussion in
237  // G4FissionProductYieldDist::G4GetFissionProduct() for more information
238  // as to why the fission events are sampled this way.
239  GenerateAlphas(Alphas);
240 
241  // Generate the first daughter product
242  FirstDaughter = new G4ReactionProduct(GetFissionProduct());
243  RemainingA_ -= FirstDaughter->GetDefinition()->GetAtomicMass();
244  RemainingZ_ -= FirstDaughter->GetDefinition()->GetAtomicNumber();
246  {
249 
250  G4cout << " -- First daughter product sampled" << G4endl;
252  G4cout << " Name: " << FirstDaughter->GetDefinition()->GetParticleName() << G4endl;
254  G4cout << " Z: " << FirstDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
256  G4cout << " A: " << FirstDaughter->GetDefinition()->GetAtomicMass() << G4endl;
258  G4cout << " Meta State: " << (FirstDaughter->GetDefinition()->GetPDGEncoding() % 10) << G4endl;
259  }
260 
261  GenerateNeutrons(Neutrons);
262 
263  // Now that all the nucleonic particles have been generated, we can
264  // calculate the composition of the second daughter product.
265  G4int NewIsotope = RemainingZ_ * 1000 + RemainingA_;
267  if(Verbosity_ & G4FFGEnumerations::DAUGHTER_INFO)
268  {
271 
272  G4cout << " -- Second daughter product sampled" << G4endl;
274  G4cout << " Name: " << SecondDaughter->GetDefinition()->GetParticleName() << G4endl;
276  G4cout << " Z: " << SecondDaughter->GetDefinition()->GetAtomicNumber() << G4endl;
278  G4cout << " A: " << SecondDaughter->GetDefinition()->GetAtomicMass() << G4endl;
280  G4cout << " Meta State: " << (SecondDaughter->GetDefinition()->GetPDGEncoding() % 10) << G4endl;
281  }
282 
283  // Calculate how much kinetic energy will be available
284  // 195 to 205 MeV are available in a fission reaction, but about 20 MeV
285  // are from delayed sources. We are concerned only with prompt sources,
286  // so sample a Gaussian distribution about 20 MeV and subtract the
287  // result from the total available energy. Also, the energy of fission
288  // neutrinos is neglected. Fission neutrinos would add ~11 MeV
289  // additional energy to the fission. (Ref 2)
290  // Finally, add in the kinetic energy contribution of the fission
291  // inducing particle, if any.
292  const G4double TotalKE = RandomEngine_->G4SampleUniform(195.0, 205.0) * MeV
293  - RandomEngine_->G4SampleGaussian(20.0, 3.0) * MeV
294  + IncidentEnergy_;
295  RemainingEnergy_ = TotalKE;
296 
297  // Calculate the energies of the alpha particles and neutrons
298  // Once again, since the alpha particles are user defined, we must
299  // sample their kinetic energy first. SampleAlphaEnergies() returns the
300  // amount of energy consumed by the alpha particles, so remove the total
301  // energy alloted to the alpha particles from the available energy
302  SampleAlphaEnergies(Alphas);
303 
304  // Second, the neutrons are sampled from the Watt fission spectrum.
305  SampleNeutronEnergies(Neutrons);
306 
307  // Calculate the total energy available to the daughter products
308  // A Gaussian distribution about the average calculated energy with
309  // a standard deviation of 1.5 MeV (Ref. 2) is used. Since the energy
310  // distribution is dependant on the alpha particle generation and the
311  // Watt fission sampling for neutrons, we only have the left-over energy
312  // to work with for the fission daughter products.
313  G4double FragmentsKE;
314  do
315  {
316  FragmentsKE = RandomEngine_->G4SampleGaussian(RemainingEnergy_, 1.5 *MeV);
317  } while(FragmentsKE > RemainingEnergy_);
318 
319  // Make sure that we don't produce any sub-gamma photons
320  if((RemainingEnergy_ - FragmentsKE) / (100 * keV) < 1.0)
321  {
322  FragmentsKE = RemainingEnergy_;
323  }
324 
325  // This energy has now been allotted to the fission fragments.
326  // Subtract FragmentsKE from the total available fission energy.
327  RemainingEnergy_ -= FragmentsKE;
328 
329  // Sample the energies of the gamma rays
330  // Assign the remainder, if any, of the energy to the gamma rays
331  SampleGammaEnergies(Gammas);
332 
333  // Prepare to balance the momenta of the system
334  // We will need these for sampling the angles and balancing the momenta
335  // of all the fission products
336  G4double NumeratorSqrt, NumeratorOther, Denominator;
337  G4double Theta, Phi, Magnitude;
338  G4ThreeVector Direction;
339  G4ParticleMomentum ResultantVector(0, 0, 0);
340 
341  if(Alphas->size() > 0)
342  {
343  // Sample the angles of the alpha particles and neutrons, then calculate
344  // the total moment contribution to the system
345  // The average angle of the alpha particles with respect to the
346  // light fragment is dependent on the ratio of the kinetic energies.
347  // This equation was determined by performing a fit on data from
348  // Ref. 3 using the website:
349  // http://soft.arquimedex.com/linear_regression.php
350  //
351  // RatioOfKE Angle (rad) Angle (degrees)
352  // 1.05 1.257 72
353  // 1.155 1.361 78
354  // 1.28 1.414 81
355  // 1.5 1.518 87
356  // 1.75 1.606 92
357  // 1.9 1.623 93
358  // 2.2 1.728 99
359  // This equation generates the angle in radians. If the RatioOfKE is
360  // greater than 2.25 the angle defaults to 1.3963 rad (100 degrees)
361  G4double MassRatio = FirstDaughter->GetDefinition()->GetPDGMass() / SecondDaughter->GetDefinition()->GetPDGMass();
362 
363  // Invert the mass ratio if the first daughter product is the lighter fragment
364  if(MassRatio < 1)
365  {
366  MassRatio = 1 / MassRatio;
367  }
368 
369  // The empirical equation is valid for mass ratios up to 2.75
370  if(MassRatio > 2.75)
371  {
372  MassRatio = 2.75;
373  }
374  const G4double MeanAlphaAngle = 0.3644 * MassRatio * MassRatio * MassRatio
375  - 1.9766 * MassRatio * MassRatio
376  + 3.8207 * MassRatio
377  - 0.9917;
378 
379  // Sample the directions of the alpha particles with respect to the
380  // light fragment. For the moment we will assume that the light
381  // fragment is traveling along the z-axis in the positive direction.
382  const G4double MeanAlphaAngleStdDev = 0.0523598776;
383  G4double PlusMinus;
384 
385  for(unsigned int i = 0; i < Alphas->size(); ++i)
386  {
387  PlusMinus = std::acos(RandomEngine_->G4SampleGaussian(0, MeanAlphaAngleStdDev)) - (pi / 2);
388  Theta = MeanAlphaAngle + PlusMinus;
389  if(Theta < 0)
390  {
391  Theta = 0.0 - Theta;
392  } else if(Theta > pi)
393  {
394  Theta = (2 * pi - Theta);
395  }
397 
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();
402  }
403  }
404 
405  // Sample the directions of the neutrons.
406  if(Neutrons->size() != 0)
407  {
408  for(unsigned int i = 0; i < Neutrons->size(); ++i)
409  {
410  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
412 
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();
417  }
418  }
419 
420  // Sample the directions of the gamma rays
421  if(Gammas->size() != 0)
422  {
423  for(unsigned int i = 0; i < Gammas->size(); ++i)
424  {
425  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
427 
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();
432  }
433  }
434 
435  // Calculate the momenta of the two daughter products
436  G4ReactionProduct* LightFragment;
437  G4ReactionProduct* HeavyFragment;
438  G4ThreeVector LightFragmentDirection;
439  G4ThreeVector HeavyFragmentDirection;
440  G4double ResultantX, ResultantY, ResultantZ;
441  ResultantX = ResultantVector.getX();
442  ResultantY = ResultantVector.getY();
443  ResultantZ = ResultantVector.getZ();
444 
445  if(FirstDaughter->GetDefinition()->GetPDGMass() < SecondDaughter->GetDefinition()->GetPDGMass())
446  {
447  LightFragment = FirstDaughter;
448  HeavyFragment = SecondDaughter;
449  } else
450  {
451  LightFragment = SecondDaughter;
452  HeavyFragment = FirstDaughter;
453  }
454  const G4double LightFragmentMass = LightFragment->GetDefinition()->GetPDGMass();
455  const G4double HeavyFragmentMass = HeavyFragment->GetDefinition()->GetPDGMass();
456 
457  LightFragmentDirection.setRThetaPhi(1.0, 0, 0);
458 
459  // Fit the momenta of the daughter products to the resultant vector of
460  // the remaining fission products. This will be done in the Cartesian
461  // coordinate system, not spherical. This is done using the following
462  // table listing the system momenta and the corresponding equations:
463  // X Y Z
464  //
465  // A 0 0 P
466  //
467  // B -R_x -R_y -P - R_z
468  //
469  // R R_x R_y R_z
470  //
471  // v = sqrt(2*m*k) -> k = v^2/(2*m)
472  // tk = k_A + k_B
473  // k_L = P^2/(2*m_L)
474  // k_H = ((-R_x)^2 + (-R_y)^2 + (-P - R_z)^2)/(2*m_H)
475  // where:
476  // P: momentum of the light daughter product
477  // R: the remaining fission products' resultant vector
478  // v: momentum
479  // m: mass
480  // k: kinetic energy
481  // tk: total kinetic energy available to the daughter products
482  //
483  // Below is the solved form for P, with the solution generated using
484  // the WolframAlpha website:
485  // http://www.wolframalpha.com/input/?i=
486  // solve+((-x)^2+%2B+(-y)^2+%2B+(-P-z)^2)%2F(2*B)+%2B+L^2%2F(2*A)+%3D+k+
487  // for+P
488  //
489  //
490  // nsqrt = sqrt(m_L*(m_L*(2*m_H*tk - R_x^2 - R_y^2) +
491  // m_H*(2*m_H*tk - R_x^2 - R_y^2 - R_z^2))
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)));
500 
501  // nother = m_L*R_z
502  NumeratorOther = LightFragmentMass * ResultantZ;
503 
504  // denom = m_L + m_H
505  Denominator = LightFragmentMass + HeavyFragmentMass;
506 
507  // P = (nsqrt + nother) / denom
508  const G4double LightFragmentMomentum = (NumeratorSqrt + NumeratorOther) / Denominator;
509  const G4double HeavyFragmentMomentum = std::sqrt(ResultantX * ResultantX
510  + ResultantY * ResultantY
511  + std::pow(LightFragmentMomentum + ResultantZ, 2));
512 
513  // Finally! We now have everything we need for the daughter products
514  LightFragment->SetMomentum(LightFragmentDirection * LightFragmentMomentum);
515  HeavyFragmentDirection.setX(-ResultantX);
516  HeavyFragmentDirection.setY(-ResultantY);
517  HeavyFragmentDirection.setZ((-LightFragmentMomentum - ResultantZ));
518  // Don't forget to normalize the vector to the unit sphere
519  HeavyFragmentDirection.setR(1.0);
520  HeavyFragment->SetMomentum(HeavyFragmentDirection * HeavyFragmentMomentum);
521 
522  if(Verbosity_ & (G4FFGEnumerations::DAUGHTER_INFO | G4FFGEnumerations::MOMENTUM_INFO))
523  {
526 
527  G4cout << " -- Daugher product momenta finalized" << G4endl;
529  }
530 
531  // Load all the particles into a contiguous set
532  //TK modifed 131108
533  //G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector(2 + Alphas->size() + Neutrons->size() + Gammas->size());
534  G4DynamicParticleVector* FissionProducts = new G4DynamicParticleVector();
535  // Load the fission fragments
536  FissionProducts->push_back(MakeG4DynamicParticle(LightFragment));
537  FissionProducts->push_back(MakeG4DynamicParticle(HeavyFragment));
538  // Load the neutrons
539  for(unsigned int i = 0; i < Neutrons->size(); i++)
540  {
541  FissionProducts->push_back(MakeG4DynamicParticle(Neutrons->at(i)));
542  }
543  // Load the gammas
544  for(unsigned int i = 0; i < Gammas->size(); i++)
545  {
546  FissionProducts->push_back(MakeG4DynamicParticle(Gammas->at(i)));
547  }
548  // Load the alphas
549  for(unsigned int i = 0; i < Alphas->size(); i++)
550  {
551  FissionProducts->push_back(MakeG4DynamicParticle(Alphas->at(i)));
552  }
553 
554  // Rotate the system to a random location so that the light fission fragment
555  // is not always traveling along the positive z-axis
556  // Sample Theta and Phi.
557  G4ThreeVector RotationAxis;
558 
559  Theta = std::acos(RandomEngine_->G4SampleUniform(-1, 1));
561  RotationAxis.setRThetaPhi(1.0, Theta, Phi);
562 
563  // We will also check the net momenta
564  ResultantVector.set(0.0, 0.0, 0.0);
565  for(unsigned int i = 0; i < FissionProducts->size(); i++)
566  {
567  Direction = FissionProducts->at(i)->GetMomentumDirection();
568  Direction.rotateUz(RotationAxis);
569  FissionProducts->at(i)->SetMomentumDirection(Direction);
570  ResultantVector += FissionProducts->at(i)->GetMomentum();
571  }
572 
573  // Warn if the sum momenta of the system is not within a reasonable
574  // tolerance
575  G4double PossibleImbalance = ResultantVector.mag();
576  if(PossibleImbalance > 0.01)
577  {
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()",
583  Temp.str().c_str(),
584  JustWarning,
585  "Results may not be valid");
586  }
587 
588  // Clean up
589  delete FirstDaughter;
590  delete SecondDaughter;
594 
596  return FissionProducts;
597 }
598 
601 {
603 
605 
607  return Product;
608 }
609 
611 G4SetAlphaProduction(G4double WhatAlphaProduction)
612 {
614 
615  AlphaProduction_ = WhatAlphaProduction;
616 
618 }
619 
621 G4SetEnergy( G4double WhatIncidentEnergy )
622 {
624 
626  {
627  IncidentEnergy_ = WhatIncidentEnergy;
628  } else
629  {
630  IncidentEnergy_ = 0 * GeV;
631  }
632 
634 }
635 
637 G4SetTernaryProbability( G4double WhatTernaryProbability )
638 {
640 
641  TernaryProbability_ = WhatTernaryProbability;
642 
644 }
645 
648 {
650 
652 
653  ENDFData_->G4SetVerbosity(Verbosity_);
654  RandomEngine_->G4SetVerbosity(Verbosity_);
655 
657 }
658 
661 {
663 
664  // This provides comfortable breathing room at 16 MeV per alpha
665  if(AlphaProduction_ > 10)
666  {
667  AlphaProduction_ = 10;
668  } else if(AlphaProduction_ < -7)
669  {
670  AlphaProduction_ = -7;
671  }
672 
674 }
675 
677 FindParticle( G4double RandomParticle )
678 {
680 
681  // Determine which energy group is currently in use
682  G4bool isExact = false;
683  G4bool lowerExists = false;
684  G4bool higherExists = false;
685  G4int energyGroup;
686  for(energyGroup = 0; energyGroup < YieldEnergyGroups_; energyGroup++)
687  {
688  if(IncidentEnergy_ == YieldEnergies_[energyGroup])
689  {
690  isExact = true;
691  break;
692  }
693 
694  if(energyGroup == 0 && IncidentEnergy_ < YieldEnergies_[energyGroup])
695  {
696  // Break if the energy is less than the lowest energy
697  higherExists = true;
698  break;
699  } else if(energyGroup == YieldEnergyGroups_ - 1)
700  {
701  // The energy is greater than any values in the yield data.
702  lowerExists = true;
703  break;
704  } else
705  {
706  // Break if the energy is less than the lowest energy
707  if(IncidentEnergy_ > YieldEnergies_[energyGroup])
708  {
709  energyGroup--;
710  lowerExists = true;
711  higherExists = true;
712  break;
713  }
714  }
715  }
716 
717  // Determine which particle it is
718  G4Ions* FoundParticle = NULL;
719  if(isExact || YieldEnergyGroups_ == 1)
720  {
721  // Determine which tree contains the random value
722  G4int tree;
723  for(tree = 0; tree < TreeCount_; tree++)
724  {
725  // Break if a tree is identified as containing the random particle
726  if(RandomParticle <= Trees_[tree].ProbabilityRangeEnd[energyGroup])
727  {
728  break;
729  }
730  }
731  ProbabilityBranch* Branch = Trees_[tree].Trunk;
732 
733  // Iteratively traverse the tree until the particle addressed by the random
734  // variable is found
735  G4bool RangeIsSmaller;
736  G4bool RangeIsGreater;
737  while((RangeIsSmaller = (RandomParticle < Branch->ProbabilityRangeBottom[energyGroup]))
738  || (RangeIsGreater = (RandomParticle > Branch->ProbabilityRangeTop[energyGroup])))
739  {
740  if(RangeIsSmaller)
741  {
742  Branch = Branch->Left;
743  } else {
744  Branch = Branch->Right;
745  }
746  }
747 
748  FoundParticle = Branch->Particle;
749  } else if(lowerExists && higherExists)
750  {
751  // We need to do some interpolation
752  FoundParticle = FindParticleInterpolation(RandomParticle, energyGroup);
753  } else
754  {
755  // We need to do some extrapolation
756  FoundParticle = FindParticleExtrapolation(RandomParticle, lowerExists);
757  }
758 
759  // Return the particle
761  return FoundParticle;
762 }
763 
766  G4bool LowerEnergyGroupExists )
767 {
769 
770  G4Ions* FoundParticle = NULL;
771  G4int NearestEnergy;
772  G4int NextNearestEnergy;
773 
774  // Check to see if we are extrapolating above or below the data set
775  if(LowerEnergyGroupExists == true)
776  {
777  NearestEnergy = YieldEnergyGroups_ - 1;
778  NextNearestEnergy = NearestEnergy - 1;
779  } else
780  {
781  NearestEnergy = 0;
782  NextNearestEnergy = 1;
783  }
784 
785  for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; Tree++)
786  {
787  FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk,
788  RandomParticle,
789  NearestEnergy,
790  NextNearestEnergy);
791  }
792 
794  return FoundParticle;
795 }
796 
799  G4int LowerEnergyGroup )
800 {
802 
803  G4Ions* FoundParticle = NULL;
804  G4int HigherEnergyGroup = LowerEnergyGroup + 1;
805 
806  for(G4int Tree = 0; Tree < TreeCount_ && FoundParticle == NULL; Tree++)
807  {
808  FoundParticle = FindParticleBranchSearch(Trees_[Tree].Trunk,
809  RandomParticle,
810  LowerEnergyGroup,
811  HigherEnergyGroup);
812  }
813 
815  return FoundParticle;
816 }
817 
820  G4double RandomParticle,
821  G4int EnergyGroup1,
822  G4int EnergyGroup2 )
823 {
825 
826  G4Ions* Particle;
827 
828  // Verify that the branch exists
829  if(Branch == NULL)
830  {
831  Particle = NULL;
832  } else if(EnergyGroup1 >= Branch->IncidentEnergiesCount
833  || EnergyGroup2 >= Branch->IncidentEnergiesCount
834  || EnergyGroup1 == EnergyGroup2
835  || Branch->IncidentEnergies[EnergyGroup1] == Branch->IncidentEnergies[EnergyGroup2])
836  {
837  // Set NULL if any invalid conditions exist
838  Particle = NULL;
839  } else
840  {
841  // Everything check out - proceed
842  G4Ions* FoundParticle = NULL;
843  G4double Intercept;
844  G4double Slope;
845  G4double RangeAtIncidentEnergy;
846  G4double Denominator = Branch->IncidentEnergies[EnergyGroup1] - Branch->IncidentEnergies[EnergyGroup2];
847 
848  // Calculate the lower probability bounds
849  Slope = (Branch->ProbabilityRangeBottom[EnergyGroup1] - Branch->ProbabilityRangeBottom[EnergyGroup2]) / Denominator;
850  Intercept = Branch->ProbabilityRangeBottom[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
851  RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
852 
853  // Go right if the particle is below the probability bounds
854  if(RandomParticle < RangeAtIncidentEnergy)
855  {
856  FoundParticle = FindParticleBranchSearch(Branch->Left,
857  RandomParticle,
858  EnergyGroup1,
859  EnergyGroup2);
860  } else
861  {
862  // Calculate the upper probability bounds
863  Slope = (Branch->ProbabilityRangeTop[EnergyGroup1] - Branch->ProbabilityRangeTop[EnergyGroup2]) / Denominator;
864  Intercept = Branch->ProbabilityRangeTop[EnergyGroup1] - Slope * Branch->IncidentEnergies[EnergyGroup1];
865  RangeAtIncidentEnergy = Slope * IncidentEnergy_ + Intercept;
866 
867  // Go left if the particle is above the probability bounds
868  if(RandomParticle > RangeAtIncidentEnergy)
869  {
870  FoundParticle = FindParticleBranchSearch(Branch->Right,
871  RandomParticle,
872  EnergyGroup1,
873  EnergyGroup2);
874  } else
875  {
876  // If the particle is bounded then we found it!
877  FoundParticle = Branch->Particle;
878  }
879  }
880 
881  Particle = FoundParticle;
882  }
883 
885  return Particle;
886 }
887 
889 GenerateAlphas( std::vector< G4ReactionProduct* >* Alphas )
890 {
892 
893  // Throw the dice to determine if ternary fission occurs
895  if(MakeAlphas)
896  {
897  G4int NumberOfAlphasToProduce;
898 
899  // Determine how many alpha particles to produce for the ternary fission
900  if(AlphaProduction_ < 0)
901  {
902  NumberOfAlphasToProduce = RandomEngine_->G4SampleIntegerGaussian(AlphaProduction_ * -1,
903  1,
905  } else
906  {
907  NumberOfAlphasToProduce = (G4int)AlphaProduction_;
908  }
909 
910  //TK modifed 131108
911  //Alphas->resize(NumberOfAlphasToProduce);
912  for(int i = 0; i < NumberOfAlphasToProduce; i++)
913  {
914  // Set the G4Ions as an alpha particle
915  Alphas->push_back(new G4ReactionProduct(AlphaDefinition_));
916 
917  // Remove 4 nucleons (2 protons and 2 neutrons) for each alpha added
918  RemainingZ_ -= 2;
919  RemainingA_ -= 4;
920  }
921  }
922 
924 }
925 
927 GenerateNeutrons( std::vector< G4ReactionProduct* >* Neutrons )
928 {
930 
931  G4int NeutronProduction;
933 
934  //TK modifed 131108
935  //Neutrons->resize(NeutronProduction);
936  for(int i = 0; i < NeutronProduction; i++)
937  {
938  // Define the fragment as a neutron
939  Neutrons->push_back(new G4ReactionProduct(NeutronDefinition_));
940 
941  // Remove 1 nucleon for each neutron added
942  RemainingA_--;
943  }
944 
946 }
947 
950  //TK modified 131108
951  //G4FFGEnumerations::MetaState MetaState )
952  G4FFGEnumerations::MetaState /*MetaState*/ )
953 {
955 
956  G4Ions* Temp;
957 
958  // Break Product down into its A and Z components
959  G4int A = Product % 1000; // Extract A
960  G4int Z = (Product - A) / 1000; // Extract Z
961 
962  // Check to see if the particle is registered using the PDG code
963  // TODO Add metastable state when supported by G4IonTable::GetIon()
964  Temp = reinterpret_cast<G4Ions*>(IonTable_->GetIon(Z, A));
965 
966  // Removed in favor of the G4IonTable::GetIon() method
967 // // Register the particle if it does not exist
968 // if(Temp == NULL)
969 // {
970 // // Define the particle properties
971 // G4String Name = MakeIsotopeName(Product, MetaState);
972 // // Calculate the rest mass using a function already in Geant4
973 // G4double Mass = G4NucleiProperties::
974 // GetNuclearMass((double)A, (double)Z );
975 // G4double Charge = Z*eplus;
976 // G4int BaryonNum = A;
977 // G4bool Stable = TRUE;
978 //
979 // // I am unsure about the following properties:
980 // // 2*Spin, Parity, C-conjugation, 2*Isospin, 2*Isospin3, G-parity.
981 // // Perhaps is would be a good idea to have a physicist familiar with
982 // // Geant4 nomenclature to review and correct these properties.
983 // Temp = new G4Ions (
984 // // Name Mass Width Charge
985 // Name, Mass, 0.0, Charge,
986 //
987 // // 2*Spin Parity C-conjugation 2*Isospin
988 // 0, 1, 0, 0,
989 //
990 // // 2*Isospin3 G-parity Type Lepton number
991 // 0, 0, "nucleus", 0,
992 //
993 // // Baryon number PDG encoding Stable Lifetime
994 // BaryonNum, PDGCode, Stable, -1,
995 //
996 // // Decay table Shortlived SubType Anti_encoding
997 // NULL, FALSE, "generic", 0,
998 //
999 // // Excitation
1000 // 0.0);
1001 // Temp->SetPDGMagneticMoment(0.0);
1002 //
1003 // // Declare that there is no anti-particle
1004 // Temp->SetAntiPDGEncoding(0);
1005 //
1006 // // Define the processes to use in transporting the particles
1007 // std::ostringstream osAdd;
1008 // osAdd << "/run/particle/addProcManager " << Name;
1009 // G4String cmdAdd = osAdd.str();
1010 //
1011 // // set /control/verbose 0
1012 // G4int tempVerboseLevel = G4UImanager::GetUIpointer()->GetVerboseLevel();
1013 // G4UImanager::GetUIpointer()->SetVerboseLevel(0);
1014 //
1015 // // issue /run/particle/addProcManage
1016 // G4UImanager::GetUIpointer()->ApplyCommand(cmdAdd);
1017 //
1018 // // retreive /control/verbose
1019 // G4UImanager::GetUIpointer()->SetVerboseLevel(tempVerboseLevel);
1020 // }
1021 
1023  return Temp;
1024 }
1025 
1028 {
1030 
1031  // Generate the file location starting in the Geant4 data directory
1032  std::ostringstream DirectoryName;
1033  DirectoryName << getenv("G4NEUTRONHPDATA") << G4FFGDefaultValues::ENDFFissionDataLocation;
1034 
1035  // Return the directory structure
1037  return DirectoryName.str();
1038 }
1039 
1043 {
1045 
1046 
1047  // Create the unique identifying name for the particle
1048  std::ostringstream FileName;
1049 
1050  // Determine if a leading 0 is needed (ZZZAAA or 0ZZAAA)
1051  if(Isotope < 100000)
1052  {
1053  FileName << "0";
1054  }
1055 
1056  // Add the name of the element and the extension
1057  FileName << MakeIsotopeName(Isotope, MetaState) << ".fpy";
1058 
1060  return FileName.str();
1061 }
1062 
1065 {
1067 
1068  G4DynamicParticle* DynamicParticle = new G4DynamicParticle(ReactionProduct->GetDefinition(), ReactionProduct->GetMomentum());
1069 
1071  return DynamicParticle;
1072 }
1073 
1077 {
1079 
1080  // Break Product down into its A and Z components
1081  G4int A = Isotope % 1000;
1082  G4int Z = (Isotope - A) / 1000;
1083 
1084  // Create the unique identifying name for the particle
1085  std::ostringstream IsotopeName;
1086 
1087  IsotopeName << Z << "_" << A;
1088 
1089  // If it is metastable then append "m" to the name
1090  if(MetaState != G4FFGEnumerations::GROUND_STATE)
1091  {
1092  IsotopeName << "m";
1093 
1094  // If it is a second isomeric state then append "2" to the name
1095  if(MetaState == G4FFGEnumerations::META_2)
1096  {
1097  IsotopeName << "2";
1098  }
1099  }
1100  // Add the name of the element and the extension
1101  IsotopeName << "_" << ElementNames_->theString[Z - 1];
1102 
1104  return IsotopeName.str();
1105 }
1106 
1108 MakeTrees( void )
1109 {
1111 
1112  // Allocate the space
1113  // We will make each tree a binary search
1114  // The maximum number of iterations required to find a single fission product
1115  // based on it's probability is defined by the following:
1116  // x = number of fission products
1117  // Trees = T(x) = ceil( ln(x) )
1118  // Rows/Tree = R(x) = ceil(( sqrt( (8 * x / T(x)) + 1) - 1) / 2)
1119  // Maximum = M(x) = T(x) + R(x)
1120  // Results: x => M(x)
1121  // 10 5
1122  // 100 10
1123  // 1000 25
1124  // 10000 54
1125  // 100000 140
1128 
1129  // Initialize the range of each node
1130  for(G4int i = 0; i < TreeCount_; i++)
1131  {
1133  Trees_[i].Trunk = NULL;
1134  Trees_[i].BranchCount = 0;
1135  Trees_[i].IsEnd = FALSE;
1136  }
1137  // Mark the last tree as the ending tree
1138  Trees_[TreeCount_ - 1].IsEnd = TRUE;
1139 
1141 }
1142 
1145 {
1147 
1148  G4int ProductCount = ENDFData_->G4GetNumberOfFissionProducts();
1149  BranchCount_ = 0;
1151 
1152  // Loop through all the products
1153  for(G4int i = 0; i < ProductCount; i++)
1154  {
1155  // Acquire the data and sort it
1157  }
1158 
1159  // Generate the true normalization factor, since round-off errors may result
1160  // in non-singular normalization of the data files. Also, reset DataTotal_
1161  // since it is used by Renormalize() to set the probability segments.
1164 
1165  // Go through all the trees one at a time
1166  for(G4int i = 0; i < TreeCount_; i++)
1167  {
1168  Renormalize(Trees_[i].Trunk);
1169  // Set the max range of the tree to DataTotal
1170  G4ArrayOps::Copy(YieldEnergyGroups_, Trees_[i].ProbabilityRangeEnd, DataTotal_);
1171  }
1172 
1174 }
1175 
1178 {
1180 
1181  // Check to see if Branch exists. Branch will be a null pointer if it
1182  // doesn't exist
1183  if(Branch != NULL)
1184  {
1185  // Call the lower branch to set the probability segment first, since it
1186  // supposed to have a lower probability segment that this node
1187  Renormalize(Branch->Left);
1188 
1189  // Set this node as the next sequential probability segment
1194 
1195  // Now call the upper branch to set those probability segments
1196  Renormalize(Branch->Right);
1197  }
1198 
1200 }
1201 
1203 SampleAlphaEnergies( std::vector< G4ReactionProduct* >* Alphas )
1204 {
1206 
1207  // The condition of sampling more energy from the fission products than is
1208  // alloted is statistically unfavorable, but it could still happen. The
1209  // do-while loop prevents such an occurrence from happening
1210  G4double MeanAlphaEnergy = 16.0;
1211  G4double TotalAlphaEnergy;
1212 
1213  do
1214  {
1215  G4double AlphaEnergy;
1216  TotalAlphaEnergy = 0;
1217 
1218  // Walk through the alpha particles one at a time and sample each's
1219  // energy
1220  for(unsigned int i = 0; i < Alphas->size(); i++)
1221  {
1222  AlphaEnergy = RandomEngine_->G4SampleGaussian(MeanAlphaEnergy,
1223  2.35,
1225  // Assign the energy to the alpha particle
1226  Alphas->at(i)->SetKineticEnergy(AlphaEnergy);
1227 
1228  // Add up the total amount of kinetic energy consumed.
1229  TotalAlphaEnergy += AlphaEnergy;
1230  }
1231 
1232  // If true, decrement the mean alpha energy by 0.1 and try again.
1233  MeanAlphaEnergy -= 0.1;
1234  } while(TotalAlphaEnergy >= RemainingEnergy_);
1235 
1236  // Subtract the total amount of energy that was assigned.
1237  RemainingEnergy_ -= TotalAlphaEnergy;
1238 
1240 }
1241 
1243 SampleGammaEnergies( std::vector< G4ReactionProduct* >* Gammas )
1244 {
1246 
1247  // Make sure that there is energy to assign to the gamma rays
1248  if(RemainingEnergy_ != 0)
1249  {
1250  G4double SampleEnergy;
1251 
1252  // Sample from RemainingEnergy until it is all gone. Also,
1253  // RemainingEnergy should not be smaller than
1254  // G4FFGDefaultValues::MeanGammaEnergy. This will prevent the
1255  // sampling of a fractional portion of the Gaussian distribution
1256  // in an attempt to find a new gamma ray energy.
1257  while(RemainingEnergy_ >= G4FFGDefaultValues::MeanGammaEnergy )
1258  {
1259  SampleEnergy = RandomEngine_->
1260  G4SampleGaussian(G4FFGDefaultValues::MeanGammaEnergy, 1.0 * MeV, G4FFGEnumerations::POSITIVE);
1261  // Make sure that we didn't sample more energy than was available
1262  if(SampleEnergy <= RemainingEnergy_)
1263  {
1264  // If this energy assignment would leave less energy than the
1265  // 'intrinsic' minimal energy of a gamma ray then just assign
1266  // all of the remaining energy
1267  if(RemainingEnergy_ - SampleEnergy < 100 * keV)
1268  {
1269  SampleEnergy = RemainingEnergy_;
1270  }
1271 
1272  // Create the new particle
1273  Gammas->push_back(new G4ReactionProduct());
1274 
1275  // Set the properties
1276  Gammas->back()->SetDefinition(GammaDefinition_);
1277  Gammas->back()->SetTotalEnergy(SampleEnergy);
1278 
1279  // Calculate how much is left
1280  RemainingEnergy_ -= SampleEnergy;
1281  }
1282  }
1283 
1284  // If there is anything left over, the energy must be above 100 keV but
1285  // less than G4FFGDefaultValues::MeanGammaEnergy. Arbitrarily assign
1286  // RemainingEnergy to a new particle
1287  if(RemainingEnergy_ > 0)
1288  {
1289  SampleEnergy = RemainingEnergy_;
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 
1302 }
1303 
1305 SampleNeutronEnergies( std::vector< G4ReactionProduct* >* Neutrons )
1306 {
1308 
1309  // The condition of sampling more energy from the fission products than is
1310  // alloted is statistically unfavorable, but it could still happen. The
1311  // do-while loop prevents such an occurrence from happening
1312  G4double TotalNeutronEnergy;
1313  G4double NeutronEnergy;
1314 
1315  // Make sure that we don't sample more energy than is available
1316  do
1317  {
1318  TotalNeutronEnergy = 0;
1319 
1320  // Walk through the neutrons one at a time and sample the energies.
1321  // The gamma rays have not yet been sampled, so the last neutron will
1322  // have a NULL value for NextFragment
1323  for(unsigned int i = 0; i < Neutrons->size(); i++)
1324  {
1325  // Assign the energy to the neutron
1326  NeutronEnergy = RandomEngine_->G4SampleWatt(Isotope_, Cause_, IncidentEnergy_);
1327  Neutrons->at(i)->SetKineticEnergy(NeutronEnergy);
1328 
1329  // Add up the total amount of kinetic energy consumed.
1330  TotalNeutronEnergy +=NeutronEnergy;
1331  }
1332  } while (TotalNeutronEnergy > RemainingEnergy_);
1333 
1334  // Subtract the total amount of energy that was assigned.
1335  RemainingEnergy_ -= TotalNeutronEnergy;
1336 
1338 }
1339 
1341 SetNubar( void )
1342 {
1344 
1345  G4int* WhichNubar;
1346  G4int* NubarWidth;
1347  G4double XFactor, BFactor;
1348 
1350  {
1351  WhichNubar = const_cast<G4int*>(&SpontaneousNubar_[0][0]);
1352  NubarWidth = const_cast<G4int*>(&SpontaneousNubarWidth_[0][0]);
1353  } else
1354  {
1355  WhichNubar = const_cast<G4int*>(&NeutronInducedNubar_[0][0]);
1356  NubarWidth = const_cast<G4int*>(&NeutronInducedNubarWidth_[0][0]);
1357  }
1358 
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)
1364  {
1365  if(*WhichNubar == Isotope_)
1366  {
1367  Nubar_ = *(WhichNubar + 1) * IncidentEnergy_ * XFactor
1368  + *(WhichNubar + 2) * BFactor;
1369 
1370  break;
1371  }
1372  WhichNubar += 3;
1373  }
1374 
1375  XFactor = std::pow((G4double)10, -6);
1376  NubarWidth_ = *(NubarWidth + 1) * XFactor;
1377  while(*WhichNubar != -1)
1378  {
1379  if(*WhichNubar == Isotope_)
1380  {
1381  NubarWidth_ = *(NubarWidth + 1) * XFactor;
1382 
1383  break;
1384  }
1385  WhichNubar += 2;
1386  }
1387 
1389 }
1390 
1393 {
1395 
1396  // Initialize the new branch
1397  ProbabilityBranch* NewBranch = new ProbabilityBranch;
1399  NewBranch->Left = NULL;
1400  NewBranch->Right = NULL;
1401  NewBranch->Particle = GetParticleDefinition(YieldData->GetProduct(), YieldData->GetMetaState());
1402  NewBranch->IncidentEnergies = new G4double[YieldEnergyGroups_];
1403  NewBranch->ProbabilityRangeTop = new G4double[YieldEnergyGroups_];
1404  NewBranch->ProbabilityRangeBottom = new G4double[YieldEnergyGroups_];
1405  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->ProbabilityRangeTop, YieldData->GetYieldProbability());
1406  G4ArrayOps::Copy(YieldEnergyGroups_, NewBranch->IncidentEnergies, YieldEnergies_);
1408 
1409  // Check to see if the this is the smallest/largest particle. First, check
1410  // to see if this is the first particle in the system
1411  if(SmallestZ_ == NULL)
1412  {
1413  SmallestZ_ = SmallestA_ = LargestZ_ = LargestA_ = NewBranch->Particle;
1414  } else
1415  {
1416  G4bool IsSmallerZ = NewBranch->Particle->GetAtomicNumber() < SmallestZ_->GetAtomicNumber();
1417  G4bool IsSmallerA = NewBranch->Particle->GetAtomicMass() < SmallestA_->GetAtomicMass();
1418  G4bool IsLargerZ = NewBranch->Particle->GetAtomicNumber() > LargestZ_->GetAtomicNumber();
1419  G4bool IsLargerA = NewBranch->Particle->GetAtomicMass() > LargestA_->GetAtomicMass();
1420 
1421  if(IsSmallerZ)
1422  {
1423  SmallestZ_ = NewBranch->Particle;
1424  }
1425 
1426  if(IsLargerZ)
1427  {
1428  LargestA_ = NewBranch->Particle;
1429  }
1430 
1431  if(IsSmallerA)
1432  {
1433  SmallestA_ = NewBranch->Particle;
1434  }
1435 
1436  if(IsLargerA)
1437  {
1438  LargestA_ = NewBranch->Particle;
1439  }
1440  }
1441 
1442  // Place the new branch
1443  // Determine which tree the new branch goes into
1444  G4int WhichTree = (G4int)floor((G4double)(BranchCount_ % TreeCount_));
1445  ProbabilityBranch** WhichBranch = &(Trees_[WhichTree].Trunk);
1446  Trees_[WhichTree].BranchCount++;
1447 
1448  // Search for the position
1449  // Determine where the branch goes
1450  G4int BranchPosition = (G4int)floor((G4double)(BranchCount_ / TreeCount_)) + 1;
1451 
1452  // Run through the tree until the end branch is reached
1453  while(BranchPosition > 1)
1454  {
1455  if(BranchPosition & 1)
1456  {
1457  // If the 1's bit is on then move to the next 'right' branch
1458  WhichBranch = &((*WhichBranch)->Right);
1459  } else
1460  {
1461  // If the 1's bit is off then move to the next 'down' branch
1462  WhichBranch = &((*WhichBranch)->Left);
1463  }
1464 
1465  BranchPosition >>= 1;
1466  }
1467 
1468  *WhichBranch = NewBranch;
1469  BranchCount_++;
1470 
1472 }
1473 
1476 {
1478 
1479  // Burn each tree, one by one
1480  G4int WhichTree = 0;
1481  while(Trees_[WhichTree].IsEnd != TRUE)
1482  {
1483  BurnTree(Trees_[WhichTree].Trunk);
1484  delete Trees_[WhichTree].Trunk;
1485  delete[] Trees_[WhichTree].ProbabilityRangeEnd;
1486  WhichTree++;
1487  }
1488 
1489  // Delete each dynamically allocated variable
1490  delete ENDFData_;
1491  delete[] Trees_;
1492  delete[] DataTotal_;
1493  delete[] MaintainNormalizedData_;
1494  delete ElementNames_;
1495  delete RandomEngine_;
1496 
1498 }
1499 
1502 {
1504 
1505  // Check to see it Branch exists. Branch will be a null pointer if it
1506  // doesn't exist
1507  if(Branch)
1508  {
1509  // Burn down before you burn up
1510  BurnTree(Branch->Left);
1511  delete Branch->Left;
1512  BurnTree(Branch->Right);
1513  delete Branch->Right;
1514 
1515  delete[] Branch->IncidentEnergies;
1516  delete[] Branch->ProbabilityRangeTop;
1517  delete[] Branch->ProbabilityRangeBottom;
1518  }
1519 
1521 }
1522 
#define G4FFG_RECURSIVE_FUNCTIONLEAVE__
void set(double x, double y, double z)
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)
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 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:449
static G4Alpha * Definition()
Definition: G4Alpha.cc:49
ProbabilityBranch * Right
void setR(double s)
G4DynamicParticleVector * G4GetFission(void)
double getY() const
void G4SetAlphaProduction(G4double WhatAlphaProduction)
G4String MakeFileName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
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
G4ParticleDefinition * GetDefinition() const
double getX() 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
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:72
#define FALSE
Definition: globals.hh:52
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
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)
void DeleteVectorOfPointers(std::vector< T > &Vector)
Definition: G4ArrayOps.hh:216
void setRThetaPhi(double r, double theta, double phi)
G4int GetAtomicMass() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:80
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 G4Neutron * Definition()
Definition: G4Neutron.cc:54
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
#define G4FFG_SPACING__
void SampleNeutronEnergies(std::vector< G4ReactionProduct * > *Neutrons)
#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
#define G4FFG_FUNCTIONENTER__
G4String MakeIsotopeName(G4int Isotope, G4FFGEnumerations::MetaState MetaState)
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)