G4PreCompoundTransitions Class Reference

#include <G4PreCompoundTransitions.hh>

Inheritance diagram for G4PreCompoundTransitions:

Public Member Functions
	G4PreCompoundTransitions ()
virtual	~G4PreCompoundTransitions ()
virtual G4double	CalculateProbability (const G4Fragment &aFragment)
virtual void	PerformTransition (G4Fragment &aFragment)
Public Member Functions inherited from G4VPreCompoundTransitions
	G4VPreCompoundTransitions ()
virtual	~G4VPreCompoundTransitions ()
G4double	GetTransitionProb1 () const
G4double	GetTransitionProb2 () const
G4double	GetTransitionProb3 () const
void	UseNGB (G4bool use)
void	UseCEMtr (G4bool use)

Additional Inherited Members
Protected Attributes inherited from G4VPreCompoundTransitions
G4bool	useNGB
G4bool	useCEMtr
G4double	TransitionProb1
G4double	TransitionProb2
G4double	TransitionProb3

Detailed Description

Definition at line 52 of file G4PreCompoundTransitions.hh.

Constructor & Destructor Documentation

G4PreCompoundTransitions::G4PreCompoundTransitions

(

)

Definition at line 57 of file G4PreCompoundTransitions.cc.

References G4PreCompoundParameters::GetFermiEnergy(), G4Pow::GetInstance(), G4PreCompoundParameters::GetLevelDensity(), G4PreCompoundParameters::GetTransitionsr0(), and G4Proton::Proton().

{
  proton = G4Proton::Proton();
  G4PreCompoundParameters param;
  FermiEnergy = param.GetFermiEnergy();
  r0 = param.GetTransitionsr0();
  aLDP = param.GetLevelDensity();
  g4pow = G4Pow::GetInstance();
}

G4PreCompoundTransitions::~G4PreCompoundTransitions ( )

virtual

Definition at line 67 of file G4PreCompoundTransitions.cc.

{}

Member Function Documentation

G4double G4PreCompoundTransitions::CalculateProbability ( const G4Fragment &  aFragment )

virtual

Implements G4VPreCompoundTransitions.

Definition at line 73 of file G4PreCompoundTransitions.cc.

References python.hepunit::c_light, python.hepunit::eV, G4UniformRand, GE, G4Fragment::GetA_asInt(), G4Fragment::GetExcitationEnergy(), G4Fragment::GetNumberOfCharged(), G4Fragment::GetNumberOfHoles(), G4Fragment::GetNumberOfParticles(), G4Fragment::GetZ_asInt(), python.hepunit::hbarc, N, python.hepunit::pi2, G4Pow::powN(), G4VPreCompoundTransitions::TransitionProb1, G4VPreCompoundTransitions::TransitionProb2, G4VPreCompoundTransitions::TransitionProb3, G4VPreCompoundTransitions::useCEMtr, G4VPreCompoundTransitions::useNGB, and test::x.

{
  // Number of holes
  G4int H = aFragment.GetNumberOfHoles();
  // Number of Particles 
  G4int P = aFragment.GetNumberOfParticles();
  // Number of Excitons 
  G4int N = P+H;
  // Nucleus 
  G4int A = aFragment.GetA_asInt();
  G4int Z = aFragment.GetZ_asInt();
  G4double U = aFragment.GetExcitationEnergy();

  //G4cout << aFragment << G4endl;
  
  if(U < 10*eV || 0==N) { return 0.0; }
  
  //J. M. Quesada (Feb. 08) new physics
  // OPT=1 Transitions are calculated according to Gudima's paper 
  //       (original in G4PreCompound from VL) 
  // OPT=2 Transitions are calculated according to Gupta's formulae
  //
  if (useCEMtr){

    // Relative Energy (T_{rel})
    G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
    
    // Sample kind of nucleon-projectile 
    G4bool ChargedNucleon(false);
    G4double chtest = 0.5;
    if (P > 0) { 
      chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P); 
    }
    if (G4UniformRand() < chtest) { ChargedNucleon = true; }
    
    // Relative Velocity: 
    // <V_{rel}>^2
    G4double RelativeVelocitySqr(0.0);
    if (ChargedNucleon) { 
      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2; 
    } else { 
      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2; 
    }
    
    // <V_{rel}>
    G4double RelativeVelocity = std::sqrt(RelativeVelocitySqr);
    
    // Proton-Proton Cross Section
    G4double ppXSection = 
      (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
      * CLHEP::millibarn;
    // Proton-Neutron Cross Section
    G4double npXSection = 
      (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
      * CLHEP::millibarn;
    
    // Averaged Cross Section: \sigma(V_{rel})
    //  G4double AveragedXSection = (ppXSection+npXSection)/2.0;
    G4double AveragedXSection(0.0);
    if (ChargedNucleon)
      {
        //JMQ: small bug fixed
        //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0);
        AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
      }
    else 
      {
        AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
        //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1);
      }
    
    // Fermi relative energy ratio
    G4double FermiRelRatio = FermiEnergy/RelativeEnergy;
    
    // This factor is introduced to take into account the Pauli principle
    G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
    if (FermiRelRatio > 0.5) {
      G4double x = 2.0 - 1.0/FermiRelRatio;
      PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
      //PauliFactor += 
      //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
    }
    // Interaction volume 
    //  G4double Vint = (4.0/3.0)
    //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
    G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
    //    G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx;
    G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
    
    // Transition probability for \Delta n = +2
    
    TransitionProb1 = AveragedXSection*PauliFactor
      *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint;

    //JMQ 281009  phenomenological factor in order to increase 
    //   equilibrium contribution
    //   G4double factor=5.0;
    //   TransitionProb1 *= factor;
    //
    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; } 
    
    // GE = g*E where E is Excitation Energy
    G4double GE = (6.0/pi2)*aLDP*A*U;
    
    //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
    G4double Fph = G4double(P*P+H*H+P-3*H)/4.0;
    
    G4bool NeverGoBack(false);
    if(useNGB) { NeverGoBack=true; }
    
    //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
    if (GE-Fph < 0.0) { NeverGoBack = true; }
    
    // F(p+1,h+1)
    G4double Fph1 = Fph + N/2.0;
    
    G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1);
    
    if (NeverGoBack)
      {
    TransitionProb2 = 0.0;
    TransitionProb3 = 0.0;
      }
    else 
      {
        // Transition probability for \Delta n = -2 (at F(p,h) = 0)
        TransitionProb2 = 
      TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph));
        if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; } 
        
        // Transition probability for \Delta n = 0 (at F(p,h) = 0)
        TransitionProb3 = TransitionProb1*(N+1)* ProbFactor  
      * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph));
        if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; }
      }
  } else {
    //JMQ: Transition probabilities from Gupta's work    
    // GE = g*E where E is Excitation Energy
    G4double GE = (6.0/pi2)*aLDP*A*U;
 
    G4double Kmfp=2.;
        
    //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
    TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1))
      /(8*Kmfp*CLHEP::c_light);
    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }

    TransitionProb2=0.;
    TransitionProb3=0.;
    
    if (!useNGB && N > 1) {
      // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);      
      TransitionProb2 = 
    3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE);      
      if (TransitionProb2 < 0.0) TransitionProb2 = 0.0; 
    }
  }
  //  G4cout<<"U = "<<U<<G4endl;
  //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
  //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2
  //   <<"  l0 ="<< TransitionProb3<<G4endl; 
  return TransitionProb1 + TransitionProb2 + TransitionProb3;
}

void G4PreCompoundTransitions::PerformTransition ( G4Fragment &  aFragment )

virtual

Implements G4VPreCompoundTransitions.

Definition at line 237 of file G4PreCompoundTransitions.cc.

References G4UniformRand, G4Fragment::GetA_asInt(), G4Fragment::GetNumberOfCharged(), G4Fragment::GetNumberOfHoles(), G4Fragment::GetNumberOfParticles(), G4Fragment::GetZ_asInt(), G4INCL::Math::max(), G4Fragment::SetNumberOfCharged(), G4Fragment::SetNumberOfHoles(), G4Fragment::SetNumberOfParticles(), G4VPreCompoundTransitions::TransitionProb1, G4VPreCompoundTransitions::TransitionProb2, and G4VPreCompoundTransitions::TransitionProb3.

{
  G4double ChosenTransition = 
    G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
  G4int deltaN = 0;
  //  G4int Nexcitons = result.GetNumberOfExcitons();
  G4int Npart     = result.GetNumberOfParticles();
  G4int Ncharged  = result.GetNumberOfCharged();
  G4int Nholes    = result.GetNumberOfHoles();
  if (ChosenTransition <= TransitionProb1) 
    {
      // Number of excitons is increased on \Delta n = +2
      deltaN = 2;
    } 
  else if (ChosenTransition <= TransitionProb1+TransitionProb2) 
    {
      // Number of excitons is increased on \Delta n = -2
      deltaN = -2;
    }

  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and  
  // in proportion to the number charges w.r.t. number of particles,  
  // PROVIDED that there are charged particles
  deltaN /= 2;

  //G4cout << "deltaN= " << deltaN << G4endl;

  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on 
  // number of charged vs. number of particles fails
  result.SetNumberOfParticles(Npart+deltaN);
  result.SetNumberOfHoles(Nholes+deltaN); 

  if(deltaN < 0) {
    if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged) 
      { 
    result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
      }

  } else if ( deltaN > 0 ) {
    // With weight Z/A, number of charged particles is increased with +1
    G4int A = result.GetA_asInt();
    G4int Z = result.GetZ_asInt();
    if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z) 
      {
    result.SetNumberOfCharged(Ncharged+deltaN);
      }
  }
  
  // Number of charged can not be greater that number of particles
  if ( Npart < Ncharged ) 
    {
      result.SetNumberOfCharged(Npart);
    }
  //G4cout << "### After transition" << G4endl;
  //G4cout << result << G4endl;
}

---

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

• G4PreCompoundTransitions.hh
• G4PreCompoundTransitions.cc