#include <G4FermiPhaseSpaceDecay.hh>

Public Member Functions
	G4FermiPhaseSpaceDecay ()

	~G4FermiPhaseSpaceDecay ()

std::vector< G4LorentzVector * > *	Decay (G4double parent_mass, const std::vector< G4double > &fragment_masses) const

Detailed Description

Definition at line 49 of file G4FermiPhaseSpaceDecay.hh.

Constructor & Destructor Documentation

◆ G4FermiPhaseSpaceDecay()

G4FermiPhaseSpaceDecay::G4FermiPhaseSpaceDecay ( )

Definition at line 46 of file G4FermiPhaseSpaceDecay.cc.

{
  g4calc = G4Pow::GetInstance();
}

◆ ~G4FermiPhaseSpaceDecay()

G4FermiPhaseSpaceDecay::~G4FermiPhaseSpaceDecay ( )

Definition at line 51 of file G4FermiPhaseSpaceDecay.cc.

52{}

Member Function Documentation

◆ Decay()

std::vector< G4LorentzVector * > * G4FermiPhaseSpaceDecay::Decay	(	G4double	parent_mass,
		const std::vector< G4double > &	fragment_masses
	)		const

Definition at line 54 of file G4FermiPhaseSpaceDecay.cc.

{
  std::size_t N = mr.size();
 
  std::vector<G4LorentzVector*>* P = 
    new std::vector<G4LorentzVector*>(N, nullptr);
 
  G4double mtot = 0.0;
  for(std::size_t k=0; k<N; ++k) { mtot += mr[k]; }
 
  G4double mu = mtot;
  G4double PFragMagCM = 0.0;
 
  // Primary mass is above the sum of mass of components
  G4double Mass = std::max(M, mtot + CLHEP::eV);
  G4double T = Mass-mtot;
 
  G4LorentzVector PFragCM(0.0,0.0,0.0,0.0);
  G4LorentzVector PRestCM(0.0,0.0,0.0,0.0);
  G4LorentzVector PRestLab(0.0,0.0,0.0,Mass);
 
  CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
 
  for (G4int k = (G4int)N-1; k>0; --k)
    {
      mu -= mr[k];
      if (k>1) { T *= BetaKopylov(k, rndmEngine); }
      else     { T = 0.0; }
      
      G4double RestMass = mu + T;
 
      PFragMagCM = PtwoBody(Mass,mr[k],RestMass);
      
      // Create a unit vector with a random direction isotropically distributed
      G4ThreeVector RandVector = PFragMagCM*G4RandomDirection();
 
      PFragCM.setVect(RandVector);
      PFragCM.setE(std::sqrt(PFragMagCM*PFragMagCM + mr[k]*mr[k]));
 
      PRestCM.setVect(-RandVector);
      PRestCM.setE(std::sqrt(PFragMagCM*PFragMagCM + RestMass*RestMass));
 
      G4ThreeVector BoostV = PRestLab.boostVector();
 
      PFragCM.boost(BoostV);
      (*P)[k] = new G4LorentzVector(PFragCM);
 
      PRestCM.boost(BoostV);
      PRestLab = PRestCM;
      
      Mass = RestMass;
    }
 
  (*P)[0] = new G4LorentzVector(PRestLab);
 
  return P;
}