ArParticleHPCaptureFS.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.          *
// ********************************************************************
//
// neutron_hp -- source file
// J.P. Wellisch, Nov-1996
// A prototype of the low energy neutron transport model.
//
// 12-April-06 Enable IC electron emissions T. Koi
// 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION flag
// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
// 101203 Bugzilla/Geant4 Problem 1155 Lack of residual in some case
// 110430 Temporary solution in the case of being MF6 final state in Capture reaction (MT102)
//
// P. Arce, June-2014 Conversion neutron_hp to particle_hp
//
// -- artg4tk includes
#include "artg4tk/lists/ArParticleHPCaptureFS.hh"

#include "Geant4/G4Fragment.hh"
#include "Geant4/G4Gamma.hh"
#include "Geant4/G4IonTable.hh"
#include "Geant4/G4Nucleus.hh"
#include "Geant4/G4ParticleHPDataUsed.hh"
#include "Geant4/G4ParticleHPManager.hh"
#include "Geant4/G4PhotonEvaporation.hh"
#include "Geant4/G4PhysicalConstants.hh"
#include "Geant4/G4ReactionProduct.hh"
#include "Geant4/G4SystemOfUnits.hh"

G4HadFinalState*
ArParticleHPCaptureFS::ApplyYourself(const G4HadProjectile& theTrack)
{

  if (theResult.Get() == NULL)
    theResult.Put(new G4HadFinalState);
  theResult.Get()->Clear();

  G4int i;

  // prepare neutron
  G4double eKinetic = theTrack.GetKineticEnergy();
  const G4HadProjectile* incidentParticle = &theTrack;
  G4ReactionProduct theNeutron(
    const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition()));
  theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect());
  theNeutron.SetKineticEnergy(eKinetic);

  // Prepare target
  G4ReactionProduct theTarget;
  G4Nucleus aNucleus;
  G4double eps = 0.0001;
  if (targetMass < 500 * MeV)
    targetMass = (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA + eps),
                                                     static_cast<G4int>(theBaseZ + eps))) /
                 G4Neutron::Neutron()->GetPDGMass();
  G4ThreeVector neutronVelocity =
    1. / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
  G4double temperature = theTrack.GetMaterial()->GetTemperature();
  theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
  theTarget.SetDefinitionAndUpdateE(
    G4IonTable::GetIonTable()->GetIon(G4int(theBaseZ), G4int(theBaseA), 0.0));

  // Put neutron in nucleus rest system
  theNeutron.Lorentz(theNeutron, theTarget);
  eKinetic = theNeutron.GetKineticEnergy();

  // Sample the photons
  G4ReactionProductVector* thePhotons = 0;
  if (useArCapGamma) {
    thePhotons = theFinalgammas.GetGammas();
  } else {
    if (HasFSData() && !G4ParticleHPManager::GetInstance()->GetUseOnlyPhotoEvaporation()) {
      // NDL has final state data
      if (hasExactMF6) {
        theMF6FinalState.SetTarget(theTarget);
        theMF6FinalState.SetProjectileRP(theNeutron);
        thePhotons = theMF6FinalState.Sample(eKinetic);
      } else {
        thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
      }
      if (thePhotons == NULL) {
        throw G4HadronicException(
          __FILE__, __LINE__, "Final state data for photon is not properly allocated");
      }
    } else {
      // NDL does not have final state data or forced to use PhotoEvaporation model
      G4ThreeVector aCMSMomentum = theNeutron.GetMomentum() + theTarget.GetMomentum();
      G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
      G4Fragment nucleus(static_cast<G4int>(theBaseA + 1), static_cast<G4int>(theBaseZ), p4);
      G4PhotonEvaporation photonEvaporation;
      // T. K. add
      photonEvaporation.SetICM(TRUE);
      G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
      G4FragmentVector::iterator it;
      thePhotons = new G4ReactionProductVector;
      for (it = products->begin(); it != products->end(); it++) {
        G4ReactionProduct* theOne = new G4ReactionProduct;
        // T. K. add
        if ((*it)->GetParticleDefinition() != 0)
          theOne->SetDefinition((*it)->GetParticleDefinition());
        else
          theOne->SetDefinition(G4Gamma::Gamma()); // this definiion will be over writen

        // T. K. comment out below line
        // theOne->SetDefinition( G4Gamma::Gamma() );
        G4IonTable* theTable = G4IonTable::GetIonTable();
        if ((*it)->GetMomentum().mag() > 10 * MeV)
          theOne->SetDefinition(
            theTable->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0));

        if ((*it)->GetExcitationEnergy() > 1.0e-2 * eV) {
          G4double ex = (*it)->GetExcitationEnergy();
          G4ReactionProduct* aPhoton = new G4ReactionProduct;
          aPhoton->SetDefinition(G4Gamma::Gamma());
          aPhoton->SetMomentum((*it)->GetMomentum().vect().unit() * ex);
          // aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum
          thePhotons->push_back(aPhoton);
        }

        theOne->SetMomentum((*it)->GetMomentum().vect() *
                            ((*it)->GetMomentum().t() - (*it)->GetExcitationEnergy()) /
                            (*it)->GetMomentum().t());
        thePhotons->push_back(theOne);
        delete *it;
      }
      delete products;
    }
  }

  // Add them to the final state
  G4int nPhotons = 0;
  nPhotons = thePhotons->size();

  if (DoNotAdjustFinalState()) {
    // Make at least one photon
    // 101203 TK
    if (nPhotons == 0) {
      G4ReactionProduct* theOne = new G4ReactionProduct;
      theOne->SetDefinition(G4Gamma::Gamma());
      // Bug #1745 DHW G4double theta = pi*G4UniformRand();
      G4double costheta = 2. * G4UniformRand() - 1.;
      G4double theta = std::acos(costheta);
      G4double phi = twopi * G4UniformRand();
      G4double sinth = std::sin(theta);
      G4ThreeVector direction(sinth * std::cos(phi), sinth * std::sin(phi), costheta);
      theOne->SetMomentum(direction);
      thePhotons->push_back(theOne);
      nPhotons++; // 0 -> 1
    }
    // One photon case: energy set to Q-value
    // 101203 TK
    // if ( nPhotons == 1 )
    if (nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0) {
      G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();

      G4double Q = G4IonTable::GetIonTable()->GetIonMass(
                     static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0) +
                   G4Neutron::Neutron()->GetPDGMass() -
                   G4IonTable::GetIonTable()->GetIonMass(
                     static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);

      thePhotons->operator[](0)->SetMomentum(Q * direction);
    }
    //
  }

  // back to lab system
  for (i = 0; i < nPhotons; i++) {
    thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1 * theTarget);
  }

  // Recoil, if only one gamma
  // if (1==nPhotons)
  if (nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0) {
    G4DynamicParticle* theOne = new G4DynamicParticle;
    G4ParticleDefinition* aRecoil = G4IonTable::GetIonTable()->GetIon(
      static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);
    theOne->SetDefinition(aRecoil);
    // Now energy;
    // Can be done slightly better @
    G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum() -
                              thePhotons->operator[](0)->GetMomentum();

    // TKDB 140520
    // G4ThreeVector theMomUnit = aMomentum.unit();
    // G4double aKinEnergy =  theTrack.GetKineticEnergy()
    //                      +theTarget.GetKineticEnergy(); // gammas come from Q-value
    // G4double theResMass = aRecoil->GetPDGMass();
    // G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
    // G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
    // G4ThreeVector theMomentum = theAbsMom*theMomUnit;
    // theOne->SetMomentum(theMomentum);

    theOne->SetMomentum(aMomentum);
    theResult.Get()->AddSecondary(theOne);
  }

  // Now fill in the gammas.
  for (i = 0; i < nPhotons; i++) {
    // back to lab system
    G4DynamicParticle* theOne = new G4DynamicParticle;
    theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
    theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
    theResult.Get()->AddSecondary(theOne);
    delete thePhotons->operator[](i);
  }
  delete thePhotons;

  // 101203TK
  G4bool residual = false;
  G4ParticleDefinition* aRecoil = G4IonTable::GetIonTable()->GetIon(
    static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA + 1), 0);
  for (G4int j = 0; j != theResult.Get()->GetNumberOfSecondaries(); j++) {
    if (theResult.Get()->GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil)
      residual = true;
  }

  if (residual == false) {
    G4int nNonZero = 0;
    G4LorentzVector p_photons(0, 0, 0, 0);
    for (G4int j = 0; j != theResult.Get()->GetNumberOfSecondaries(); j++) {
      p_photons += theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum();
      // To many 0 momentum photons -> Check PhotonDist
      if (theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0)
        nNonZero++;
    }

    // Can we include kinetic energy here?
    G4double deltaE = (theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy()) -
                      (p_photons.e() + aRecoil->GetPDGMass());

    // Add photons
    if (nPhotons - nNonZero > 0) {
      // G4cout << "TKDB ArParticleHPCaptureFS::ApplyYourself we will create additional " <<
      // nPhotons - nNonZero << " photons" << G4endl;
      std::vector<G4double> vRand;
      vRand.push_back(0.0);
      for (G4int j = 0; j != nPhotons - nNonZero - 1; j++) {
        vRand.push_back(G4UniformRand());
      }
      vRand.push_back(1.0);
      std::sort(vRand.begin(), vRand.end());

      std::vector<G4double> vEPhoton;
      for (G4int j = 0; j < (G4int)vRand.size() - 1; j++) {
        vEPhoton.push_back(deltaE * (vRand[j + 1] - vRand[j]));
      }
      std::sort(vEPhoton.begin(), vEPhoton.end());

      for (G4int j = 0; j < nPhotons - nNonZero - 1; j++) {
        // Isotopic in LAB OK?
        // Bug # 1745 DHW G4double theta = pi*G4UniformRand();
        G4double costheta = 2. * G4UniformRand() - 1.;
        G4double theta = std::acos(costheta);
        G4double phi = twopi * G4UniformRand();
        G4double sinth = std::sin(theta);
        G4double en = vEPhoton[j];
        G4ThreeVector tempVector(
          en * sinth * std::cos(phi), en * sinth * std::sin(phi), en * costheta);

        p_photons += G4LorentzVector(tempVector, tempVector.mag());
        G4DynamicParticle* theOne = new G4DynamicParticle;
        theOne->SetDefinition(G4Gamma::Gamma());
        theOne->SetMomentum(tempVector);
        theResult.Get()->AddSecondary(theOne);
      }

      //        Add last photon
      G4DynamicParticle* theOne = new G4DynamicParticle;
      theOne->SetDefinition(G4Gamma::Gamma());
      //        For better momentum conservation
      G4ThreeVector lastPhoton = -p_photons.vect().unit() * vEPhoton.back();
      p_photons += G4LorentzVector(lastPhoton, lastPhoton.mag());
      theOne->SetMomentum(lastPhoton);
      theResult.Get()->AddSecondary(theOne);
    }

    // Add residual
    G4DynamicParticle* theOne = new G4DynamicParticle;
    G4ThreeVector aMomentum =
      theTrack.Get4Momentum().vect() + theTarget.GetMomentum() - p_photons.vect();
    theOne->SetDefinition(aRecoil);
    theOne->SetMomentum(aMomentum);
    theResult.Get()->AddSecondary(theOne);
  }
  // 101203TK END

  // clean up the primary neutron
  theResult.Get()->SetStatusChange(stopAndKill);

  //<--for(int i=0;i<100;i++){
  //<--    std::cout<<"I am ArParticleHPCaptureFS"<<std::endl;
  //<--}

  return theResult.Get();
}

#include <sstream>
void
ArParticleHPCaptureFS::Init(G4double A,
                            G4double Z,
                            G4int M,
                            G4String& dirName,
                            G4String&,
                            G4ParticleDefinition*)
{

  // TK110430 BEGIN
  std::stringstream ss;
  ss << static_cast<G4int>(Z);
  G4String sZ;
  ss >> sZ;
  ss.clear();
  ss << static_cast<G4int>(A);
  G4String sA;
  ss >> sA;

  ss.clear();
  G4String sM;
  if (M > 0) {
    ss << "m";
    ss << M;
    ss >> sM;
    ss.clear();
  }

  G4String element_name = theNames.GetName(static_cast<G4int>(Z) - 1);
  G4String filenameMF6 = dirName + "/FSMF6/" + sZ + "_" + sA + sM + "_" + element_name;
  // std::ifstream dummyIFS(filenameMF6, std::ios::in);
  // if ( dummyIFS.good() == true ) hasExactMF6=true;
  std::istringstream theData(std::ios::in);
  G4ParticleHPManager::GetInstance()->GetDataStream(filenameMF6, theData);

  // TK110430 Only use MF6MT102 which has exactly same A and Z
  // Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0
  if (theData.good() == true) {
    hasExactMF6 = true;
    theMF6FinalState.Init(theData);
    // theData.close();
    return;
  }
  // TK110430 END

  G4String tString = "/FS";
  G4bool dbool;
  G4ParticleHPDataUsed aFile =
    theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);

  G4String filename = aFile.GetName();
  SetAZMs(A, Z, M, aFile);
  // theBaseA = A;
  // theBaseZ = G4int(Z+.5);
  if (!dbool || (Z < 2.5 && (std::abs(theBaseZ - Z) > 0.0001 || std::abs(theBaseA - A) > 0.0001))) {
    hasAnyData = false;
    hasFSData = false;
    hasXsec = false;
    return;
  }
  // std::ifstream theData(filename, std::ios::in);
  // std::istringstream theData(std::ios::in);
  theData.clear();
  G4ParticleHPManager::GetInstance()->GetDataStream(filename, theData);
  hasFSData = theFinalStatePhotons.InitMean(theData);
  if (hasFSData) {
    targetMass = theFinalStatePhotons.GetTargetMass();
    theFinalStatePhotons.InitAngular(theData);
    theFinalStatePhotons.InitEnergies(theData);
  }
  // theData.close();
}