G4FissionFragmentGenerator.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
/*
 * File:   G4FissionFragmentGenerator.cc
 * Author: B. Wendt (wendbryc@isu.edu)
 *
 * Created on May 11, 2011, 12:04 PM
 */
 
#include "G4FFGDebuggingMacros.hh"
#include "G4FFGDefaultValues.hh"
#include "G4HadFinalState.hh"
#include "G4Ions.hh"
#include "G4Neutron.hh"
#include "globals.hh"
 
#include <iostream>
#include <vector>
// Use a few select constant of CLHEP namespace
using CLHEP::eV;
using CLHEP::GeV;
using CLHEP::keV;
using CLHEP::MeV;
#include "G4FFGEnumerations.hh"
#include "G4FPYBiasedLightFragmentDist.hh"
#include "G4FPYNormalFragmentDist.hh"
#include "G4FissionFragmentGenerator.hh"
#include "G4TableTemplate.hh"
 
G4FissionFragmentGenerator::G4FissionFragmentGenerator()
{
  // Set the default verbosity
  Verbosity_ = G4FFGDefaultValues::Verbosity;
 
  // Initialize the class
  Initialize();
}
 
G4FissionFragmentGenerator::G4FissionFragmentGenerator(G4int Verbosity)
{
  // Set the verbosity
  Verbosity_ = Verbosity;
 
  // Initialize the class
  Initialize();
}
 
void G4FissionFragmentGenerator::Initialize()
{
  G4FFG_FUNCTIONENTER__
 
  // Initialize the class descriptor variables to the default values. These
  // will be used unless the user redefines them.
  Isotope_ = G4FFGDefaultValues::Isotope;
  MetaState_ = G4FFGDefaultValues::MetaState;
  Cause_ = G4FFGDefaultValues::FissionCause;
  IncidentEnergy_ = G4FFGDefaultValues::ThermalNeutronEnergy;
  YieldType_ = G4FFGDefaultValues::YieldType;
  TernaryProbability_ = G4FFGDefaultValues::TernaryProbability;
  AlphaProduction_ = G4FFGDefaultValues::AlphaProduction;
  SamplingScheme_ = G4FFGDefaultValues::SamplingScheme;
 
  // No data class has been created yet
  YieldData_ = nullptr;
  IsReconstructionNeeded_ = TRUE;
 
  G4FFG_FUNCTIONLEAVE__
}
 
G4DynamicParticleVector* G4FissionFragmentGenerator::G4GenerateFission()
{
  G4FFG_FUNCTIONENTER__
 
  const G4HadProjectile Projectile(G4DynamicParticle(G4Neutron::Neutron(), G4ThreeVector(0, 0, 0),
                                                     G4FFGDefaultValues::ThermalNeutronEnergy));
 
  // Call the overloaded function and generate 1 fission
  std::vector<G4DynamicParticleVector*> FissionEvent = G4GenerateFission(1, Projectile);
  G4DynamicParticleVector* Container = FissionEvent[0];
 
  G4FFG_FUNCTIONLEAVE__
  return Container;
}
 
G4DynamicParticleVector*
G4FissionFragmentGenerator::G4GenerateFission(const G4HadProjectile& Projectile)
{
  G4FFG_FUNCTIONENTER__
 
  // Call the overloaded function and generate 1 fission
  std::vector<G4DynamicParticleVector*> FissionEvent = G4GenerateFission(1, Projectile);
  G4DynamicParticleVector* const Container = FissionEvent[0];
 
  G4FFG_FUNCTIONLEAVE__
  return Container;
}
 
const std::vector<G4DynamicParticleVector*>
G4FissionFragmentGenerator::G4GenerateFission(G4long NumberOfFissions,
                                              const G4HadProjectile& Projectile)
{
  G4FFG_FUNCTIONENTER__
 
  // TK Modified 131107
  // std::vector< G4DynamicParticleVector* > FissionEvents(NumberOfFissions);
  std::vector<G4DynamicParticleVector*> FissionEvents(0);
 
  if (Projectile.GetDefinition() == G4Neutron::Neutron()) {
    if (static_cast<int>(IsReconstructionNeeded_) == TRUE) {
      // TODO Eliminate potential need for restructuring during run phase
      // InitializeFissionProductYieldClass();
    }
 
    if (nullptr != YieldData_) {
      for (G4long i = 0; i < NumberOfFissions; i++) {
    FissionEvents.push_back(YieldData_->G4GetFission());
    // FIXME Use particle momentum in balance equation
    // FissionEvents.push_back(YieldData_->G4GetFission(Projectile.Get4Momentum()));
      }
    }
  }
  else {
    FissionEvents.push_back(nullptr);
  }
 
  G4FFG_FUNCTIONLEAVE__
  return FissionEvents;
}
 
G4Ions* G4FissionFragmentGenerator::G4GenerateFissionProduct()
{
  G4FFG_FUNCTIONENTER__
 
  if (static_cast<int>(IsReconstructionNeeded_) == TRUE) {
    // TODO Eliminate potential need for restructuring during run phase
    // InitializeFissionProductYieldClass();
  }
 
  G4Ions* Product = YieldData_->G4GetFissionProduct();
 
  G4FFG_FUNCTIONLEAVE__
  return Product;
}
 
G4double G4FissionFragmentGenerator::G4GetAlphaProduction()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return AlphaProduction_;
}
 
G4double G4FissionFragmentGenerator::G4GetTernaryProbability()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return TernaryProbability_;
}
 
G4FFGEnumerations::FissionCause G4FissionFragmentGenerator::G4GetCause()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return Cause_;
}
 
G4double G4FissionFragmentGenerator::G4GetIncidentEnergy()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return IncidentEnergy_;
}
 
G4int G4FissionFragmentGenerator::G4GetIsotope()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return Isotope_;
}
 
G4FFGEnumerations::MetaState G4FissionFragmentGenerator::G4GetMetaState()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return MetaState_;
}
 
G4FFGEnumerations::FissionSamplingScheme G4FissionFragmentGenerator::G4GetSamplingScheme()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return SamplingScheme_;
}
 
G4FFGEnumerations::YieldType G4FissionFragmentGenerator::G4GetYieldType()
{
  G4FFG_FUNCTIONENTER__
 
  G4FFG_FUNCTIONLEAVE__
  return YieldType_;
}
 
G4int G4FissionFragmentGenerator::G4MakeIsotopeCode(G4int Z, G4int A, G4int M)
{
  // Sanity check;
  A %= 1000;
  Z %= 1000;
  M %= 10;
 
  return (A + Z * 1000) * 10 + M;
}
 
void G4FissionFragmentGenerator::G4SetAlphaProduction(G4double WhatAlphaProduction)
{
  G4FFG_FUNCTIONENTER__
 
  AlphaProduction_ = WhatAlphaProduction;
  if (YieldData_ != nullptr) {
    YieldData_->G4SetAlphaProduction(AlphaProduction_);
  }
 
  if ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) {
    G4FFG_SPACING__
    G4FFG_LOCATION__
 
    G4cout << " -- Alpha production set to " << AlphaProduction_ << G4endl;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetTernaryProbability(G4double WhatTernaryProbability)
{
  G4FFG_FUNCTIONENTER__
 
  TernaryProbability_ = WhatTernaryProbability;
  if (YieldData_ != nullptr) {
    YieldData_->G4SetTernaryProbability(TernaryProbability_);
  }
 
  if ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) {
    G4FFG_SPACING__
    G4FFG_LOCATION__
 
    G4cout << " -- Ternary fission probability set to " << TernaryProbability_ << G4endl;
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetCause(G4FFGEnumerations::FissionCause WhichCause)
{
  G4FFG_FUNCTIONENTER__
 
  G4bool IsValidCause = (WhichCause == G4FFGEnumerations::SPONTANEOUS
                         || WhichCause == G4FFGEnumerations::NEUTRON_INDUCED);
  G4bool IsSameCause = (Cause_ == WhichCause);
 
  if (!IsSameCause && IsValidCause) {
    Cause_ = WhichCause;
    if (Cause_ == G4FFGEnumerations::SPONTANEOUS) {
      IncidentEnergy_ = 0;
    }
    IsReconstructionNeeded_ = TRUE;
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    G4String CauseString;
    switch (WhichCause) {
      case G4FFGEnumerations::SPONTANEOUS:
        CauseString = "SPONTANEOUS";
        break;
      case G4FFGEnumerations::NEUTRON_INDUCED:
        CauseString = "NEUTRON_INDUCED";
        break;
      case G4FFGEnumerations::PROTON_INDUCED:
        CauseString = "PROTON_INDUCED";
        break;
      case G4FFGEnumerations::GAMMA_INDUCED:
        CauseString = "GAMMA_INDUCED";
        break;
    }
 
    if ((Verbosity_ & G4FFGEnumerations::WARNING) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      if (IsValidCause) {
        if (IsSameCause && YieldData_ != nullptr) {
          G4cout << " -- Already set to use " << CauseString
                 << " as the fission cause. Yield data class will not be reconstructed." << G4endl;
        }
        else if (YieldData_ == nullptr) {
          G4cout << " -- Yield data class not yet constructed. " << CauseString
                 << " will be applied when it is constructed." << G4endl;
        }
      }
      else {
        G4cout << " -- Invalid cause of fission" << G4endl;
      }
    }
 
    if (((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) && IsValidCause) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Fission cause set to " << CauseString << "." << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetIncidentEnergy(G4double WhatIncidentEnergy)
{
  G4FFG_FUNCTIONENTER__
 
  if (Cause_ != G4FFGEnumerations::SPONTANEOUS) {
    IncidentEnergy_ = WhatIncidentEnergy;
    if (YieldData_ != nullptr) {
      YieldData_->G4SetEnergy(IncidentEnergy_);
    }
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    std::ostringstream EnergyString;
    if (IncidentEnergy_ / GeV > 1) {
      EnergyString << IncidentEnergy_ / GeV << " GeV";
    }
    else if (IncidentEnergy_ / MeV > 1) {
      EnergyString << IncidentEnergy_ / MeV << " MeV";
    }
    else if (IncidentEnergy_ / keV > 1) {
      EnergyString << IncidentEnergy_ / keV << " keV";
    }
    else {
      EnergyString << IncidentEnergy_ / eV << " eV";
    }
 
    if (((Verbosity_ & G4FFGEnumerations::ENERGY_INFO) != 0)
        || ((Verbosity_ & G4FFGEnumerations::WARNING) != 0))
    {
      if (Cause_ == G4FFGEnumerations::SPONTANEOUS && IncidentEnergy_ != 0) {
        G4FFG_SPACING__
        G4FFG_LOCATION__
 
        G4cout << " -- Cannot set a non-zero energy for spontaneous fission" << G4endl;
      }
      else if (YieldData_ == nullptr) {
        G4FFG_SPACING__
        G4FFG_LOCATION__
 
        G4cout << " -- Yield data class not yet constructed. " << EnergyString.str()
               << " will be applied when it is constructed." << G4endl;
      }
    }
 
    if (((Verbosity_ & G4FFGEnumerations::ENERGY_INFO) != 0)
        || ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0))
    {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Incident neutron energy set to " << EnergyString.str() << "." << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetIsotope(G4int WhichIsotope)
{
  G4FFG_FUNCTIONENTER__
 
  G4bool IsSameIsotope = (Isotope_ == WhichIsotope);
 
  if (!IsSameIsotope) {
    Isotope_ = WhichIsotope;
    IsReconstructionNeeded_ = TRUE;
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    if ((Verbosity_ & G4FFGEnumerations::WARNING) != 0) {
      if (IsSameIsotope && YieldData_ != nullptr) {
        G4FFG_SPACING__
        G4FFG_LOCATION__
 
        G4cout << " -- Isotope " << Isotope_
               << " already in use. Yield data class will not be reconstructed." << G4endl;
      }
      else if (YieldData_ == nullptr) {
        G4FFG_SPACING__
        G4FFG_LOCATION__
 
        G4cout << " -- Yield data class not yet constructed. The isotope will be set to "
               << Isotope_ << " when it is constructed." << G4endl;
      }
    }
 
    if ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Isotope set to " << Isotope_ << "." << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetMetaState(G4FFGEnumerations::MetaState WhichMetaState)
{
  G4FFG_FUNCTIONENTER__
 
  G4bool IsValidMetaState = (WhichMetaState >= G4FFGEnumerations::MetaStateFirst
                             && WhichMetaState <= G4FFGEnumerations::MetaStateLast);
  G4bool IsSameMetaState = (MetaState_ == WhichMetaState);
 
  if (!IsSameMetaState && IsValidMetaState) {
    MetaState_ = WhichMetaState;
    IsReconstructionNeeded_ = TRUE;
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    G4String MetaName;
    switch (MetaState_) {
      case G4FFGEnumerations::GROUND_STATE:
        MetaName = "GROUND_STATE";
        break;
 
      case G4FFGEnumerations::META_1:
        MetaName = "META_1";
        break;
 
      case G4FFGEnumerations::META_2:
        MetaName = "META_2";
        break;
    }
 
    if ((Verbosity_ & G4FFGEnumerations::WARNING) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      std::ostringstream Temp;
      if (IsValidMetaState) {
        if (IsSameMetaState && YieldData_ != nullptr) {
          G4cout << " -- Already set to use " << MetaName
                 << " as the metastable state. Yield data class will not be reconstructed"
                 << G4endl;
        }
        else if (YieldData_ == nullptr) {
          G4cout << " -- Yield data class not yet constructed. " << MetaName
                 << " will be applied when it is constructed." << G4endl;
        }
      }
      else {
        G4cout << " -- Invalid metastable state." << G4endl;
      }
    }
 
    if (((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) && IsValidMetaState) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Metastable state set to " << MetaName << "." << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator ::G4SetSamplingScheme(
  G4FFGEnumerations::FissionSamplingScheme NewScheme)
{
  G4FFG_FUNCTIONENTER__
 
  G4bool IsValidScheme = (NewScheme >= G4FFGEnumerations::FissionSamplingSchemeFirst
                          && NewScheme <= G4FFGEnumerations::FissionSamplingSchemeLast);
  G4bool IsSameScheme = (NewScheme == SamplingScheme_);
 
  if (!IsSameScheme && IsValidScheme) {
    SamplingScheme_ = NewScheme;
    IsReconstructionNeeded_ = TRUE;
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    G4String SchemeString;
    switch (SamplingScheme_) {
      case G4FFGEnumerations::NORMAL:
        SchemeString = "NORMAL";
        break;
 
      case G4FFGEnumerations::LIGHT_FRAGMENT:
        SchemeString = "LIGHT_FRAGMENT";
        break;
 
      default:
        SchemeString = "UNSUPPORTED";
        break;
    }
 
    if ((Verbosity_ & G4FFGEnumerations::WARNING) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      if (IsValidScheme) {
        if (IsSameScheme && YieldData_ != nullptr) {
          G4cout << " -- Already set to use " << SchemeString
                 << " as the sampling scheme. Yield data class will not be reconstructed."
                 << G4endl;
        }
        else if (YieldData_ == nullptr) {
          G4cout << " -- Yield data class not yet constructed. " << SchemeString
                 << " will be applied when it is constructed." << G4endl;
        }
      }
      else {
        G4cout << " -- Invalid sampling scheme." << G4endl;
      }
    }
 
    if (((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) && IsValidScheme) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Sampling scheme set to " << SchemeString << "." << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetYieldType(G4FFGEnumerations::YieldType WhichYieldType)
{
  G4FFG_FUNCTIONENTER__
 
  G4bool IsValidYieldType = (WhichYieldType == G4FFGEnumerations::INDEPENDENT
                             || WhichYieldType == G4FFGEnumerations::CUMULATIVE);
  G4bool IsSameYieldType = (YieldType_ == WhichYieldType);
 
  if (!IsSameYieldType && IsValidYieldType) {
    YieldType_ = WhichYieldType;
    IsReconstructionNeeded_ = TRUE;
  }
 
  if (Verbosity_ != G4FFGEnumerations::SILENT) {
    G4String YieldString;
    switch ((int)YieldType_) {
      case G4FFGEnumerations::INDEPENDENT:
        YieldString = "INDEPENDENT";
        break;
 
      case G4FFGEnumerations::SPONTANEOUS:
        YieldString = "SPONTANEOUS";
        break;
 
      default:
        YieldString = "UNSUPPORTED";
        break;
    }
 
    if ((Verbosity_ & G4FFGEnumerations::WARNING) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      if (IsValidYieldType) {
        if (IsSameYieldType && YieldData_ != nullptr) {
        }
        else if (YieldData_ == nullptr) {
          G4cout << " -- Yield data class not yet constructed. Yield type " << YieldString
                 << " will be applied when it is constructed." << G4endl;
        }
      }
      else {
        G4cout << " -- Invalid yield type." << G4endl;
      }
    }
 
    if (((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) && IsValidYieldType) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Yield type set to " << YieldString << G4endl;
    }
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
void G4FissionFragmentGenerator::G4SetVerbosity(G4int Verbosity)
{
  G4FFG_FUNCTIONENTER__
 
  Verbosity_ = Verbosity;
 
  if (YieldData_ != nullptr) {
    YieldData_->G4SetVerbosity(Verbosity_);
  }
 
  G4FFG_FUNCTIONLEAVE__
}
 
bool G4FissionFragmentGenerator::InitializeFissionProductYieldClass(std::istringstream& dataStream)
{
  G4FFG_FUNCTIONENTER__
 
  if (YieldData_ != nullptr) {
    delete YieldData_;
 
    if ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Old yield data class deleted." << G4endl;
    }
  }
 
  try {
    if (SamplingScheme_ == G4FFGEnumerations::NORMAL) {
      YieldData_ = new G4FPYNormalFragmentDist(Isotope_, MetaState_, Cause_, YieldType_, Verbosity_,
                                               dataStream);
    }
    else {
      YieldData_ = new G4FPYBiasedLightFragmentDist(Isotope_, MetaState_, Cause_, YieldType_,
                                                    Verbosity_, dataStream);
    }
 
    if (AlphaProduction_ != 0 && TernaryProbability_ != 0) {
      YieldData_->G4SetTernaryProbability(TernaryProbability_);
      YieldData_->G4SetAlphaProduction(AlphaProduction_);
    }
 
    if ((Verbosity_ & G4FFGEnumerations::UPDATES) != 0) {
      G4FFG_SPACING__
      G4FFG_LOCATION__
 
      G4cout << " -- Yield data class constructed with defined values." << G4endl;
    }
  }
  catch (std::exception& e) {
    YieldData_ = nullptr;
  }
 
  IsReconstructionNeeded_ = FALSE;
 
  G4FFG_FUNCTIONLEAVE__
  return YieldData_ != nullptr;
}
 
G4FissionFragmentGenerator::~G4FissionFragmentGenerator()
{
  G4FFG_FUNCTIONENTER__
 
  delete YieldData_;
 
  G4FFG_FUNCTIONLEAVE__
}