G4ParticleHPCaptureFS Class Reference

#include <G4ParticleHPCaptureFS.hh>

Inheritance diagram for G4ParticleHPCaptureFS:

Public Member Functions
	G4ParticleHPCaptureFS ()
	~G4ParticleHPCaptureFS ()
void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
G4HadFinalState *	ApplyYourself (const G4HadProjectile &theTrack)
G4ParticleHPFinalState *	New ()
Public Member Functions inherited from G4ParticleHPFinalState
	G4ParticleHPFinalState ()
virtual	~G4ParticleHPFinalState ()
void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
virtual void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)=0
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &)
virtual G4ParticleHPFinalState *	New ()=0
G4bool	HasXsec ()
G4bool	HasFSData ()
G4bool	HasAnyData ()
virtual G4double	GetXsec (G4double)
virtual G4ParticleHPVector *	GetXsec ()
void	SetA_Z (G4double anA, G4double aZ, G4int aM=0)
G4double	GetZ ()
G4double	GetN ()
G4double	GetA ()
G4int	GetM ()
void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
void	SetProjectile (G4ParticleDefinition *projectile)

Additional Inherited Members
Protected Member Functions inherited from G4ParticleHPFinalState
void	adjust_final_state (G4LorentzVector)
G4bool	DoNotAdjustFinalState ()
Protected Attributes inherited from G4ParticleHPFinalState
G4bool	hasXsec
G4bool	hasFSData
G4bool	hasAnyData
G4ParticleHPNames	theNames
G4Cache< G4HadFinalState * >	theResult
G4ParticleDefinition *	theProjectile
G4double	theBaseA
G4double	theBaseZ
G4int	theBaseM
G4int	theNDLDataZ
G4int	theNDLDataA
G4int	theNDLDataM

Detailed Description

Definition at line 42 of file G4ParticleHPCaptureFS.hh.

Constructor & Destructor Documentation

G4ParticleHPCaptureFS()

G4ParticleHPCaptureFS::G4ParticleHPCaptureFS ( )

inline

Definition at line 46 of file G4ParticleHPCaptureFS.hh.

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

Referenced by New().

~G4ParticleHPCaptureFS()

G4ParticleHPCaptureFS::~G4ParticleHPCaptureFS ( )

inline

Definition at line 53 of file G4ParticleHPCaptureFS.hh.

  {
  }

Member Function Documentation

ApplyYourself()

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

virtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 50 of file G4ParticleHPCaptureFS.cc.

  {
 
   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 ( 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 ( std::size_t 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 ( std::size_t 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 G4ParticleHPCaptureFS::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);
    return theResult.Get();
  }

Init()

void G4ParticleHPCaptureFS::Init ( G4double  A, G4double  Z, G4int  M, G4String &  dirName, G4String &  aFSType, G4ParticleDefinition *    )

virtual

Implements G4ParticleHPFinalState.

Definition at line 312 of file G4ParticleHPCaptureFS.cc.

  {
 
     //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();
  }

New()

G4ParticleHPFinalState * G4ParticleHPCaptureFS::New ( )

inlinevirtual

Implements G4ParticleHPFinalState.

Definition at line 59 of file G4ParticleHPCaptureFS.hh.

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

---

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

• G4ParticleHPCaptureFS.hh
• G4ParticleHPCaptureFS.cc