G4NeutronHPCaptureFS Class Reference

#include <G4NeutronHPCaptureFS.hh>

Inheritance diagram for G4NeutronHPCaptureFS:

Public Member Functions
	G4NeutronHPCaptureFS ()
	~G4NeutronHPCaptureFS ()
void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType)
G4HadFinalState *	ApplyYourself (const G4HadProjectile &theTrack)
G4NeutronHPFinalState *	New ()
Public Member Functions inherited from G4NeutronHPFinalState
	G4NeutronHPFinalState ()
virtual	~G4NeutronHPFinalState ()
void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType)
G4bool	HasXsec ()
G4bool	HasFSData ()
G4bool	HasAnyData ()
virtual G4double	GetXsec (G4double)
virtual G4NeutronHPVector *	GetXsec ()
void	SetA_Z (G4double anA, G4double aZ, G4int aM=0)
G4double	GetZ ()
G4double	GetN ()
G4int	GetM ()

Additional Inherited Members
Protected Member Functions inherited from G4NeutronHPFinalState
void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4NeutronHPDataUsed used)
void	adjust_final_state (G4LorentzVector)
G4bool	DoNotAdjustFinalState ()
Protected Attributes inherited from G4NeutronHPFinalState
G4bool	hasXsec
G4bool	hasFSData
G4bool	hasAnyData
G4NeutronHPNames	theNames
G4HadFinalState	theResult
G4double	theBaseA
G4double	theBaseZ
G4int	theBaseM
G4int	theNDLDataZ
G4int	theNDLDataA
G4int	theNDLDataM

Detailed Description

Definition at line 40 of file G4NeutronHPCaptureFS.hh.

Constructor & Destructor Documentation

G4NeutronHPCaptureFS::G4NeutronHPCaptureFS ( )

inline

Definition at line 44 of file G4NeutronHPCaptureFS.hh.

References G4NeutronHPFinalState::hasXsec.

Referenced by New().

  {
    hasXsec = false; 
    hasExactMF6 = false;
    targetMass = 0;
  }

G4NeutronHPCaptureFS::~G4NeutronHPCaptureFS ( )

inline

Definition at line 51 of file G4NeutronHPCaptureFS.hh.

  {
  }

Member Function Documentation

G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself ( const G4HadProjectile &  theTrack )

virtual

Reimplemented from G4NeutronHPFinalState.

Definition at line 48 of file G4NeutronHPCaptureFS.cc.

References G4HadFinalState::AddSecondary(), G4PhotonEvaporation::BreakItUp(), G4HadFinalState::Clear(), G4NeutronHPFinalState::DoNotAdjustFinalState(), CLHEP::HepLorentzVector::e(), python.hepunit::eV, G4UniformRand, G4Gamma::Gamma(), G4HadProjectile::Get4Momentum(), G4DynamicParticle::Get4Momentum(), G4Nucleus::GetBiasedThermalNucleus(), G4HadProjectile::GetDefinition(), G4DynamicParticle::GetDefinition(), G4IonTable::GetIon(), G4IonTable::GetIonMass(), G4IonTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4NucleiProperties::GetNuclearMass(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4ParticleDefinition::GetPDGMass(), G4NeutronHPPhotonDist::GetPhotons(), G4HadFinalState::GetSecondary(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4NeutronHPFinalState::HasFSData(), G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag(), python.hepunit::MeV, G4Neutron::Neutron(), python.hepunit::pi, G4NeutronHPEnAngCorrelation::Sample(), G4ReactionProduct::SetDefinition(), G4DynamicParticle::SetDefinition(), G4PhotonEvaporation::SetICM(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), G4DynamicParticle::SetMomentum(), G4NeutronHPEnAngCorrelation::SetNeutron(), G4HadFinalState::SetStatusChange(), G4NeutronHPEnAngCorrelation::SetTarget(), sort(), stopAndKill, G4NeutronHPFinalState::theBaseA, G4NeutronHPFinalState::theBaseZ, G4NeutronHPFinalState::theResult, theTarget, TRUE, python.hepunit::twopi, CLHEP::Hep3Vector::unit(), and CLHEP::HepLorentzVector::vect().

  {

    G4int i;
    theResult.Clear();
// 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);

// go to nucleus rest system
    theNeutron.Lorentz(theNeutron, -1*theTarget);
    eKinetic = theNeutron.GetKineticEnergy();

// dice the photons

    G4ReactionProductVector * thePhotons = 0;
    if ( HasFSData() && !getenv ( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) ) 
    { 
       //TK110430
       if ( hasExactMF6 )
       {
          theMF6FinalState.SetTarget(theTarget);
          theMF6FinalState.SetNeutron(theNeutron);
          thePhotons = theMF6FinalState.Sample( eKinetic );
       }
       else
          thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
    }
    else
    {
      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->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)) );
        if( (*it)->GetMomentum().mag() > 10*MeV ) theOne->SetDefinition( theTable->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0 ) );

        //if ( (*i)->GetExcitationEnergy() > 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() ) ;
        //theOne->SetTotalEnergy( (*i)->GetMomentum().t() - (*i)->GetExcitationEnergy() ); //will be calculated from momentum 
        thePhotons->push_back(theOne);
        delete *it;
      } 
      delete products;
    }

// add them to the final state

    G4int nPhotons = 0;
    if(thePhotons!=0) nPhotons=thePhotons->size();

///*
   if ( DoNotAdjustFinalState() ) {
//Make at least one photon  
//101203 TK
    if ( nPhotons == 0 )
    {
       G4ReactionProduct * theOne = new G4ReactionProduct;
       theOne->SetDefinition( G4Gamma::Gamma() ); 
       G4double theta = pi*G4UniformRand();
       G4double phi = twopi*G4UniformRand();
       G4double sinth = std::sin(theta);
       G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) );
       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 = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass()
       //  - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass();
       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)), theTarget);
    }
    
    // Recoil, if only one gamma
    //if (1==nPhotons)
    if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
    {
       G4DynamicParticle * theOne = new G4DynamicParticle;
       //G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
       //                                 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
       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();

       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);
       theResult.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.AddSecondary(theOne);
      delete thePhotons->operator[](i);
    }
    delete thePhotons; 

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

    if ( residual == false )
    {
       //G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
       //                                 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
       G4int nNonZero = 0;
       G4LorentzVector p_photons(0,0,0,0);
       for ( G4int j = 0 ; j != theResult.GetNumberOfSecondaries() ; j++ )
       {
          p_photons += theResult.GetSecondary(j)->GetParticle()->Get4Momentum();
          // To many 0 momentum photons -> Check PhotonDist 
          if ( theResult.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 G4NeutronHPCaptureFS::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?
             G4double theta = pi*G4UniformRand();
             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*std::cos(theta) );
              
             p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
             G4DynamicParticle * theOne = new G4DynamicParticle;
             theOne->SetDefinition( G4Gamma::Gamma() );
             theOne->SetMomentum( tempVector );
             theResult.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.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.AddSecondary(theOne);

    }
//101203TK END

// clean up the primary neutron
    theResult.SetStatusChange(stopAndKill);
    return &theResult;
  }

void G4NeutronHPCaptureFS::Init ( G4double  A, G4double  Z, G4int  M, G4String &  dirName, G4String &  aFSType  )

virtual

Implements G4NeutronHPFinalState.

Definition at line 309 of file G4NeutronHPCaptureFS.cc.

References G4NeutronHPManager::GetDataStream(), G4NeutronHPManager::GetInstance(), G4NeutronHPNames::GetName(), G4NeutronHPDataUsed::GetName(), G4NeutronHPPhotonDist::GetTargetMass(), G4NeutronHPFinalState::hasAnyData, G4NeutronHPFinalState::hasFSData, G4NeutronHPFinalState::hasXsec, G4NeutronHPEnAngCorrelation::Init(), G4NeutronHPPhotonDist::InitAngular(), G4NeutronHPPhotonDist::InitEnergies(), G4NeutronHPPhotonDist::InitMean(), G4NeutronHPFinalState::SetAZMs(), G4NeutronHPFinalState::theBaseA, and G4NeutronHPFinalState::theBaseZ.

  {

     //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);
   G4NeutronHPManager::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;
    G4NeutronHPDataUsed 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();
   G4NeutronHPManager::GetInstance()->GetDataStream(filename,theData);
    hasFSData = theFinalStatePhotons.InitMean(theData); 
    if(hasFSData)
    {
      targetMass = theFinalStatePhotons.GetTargetMass();
      theFinalStatePhotons.InitAngular(theData); 
      theFinalStatePhotons.InitEnergies(theData); 
    }
    //theData.close();
  }

G4NeutronHPFinalState* G4NeutronHPCaptureFS::New ( )

inlinevirtual

Implements G4NeutronHPFinalState.

Definition at line 57 of file G4NeutronHPCaptureFS.hh.

References G4NeutronHPCaptureFS().

  {
   G4NeutronHPCaptureFS * theNew = new G4NeutronHPCaptureFS;
   return theNew;
  }

---

The documentation for this class was generated from the following files:

• G4NeutronHPCaptureFS.hh
• G4NeutronHPCaptureFS.cc