#include <G4GDecay3.hh>

Public Member Functions
	G4GDecay3 ()

	G4GDecay3 (const G4double &pMass, const G4double &dMass0, const G4double &dMass1, const G4double &dMass2)

std::vector< G4ThreeVector >	GetThreeBodyMomenta ()

	~G4GDecay3 ()

Private Member Functions
G4bool	CalculateMomentumMagnitudes ()

Private Attributes
G4int	loopMax

G4double	mDaughter0

G4double	mDaughter1

G4double	mDaughter2

G4double	parentMass

G4double	pDaughter0

G4double	pDaughter1

G4double	pDaughter2

Detailed Description

Definition at line 42 of file G4GDecay3.hh.

Constructor & Destructor Documentation

◆ G4GDecay3() [1/2]

G4GDecay3::G4GDecay3 ( )

inline

Definition at line 45 of file G4GDecay3.hh.

45{;}

◆ G4GDecay3() [2/2]

G4GDecay3::G4GDecay3	(	const G4double &	pMass,
		const G4double &	dMass0,
		const G4double &	dMass1,
		const G4double &	dMass2
	)

Definition at line 42 of file G4GDecay3.cc.

 : loopMax(100), mDaughter0(dMass0), mDaughter1(dMass1),
   mDaughter2(dMass2), pDaughter0(0.), pDaughter1(0.), pDaughter2(0.) 
{
  parentMass = std::max(pMass, mDaughter0 + mDaughter1 + mDaughter2 + CLHEP::keV);
}

References CLHEP::keV, G4INCL::Math::max(), mDaughter0, mDaughter1, mDaughter2, and parentMass.

◆ ~G4GDecay3()

G4GDecay3::~G4GDecay3 ( )

inline

Definition at line 48 of file G4GDecay3.hh.

48{;}

Member Function Documentation

◆ CalculateMomentumMagnitudes()

G4bool G4GDecay3::CalculateMomentumMagnitudes ( )

private

Definition at line 50 of file G4GDecay3.cc.

{
  G4double rndm;
  G4double rndm1;
  G4double rndm2;
 
  G4double momentummax;
  G4double momentumsum;
  G4double energy;
 
  G4double availableE = parentMass - mDaughter0 - mDaughter1 - mDaughter2;
  do {              /* Loop checking 08.06.2015 MHK */
    rndm1 = G4UniformRand();
    rndm2 = G4UniformRand();
    if (rndm2 > rndm1) {
      // keep randoms in descending order 
      rndm = rndm1;
      rndm1 = rndm2;
      rndm2 = rndm;
    }
    momentummax = 0.0;
    momentumsum = 0.0;
 
    // daughter 0
    energy = rndm2*availableE;
    pDaughter0 = std::sqrt(energy*energy + 2.0*energy*mDaughter0);
    if (pDaughter0 > momentummax) momentummax = pDaughter0;
    momentumsum += pDaughter0;
 
    // daughter 1
    energy = (1.-rndm1)*availableE;
    pDaughter1 = std::sqrt(energy*energy + 2.0*energy*mDaughter1);
    if (pDaughter1 > momentummax) momentummax = pDaughter1;
    momentumsum += pDaughter1;
 
    // daughter 2
    energy = (rndm1-rndm2)*availableE;
    pDaughter2 = std::sqrt(energy*energy + 2.0*energy*mDaughter2);
    if (pDaughter2 > momentummax) momentummax = pDaughter2;
    momentumsum += pDaughter2;
  } while (momentummax > momentumsum - momentummax);
 
  return true;
}

References G4INCL::KinematicsUtils::energy(), G4UniformRand, mDaughter0, mDaughter1, mDaughter2, parentMass, pDaughter0, pDaughter1, and pDaughter2.

Referenced by GetThreeBodyMomenta().

◆ GetThreeBodyMomenta()

std::vector< G4ThreeVector > G4GDecay3::GetThreeBodyMomenta ( )

Definition at line 96 of file G4GDecay3.cc.

{
 
  std::vector<G4ThreeVector> pVect;
 
  if (CalculateMomentumMagnitudes() ) {
    
    // Calculate directions
    G4double costheta = 2.*G4UniformRand()-1.;
    G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
    G4double phi = twopi*G4UniformRand();
    G4double sinphi = std::sin(phi);
    G4double cosphi = std::cos(phi);
    G4ThreeVector direction0(sintheta*cosphi, sintheta*sinphi, costheta);
 
    G4double costhetan = (pDaughter1*pDaughter1 - pDaughter2*pDaughter2
                          - pDaughter0*pDaughter0)/(2.0*pDaughter2*pDaughter0);
    G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
    G4double phin = twopi*G4UniformRand();
    G4double sinphin = std::sin(phin);
    G4double cosphin = std::cos(phin);
    G4ThreeVector direction2;
    direction2.setX(sinthetan*cosphin*costheta*cosphi -
                    sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
    direction2.setY(sinthetan*cosphin*costheta*sinphi +
                    sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
    direction2.setZ(-sinthetan*cosphin*sintheta + costhetan*costheta);
 
    // Return momentum vectors
    pVect.push_back(pDaughter0*direction0);
    pVect.push_back(-direction0*pDaughter0 - direction2*pDaughter2);
    pVect.push_back(pDaughter2*direction2);
 
  } else {
    G4cerr << "G4GDecay3::GetThreeBodyMomenta: " << loopMax
           << " or more loops in momentum magnitude calculation " << G4endl;
  }
 
  return pVect;
}