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