G4INCL::PhaseSpaceDecay Namespace Reference

Functions
void	decay (G4double initialMass, const ThreeVector &theBoostVector, ParticleList &particles)
	Generate decay momenta according to a uniform phase-space model. More...

Function Documentation

void G4INCL::PhaseSpaceDecay::decay	(	G4double	initialMass,
		const ThreeVector &	theBoostVector,
		ParticleList &	particles
	)

Generate decay momenta according to a uniform phase-space model.

This function will assign momenta to the particles in the list that is passed as an argument.

Parameters

initialMass	mass of the decaying system
theBoostVector	boost vector of the decaying system
particles	list of decay particles

Definition at line 62 of file G4INCLPhaseSpaceDecay.cc.

References G4INCL::Particle::adjustEnergyFromMomentum(), G4INCL::Particle::boost(), G4INCL::Particle::boostVector(), G4INCL::Particle::getMass(), G4INCL::KinematicsUtils::momentumInCM(), N, G4INCL::Random::normVector(), G4INCL::Particle::setMass(), and G4INCL::Particle::setMomentum().

                                                                                                       {

      size_t N = particles.size();
      std::vector<G4double> masses(N);
      std::vector<G4double> sumMasses(N);
      std::transform(particles.begin(), particles.end(), masses.begin(), std::mem_fun(&Particle::getMass));
      std::partial_sum(masses.begin(), masses.end(), sumMasses.begin());

      G4double PFragMagCM = 0.0;
      G4double T = initialMass-sumMasses.back();
// assert(T>-1.e-5);
      if(T<0.)
        T=0.;

      ThreeVector PFragCM, PRestCM, PRestLab;

      // The first particle in the list will pick up all the recoil
      Particle *restParticle = particles.front();
      restParticle->setMass(initialMass);
      restParticle->adjustEnergyFromMomentum();

      ThreeVector boostV;

      G4int k=N-1;
      for (ParticleList::reverse_iterator p=particles.rbegin(); k>0; ++p, --k) {
        const G4double mu = sumMasses[k-1];
        T *= (k>1) ? betaKopylov(k) : 0.;

        const G4double restMass = mu + T;

        PFragMagCM = KinematicsUtils::momentumInCM(restParticle->getMass(), masses[k], restMass);
        PFragCM = Random::normVector(PFragMagCM);
        (*p)->setMomentum(PFragCM);
        (*p)->adjustEnergyFromMomentum();
        restParticle->setMass(restMass);
        restParticle->setMomentum(-PFragCM);
        restParticle->adjustEnergyFromMomentum();

        (*p)->boost(boostV);
        (*p)->boost(theBoostVector);
        restParticle->boost(boostV);

        boostV = -restParticle->boostVector();
      }
      restParticle->setMass(masses[0]);
      restParticle->adjustEnergyFromMomentum();
      restParticle->boost(theBoostVector);
    }