G4NRESP71M03 Class Reference

#include <G4NRESP71M03.hh>

Public Member Functions
	G4NRESP71M03 ()
	~G4NRESP71M03 ()
void	DKINMA (G4ReactionProduct *p1, G4ReactionProduct *p2, G4ReactionProduct *p3, G4ReactionProduct *p4, const G4double Q, const G4double costhcm3)
G4int	ApplyMechanismI_NBeA2A (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
G4int	ApplyMechanismII_ACN2A (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
G4int	ApplyMechanismABE (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds)

Detailed Description

Definition at line 36 of file G4NRESP71M03.hh.

Constructor & Destructor Documentation

G4NRESP71M03()

G4NRESP71M03::G4NRESP71M03 ( )

inline

Definition at line 39 of file G4NRESP71M03.hh.

{ ; };

~G4NRESP71M03()

G4NRESP71M03::~G4NRESP71M03 ( )

inline

Definition at line 40 of file G4NRESP71M03.hh.

{ ; };

Member Function Documentation

ApplyMechanismABE()

G4int G4NRESP71M03::ApplyMechanismABE	(	G4ReactionProduct &	neut,
		G4ReactionProduct &	carb,
		G4ReactionProduct *	theProds )

Definition at line 350 of file G4NRESP71M03.cc.

{
  G4double CosThetaCMAlpha = 0.;  // Cosine of the angle of emission of the alpha particle (theta).
 
  G4double Kn = neut.GetKineticEnergy();  // Neutron energy in the center of mass system.
  if (Kn > 5.7 * MeV) {
    // Sorting.
    for (G4int i = 1; i < ndist; i++) {
      if (BEN2[i] >= Kn / keV) {
        // Ok. The neutron energy is between values i-1 and i of BEN2. Taking them.
        G4double BE1 = BEN2[i - 1];
        G4double BE2 = BEN2[i];
 
        // Performing energy and angle interpolation.
 
        G4double FRA =
          G4UniformRand()
          * 49.99999999;  // Sorting the bin of the cumulative probability FRA (Rho). There are 51
                          // values of Rho from 0 to 1 (0; 0.02; 0.04 ... 1).
        G4double DJTETA =
          FRA
          - G4int(FRA);  // Distance in bin units (DJTETA) from the low edge of the bin with Rho.
        G4int JTETA =
          G4int(FRA) + 1;  // Getting the bin (JTETA) next to the bin with the value of Rho.
 
        // Calculating the angle from the cumulative distribution at energy i.
 
        G4double B11 = B2[i - 1][JTETA - 1];
        G4double B12 = B2[i - 1][JTETA];
 
        G4double TCM1 = B11 + (B12 - B11) * DJTETA;  // Angle at energy i.
 
        // Calculating the angle from the cumulative distribution at energy i-1.
 
        G4double B21 = B2[i][JTETA - 1];
        G4double B22 = B2[i][JTETA];
 
        G4double TCM2 = B21 + (B22 - B21) * DJTETA;  // Angle at energy i-1.
 
        // Interpolating the angle between energy values i and i-1.
        G4double TCM = (TCM1 + (TCM2 - TCM1) * (Kn / keV - BE1) / (BE2 - BE1)) * 1.E-4;
        CosThetaCMAlpha = std::cos(TCM);
 
        break;
      }
    }
  }
  else {
    // Isotropic distribution in CM.
    CosThetaCMAlpha = 1. - 2. * G4UniformRand();
  }
 
  // N+12C --> A+9BE
  theProds[0].SetDefinition(G4Alpha::Alpha());
  theProds[1].SetDefinition(G4IonTable::GetIonTable()->GetIon(4, 9, 0.));
 
  DKINMA(&neut, &carb, &(theProds[0]), &(theProds[1]), -5.71 * MeV, CosThetaCMAlpha);
 
  return 0;
}

ApplyMechanismI_NBeA2A()

G4int G4NRESP71M03::ApplyMechanismI_NBeA2A	(	G4ReactionProduct &	neut,
		G4ReactionProduct &	carb,
		G4ReactionProduct *	theProds,
		const G4double	QI )

Definition at line 291 of file G4NRESP71M03.cc.

{
  // N+12C --> A+9BE*
  G4ReactionProduct p4;
 
  theProds[0].SetDefinition(G4Alpha::Alpha());
 
  DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2. * G4UniformRand() - 1.);
 
  // 9BE* --> N+8BE
  G4ReactionProduct p1(p4);
 
  theProds[1].SetDefinition(G4Neutron::Neutron());
 
  DKINMA(&p1, nullptr, &(theProds[1]), &p4, -QI - 7.369, 2. * G4UniformRand() - 1.);
 
  // 8BE --> 2*A
  p1 = p4;
 
  theProds[2].SetDefinition(G4Alpha::Alpha());
  theProds[3].SetDefinition(G4Alpha::Alpha());
 
  DKINMA(&p1, nullptr, &(theProds[2]), &(theProds[3]), 0.09538798439007223351,
         2. * G4UniformRand() - 1.);
 
  return 0;
}

ApplyMechanismII_ACN2A()

G4int G4NRESP71M03::ApplyMechanismII_ACN2A	(	G4ReactionProduct &	neut,
		G4ReactionProduct &	carb,
		G4ReactionProduct *	theProds,
		const G4double	QI )

Definition at line 321 of file G4NRESP71M03.cc.

{
  // 12C(N,N')12C'
  G4ReactionProduct p4;
 
  theProds[0].SetDefinition(G4Neutron::Neutron());
 
  DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2. * G4UniformRand() - 1.);
 
  // 12C' --> ALPHA+8BE'
  G4ReactionProduct p1(p4);
 
  theProds[1].SetDefinition(G4Alpha::Alpha());
 
  DKINMA(&p1, nullptr, &(theProds[1]), &p4, -QI - 7.369, 2. * G4UniformRand() - 1.);
 
  // 8BE --> 2*ALPHA
  p1 = p4;
 
  theProds[2].SetDefinition(G4Alpha::Alpha());
  theProds[3].SetDefinition(G4Alpha::Alpha());
 
  DKINMA(&p1, nullptr, &(theProds[2]), &(theProds[3]), 0.09538798439007223351,
         2. * G4UniformRand() - 1.);
 
  return 0;
}

DKINMA()

void G4NRESP71M03::DKINMA	(	G4ReactionProduct *	p1,
		G4ReactionProduct *	p2,
		G4ReactionProduct *	p3,
		G4ReactionProduct *	p4,
		const G4double	Q,
		const G4double	costhcm3 )

Definition at line 206 of file G4NRESP71M03.cc.

{
  G4ReactionProduct CM;
  G4double TotalEnergyCM;
 
  if (p2 != nullptr)  // If it is not a decay reaction...
  {
    // Calculating (total momentum, energy and mass) of the center of mass.
    const G4ThreeVector TotalMomentumLAB = p1->GetMomentum() + p2->GetMomentum();
    CM.SetMomentum(TotalMomentumLAB);
 
    const G4double TotalEnergyLAB = p1->GetTotalEnergy() + p2->GetTotalEnergy();
    CM.SetTotalEnergy(TotalEnergyLAB);
 
    CM.SetMass(std::sqrt(TotalEnergyLAB * TotalEnergyLAB - TotalMomentumLAB * TotalMomentumLAB));
 
    // Transforming primary particles from laboratory to center of mass system.
    p1->Lorentz(*p1, CM);
    p2->Lorentz(*p2, CM);
 
    TotalEnergyCM = p1->GetTotalEnergy() + p2->GetTotalEnergy();
 
    const G4double m4 =
      (p1->GetMass() + p2->GetMass())
      - (p3->GetMass()
         + Q);  // Mass of the residual nucleus in the excited state (not in the ground state).
    p4->SetMass(m4);
  }
  else  // If it is a decay reaction...
  {
    const G4ThreeVector TotalMomentumLAB = p1->GetMomentum();
    CM.SetMomentum(TotalMomentumLAB);
 
    const G4double TotalEnergyLAB = p1->GetTotalEnergy();
    CM.SetTotalEnergy(TotalEnergyLAB);
 
    CM.SetMass(std::sqrt(TotalEnergyLAB * TotalEnergyLAB - TotalMomentumLAB * TotalMomentumLAB));
 
    // Transforming primary particles from laboratory to center of mass system (not really necessary
    // in this case).
    p1->Lorentz(*p1, CM);
 
    const G4double m4 =
      p1->GetMass()
      - (p3->GetMass()
         + Q);  // Mass of the residual nucleus in the excited state (not in the ground state).
    p4->SetMass(m4);
 
    TotalEnergyCM = p1->GetTotalEnergy();
  }
 
  // Calculating momentum and total energy of the reaction products.
 
  const G4ThreeVector p1unit = p1->GetMomentum().unit();
 
  G4RotationMatrix rot(std::acos(p1unit * G4ThreeVector(0., 1., 0.)),
                       std::acos(p1unit * G4ThreeVector(0., 0., 1.)), 0.);
  rot.inverse();
 
  const G4double theta = std::acos(costhcm3);
  const G4double phi = twopi * G4UniformRand();
 
  const G4double Energy3CM =
    (std::pow(TotalEnergyCM, 2.) + std::pow(p3->GetMass(), 2.) - std::pow(p4->GetMass(), 2.))
    / (2. * TotalEnergyCM);
  p3->SetTotalEnergy(Energy3CM);
 
  const G4double Momentum3CM = std::sqrt(std::pow(Energy3CM, 2.) - std::pow(p3->GetMass(), 2.));
  p3->SetMomentum(rot
                  * G4ThreeVector(Momentum3CM * std::sin(theta) * std::cos(phi),
                                  Momentum3CM * std::sin(theta) * std::sin(phi),
                                  Momentum3CM * costhcm3));
 
  const G4double Energy4CM = TotalEnergyCM - Energy3CM;
  p4->SetTotalEnergy(Energy4CM);
 
  const G4double Momentum4CM = std::sqrt(std::pow(Energy4CM, 2.) - std::pow(p4->GetMass(), 2.));
  p4->SetMomentum(-Momentum4CM * p3->GetMomentum().unit());
 
  // Transforming reaction products from center of mass to laboratory system.
  p3->Lorentz(*p3, -1. * CM);
  p4->Lorentz(*p4, -1. * CM);
}

Referenced by ApplyMechanismABE(), ApplyMechanismI_NBeA2A(), and ApplyMechanismII_ACN2A().

---

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

• G4NRESP71M03.hh
• G4NRESP71M03.cc