G4ParticleHPInelasticCompFS Class Reference

#include <G4ParticleHPInelasticCompFS.hh>

Inheritance diagram for G4ParticleHPInelasticCompFS:

Public Member Functions
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &theTrack)=0
void	CompositeApply (const G4HadProjectile &theTrack, G4ParticleDefinition *aHadron)
	G4ParticleHPInelasticCompFS ()
G4double	GetA ()
G4int	GetM ()
G4double	GetN ()
virtual G4ParticleHPVector *	GetXsec ()
virtual G4double	GetXsec (G4double anEnergy)
G4double	GetZ ()
G4bool	HasAnyData ()
G4bool	HasFSData ()
G4bool	HasXsec ()
void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aSFType, G4ParticleDefinition *)
void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
void	InitDistributionInitialState (G4ReactionProduct &anIncidentPart, G4ReactionProduct &aTarget, G4int it)
void	InitGammas (G4double AR, G4double ZR)
virtual G4ParticleHPFinalState *	New ()=0
G4int	SelectExitChannel (G4double eKinetic)
void	SetA_Z (G4double anA, G4double aZ, G4int aM=0)
void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
void	SetProjectile (G4ParticleDefinition *projectile)
virtual	~G4ParticleHPInelasticCompFS ()

Protected Member Functions
void	adjust_final_state (G4LorentzVector)

Protected Attributes
G4String	gammaPath
G4bool	hasAnyData
G4bool	hasFSData
G4bool	hasXsec
std::vector< G4int >	LR
std::vector< G4double >	QI
G4int	secID
G4ParticleHPAngular *	theAngularDistribution [51]
G4double	theBaseA
G4int	theBaseM
G4double	theBaseZ
G4ParticleHPEnAngCorrelation *	theEnergyAngData [51]
G4ParticleHPEnergyDistribution *	theEnergyDistribution [51]
G4ParticleHPPhotonDist *	theFinalStatePhotons [51]
G4ParticleHPDeExGammas	theGammas
G4ParticleHPNames	theNames
G4int	theNDLDataA
G4int	theNDLDataM
G4int	theNDLDataZ
G4ParticleDefinition *	theProjectile
G4Cache< G4HadFinalState * >	theResult
G4ParticleHPVector *	theXsection [51]

Private Member Functions
void	two_body_reaction (G4ReactionProduct *proj, G4ReactionProduct *targ, G4ReactionProduct *product, G4double resExcitationEnergy)
G4bool	use_nresp71_model (const G4ParticleDefinition *aDefinition, const G4int it, const G4ReactionProduct &theTarget, G4ReactionProduct &boosted)

Private Attributes
G4NRESP71M03	nresp71_model

Detailed Description

Definition at line 50 of file G4ParticleHPInelasticCompFS.hh.

Constructor & Destructor Documentation

G4ParticleHPInelasticCompFS()

G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS ( )

inline

Definition at line 54 of file G4ParticleHPInelasticCompFS.hh.

    {
      QI.resize(51);
      LR.resize(51);
      for(G4int i=0; i<51; i++) {
        hasXsec = true; 
        theXsection[i] = 0;
        theEnergyDistribution[i] = 0;
        theAngularDistribution[i] = 0;
        theEnergyAngData[i] = 0;
        theFinalStatePhotons[i] = 0;
        QI[i] = 0.0;
        LR[i] = 0;
      }
    }

References G4ParticleHPFinalState::hasXsec, LR, QI, theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

~G4ParticleHPInelasticCompFS()

virtual G4ParticleHPInelasticCompFS::~G4ParticleHPInelasticCompFS ( )

inlinevirtual

Definition at line 70 of file G4ParticleHPInelasticCompFS.hh.

    {
      for(G4int i=0; i<51; i++) {
        if (theXsection[i] != 0) delete theXsection[i];
        if (theEnergyDistribution[i] != 0) delete theEnergyDistribution[i];
        if (theAngularDistribution[i] != 0) delete theAngularDistribution[i];
        if (theEnergyAngData[i] != 0) delete theEnergyAngData[i];
        if (theFinalStatePhotons[i] != 0) delete theFinalStatePhotons[i];
      }
    }

References theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, and theXsection.

Member Function Documentation

adjust_final_state()

void G4ParticleHPFinalState::adjust_final_state ( G4LorentzVector  init_4p_lab )

protectedinherited

Definition at line 47 of file G4ParticleHPFinalState.cc.

{
 
   G4double minimum_energy = 1*keV;
 
   if ( G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState() ) return;
 
   G4int nSecondaries = theResult.Get()->GetNumberOfSecondaries();
 
   G4int sum_Z = 0;
   G4int sum_A = 0;
   G4int max_SecZ = 0;
   G4int max_SecA = 0;
   G4int imaxA = -1;
   for ( int i = 0 ; i < nSecondaries ; i++ )
   {
      //G4cout << "G4ParticleHPFinalState::adjust_final_state theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() = " << theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() << G4endl;
      sum_Z += theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber();
      max_SecZ = std::max ( max_SecZ , theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
      sum_A += theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass();
      max_SecA = std::max ( max_SecA , theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
      if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))    G4cout << "G4ParticleHPFinalState::adjust_final_stat SECO " << i << " " <<theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() << G4endl;
#endif
 
   }
 
   G4ParticleDefinition* resi_pd = 0;
 
   G4double baseZNew = theBaseZ;
   G4double baseANew = theBaseA;
   if( theProjectile == G4Neutron::Neutron() ) {
     baseANew ++;
   } else if( theProjectile == G4Proton::Proton() ) {
     baseZNew ++;
     baseANew ++;
   } else if( theProjectile == G4Deuteron::Deuteron() ) {
     baseZNew ++;
     baseANew += 2;
   } else if( theProjectile == G4Triton::Triton() ) {
     baseZNew ++;
     baseANew += 3;
   } else if( theProjectile == G4He3::He3() ) {
     baseZNew += 2;
     baseANew += 3;
   } else if( theProjectile == G4Alpha::Alpha() ) {
     baseZNew += 2;
     baseANew += 4;
   }
 
#ifdef G4PHPDEBUG
   if( std::getenv("G4ParticleHPDebug")) G4cout  << "G4ParticleHPFinalState::adjust_final_stat  BaseZ " << baseZNew << " BaseA " << baseANew << " sum_Z " << sum_Z << " sum_A " << sum_A << G4endl;
#endif
 
   G4bool needOneMoreSec = false;
   G4ParticleDefinition* oneMoreSec_pd = 0;
   if ( (int)(baseZNew - sum_Z) == 0 && (int)(baseANew - sum_A) == 0 )
   {
      //All secondaries are already created;
      resi_pd = theResult.Get()->GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
   }
   else
   {
      if ( max_SecA >= int(baseANew - sum_A) ) 
      {
         //Most heavy secondary is interpreted as residual
         resi_pd = theResult.Get()->GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
         needOneMoreSec = true;
      }
      else 
      {
         //creation of residual is requierd 
    resi_pd = G4IonTable::GetIonTable()->GetIon ( int(baseZNew - sum_Z) , (int)(baseANew - sum_A) , 0.0 );
      }
 
      if ( needOneMoreSec )
      {
         if ( int(baseZNew - sum_Z) == 0 && (int)(baseANew - sum_A) > 0 )
         {
            //In this case, one neutron is added to secondaries 
            if ( int(baseANew - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
            oneMoreSec_pd = G4Neutron::Neutron();
         }
         else 
         {
#ifdef G4PHPDEBUG
       if( std::getenv("G4ParticleHPDebug")) G4cout << this << "G4ParticleHPFinalState oneMoreSec_pd Z " << baseZNew << " - " << sum_Z << " A " << baseANew << " - " << sum_A << " projectile " << theProjectile->GetParticleName() << G4endl;
#endif
       oneMoreSec_pd = G4IonTable::GetIonTable()->GetIon ( int(baseZNew - sum_Z) , (int)(baseANew - sum_A) , 0.0 );
       if( !oneMoreSec_pd ) {
         G4cerr << this << "G4ParticleHPFinalState oneMoreSec_pd Z " << baseZNew << " - " << sum_Z << " A " << baseANew << " - " << sum_A << " projectile " << theProjectile->GetParticleName() << G4endl;
         G4Exception("G4ParticleHPFinalState:adjust_final_state",
             "Warning",
             JustWarning,
             "No adjustment will be done!");
         return;
       }
         }
      }
  
      if ( resi_pd == 0 )
      {
         // theNDLDataZ,A has the Z and A of used NDL file
    G4double ndlZNew = theNDLDataZ;
    G4double ndlANew = theNDLDataA;
    if( theProjectile == G4Neutron::Neutron() ) {
      ndlANew ++;
    } else if( theProjectile == G4Proton::Proton() ) {
      ndlZNew ++;
      ndlANew ++;
    } else if( theProjectile == G4Deuteron::Deuteron() ) {
      ndlZNew ++;
      ndlANew += 2;
    } else if( theProjectile == G4Triton::Triton() ) {
      ndlZNew ++;
      ndlANew += 3;
    } else if( theProjectile == G4He3::He3() ) {
      ndlZNew += 2;
      ndlANew += 3;
    } else if( theProjectile == G4Alpha::Alpha() ) {
      ndlZNew += 2;
      ndlANew += 4;
    }
        // theNDLDataZ,A has the Z and A of used NDL file
        if ( (int)(ndlZNew - sum_Z) == 0 && (int)(ndlANew - sum_A) == 0 )
        {
           G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
           G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
           resi_pd = G4IonTable::GetIonTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
       if( !resi_pd ) {
         G4cerr << "G4ParticleHPFinalState resi_pd Z " << max_SecZ << " - " << dif_Z << " A " << max_SecA << " - " << dif_A << " projectile " << theProjectile->GetParticleName() << G4endl;
         G4Exception("G4ParticleHPFinalState:adjust_final_state",
             "Warning",
             JustWarning,
             "No adjustment will be done!");
         return;
       }
 
           for ( int i = 0 ; i < nSecondaries ; i++ )
           {
              if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ  
                && theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
              {
                 G4ThreeVector p = theResult.Get()->GetSecondary( i )->GetParticle()->GetMomentum();
                 p = p * resi_pd->GetPDGMass()/ G4IonTable::GetIonTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass(); 
                 theResult.Get()->GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
                 theResult.Get()->GetSecondary( i )->GetParticle()->SetMomentum( p );
              } 
           }
        }
      }
   }
 
 
   G4LorentzVector secs_4p_lab( 0.0 );
 
   G4int n_sec = theResult.Get()->GetNumberOfSecondaries();
   G4double fast = 0;
   G4double slow = 1;
   G4int ifast = 0;
   G4int islow = 0;
   G4int ires = -1;
 
   for ( G4int i = 0 ; i < n_sec ; i++ )
   {
 
      //G4cout << "HP_DB " << i 
      //       << " " << theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() 
      //       << " 4p " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum() 
      //       << " ke " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e() - theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetPDGMass() 
      //       << G4endl;
 
      secs_4p_lab += theResult.Get()->GetSecondary( i )->GetParticle()->Get4Momentum();
 
      G4double beta = 0;
      if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition() != G4Gamma::Gamma() )
      {
         beta = theResult.Get()->GetSecondary( i )->GetParticle()->Get4Momentum().beta(); 
      }
      else
      {
         beta = 1;
      }
 
      if ( theResult.Get()->GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i; 
 
      if ( slow > beta && beta != 0 ) 
      {
         slow = beta; 
         islow = i;
      }
 
      if ( fast <= beta ) 
      {
         if ( fast != 1 )
         {
            fast = beta; 
            ifast = i;
         } 
         else
         {
//          fast is already photon then check E
            G4double e = theResult.Get()->GetSecondary( i )->GetParticle()->Get4Momentum().e();
            if ( e > theResult.Get()->GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )   
            {
//             among photons, the highest E becomes the fastest 
               ifast = i; 
            }
         } 
      }
   }
 
 
   G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
 
   //G4cout << "HP_DB dif_4p " << init_4p_lab - secs_4p_lab << G4endl;
   //G4cout << "HP_DB dif_3p mag " << ( dif_4p.v() ).mag() << G4endl;
   //G4cout << "HP_DB dif_e " << dif_4p.e() - ( dif_4p.v() ).mag()<< G4endl;
 
   G4LorentzVector p4(0);
   if ( ires == -1 ) 
   {
//    Create and Add Residual Nucleus 
      ires = nSecondaries;
      nSecondaries += 1;
 
      G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );    
      theResult.Get()->AddSecondary ( res, secID );    
 
      p4 = res->Get4Momentum(); 
      if ( slow > p4.beta() ) 
      {
         slow = p4.beta(); 
         islow = ires;
      }
      dif_4p = init_4p_lab - ( secs_4p_lab + p4 );  
   }
 
   if ( needOneMoreSec && oneMoreSec_pd)
     //
     // fca: this is not a fix, this is a crash avoidance...
     // fca: the baryon number is still wrong, most probably because it
     // fca: should have been decreased, but since we could not create a particle
     // fca: we just do not add it
     //
   {
      nSecondaries += 1;
      G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );    
      theResult.Get()->AddSecondary ( one, secID );    
      p4 = one->Get4Momentum(); 
      if ( slow > p4.beta() ) 
      {
         slow = p4.beta(); 
         islow = nSecondaries-1; //Because the first is 0th, so the last becomes "nSecondaries-1"
      }
      dif_4p = init_4p_lab - ( secs_4p_lab + p4 );  
   }
 
   //Which is bigger dif_p or dif_e
 
   if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
   {
 
      // Adjust p
      //if ( dif_4p.v().mag() < 1*MeV )
      if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
      {
 
         nSecondaries += 1;
         theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ), secID );    
 
      }
      else
      {
         //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
      }
 
   }
   else
   {
 
      // dif_p > dif_e 
      // at first momentum 
      // Move residual momentum
 
      p4 = theResult.Get()->GetSecondary( ires )->GetParticle()->Get4Momentum();
      theResult.Get()->GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() ); 
      dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.Get()->GetSecondary( ires )->GetParticle()->Get4Momentum()  );  
 
      //G4cout << "HP_DB new residual kinetic energy " << theResult.GetSecondary( ires )->GetParticle()->GetKineticEnergy() << G4endl;
 
   }
 
   G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
   //G4cout << "HP_DB dif_e " << dif_e << G4endl;
 
   if ( dif_e > 0 )
   {
 
//    create 2 gamma 
 
      nSecondaries += 2;
      G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
      
      if ( minimum_energy < e1 )  
      {
         G4double costh = 2.*G4UniformRand()-1.;
         G4double phi = twopi*G4UniformRand();
         G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi), 
                            std::sin(std::acos(costh))*std::sin(phi),
                            costh);
         theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ),  secID );    
         theResult.Get()->AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ), secID );    
      }
      else
      {
         //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
      }
 
   }
   else    //dif_e < 0
   {
 
//    At first reduce KE of the fastest secondary;
      G4double ke0 = theResult.Get()->GetSecondary( ifast )->GetParticle()->GetKineticEnergy();
      G4ThreeVector p0 = theResult.Get()->GetSecondary( ifast )->GetParticle()->GetMomentum();
      G4ThreeVector dir = ( theResult.Get()->GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
 
      //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
 
      if ( ke0 + dif_e > 0 )
      {
         theResult.Get()->GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e ); 
         G4ThreeVector dp = p0 - theResult.Get()->GetSecondary( ifast )->GetParticle()->GetMomentum();  
 
         G4ThreeVector p = theResult.Get()->GetSecondary( islow )->GetParticle()->GetMomentum();
         //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir ); 
         theResult.Get()->GetSecondary( islow )->GetParticle()->SetMomentum( p + dp ); 
      }
      else
      {
         //G4cout << "HP_DB Difference in dif_e too large ( <0MeV ) to adjust, so that give up tuning" << G4endl;
      }
 
   }
 
}

References G4HadFinalState::AddSecondary(), G4Alpha::Alpha(), CLHEP::HepLorentzVector::beta(), anonymous_namespace{G4PionRadiativeDecayChannel.cc}::beta, G4Deuteron::Deuteron(), CLHEP::HepLorentzVector::e(), e1, G4cerr, G4cout, G4endl, G4Exception(), G4UniformRand, G4Gamma::Gamma(), G4Cache< VALTYPE >::Get(), G4DynamicParticle::Get4Momentum(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4DynamicParticle::GetDefinition(), G4ParticleHPManager::GetInstance(), G4IonTable::GetIon(), G4IonTable::GetIonTable(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentum(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGMass(), G4HadFinalState::GetSecondary(), G4He3::He3(), JustWarning, keV, CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), MeV, G4Neutron::Neutron(), G4Proton::Proton(), G4ParticleHPFinalState::secID, G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentum(), G4ParticleHPFinalState::theBaseA, G4ParticleHPFinalState::theBaseZ, G4ParticleHPFinalState::theNDLDataA, G4ParticleHPFinalState::theNDLDataZ, G4ParticleHPFinalState::theProjectile, G4ParticleHPFinalState::theResult, G4Triton::Triton(), twopi, and CLHEP::HepLorentzVector::v().

Referenced by G4ParticleHPInelasticBaseFS::BaseApply(), and CompositeApply().

ApplyYourself()

virtual G4HadFinalState * G4ParticleHPInelasticCompFS::ApplyYourself ( const G4HadProjectile &  theTrack )

pure virtual

Reimplemented from G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

CompositeApply()

void G4ParticleHPInelasticCompFS::CompositeApply	(	const G4HadProjectile &	theTrack,
		G4ParticleDefinition *	aHadron
	)

Definition at line 235 of file G4ParticleHPInelasticCompFS.cc.

{
 
// prepare neutron
    if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
    theResult.Get()->Clear();
    G4double eKinetic = theTrack.GetKineticEnergy();
    const G4HadProjectile *hadProjectile = &theTrack;
    G4ReactionProduct incidReactionProduct( const_cast<G4ParticleDefinition *>(hadProjectile->GetDefinition()) ); // incidReactionProduct
    incidReactionProduct.SetMomentum( hadProjectile->Get4Momentum().vect() );
    incidReactionProduct.SetKineticEnergy( eKinetic );
 
// prepare target
    G4int i;
    for(i=0; i<50; i++)
    { if(theXsection[i] != 0) { break; } } 
 
    G4double targetMass=0;
    G4double eps = 0.0001;
    targetMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps));
#ifdef G4PHPDEBUG
    if( std::getenv("G4ParticleHPDebug"))  G4cout <<this <<" G4ParticleHPInelasticCompFS::CompositeApply A " <<theBaseA <<" Z " <<theBaseZ <<" incident " <<hadProjectile->GetDefinition()->GetParticleName() <<G4endl;
#endif
//    if(theEnergyAngData[i]!=0)
//        targetMass = theEnergyAngData[i]->GetTargetMass();
//    else if(theAngularDistribution[i]!=0)
//        targetMass = theAngularDistribution[i]->GetTargetMass();
//    else if(theFinalStatePhotons[50]!=0)
//        targetMass = theFinalStatePhotons[50]->GetTargetMass();
    G4ReactionProduct theTarget; 
    G4Nucleus aNucleus;
    //G4ThreeVector neuVelo = (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum();
    //theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() , neuVelo, theTrack.GetMaterial()->GetTemperature());
    //G4Nucleus::GetBiasedThermalNucleus requests normalization of mass and velocity in neutron mass
    G4ThreeVector neuVelo = ( 1./G4Neutron::Neutron()->GetPDGMass() )*incidReactionProduct.GetMomentum();
    theTarget = aNucleus.GetBiasedThermalNucleus( targetMass/G4Neutron::Neutron()->GetPDGMass()
                                                , neuVelo, theTrack.GetMaterial()->GetTemperature() );
    
    theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) );  //XX
 
// prepare the residual mass
    G4double residualMass=0;
    G4double residualZ = theBaseZ + theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge();
    G4double residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
    residualMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps));
 
// prepare energy in target rest frame
    G4ReactionProduct boosted;
    boosted.Lorentz(incidReactionProduct, theTarget);
    eKinetic = boosted.GetKineticEnergy();
//    G4double momentumInCMS = boosted.GetTotalMomentum();
  
// select exit channel for composite FS class.
    G4int it = SelectExitChannel( eKinetic );
  
    //E. Mendoza (2018) -- to use JENDL/AN-2005
    if(theEnergyDistribution[it]==0 && theAngularDistribution[it]==0 && theEnergyAngData[it]==0){
      if(theEnergyDistribution[50]!=0 || theAngularDistribution[50]!=0 || theEnergyAngData[50]!=0){
        it=50;
      }
    }
 
    // set target and neutron in the relevant exit channel
    InitDistributionInitialState(incidReactionProduct, theTarget, it);    
 
   //---------------------------------------------------------------------//
   //Hook for NRESP71MODEL
   if ( G4ParticleHPManager::GetInstance()->GetUseNRESP71Model() && eKinetic<20*MeV) {
      if ( (G4int)(theBaseZ+0.1) == 6 ) // If the reaction is with Carbon...
      {
         if ( theProjectile == G4Neutron::Definition() ) {
            if ( use_nresp71_model( aDefinition , it , theTarget , boosted ) ) return;
         }
      }
   }
   //---------------------------------------------------------------------//
 
    G4ReactionProductVector * thePhotons = 0;
    G4ReactionProductVector * theParticles = 0;
    G4ReactionProduct aHadron;
    aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@    
    G4double availableEnergy = incidReactionProduct.GetKineticEnergy() + incidReactionProduct.GetMass() - aHadron.GetMass() +
                             (targetMass - residualMass);
//080730c
    if ( availableEnergy < 0 )
    {
       //G4cout << "080730c Adjust availavleEnergy " << G4endl; 
       availableEnergy = 0; 
    }
    G4int nothingWasKnownOnHadron = 0;
    G4int dummy;
    G4double eGamm = 0;
    G4int iLevel=it-1;
 
//  TK without photon has it = 0
    if( 50 == it ) 
    {
 
//    TK Excitation level is not determined
      iLevel=-1;
      aHadron.SetKineticEnergy(availableEnergy*residualMass/
                               (aHadron.GetMass()+residualMass));
 
      //aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*
      //                  std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
      //                            aHadron.GetMass()*aHadron.GetMass()));
 
      //TK add safty 100909
      G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy() - aHadron.GetMass()*aHadron.GetMass() );
      G4double p = 0.0;
      if ( p2 > 0.0 ) p = std::sqrt( p ); 
 
      aHadron.SetMomentum(incidReactionProduct.GetMomentum()*(1./incidReactionProduct.GetTotalMomentum())*p );
 
    }
    else
    {
      while ( iLevel!=-1 && theGammas.GetLevel(iLevel) == 0 ) { iLevel--; } // Loop checking, 11.05.2015, T. Koi
    }
 
 
    if ( theAngularDistribution[it] != 0 ) // MF4
    {
      if(theEnergyDistribution[it]!=0) // MF5
      {
    //************************************************************
    /*
        aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
        G4double eSecN = aHadron.GetKineticEnergy();
    */
    //************************************************************
    //EMendoza --> maximum allowable energy should be taken into account.
        G4double dqi = 0.0;
        if ( QI[it] < 0 || 849 < QI[it] ) dqi = QI[it]; //For backword compatibility QI introduced since G4NDL3.15
    G4double MaxEne=eKinetic+dqi;
    G4double eSecN=0.;
 
        G4int icounter=0;
        G4int icounter_max=1024;
        do {
           icounter++;
           if ( icounter > icounter_max ) {
          G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
              break;
           }
       eSecN=theEnergyDistribution[it]->Sample(eKinetic, dummy);
    }while(eSecN>MaxEne); // Loop checking, 11.05.2015, T. Koi
    aHadron.SetKineticEnergy(eSecN);
    //************************************************************
        eGamm = eKinetic-eSecN;
        for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
        {
          if(theGammas.GetLevelEnergy(iLevel)<eGamm) break;
        }
        G4double random = 2*G4UniformRand();
        iLevel+=G4int(random);
        if(iLevel>theGammas.GetNumberOfLevels()-1)iLevel = theGammas.GetNumberOfLevels()-1;
      }
      else
      {
        G4double eExcitation = 0;
        if(iLevel>=0) eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
        while (eKinetic-eExcitation < 0 && iLevel>0) // Loop checking, 11.05.2015, T. Koi
    {
      iLevel--;
      eExcitation = theGammas.GetLevel(iLevel)->GetLevelEnergy();    
    }
        //110610TK BEGIN
        //Use QI value for calculating excitation energy of residual.
        G4bool useQI=false;
        G4double dqi = QI[it]; 
        if (dqi < 0 || 849 < dqi) useQI = true; // Former libraries do not have values in this range
 
        if (useQI) {
          // QI introudced since G4NDL3.15
          // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
          // eExcitation = QM-QI[it];
          // eExcitation = QI[0] - QI[it];   // Bug fix #1838
          // if(eExcitation < 20*CLHEP::keV) eExcitation = 0;
 
          eExcitation = std::max(0.,QI[0] - QI[it]);  // Bug fix 2333
 
          // Re-evluate iLevel based on this eExcitation
          iLevel = 0;
          G4bool find = false;
          G4int imaxEx = 0;
          G4double level_tolerance = 1.0*CLHEP::keV;
 
          while( theGammas.GetLevel(iLevel+1) != 0 ) { // Loop checking, 11.05.2015, T. Koi
            G4double maxEx = 0.0;
            if (maxEx < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) {
              maxEx = theGammas.GetLevel(iLevel)->GetLevelEnergy();
              imaxEx = iLevel;
            }
 
            // Fix bug 1789 DHW - first if-branch added because gamma data come from ENSDF
            //                    and do not necessarily match the excitations used in ENDF-B.VII
            //                    Compromise solution: use 1 keV tolerance suggested by T. Koi
            if (std::abs(eExcitation - theGammas.GetLevel(iLevel)->GetLevelEnergy() ) < level_tolerance) {
              find = true;
              break;
 
            } else if (eExcitation < theGammas.GetLevel(iLevel)->GetLevelEnergy() ) {
              find = true;
              iLevel--;
              // very small eExcitation, iLevel becomes -1, this is protected below
              if (theTrack.GetDefinition() == aDefinition) { // this line added as part of fix #1838
                if (iLevel == -1) iLevel = 0;
              }
              break;
            }
            iLevel++;
          }
 
          // If proper level cannot be found, use the maximum level
          if (!find) iLevel = imaxEx;
        }
    
    if(std::getenv("G4ParticleHPDebug") && eKinetic-eExcitation < 0) 
    {
      throw G4HadronicException(__FILE__, __LINE__, "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
    }
    if(eKinetic-eExcitation < 0) eExcitation = 0;
    if(iLevel!= -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
    
      }
      theAngularDistribution[it]->SampleAndUpdate(aHadron);
 
      if( theFinalStatePhotons[it] == 0 )
      {
        //G4cout << "110610 USE Gamma Level" << G4endl;
// TK comment Most n,n* eneter to this  
    thePhotons = theGammas.GetDecayGammas(iLevel);
    eGamm -= theGammas.GetLevelEnergy(iLevel);
    if(eGamm>0) // @ ok for now, but really needs an efficient way of correllated sampling @
    {
          G4ReactionProduct * theRestEnergy = new G4ReactionProduct;
          theRestEnergy->SetDefinition(G4Gamma::Gamma());
          theRestEnergy->SetKineticEnergy(eGamm);
          G4double costh = 2.*G4UniformRand()-1.;
          G4double phi = CLHEP::twopi*G4UniformRand();
          theRestEnergy->SetMomentum(eGamm*std::sin(std::acos(costh))*std::cos(phi), 
                                     eGamm*std::sin(std::acos(costh))*std::sin(phi),
                                     eGamm*costh);
          if(thePhotons == 0) { thePhotons = new G4ReactionProductVector; }
          thePhotons->push_back(theRestEnergy);
    }
      }
    }
    else if(theEnergyAngData[it] != 0) // MF6  
    {
 
      theParticles = theEnergyAngData[it]->Sample(eKinetic);
 
      //141017 Fix BEGIN
      //Adjust A and Z in the case of miss much between selected data and target nucleus 
      if ( theParticles != NULL ) {
         G4int sumA = 0;
         G4int sumZ = 0;
         G4int maxA = 0;
         G4int jAtMaxA = 0;
         for ( G4int j = 0 ; j != (G4int)theParticles->size() ; j++ ) {
            if ( theParticles->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
               maxA = theParticles->at(j)->GetDefinition()->GetBaryonNumber(); 
               jAtMaxA = j; 
            }
            sumA += theParticles->at(j)->GetDefinition()->GetBaryonNumber();
            sumZ += G4int( theParticles->at(j)->GetDefinition()->GetPDGCharge() + eps );
         }
         G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
         G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
         if ( dA < 0 || dZ < 0 ) {
            G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
            G4int newZ = G4int( theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
            G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
            theParticles->at( jAtMaxA )->SetDefinition( pd );
         }
      }
      //141017 Fix END
 
    }
    else
    {
      // @@@ what to do, if we have photon data, but no info on the hadron itself
      nothingWasKnownOnHadron = 1;
    }
 
    //G4cout << "theFinalStatePhotons it " << it << G4endl;
    //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
    //G4cout << "theFinalStatePhotons it " << it << G4endl;
    //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
    //G4cout << "thePhotons " << thePhotons << G4endl;
 
    if ( theFinalStatePhotons[it] != 0 ) 
    {
       // the photon distributions are in the Nucleus rest frame.
       // TK residual rest frame
      G4ReactionProduct boosted_tmp;
      boosted_tmp.Lorentz(incidReactionProduct, theTarget);
      G4double anEnergy = boosted_tmp.GetKineticEnergy();
      thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
      G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
      G4double testEnergy = 0;
      if(thePhotons!=0 && thePhotons->size()!=0)
      { aBaseEnergy-=thePhotons->operator[](0)->GetTotalEnergy(); }
      if(theFinalStatePhotons[it]->NeedsCascade())
      {
    while(aBaseEnergy>0.01*CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
        {
          // cascade down the levels
      G4bool foundMatchingLevel = false;
          G4int closest = 2;
      G4double deltaEold = -1;
      for(G4int j=1; j<it; j++)
          {
            if(theFinalStatePhotons[j]!=0) 
            {
              testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
            }
            else
            {
              testEnergy = 0;
            }
        G4double deltaE = std::abs(testEnergy-aBaseEnergy);
            if(deltaE<0.1*CLHEP::keV)
            {
              G4ReactionProductVector * theNext = 
            theFinalStatePhotons[j]->GetPhotons(anEnergy);
              if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
              aBaseEnergy = testEnergy-theNext->operator[](0)->GetTotalEnergy();
              delete theNext;
          foundMatchingLevel = true;
              break; // ===>
            }
        if(theFinalStatePhotons[j]!=0 && ( deltaE<deltaEold||deltaEold<0.) )
        {
          closest = j;
          deltaEold = deltaE;     
        }
          } // <=== the break goes here.
      if(!foundMatchingLevel)
      {
            G4ReactionProductVector * theNext = 
               theFinalStatePhotons[closest]->GetPhotons(anEnergy);
            if ( thePhotons != NULL ) thePhotons->push_back(theNext->operator[](0));
            aBaseEnergy = aBaseEnergy-theNext->operator[](0)->GetTotalEnergy();
            delete theNext;
      }
        } 
      }
    }
    unsigned int i0;
    if(thePhotons!=0)
    {
      for(i0=0; i0<thePhotons->size(); i0++)
      {
    // back to lab
    thePhotons->operator[](i0)->Lorentz(*(thePhotons->operator[](i0)), -1.*theTarget);
      }
    }
    //G4cout << "nothingWasKnownOnHadron " << nothingWasKnownOnHadron << G4endl;
    if (nothingWasKnownOnHadron)
    {
//    In this case, hadron should be isotropic in CM
// Next 12 lines are Emilio's replacement 
      // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
      // G4double eExcitation = QM-QI[it];
      // G4double eExcitation = QI[0] - QI[it];  // Fix of bug #1838
      // if(eExcitation<20*CLHEP::keV){eExcitation=0;}
 
      G4double eExcitation = std::max(0.,QI[0] - QI[it]);  // Fix of bug #2333
 
      two_body_reaction(&incidReactionProduct,&theTarget,&aHadron,eExcitation);
      if(thePhotons==0 && eExcitation>0){
        for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
        {
          if(theGammas.GetLevelEnergy(iLevel)<eExcitation+5*keV) break; // 5 keV tolerance
        }
        thePhotons = theGammas.GetDecayGammas(iLevel);
      }
    }
// Emilio's replacement done
/*
// This code replaced by Emilio (previous 12 lines) 
//    mu and p should be correlated
//
      //isotropic distribution in CM 
      G4double mu = 1.0 - 2.*G4UniformRand();
 
      // Need momenta in target rest frame
      G4LorentzVector target_in_LAB ( theTarget.GetMomentum() , theTarget.GetTotalEnergy() );
      G4ThreeVector boostToTargetRest = -target_in_LAB.boostVector();
      G4LorentzVector proj_in_LAB = hadProjectile->Get4Momentum();
 
      G4DynamicParticle* proj = new G4DynamicParticle(theProjectile, proj_in_LAB.boost(boostToTargetRest) ); 
//      G4DynamicParticle* targ = 
//        new G4DynamicParticle(G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, totalPhotonEnergy), G4ThreeVector(0) );
//    Fix bug 2166 (A. Zontikov): replace above two lines with next three lines
      G4double excitationEnergy = theFinalStatePhotons[it] ? theFinalStatePhotons[it]->GetLevelEnergy() : 0.0;
      G4DynamicParticle* targ = 
        new G4DynamicParticle(G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, excitationEnergy), G4ThreeVector(0) );
      G4DynamicParticle* hadron = 
        new G4DynamicParticle(aHadron.GetDefinition(), G4ThreeVector(0) );  // Will fill in the momentum
 
      two_body_reaction ( proj , targ , hadron , mu );
 
      G4LorentzVector hadron_in_trag_rest = hadron->Get4Momentum();
      G4LorentzVector hadron_in_LAB = hadron_in_trag_rest.boost ( -boostToTargetRest );
      aHadron.SetMomentum( hadron_in_LAB.v() );
      aHadron.SetKineticEnergy ( hadron_in_LAB.e() - hadron_in_LAB.m() );
 
      delete proj;
      delete targ; 
      delete hadron;
 
    }
*/
 
// fill the result
// Beware - the recoil is not necessarily in the particles...
// Can be calculated from momentum conservation?
// The idea is that the particles ar emitted forst, and the gammas only once the
// recoil is on the residual; assumption is that gammas do not contribute to 
// the recoil.
// This needs more design @@@
 
    G4int nSecondaries = 2; // the hadron and the recoil
    G4bool needsSeparateRecoil = false;
    G4int totalBaryonNumber = 0;
    G4int totalCharge = 0;
    G4ThreeVector totalMomentum(0);
    if(theParticles != 0) 
    {
      nSecondaries = theParticles->size();
      const G4ParticleDefinition * aDef;
      unsigned int ii0;
      for(ii0=0; ii0<theParticles->size(); ii0++)
      {
        aDef = theParticles->operator[](ii0)->GetDefinition();
    totalBaryonNumber+=aDef->GetBaryonNumber();
    totalCharge+=G4int(aDef->GetPDGCharge()+eps);
        totalMomentum += theParticles->operator[](ii0)->GetMomentum();
      } 
      if(totalBaryonNumber!=G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber())) 
      {
        needsSeparateRecoil = true;
    nSecondaries++;
    residualA = G4int(theBaseA+eps+hadProjectile->GetDefinition()->GetBaryonNumber()
                      -totalBaryonNumber);
    residualZ = G4int(theBaseZ+eps+hadProjectile->GetDefinition()->GetPDGCharge()
                      -totalCharge);
      }
    }
    
    G4int nPhotons = 0;
    if(thePhotons!=0) { nPhotons = thePhotons->size(); }
    nSecondaries += nPhotons;
        
    G4DynamicParticle * theSec;
    
    if( theParticles==0 )
    {
      theSec = new G4DynamicParticle;   
      theSec->SetDefinition(aHadron.GetDefinition());
      theSec->SetMomentum(aHadron.GetMomentum());
      theResult.Get()->AddSecondary(theSec, secID);    
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply  add secondary1 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
      
      aHadron.Lorentz(aHadron, theTarget);
      G4ReactionProduct theResidual;   
      theResidual.SetDefinition(G4IonTable::GetIonTable()
                ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
      theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
      
      //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
      //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
      G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
      theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
      
      theResidual.Lorentz(theResidual, -1.*theTarget);
      G4ThreeVector totalPhotonMomentum(0,0,0);
      if(thePhotons!=0)
    {
          for(i=0; i<nPhotons; i++)
        {
          totalPhotonMomentum += thePhotons->operator[](i)->GetMomentum();
        }
    }
      theSec = new G4DynamicParticle;   
      theSec->SetDefinition(theResidual.GetDefinition());
      theSec->SetMomentum(theResidual.GetMomentum()-totalPhotonMomentum);
      theResult.Get()->AddSecondary(theSec, secID);    
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
    }
    else
    {
      for(i0=0; i0<theParticles->size(); i0++)
      {
        theSec = new G4DynamicParticle; 
        theSec->SetDefinition(theParticles->operator[](i0)->GetDefinition());
        theSec->SetMomentum(theParticles->operator[](i0)->GetMomentum());
        theResult.Get()->AddSecondary(theSec, secID); 
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
        delete theParticles->operator[](i0); 
      } 
      delete theParticles;
      if(needsSeparateRecoil && residualZ!=0)
      {
        G4ReactionProduct theResidual;   
        theResidual.SetDefinition(G4IonTable::GetIonTable()
                              ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
        G4double resiualKineticEnergy  = theResidual.GetMass()*theResidual.GetMass();
                 resiualKineticEnergy += totalMomentum*totalMomentum;
             resiualKineticEnergy  = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
//        cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
    theResidual.SetKineticEnergy(resiualKineticEnergy);
 
        //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
        //theResidual.SetMomentum(-1.*totalMomentum);
    //G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
        //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
//080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
        theResidual.SetMomentum( incidReactionProduct.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
 
        theSec = new G4DynamicParticle;   
        theSec->SetDefinition(theResidual.GetDefinition());
        theSec->SetMomentum(theResidual.GetMomentum());
        theResult.Get()->AddSecondary(theSec, secID);  
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
 
      }  
    }
    if(thePhotons!=0)
    {
      for(i=0; i<nPhotons; i++)
      {
        theSec = new G4DynamicParticle;    
        //Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25, 2009 
        //theSec->SetDefinition(G4Gamma::Gamma());
        theSec->SetDefinition( thePhotons->operator[](i)->GetDefinition() );
        //But never cause real effect at least with G4NDL3.13 TK
        theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
        theResult.Get()->AddSecondary(theSec, secID); 
#ifdef G4PHPDEBUG
      if( std::getenv("G4ParticleHPDebug"))  G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
#endif
 
        delete thePhotons->operator[](i);
      }
// some garbage collection
      delete thePhotons;
    }
 
//080721 
   G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
   G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
   G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
   G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
   adjust_final_state ( init_4p_lab ); 
 
// clean up the primary neutron
    theResult.Get()->SetStatusChange( stopAndKill );
}

References G4HadFinalState::AddSecondary(), G4ParticleHPFinalState::adjust_final_state(), G4HadFinalState::Clear(), G4Neutron::Definition(), eps, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4Cache< VALTYPE >::Get(), G4HadProjectile::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4Nucleus::GetBiasedThermalNucleus(), G4ParticleHPDeExGammas::GetDecayGammas(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4ParticleHPManager::GetInstance(), G4IonTable::GetIon(), G4IonTable::GetIonTable(), G4DynamicParticle::GetKineticEnergy(), G4HadProjectile::GetKineticEnergy(), G4ReactionProduct::GetKineticEnergy(), G4ParticleHPDeExGammas::GetLevel(), G4ParticleHPLevel::GetLevelEnergy(), G4ParticleHPPhotonDist::GetLevelEnergy(), G4ParticleHPDeExGammas::GetLevelEnergy(), G4ReactionProduct::GetMass(), G4HadProjectile::GetMaterial(), G4ReactionProduct::GetMomentum(), G4NucleiProperties::GetNuclearMass(), G4ParticleHPDeExGammas::GetNumberOfLevels(), G4HadFinalState::GetNumberOfSecondaries(), G4DynamicParticle::GetParticleDefinition(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ParticleHPPhotonDist::GetPhotons(), G4Material::GetTemperature(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), InitDistributionInitialState(), CLHEP::keV, keV, G4ReactionProduct::Lorentz(), CLHEP::Hep3Vector::mag2(), G4INCL::Math::max(), maxA, MeV, G4Neutron::Neutron(), G4Cache< VALTYPE >::Put(), QI, G4ParticleHPEnAngCorrelation::Sample(), G4ParticleHPEnergyDistribution::Sample(), G4ParticleHPAngular::SampleAndUpdate(), G4ParticleHPFinalState::secID, SelectExitChannel(), G4DynamicParticle::SetDefinition(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetStatusChange(), stopAndKill, theAngularDistribution, G4ParticleHPFinalState::theBaseA, G4ParticleHPFinalState::theBaseZ, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, theGammas, G4ParticleHPFinalState::theProjectile, G4ParticleHPFinalState::theResult, theXsection, two_body_reaction(), CLHEP::twopi, use_nresp71_model(), and CLHEP::HepLorentzVector::vect().

Referenced by G4ParticleHPAInelasticFS::ApplyYourself(), G4ParticleHPDInelasticFS::ApplyYourself(), G4ParticleHPHe3InelasticFS::ApplyYourself(), G4ParticleHPNInelasticFS::ApplyYourself(), G4ParticleHPPInelasticFS::ApplyYourself(), and G4ParticleHPTInelasticFS::ApplyYourself().

GetA()

G4double G4ParticleHPFinalState::GetA ( void  )

inlineinherited

Definition at line 104 of file G4ParticleHPFinalState.hh.

{ return theBaseA; }

References G4ParticleHPFinalState::theBaseA.

GetM()

G4int G4ParticleHPFinalState::GetM ( )

inlineinherited

Definition at line 105 of file G4ParticleHPFinalState.hh.

{ return theBaseM; }

References G4ParticleHPFinalState::theBaseM.

Referenced by G4ParticleHPChannel::GetM().

GetN()

G4double G4ParticleHPFinalState::GetN ( )

inlineinherited

Definition at line 103 of file G4ParticleHPFinalState.hh.

{ return theBaseA; }

References G4ParticleHPFinalState::theBaseA.

Referenced by G4ParticleHPChannel::ApplyYourself(), and G4ParticleHPChannel::GetN().

GetXsec() [1/2]

virtual G4ParticleHPVector * G4ParticleHPInelasticCompFS::GetXsec ( )

inlinevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 95 of file G4ParticleHPInelasticCompFS.hh.

{ return theXsection[50]; }

References theXsection.

Referenced by SelectExitChannel().

GetXsec() [2/2]

virtual G4double G4ParticleHPInelasticCompFS::GetXsec ( G4double  anEnergy )

inlinevirtual

Reimplemented from G4ParticleHPFinalState.

Definition at line 90 of file G4ParticleHPInelasticCompFS.hh.

    {
      return std::max(0., theXsection[50]->GetY(anEnergy));
    }

References G4INCL::Math::max(), and theXsection.

GetZ()

G4double G4ParticleHPFinalState::GetZ ( void  )

inlineinherited

Definition at line 102 of file G4ParticleHPFinalState.hh.

{ return theBaseZ; }

References G4ParticleHPFinalState::theBaseZ.

Referenced by G4ParticleHPChannel::ApplyYourself(), and G4ParticleHPChannel::GetZ().

HasAnyData()

G4bool G4ParticleHPFinalState::HasAnyData ( )

inlineinherited

Definition at line 96 of file G4ParticleHPFinalState.hh.

{ return hasAnyData; }

References G4ParticleHPFinalState::hasAnyData.

Referenced by G4ParticleHPChannel::HasAnyData().

HasFSData()

G4bool G4ParticleHPFinalState::HasFSData ( )

inlineinherited

Definition at line 95 of file G4ParticleHPFinalState.hh.

{ return hasFSData; }

References G4ParticleHPFinalState::hasFSData.

Referenced by G4ParticleHPCaptureFS::ApplyYourself(), G4ParticleHPFissionBaseFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), G4ParticleHPChannel::HasFSData(), and G4ParticleHPFissionFS::Init().

HasXsec()

G4bool G4ParticleHPFinalState::HasXsec ( )

inlineinherited

Definition at line 94 of file G4ParticleHPFinalState.hh.

{ return hasXsec; }

References G4ParticleHPFinalState::hasXsec.

Referenced by G4ParticleHPChannel::DumpInfo().

Init() [1/2]

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

virtual

Implements G4ParticleHPFinalState.

Reimplemented in G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

Definition at line 86 of file G4ParticleHPInelasticCompFS.cc.

{
  gammaPath = "/Inelastic/Gammas/";  //only in neutron data base 
    if(!std::getenv("G4NEUTRONHPDATA")) 
       throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files where Inelastic/Gammas data is found.");
  G4String tBase = std::getenv("G4NEUTRONHPDATA");
  gammaPath = tBase+gammaPath;
  G4String tString = dirName;
  G4bool dbool;
  G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, aFSType, dbool);
  G4String filename = aFile.GetName();
#ifdef G4PHPDEBUG
  if( std::getenv("G4ParticleHPDebug") ) G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
#endif
 
  SetAZMs( A, Z, M, aFile ); 
  //theBaseA = aFile.GetA();
  //theBaseZ = aFile.GetZ();
  //theNDLDataA = (int)aFile.GetA();
  //theNDLDataZ = aFile.GetZ();
  //if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
  if ( !dbool || ( Z<2.5 && ( std::abs(theNDLDataZ - Z)>0.0001 || std::abs(theNDLDataA - A)>0.0001)) )
  {
#ifdef G4PHPDEBUG
    if(std::getenv("G4ParticleHPDebug_NamesLogging")) G4cout << "Skipped = "<< filename <<" "<<A<<" "<<Z<<G4endl;
#endif
    hasAnyData = false;
    hasFSData = false; 
    hasXsec = false;
    return;
  }
   //  theBaseA = A;
   //  theBaseZ = G4int(Z+.5);
//std::ifstream theData(filename, std::ios::in);
   std::istringstream theData(std::ios::in);
   G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
  if(!theData) //"!" is a operator of ios
  {
    hasAnyData = false;
    hasFSData = false; 
    hasXsec = false;
    //    theData.close();
    return;
  }
  // here we go
  G4int infoType, dataType, dummy;
  G4int sfType, it;
  hasFSData = false; 
  while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
  {
    hasFSData = true; 
    theData >> dataType;
    theData >> sfType >> dummy;
    it = 50;
    if(sfType>=600||(sfType<100&&sfType>=50)) it = sfType%50;
    if(dataType==3) 
    {
      //theData >> dummy >> dummy;
      //TK110430
      // QI and LR introudced since G4NDL3.15
      G4double dqi;
      G4int ilr;
      theData >> dqi >> ilr;
 
      QI[ it ] = dqi*CLHEP::eV;
      LR[ it ] = ilr;
      theXsection[it] = new G4ParticleHPVector;
      G4int total;
      theData >> total;
      theXsection[it]->Init(theData, total, CLHEP::eV);
      //std::cout << theXsection[it]->GetXsec(1*MeV) << std::endl;
    }
    else if(dataType==4)
    {
      theAngularDistribution[it] = new G4ParticleHPAngular;
      theAngularDistribution[it]->Init(theData);
    }
    else if(dataType==5)
    {
      theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
      theEnergyDistribution[it]->Init(theData); 
    }
    else if(dataType==6)
    {
      theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
      //      G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl; //GDEB
      theEnergyAngData[it]->Init(theData);
    }
    else if(dataType==12)
    {
      theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
      theFinalStatePhotons[it]->InitMean(theData);
    }
    else if(dataType==13)
    {
      theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
      theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]);
    }
    else if(dataType==14)
    {
      theFinalStatePhotons[it]->InitAngular(theData);
    }
    else if(dataType==15)
    {
      theFinalStatePhotons[it]->InitEnergies(theData);
    }
    else
    {
      throw G4HadronicException(__FILE__, __LINE__, "Data-type unknown to G4ParticleHPInelasticCompFS");
    }
  }
  //  theData.close();
}

References A, CLHEP::eV, G4cout, G4endl, gammaPath, G4ParticleHPManager::GetDataStream(), G4ParticleHPManager::GetInstance(), G4ParticleHPDataUsed::GetName(), G4ParticleHPNames::GetName(), G4ParticleHPFinalState::hasAnyData, G4ParticleHPFinalState::hasFSData, G4ParticleHPFinalState::hasXsec, G4ParticleHPAngular::Init(), G4ParticleHPEnAngCorrelation::Init(), G4ParticleHPVector::Init(), G4ParticleHPEnergyDistribution::Init(), G4ParticleHPPhotonDist::InitAngular(), G4ParticleHPPhotonDist::InitEnergies(), G4ParticleHPPhotonDist::InitMean(), G4ParticleHPPhotonDist::InitPartials(), LR, M, QI, G4ParticleHPFinalState::SetAZMs(), theAngularDistribution, theEnergyAngData, theEnergyDistribution, theFinalStatePhotons, G4ParticleHPFinalState::theNames, G4ParticleHPFinalState::theNDLDataA, G4ParticleHPFinalState::theNDLDataZ, G4ParticleHPFinalState::theProjectile, theXsection, G4INCL::CrossSections::total(), and Z.

Referenced by G4ParticleHPAInelasticFS::Init(), G4ParticleHPDInelasticFS::Init(), G4ParticleHPHe3InelasticFS::Init(), G4ParticleHPNInelasticFS::Init(), G4ParticleHPPInelasticFS::Init(), and G4ParticleHPTInelasticFS::Init().

Init() [2/2]

void G4ParticleHPFinalState::Init ( G4double  A, G4double  Z, G4String &  dirName, G4String &  aFSType, G4ParticleDefinition *  projectile  )

inlineinherited

Definition at line 76 of file G4ParticleHPFinalState.hh.

  {
    G4int M = 0;
    Init ( A, Z, M, dirName, aFSType,const_cast<G4ParticleDefinition*>(projectile));
  }

References A, G4ParticleHPFinalState::Init(), M, and Z.

Referenced by G4ParticleHPFinalState::Init(), and G4ParticleHPChannel::UpdateData().

InitDistributionInitialState()

void G4ParticleHPInelasticCompFS::InitDistributionInitialState ( G4ReactionProduct &  anIncidentPart, G4ReactionProduct &  aTarget, G4int  it  )

inline

Definition at line 102 of file G4ParticleHPInelasticCompFS.hh.

    {
      if (theAngularDistribution[it] != 0) {
        theAngularDistribution[it]->SetTarget(aTarget);
        theAngularDistribution[it]->SetProjectileRP(anIncidentPart);
      }
 
      if (theEnergyAngData[it] != 0) {
        theEnergyAngData[it]->SetTarget(aTarget);
        theEnergyAngData[it]->SetProjectileRP(anIncidentPart);
      }
    }

References G4ParticleHPAngular::SetProjectileRP(), G4ParticleHPEnAngCorrelation::SetProjectileRP(), G4ParticleHPAngular::SetTarget(), G4ParticleHPEnAngCorrelation::SetTarget(), theAngularDistribution, and theEnergyAngData.

Referenced by CompositeApply().

InitGammas()

void G4ParticleHPInelasticCompFS::InitGammas	(	G4double	AR,
		G4double	ZR
	)

Definition at line 64 of file G4ParticleHPInelasticCompFS.cc.

{
  //   char the[100] = {""};
  //   std::ostrstream ost(the, 100, std::ios::out);
  //   ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
  //   G4String * aName = new G4String(the);
  //   std::ifstream from(*aName, std::ios::in);
 
   std::ostringstream ost;
   ost <<gammaPath<<"z"<<ZR<<".a"<<AR;
   G4String aName = ost.str();
   std::ifstream from(aName, std::ios::in);
 
   if(!from) return; // no data found for this isotope
   //   std::ifstream theGammaData(*aName, std::ios::in);
   std::ifstream theGammaData(aName, std::ios::in);
    
   theGammas.Init(theGammaData);
   //   delete aName;
 
}

References gammaPath, G4ParticleHPDeExGammas::Init(), and theGammas.

Referenced by G4ParticleHPAInelasticFS::Init(), G4ParticleHPDInelasticFS::Init(), G4ParticleHPHe3InelasticFS::Init(), G4ParticleHPNInelasticFS::Init(), G4ParticleHPPInelasticFS::Init(), and G4ParticleHPTInelasticFS::Init().

New()

virtual G4ParticleHPFinalState * G4ParticleHPInelasticCompFS::New ( )

pure virtual

Implements G4ParticleHPFinalState.

Implemented in G4ParticleHPAInelasticFS, G4ParticleHPDInelasticFS, G4ParticleHPHe3InelasticFS, G4ParticleHPNInelasticFS, G4ParticleHPPInelasticFS, and G4ParticleHPTInelasticFS.

SelectExitChannel()

G4int G4ParticleHPInelasticCompFS::SelectExitChannel

(

G4double

eKinetic

)

Definition at line 200 of file G4ParticleHPInelasticCompFS.cc.

{
 
// it = 0 has without Photon
  G4double running[50];
  running[0] = 0;
  unsigned int i;
  for(i=0; i<50; i++)
  {
    if(i!=0) running[i]=running[i-1];
    if(theXsection[i] != 0) 
    {
      running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
    }
  }
  G4double random = G4UniformRand();
  G4double sum = running[49];
  G4int it = 50;
  if(0!=sum)
  {
    G4int i0;
    for(i0=0; i0<50; i0++)
    {
      it = i0;
      //       G4cout << " SelectExitChannel " << it << " " << random << " " << running[i0]/sum << " " << running[i0] << G4endl; //GDEB
      if(random < running[i0]/sum) break;
    }
  }
//debug:  it = 1;
//  G4cout << " SelectExitChannel " << it << " " << sum << G4endl; //GDEB
  return it;
}

References G4UniformRand, GetXsec(), G4INCL::Math::max(), and theXsection.

Referenced by CompositeApply().

SetA_Z()

void G4ParticleHPFinalState::SetA_Z ( G4double  anA, G4double  aZ, G4int  aM = 0  )

inlineinherited

Definition at line 101 of file G4ParticleHPFinalState.hh.

{theBaseA = anA; theBaseZ = aZ; theBaseM=aM; }

References G4ParticleHPFinalState::theBaseA, G4ParticleHPFinalState::theBaseM, and G4ParticleHPFinalState::theBaseZ.

Referenced by G4ParticleHPChannel::Register().

SetAZMs()

void G4ParticleHPFinalState::SetAZMs ( G4double  anA, G4double  aZ, G4int  aM, G4ParticleHPDataUsed  used  )

inlineinherited

Definition at line 107 of file G4ParticleHPFinalState.hh.

  {
    theBaseA = anA; theBaseZ = aZ; theBaseM=aM; 
    theNDLDataA=(G4int)used.GetA();
    theNDLDataZ=(G4int)used.GetZ();
    theNDLDataM=used.GetM();
  }

References G4ParticleHPDataUsed::GetA(), G4ParticleHPDataUsed::GetM(), G4ParticleHPDataUsed::GetZ(), G4ParticleHPFinalState::theBaseA, G4ParticleHPFinalState::theBaseM, G4ParticleHPFinalState::theBaseZ, G4ParticleHPFinalState::theNDLDataA, G4ParticleHPFinalState::theNDLDataM, and G4ParticleHPFinalState::theNDLDataZ.

Referenced by G4ParticleHPCaptureFS::Init(), G4ParticleHPElasticFS::Init(), G4ParticleHPFSFissionFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

SetProjectile()

void G4ParticleHPFinalState::SetProjectile ( G4ParticleDefinition *  projectile )

inlineinherited

Definition at line 115 of file G4ParticleHPFinalState.hh.

  {
    theProjectile = projectile;
  }

References G4ParticleHPFinalState::theProjectile.

Referenced by G4ParticleHPChannel::Register().

two_body_reaction()

void G4ParticleHPInelasticCompFS::two_body_reaction ( G4ReactionProduct *  proj, G4ReactionProduct *  targ, G4ReactionProduct *  product, G4double  resExcitationEnergy  )

private

Definition at line 816 of file G4ParticleHPInelasticCompFS.cc.

{
  //CMS system:
  G4ReactionProduct theCMS= *proj+ *targ;
 
  //Residual definition:
  G4int resZ=(G4int)(proj->GetDefinition()->GetPDGCharge()+targ->GetDefinition()->GetPDGCharge()-product->GetDefinition()->GetPDGCharge()+0.1);
  G4int resA=proj->GetDefinition()->GetBaryonNumber()+targ->GetDefinition()->GetBaryonNumber()-product->GetDefinition()->GetBaryonNumber();
  G4ReactionProduct theResidual;
  theResidual.SetDefinition(G4IonTable::GetIonTable()->GetIon(resZ,resA,0.0));
 
  //CMS system:
  G4ReactionProduct theCMSproj;
  G4ReactionProduct theCMStarg;
  theCMSproj.Lorentz(*proj,theCMS);
  theCMStarg.Lorentz(*targ,theCMS);
  //final Momentum in the CMS:
  G4double totE=std::sqrt(theCMSproj.GetMass()*theCMSproj.GetMass()+theCMSproj.GetTotalMomentum()*theCMSproj.GetTotalMomentum())+std::sqrt(theCMStarg.GetMass()*theCMStarg.GetMass()+theCMStarg.GetTotalMomentum()*theCMStarg.GetTotalMomentum());
  G4double prodmass=product->GetMass();
  G4double resmass=theResidual.GetMass()+resExcitationEnergy;
  G4double fmomsquared=1./4./totE/totE*(totE*totE-(prodmass-resmass)*(prodmass-resmass))*(totE*totE-(prodmass+resmass)*(prodmass+resmass));
  G4double fmom=0;
  if(fmomsquared>0){
    fmom=std::sqrt(fmomsquared);
  }
 
  //random (isotropic direction):
  G4double cosTh = 2.*G4UniformRand()-1.;
  G4double phi = CLHEP::twopi*G4UniformRand();
  G4double theta = std::acos(cosTh);
  G4double sinth = std::sin(theta);
  product->SetMomentum(fmom*sinth*std::cos(phi),fmom*sinth*std::sin(phi),fmom*cosTh); //CMS
  product->SetTotalEnergy(std::sqrt(prodmass*prodmass+fmom*fmom)); //CMS
  //Back to the LAB system:
  product->Lorentz(*product,-1.*theCMS);
 
}

References G4UniformRand, G4ParticleDefinition::GetBaryonNumber(), G4ReactionProduct::GetDefinition(), G4IonTable::GetIon(), G4IonTable::GetIonTable(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGCharge(), G4ReactionProduct::GetTotalMomentum(), G4ReactionProduct::Lorentz(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTotalEnergy(), and CLHEP::twopi.

Referenced by CompositeApply().

use_nresp71_model()

G4bool G4ParticleHPInelasticCompFS::use_nresp71_model ( const G4ParticleDefinition *  aDefinition, const G4int  it, const G4ReactionProduct &  theTarget, G4ReactionProduct &  boosted  )

private

Definition at line 858 of file G4ParticleHPInelasticCompFS.cc.

{
    if ( aDefinition == G4Neutron::Definition() ) // If the outgoing particle is a neutron...
        {
        // LR: flag LR in ENDF. It indicates whether there is breakup of the residual nucleus or not.
        // it: exit channel (index of the carbon excited state)
 
        //if ( (G4int)(theBaseZ+0.1) == 6 ) G4cout << "LR[" << it << "] = " << LR[it] << G4endl;
 
        // Added by A. R. García (CIEMAT) to include the physics of C(N,N'3A) reactions from NRESP71.
 
        if ( LR[it] > 0 ) // If there is breakup of the residual nucleus LR(flag LR in ENDF)>0 (i.e. Z=6 MT=52-91 (it=MT-50)).
            {
            // Defining carbon as the target in the reference frame at rest.
            G4ReactionProduct theCarbon(theTarget);
 
            theCarbon.SetMomentum(G4ThreeVector());
            theCarbon.SetKineticEnergy(0.);
 
            // Creating four reaction products.
            G4ReactionProduct theProds[4];
 
            // Applying C(N,N'3A) reaction mechanisms in the target rest frame.
            if ( it == 41 )
                {
                // QI=QM=-7.275 MeV for C-0(N,N')C-C(3A) in ENDF/B-VII.1. This is not the value of the QI of the first step according to the model. So we don't take it. Instead, we set the one we have calculated: QI=(mn+m12C)-(ma+m9Be+Ex9Be)=-8.130 MeV.
                nresp71_model.ApplyMechanismI_NBeA2A(boosted, theCarbon, theProds, -8.130/*QI[it]*/); // N+C --> A[0]+9BE* | 9BE* --> N[1]+8BE | 8BE --> 2*A[2,3].
                //printf("- QI=%f\n", QI[it]);
                }
            else
                {
                nresp71_model.ApplyMechanismII_ACN2A(boosted, theCarbon, theProds, QI[it]); // N+C --> N'[0]+C* | C* --> A[1]+8BE | 8BE --> 2*A[2,3].
                }
 
            //printf("it=%d   qi=%f  \n", it, QI[it]);
 
            // Returning to the reference frame where the target was in motion.
            for ( G4int j=0; j<4; j++ )
                {
                theProds[j].Lorentz(theProds[j], -1.*theTarget);
                theResult.Get()->AddSecondary(new G4DynamicParticle(theProds[j].GetDefinition(), theProds[j].GetMomentum()), secID);
                }
 
            /*G4double EN0 = theNeutron.GetKineticEnergy();
            G4double EN1 = theProds[0].GetKineticEnergy();
 
            G4double EA1 = theProds[1].GetKineticEnergy();
            G4double EA2 = theProds[2].GetKineticEnergy();
            G4double EA3 = theProds[3].GetKineticEnergy();
 
            printf("Q=%f\n", EN1+EA1+EA2+EA3-EN0);*/
 
            // Killing the primary neutron.
            theResult.Get()->SetStatusChange(stopAndKill);
 
            return true;
            }
        }
    else if ( aDefinition == G4Alpha::Definition() ) // If the outgoing particle is an alpha, ...
        {
        // Added by A. R. García (CIEMAT) to include the physics of C(N,A)9BE reactions from NRESP71.
 
        if ( LR[it] == 0 ) // If Z=6, an alpha particle is emitted and there is no breakup of the residual nucleus LR(flag LR in ENDF)==0.
            {
            // Defining carbon as the target in the reference frame at rest.
            G4ReactionProduct theCarbon(theTarget);
 
            theCarbon.SetMomentum(G4ThreeVector());
            theCarbon.SetKineticEnergy(0.);
 
            // Creating four reaction products.
            G4ReactionProduct theProds[2];
 
            // Applying C(N,A)9BE reaction mechanism.
            nresp71_model.ApplyMechanismABE(boosted, theCarbon, theProds); // N+C --> A[0]+9BE[1].
 
            //G4DynamicParticle *theSec;
            for ( G4int j=0; j<2; j++ )
                {
                // Returning to the system of reference where the target was in motion.
                theProds[j].Lorentz(theProds[j], -1.*theTarget);
                theResult.Get()->AddSecondary(new G4DynamicParticle(theProds[j].GetDefinition(), theProds[j].GetMomentum()), secID);
                }
 
            // Killing the primary neutron.
            theResult.Get()->SetStatusChange(stopAndKill);;
 
            return true;
            }
        else
            {
            G4Exception("G4ParticleHPInelasticCompFS::CompositeApply()", "G4ParticleInelasticCompFS.cc", FatalException, "Alpha production with LR!=0.");
            }
    }
 
   return false;
}

References G4HadFinalState::AddSecondary(), G4NRESP71M03::ApplyMechanismABE(), G4NRESP71M03::ApplyMechanismI_NBeA2A(), G4NRESP71M03::ApplyMechanismII_ACN2A(), G4Neutron::Definition(), G4Alpha::Definition(), FatalException, G4Exception(), G4Cache< VALTYPE >::Get(), G4ReactionProduct::Lorentz(), LR, nresp71_model, QI, G4ParticleHPFinalState::secID, G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMomentum(), G4HadFinalState::SetStatusChange(), stopAndKill, and G4ParticleHPFinalState::theResult.

Referenced by CompositeApply().

Field Documentation

gammaPath

G4String G4ParticleHPInelasticCompFS::gammaPath

protected

Definition at line 127 of file G4ParticleHPInelasticCompFS.hh.

Referenced by Init(), and InitGammas().

hasAnyData

G4bool G4ParticleHPFinalState::hasAnyData

protectedinherited

Definition at line 126 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFinalState::HasAnyData(), G4FissionLibrary::Init(), G4ParticleHPCaptureFS::Init(), G4ParticleHPElasticFS::Init(), G4ParticleHPFFFissionFS::Init(), G4ParticleHPFSFissionFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

hasFSData

G4bool G4ParticleHPFinalState::hasFSData

protectedinherited

Definition at line 125 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFinalState::HasFSData(), G4FissionLibrary::Init(), G4ParticleHPCaptureFS::Init(), G4ParticleHPElasticFS::Init(), G4ParticleHPFFFissionFS::Init(), G4ParticleHPFSFissionFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

hasXsec

G4bool G4ParticleHPFinalState::hasXsec

protectedinherited

Definition at line 124 of file G4ParticleHPFinalState.hh.

Referenced by G4FissionLibrary::G4FissionLibrary(), G4ParticleHPCaptureFS::G4ParticleHPCaptureFS(), G4ParticleHPElasticFS::G4ParticleHPElasticFS(), G4ParticleHPFCFissionFS::G4ParticleHPFCFissionFS(), G4ParticleHPFFFissionFS::G4ParticleHPFFFissionFS(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFissionBaseFS::G4ParticleHPFissionBaseFS(), G4ParticleHPFissionFS::G4ParticleHPFissionFS(), G4ParticleHPFSFissionFS::G4ParticleHPFSFissionFS(), G4ParticleHPInelasticBaseFS::G4ParticleHPInelasticBaseFS(), G4ParticleHPInelasticCompFS(), G4ParticleHPLCFissionFS::G4ParticleHPLCFissionFS(), G4ParticleHPSCFissionFS::G4ParticleHPSCFissionFS(), G4ParticleHPTCFissionFS::G4ParticleHPTCFissionFS(), G4ParticleHPFinalState::HasXsec(), G4FissionLibrary::Init(), G4ParticleHPCaptureFS::Init(), G4ParticleHPElasticFS::Init(), G4ParticleHPFFFissionFS::Init(), G4ParticleHPFSFissionFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

LR

std::vector<G4int> G4ParticleHPInelasticCompFS::LR

protected

Definition at line 131 of file G4ParticleHPInelasticCompFS.hh.

Referenced by G4ParticleHPInelasticCompFS(), Init(), and use_nresp71_model().

nresp71_model

G4NRESP71M03 G4ParticleHPInelasticCompFS::nresp71_model

private

Definition at line 138 of file G4ParticleHPInelasticCompFS.hh.

Referenced by use_nresp71_model().

QI

std::vector<G4double> G4ParticleHPInelasticCompFS::QI

protected

Definition at line 130 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), Init(), and use_nresp71_model().

secID

G4int G4ParticleHPFinalState::secID

protectedinherited

Definition at line 140 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPCaptureFS::ApplyYourself(), G4ParticleHPElasticFS::ApplyYourself(), G4ParticleHPFissionFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), CompositeApply(), G4ParticleHP2AInelasticFS::G4ParticleHP2AInelasticFS(), G4ParticleHP2N2AInelasticFS::G4ParticleHP2N2AInelasticFS(), G4ParticleHP2NAInelasticFS::G4ParticleHP2NAInelasticFS(), G4ParticleHP2NDInelasticFS::G4ParticleHP2NDInelasticFS(), G4ParticleHP2NInelasticFS::G4ParticleHP2NInelasticFS(), G4ParticleHP2NPInelasticFS::G4ParticleHP2NPInelasticFS(), G4ParticleHP2PInelasticFS::G4ParticleHP2PInelasticFS(), G4ParticleHP3AInelasticFS::G4ParticleHP3AInelasticFS(), G4ParticleHP3NAInelasticFS::G4ParticleHP3NAInelasticFS(), G4ParticleHP3NInelasticFS::G4ParticleHP3NInelasticFS(), G4ParticleHP3NPInelasticFS::G4ParticleHP3NPInelasticFS(), G4ParticleHP4NInelasticFS::G4ParticleHP4NInelasticFS(), G4ParticleHPAInelasticFS::G4ParticleHPAInelasticFS(), G4ParticleHPCaptureFS::G4ParticleHPCaptureFS(), G4ParticleHPD2AInelasticFS::G4ParticleHPD2AInelasticFS(), G4ParticleHPDAInelasticFS::G4ParticleHPDAInelasticFS(), G4ParticleHPDInelasticFS::G4ParticleHPDInelasticFS(), G4ParticleHPElasticFS::G4ParticleHPElasticFS(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFissionFS::G4ParticleHPFissionFS(), G4ParticleHPHe3InelasticFS::G4ParticleHPHe3InelasticFS(), G4ParticleHPN2AInelasticFS::G4ParticleHPN2AInelasticFS(), G4ParticleHPN2PInelasticFS::G4ParticleHPN2PInelasticFS(), G4ParticleHPN3AInelasticFS::G4ParticleHPN3AInelasticFS(), G4ParticleHPNAInelasticFS::G4ParticleHPNAInelasticFS(), G4ParticleHPND2AInelasticFS::G4ParticleHPND2AInelasticFS(), G4ParticleHPNDInelasticFS::G4ParticleHPNDInelasticFS(), G4ParticleHPNHe3InelasticFS::G4ParticleHPNHe3InelasticFS(), G4ParticleHPNInelasticFS::G4ParticleHPNInelasticFS(), G4ParticleHPNPAInelasticFS::G4ParticleHPNPAInelasticFS(), G4ParticleHPNPInelasticFS::G4ParticleHPNPInelasticFS(), G4ParticleHPNT2AInelasticFS::G4ParticleHPNT2AInelasticFS(), G4ParticleHPNTInelasticFS::G4ParticleHPNTInelasticFS(), G4ParticleHPNXInelasticFS::G4ParticleHPNXInelasticFS(), G4ParticleHPPAInelasticFS::G4ParticleHPPAInelasticFS(), G4ParticleHPPDInelasticFS::G4ParticleHPPDInelasticFS(), G4ParticleHPPInelasticFS::G4ParticleHPPInelasticFS(), G4ParticleHPPTInelasticFS::G4ParticleHPPTInelasticFS(), G4ParticleHPT2AInelasticFS::G4ParticleHPT2AInelasticFS(), G4ParticleHPTInelasticFS::G4ParticleHPTInelasticFS(), and use_nresp71_model().

theAngularDistribution

G4ParticleHPAngular* G4ParticleHPInelasticCompFS::theAngularDistribution[51]

protected

Definition at line 121 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4ParticleHPInelasticCompFS().

theBaseA

G4double G4ParticleHPFinalState::theBaseA

protectedinherited

Definition at line 132 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPCaptureFS::ApplyYourself(), G4ParticleHPElasticFS::ApplyYourself(), G4ParticleHPFissionFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), CompositeApply(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFinalState::GetA(), G4ParticleHPFinalState::GetN(), G4ParticleHPCaptureFS::Init(), G4ParticleHPFFFissionFS::Init(), G4ParticleHPInelasticBaseFS::InitGammas(), G4ParticleHPFinalState::SetA_Z(), and G4ParticleHPFinalState::SetAZMs().

theBaseM

G4int G4ParticleHPFinalState::theBaseM

protectedinherited

Definition at line 134 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFinalState::GetM(), G4ParticleHPFinalState::SetA_Z(), and G4ParticleHPFinalState::SetAZMs().

theBaseZ

G4double G4ParticleHPFinalState::theBaseZ

protectedinherited

Definition at line 133 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPCaptureFS::ApplyYourself(), G4ParticleHPElasticFS::ApplyYourself(), G4ParticleHPFissionFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), CompositeApply(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPFinalState::GetZ(), G4ParticleHPCaptureFS::Init(), G4ParticleHPFFFissionFS::Init(), G4ParticleHPInelasticBaseFS::InitGammas(), G4ParticleHPFinalState::SetA_Z(), and G4ParticleHPFinalState::SetAZMs().

theEnergyAngData

G4ParticleHPEnAngCorrelation* G4ParticleHPInelasticCompFS::theEnergyAngData[51]

protected

Definition at line 122 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), Init(), InitDistributionInitialState(), and ~G4ParticleHPInelasticCompFS().

theEnergyDistribution

G4ParticleHPEnergyDistribution* G4ParticleHPInelasticCompFS::theEnergyDistribution[51]

protected

Definition at line 120 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), Init(), and ~G4ParticleHPInelasticCompFS().

theFinalStatePhotons

G4ParticleHPPhotonDist* G4ParticleHPInelasticCompFS::theFinalStatePhotons[51]

protected

Definition at line 124 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), Init(), and ~G4ParticleHPInelasticCompFS().

theGammas

G4ParticleHPDeExGammas G4ParticleHPInelasticCompFS::theGammas

protected

Definition at line 126 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), and InitGammas().

theNames

G4ParticleHPNames G4ParticleHPFinalState::theNames

protectedinherited

Definition at line 127 of file G4ParticleHPFinalState.hh.

Referenced by G4FissionLibrary::Init(), G4ParticleHPElasticFS::Init(), G4ParticleHPFFFissionFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

theNDLDataA

G4int G4ParticleHPFinalState::theNDLDataA

protectedinherited

Definition at line 137 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPElasticFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), G4ParticleHPInelasticBaseFS::Init(), and G4ParticleHPFinalState::SetAZMs().

theNDLDataM

G4int G4ParticleHPFinalState::theNDLDataM

protectedinherited

Definition at line 138 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPElasticFS::Init(), and G4ParticleHPFinalState::SetAZMs().

theNDLDataZ

G4int G4ParticleHPFinalState::theNDLDataZ

protectedinherited

Definition at line 136 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPFinalState::G4ParticleHPFinalState(), G4ParticleHPElasticFS::Init(), Init(), G4ParticleHPFissionBaseFS::Init(), G4ParticleHPInelasticBaseFS::Init(), and G4ParticleHPFinalState::SetAZMs().

theProjectile

G4ParticleDefinition* G4ParticleHPFinalState::theProjectile

protectedinherited

Definition at line 130 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4ParticleHPInelasticBaseFS::BaseApply(), CompositeApply(), G4ParticleHPFinalState::G4ParticleHPFinalState(), Init(), G4ParticleHPInelasticBaseFS::Init(), and G4ParticleHPFinalState::SetProjectile().

theResult

G4Cache< G4HadFinalState* > G4ParticleHPFinalState::theResult

protectedinherited

Definition at line 129 of file G4ParticleHPFinalState.hh.

Referenced by G4ParticleHPFinalState::adjust_final_state(), G4FissionLibrary::ApplyYourself(), G4ParticleHP2AInelasticFS::ApplyYourself(), G4ParticleHP2N2AInelasticFS::ApplyYourself(), G4ParticleHP2NAInelasticFS::ApplyYourself(), G4ParticleHP2NDInelasticFS::ApplyYourself(), G4ParticleHP2NInelasticFS::ApplyYourself(), G4ParticleHP2NPInelasticFS::ApplyYourself(), G4ParticleHP2PInelasticFS::ApplyYourself(), G4ParticleHP3AInelasticFS::ApplyYourself(), G4ParticleHP3NAInelasticFS::ApplyYourself(), G4ParticleHP3NInelasticFS::ApplyYourself(), G4ParticleHP3NPInelasticFS::ApplyYourself(), G4ParticleHP4NInelasticFS::ApplyYourself(), G4ParticleHPAInelasticFS::ApplyYourself(), G4ParticleHPCaptureFS::ApplyYourself(), G4ParticleHPD2AInelasticFS::ApplyYourself(), G4ParticleHPDAInelasticFS::ApplyYourself(), G4ParticleHPDInelasticFS::ApplyYourself(), G4ParticleHPElasticFS::ApplyYourself(), G4ParticleHPFissionFS::ApplyYourself(), G4ParticleHPHe3InelasticFS::ApplyYourself(), G4ParticleHPN2AInelasticFS::ApplyYourself(), G4ParticleHPN2PInelasticFS::ApplyYourself(), G4ParticleHPN3AInelasticFS::ApplyYourself(), G4ParticleHPNAInelasticFS::ApplyYourself(), G4ParticleHPND2AInelasticFS::ApplyYourself(), G4ParticleHPNDInelasticFS::ApplyYourself(), G4ParticleHPNHe3InelasticFS::ApplyYourself(), G4ParticleHPNInelasticFS::ApplyYourself(), G4ParticleHPNPAInelasticFS::ApplyYourself(), G4ParticleHPNPInelasticFS::ApplyYourself(), G4ParticleHPNT2AInelasticFS::ApplyYourself(), G4ParticleHPNTInelasticFS::ApplyYourself(), G4ParticleHPNXInelasticFS::ApplyYourself(), G4ParticleHPPAInelasticFS::ApplyYourself(), G4ParticleHPPDInelasticFS::ApplyYourself(), G4ParticleHPPInelasticFS::ApplyYourself(), G4ParticleHPPTInelasticFS::ApplyYourself(), G4ParticleHPT2AInelasticFS::ApplyYourself(), G4ParticleHPTInelasticFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), CompositeApply(), G4ParticleHPFinalState::G4ParticleHPFinalState(), use_nresp71_model(), and G4ParticleHPFinalState::~G4ParticleHPFinalState().

theXsection

G4ParticleHPVector* G4ParticleHPInelasticCompFS::theXsection[51]

protected

Definition at line 119 of file G4ParticleHPInelasticCompFS.hh.

Referenced by CompositeApply(), G4ParticleHPInelasticCompFS(), GetXsec(), Init(), SelectExitChannel(), and ~G4ParticleHPInelasticCompFS().

---

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

• source/processes/hadronic/models/particle_hp/include/G4ParticleHPInelasticCompFS.hh
• source/processes/hadronic/models/particle_hp/src/G4ParticleHPInelasticCompFS.cc