#include <G4NeutronHPFinalState.hh>

Inheritance diagram for G4NeutronHPFinalState:

Public Member Functions
	G4NeutronHPFinalState ()

virtual	~G4NeutronHPFinalState ()

void	Init (G4double A, G4double Z, G4String &dirName, G4String &aFSType)

virtual void	Init (G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType)=0

virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &)

virtual G4NeutronHPFinalState *	New ()=0

G4bool	HasXsec ()

G4bool	HasFSData ()

G4bool	HasAnyData ()

virtual G4double	GetXsec (G4double)

virtual G4NeutronHPVector *	GetXsec ()

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

G4double	GetZ ()

G4double	GetN ()

G4int	GetM ()

Protected Member Functions
void	SetAZMs (G4double anA, G4double aZ, G4int aM, G4NeutronHPDataUsed used)

void	adjust_final_state (G4LorentzVector)

G4bool	DoNotAdjustFinalState ()

Protected Attributes
G4bool	hasXsec

G4bool	hasFSData

G4bool	hasAnyData

G4NeutronHPNames	theNames

G4HadFinalState	theResult

G4double	theBaseA

G4double	theBaseZ

G4int	theBaseM

G4int	theNDLDataZ

G4int	theNDLDataA

G4int	theNDLDataM

Detailed Description

Definition at line 41 of file G4NeutronHPFinalState.hh.

Constructor & Destructor Documentation

G4NeutronHPFinalState::G4NeutronHPFinalState ( )

inline

Definition at line 45 of file G4NeutronHPFinalState.hh.

References hasAnyData, hasFSData, hasXsec, theBaseA, theBaseM, theBaseZ, theNDLDataA, and theNDLDataZ.

   { 
     hasFSData = true; 
     hasXsec = true;
     hasAnyData = true;
     theBaseZ = 0;
     theBaseA = 0;
     theBaseM = 0;
 
     theNDLDataZ = 0;
     theNDLDataA = 0;
 
      adjustResult = true;
      if ( getenv( "G4NEUTRONHP_DO_NOT_ADJUST_FINAL_STATE" ) ) adjustResult = false;
 
   };

virtual G4NeutronHPFinalState::~G4NeutronHPFinalState ( )

inlinevirtual

Definition at line 62 of file G4NeutronHPFinalState.hh.

62 {};

Member Function Documentation

void G4NeutronHPFinalState::adjust_final_state ( G4LorentzVector init_4p_lab )

protected

Definition at line 38 of file G4NeutronHPFinalState.cc.

References G4HadFinalState::AddSecondary(), CLHEP::HepLorentzVector::beta(), CLHEP::HepLorentzVector::e(), G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::Get4Momentum(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4DynamicParticle::GetDefinition(), G4IonTable::GetIon(), G4IonTable::GetIonTable(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentum(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4ParticleDefinition::GetPDGMass(), G4HadFinalState::GetSecondary(), int(), python.hepunit::keV, CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), python.hepunit::MeV, G4Neutron::Neutron(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentum(), theBaseA, theBaseZ, theNDLDataA, theNDLDataZ, theResult, python.hepunit::twopi, and CLHEP::HepLorentzVector::v().

Referenced by G4NeutronHPInelasticBaseFS::BaseApply(), and G4NeutronHPInelasticCompFS::CompositeApply().

 {
 
    G4double minimum_energy = 1*keV;
 
    if ( adjustResult != true ) return;
 
    G4int nSecondaries = theResult.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++ )
    {
       sum_Z += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber();
       max_SecZ = std::max ( max_SecZ , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
       sum_A += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass();
       max_SecA = std::max ( max_SecA , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
       if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
    }
 
    G4ParticleDefinition* resi_pd = NULL;
    G4bool needOneMoreSec = false;
    G4ParticleDefinition* oneMoreSec_pd = NULL;
    if ( (int)(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) == 0 )
    {
       //All secondaries are already created;
       resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
    }
    else
    {
       if ( max_SecA > int(theBaseA + 1 - sum_A) ) 
       {
          //Most heavy secondary is interpreted as residual
          resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
          needOneMoreSec = true;
       }
       else 
       {
          //creation of residual is requierd 
          resi_pd = G4IonTable::GetIonTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
       }
 
       if ( needOneMoreSec )
       {
          if ( int(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) > 0 )
          {
             if ( int(theBaseA + 1 - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
             oneMoreSec_pd = G4Neutron::Neutron();
          }
          else 
          {
             oneMoreSec_pd = G4IonTable::GetIonTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
          }
       }
   
       if ( resi_pd == NULL )
       {
          // theNDLDataZ,A has the Z and A of used NDL file
          if ( (int)(theNDLDataZ - sum_Z) == 0 && (int)(theNDLDataA + 1 - 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 );
             for ( int i = 0 ; i < nSecondaries ; i++ )
             {
                if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ  
                  && theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
                {
                   G4ThreeVector p = theResult.GetSecondary( i )->GetParticle()->GetMomentum();
                   p = p * resi_pd->GetPDGMass()/ G4IonTable::GetIonTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass(); 
                   theResult.GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
                   theResult.GetSecondary( i )->GetParticle()->SetMomentum( p );
                } 
             }
          }
       }
    }
 
 
    G4LorentzVector secs_4p_lab( 0.0 );
 
    G4int n_sec = theResult.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.GetSecondary( i )->GetParticle()->Get4Momentum();
 
       G4double beta = 0;
       if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() != G4Gamma::Gamma() )
       {
          beta = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().beta(); 
       }
       else
       {
          beta = 1;
       }
 
       if ( theResult.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.GetSecondary( i )->GetParticle()->Get4Momentum().e();
             if ( e > theResult.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.AddSecondary ( res );    
 
       p4 = res->Get4Momentum(); 
       if ( slow > p4.beta() ) 
       {
          slow = p4.beta(); 
          islow = ires;
       }
       dif_4p = init_4p_lab - ( secs_4p_lab + p4 );  
    }
 
    if ( needOneMoreSec )
    {
       nSecondaries += 1;
       G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );    
       theResult.AddSecondary ( one );    
       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.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );    
 
       }
       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.GetSecondary( ires )->GetParticle()->Get4Momentum();
       theResult.GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() ); 
       dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.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.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );    
          theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );    
       }
       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.GetSecondary( ifast )->GetParticle()->GetKineticEnergy();
       G4ThreeVector p0 = theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
       G4ThreeVector dir = ( theResult.GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
 
       //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
 
       if ( ke0 + dif_e > 0 )
       {
          theResult.GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e ); 
          G4ThreeVector dp = p0 - theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();  
 
          G4ThreeVector p = theResult.GetSecondary( islow )->GetParticle()->GetMomentum();
          //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir ); 
          theResult.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;
       }
 
    }
 
 }

virtual G4HadFinalState* G4NeutronHPFinalState::ApplyYourself ( const G4HadProjectile & )

inlinevirtual

Definition at line 67 of file G4NeutronHPFinalState.hh.

Referenced by G4NeutronHPChannel::ApplyYourself().

   {
     throw G4HadronicException(__FILE__, __LINE__, "G4HadFinalState * ApplyYourself(const G4HadProjectile & theTrack) needs implementation");
     return 0;
   };