#include <G4CascadeFinalStateAlgorithm.hh>

Inheritance diagram for G4CascadeFinalStateAlgorithm:

Public Member Functions
	G4CascadeFinalStateAlgorithm ()

virtual	~G4CascadeFinalStateAlgorithm ()

virtual void	SetVerboseLevel (G4int verbose)

void	Configure (G4InuclElementaryParticle bullet, G4InuclElementaryParticle target, const std::vector< G4int > &particle_kinds)

Public Member Functions inherited from G4VHadDecayAlgorithm
	G4VHadDecayAlgorithm (const G4String &algName, G4int verbose=0)

virtual	~G4VHadDecayAlgorithm ()

void	Generate (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

G4int	GetVerboseLevel () const

const G4String &	GetName () const

Protected Member Functions
virtual void	GenerateTwoBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

virtual void	GenerateMultiBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

void	SaveKinematics (G4InuclElementaryParticle bullet, G4InuclElementaryParticle target)

void	ChooseGenerators (G4int is, G4int fs)

void	FillMagnitudes (G4double initialMass, const std::vector< G4double > &masses)

G4bool	satisfyTriangle (const std::vector< G4double > &pmod) const

void	FillDirections (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

void	FillDirThreeBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

void	FillDirManyBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

G4double	GenerateCosTheta (G4int ptype, G4double pmod) const

void	FillUsingKopylov (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)

G4double	BetaKopylov (G4int K) const

Protected Member Functions inherited from G4VHadDecayAlgorithm
virtual G4bool	IsDecayAllowed (G4double initialMass, const std::vector< G4double > &masses) const

G4double	TwoBodyMomentum (G4double M0, G4double M1, G4double M2) const

G4double	UniformTheta () const

G4double	UniformPhi () const

void	PrintVector (const std::vector< G4double > &v, const G4String &name, std::ostream &os) const

Detailed Description

Definition at line 48 of file G4CascadeFinalStateAlgorithm.hh.

Constructor & Destructor Documentation

◆ G4CascadeFinalStateAlgorithm()

G4CascadeFinalStateAlgorithm::G4CascadeFinalStateAlgorithm ( )

Definition at line 75 of file G4CascadeFinalStateAlgorithm.cc.

76 : G4VHadDecayAlgorithm("G4CascadeFinalStateAlgorithm"),

77 momDist(0), angDist(0), multiplicity(0), bullet_ekin(0.) {;}

G4VHadDecayAlgorithm::G4VHadDecayAlgorithm

G4VHadDecayAlgorithm(const G4String &algName, G4int verbose=0)

Definition G4VHadDecayAlgorithm.hh:44

◆ ~G4CascadeFinalStateAlgorithm()

G4CascadeFinalStateAlgorithm::~G4CascadeFinalStateAlgorithm ( )

virtual

Definition at line 79 of file G4CascadeFinalStateAlgorithm.cc.

79{;}

Member Function Documentation

◆ BetaKopylov()

G4double G4CascadeFinalStateAlgorithm::BetaKopylov ( G4int K ) const

protected

Definition at line 507 of file G4CascadeFinalStateAlgorithm.cc.

                                                                {
  G4Pow* g4pow = G4Pow::GetInstance();
 
  G4int N = 3*K - 5;
  G4double xN = G4double(N);
  G4double Fmax = std::sqrt(g4pow->powN(xN/(xN+1.),N)/(xN+1.)); 
 
  G4double F, chi;
  do {                  /* Loop checking 08.06.2015 MHK */
    chi = G4UniformRand();
    F = std::sqrt(g4pow->powN(chi,N)*(1.-chi));      
  } while ( Fmax*G4UniformRand() > F);
  return chi;
}

Referenced by FillUsingKopylov().

◆ ChooseGenerators()

void G4CascadeFinalStateAlgorithm::ChooseGenerators	(	G4int	is,
		G4int	fs )

protected

Definition at line 136 of file G4CascadeFinalStateAlgorithm.cc.

                                                                      {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::ChooseGenerators"
       << " is " << is << " fs " << fs << G4endl;
 
  // Get generators for momentum and angle
  if (G4CascadeParameters::usePhaseSpace()) momDist = 0;
  else momDist = G4MultiBodyMomentumDist::GetDist(is, multiplicity);
 
  if (fs > 0 && multiplicity == 2) {
    G4int kw = (fs==is) ? 1 : 2;
    angDist = G4TwoBodyAngularDist::GetDist(is, fs, kw);
  } else if (multiplicity == 3) {
    angDist = G4TwoBodyAngularDist::GetDist(is);
  } else {
    angDist = 0;
  }
 
  if (GetVerboseLevel()>1) {
    G4cout << " " << (momDist?momDist->GetName().c_str():"") << " "
       << (angDist?angDist->GetName().c_str():"") << G4endl;
  }
}

Referenced by Configure().

◆ Configure()

void G4CascadeFinalStateAlgorithm::Configure	(	G4InuclElementaryParticle *	bullet,
		G4InuclElementaryParticle *	target,
		const std::vector< G4int > &	particle_kinds )

Definition at line 91 of file G4CascadeFinalStateAlgorithm.cc.

                                              {
  if (GetVerboseLevel()>1)
    G4cout << " >>> " << GetName() << "::Configure" << G4endl;
 
  // Identify initial and final state (if two-body) for algorithm selection
  multiplicity = (G4int)particle_kinds.size();
  G4int is = bullet->type() * target->type();
  G4int fs = (multiplicity==2) ? particle_kinds[0]*particle_kinds[1] : 0;
 
  ChooseGenerators(is, fs);
 
  // Save kinematics for use with distributions
  SaveKinematics(bullet, target);
 
  // Save particle types for use with distributions
  kinds = particle_kinds;
}

Referenced by G4CascadeFinalStateGenerator::Configure().

◆ FillDirections()

void G4CascadeFinalStateAlgorithm::FillDirections	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protected

Definition at line 312 of file G4CascadeFinalStateAlgorithm.cc.

                                                   {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::FillDirections" << G4endl;
 
  finalState.clear();           // Initialization and sanity check
  if ((G4int)modules.size() != multiplicity) return;
 
  // Different order of processing for three vs. N body
  if (multiplicity == 3)
    FillDirThreeBody(initialMass, masses, finalState);
  else
    FillDirManyBody(initialMass, masses, finalState);
}

Referenced by GenerateMultiBody().

◆ FillDirManyBody()

void G4CascadeFinalStateAlgorithm::FillDirManyBody	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protected

Definition at line 361 of file G4CascadeFinalStateAlgorithm.cc.

                                                {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::FillDirManyBody" << G4endl;
 
  // Fill all but the last two particles randomly
  G4double costh = 0.;
 
  finalState.resize(multiplicity);
 
  for (G4int i=0; i<multiplicity-2; i++) {
    costh = GenerateCosTheta(kinds[i], modules[i]);
    finalState[i] = generateWithFixedTheta(costh, modules[i], masses[i]);
    finalState[i] = toSCM.rotate(finalState[i]);    // Align target axis
  }
 
  // Total momentum so far, to compute recoil of last two particles
  G4LorentzVector psum =
    std::accumulate(finalState.begin(), finalState.end()-2, G4LorentzVector());
  G4double pmod = psum.rho();
 
  costh = -0.5 * (pmod*pmod +
          modules[multiplicity-2]*modules[multiplicity-2] -
          modules[multiplicity-1]*modules[multiplicity-1])
    / pmod / modules[multiplicity-2];
 
  if (GetVerboseLevel() > 2) G4cout << " ct last " << costh << G4endl;
 
  if (std::fabs(costh) >= maxCosTheta) {  // Bad kinematics; abort generation
    finalState.clear();
    return;
  }
 
  // Report success
  if (GetVerboseLevel()>2) G4cout << " ok for mult " << multiplicity << G4endl;
 
  // First particle is at fixed angle to recoil system
  finalState[multiplicity-2] =
    generateWithFixedTheta(costh, modules[multiplicity-2],
               masses[multiplicity-2]);
  finalState[multiplicity-2] = toSCM.rotate(psum, finalState[multiplicity-2]);
 
  // Remaining particle is constrained to recoil from entire rest of system
  finalState[multiplicity-1].set(0.,0.,0.,initialMass);
  finalState[multiplicity-1] -= psum + finalState[multiplicity-2];
}

Referenced by FillDirections().

◆ FillDirThreeBody()

void G4CascadeFinalStateAlgorithm::FillDirThreeBody	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protected

Definition at line 328 of file G4CascadeFinalStateAlgorithm.cc.

                                                 {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::FillDirThreeBody" << G4endl;
 
  finalState.resize(3);
 
  G4double costh = GenerateCosTheta(kinds[2], modules[2]);
  finalState[2] = generateWithFixedTheta(costh, modules[2], masses[2]);
  finalState[2] = toSCM.rotate(finalState[2]);  // Align target axis
 
  // Generate direction of first particle
  costh = -0.5 * (modules[2]*modules[2] + modules[0]*modules[0] -
          modules[1]*modules[1]) / modules[2] / modules[0];
 
  if (std::fabs(costh) >= maxCosTheta) {  // Bad kinematics; abort generation
    finalState.clear();
    return;
  }
 
  // Report success
  if (GetVerboseLevel()>2) G4cout << " ok for mult 3" << G4endl;
 
  // First particle is at fixed angle to recoil system
  finalState[0] = generateWithFixedTheta(costh, modules[0], masses[0]);
  finalState[0] = toSCM.rotate(finalState[2], finalState[0]);
 
  // Remaining particle is constrained to recoil from entire rest of system
  finalState[1].set(0.,0.,0.,initialMass);
  finalState[1] -= finalState[0] + finalState[2];
}

Referenced by FillDirections().

◆ FillMagnitudes()

void G4CascadeFinalStateAlgorithm::FillMagnitudes	(	G4double	initialMass,
		const std::vector< G4double > &	masses )

protected

Definition at line 227 of file G4CascadeFinalStateAlgorithm.cc.

                                                                        {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::FillMagnitudes" << G4endl;
 
  modules.clear();      // Initialization and sanity checks
  if (!momDist) return;
 
  modules.resize(multiplicity,0.);  // Pre-allocate to avoid resizing
 
  G4double mass_last = masses.back();
  G4double pmod = 0.;
 
  if (GetVerboseLevel() > 3){
    G4cout << " knd_last " << kinds.back() << " mass_last " 
           << mass_last << G4endl;
  }
 
  G4int itry = -1;
  while (++itry < itry_max) {       /* Loop checking 08.06.2015 MHK */
    if (GetVerboseLevel() > 3) {
      G4cout << " itry in fillMagnitudes " << itry << G4endl;
    }
 
    G4double eleft = initialMass;
 
    G4int i;    // For access outside of loop
    for (i=0; i < multiplicity-1; i++) {
      pmod = momDist->GetMomentum(kinds[i], bullet_ekin);
 
      if (pmod < small) break;
      eleft -= std::sqrt(pmod*pmod + masses[i]*masses[i]);
 
      if (GetVerboseLevel() > 3) {
    G4cout << " kp " << kinds[i] << " pmod " << pmod
           << " mass2 " << masses[i]*masses[i] << " eleft " << eleft
           << "\n x1 " << eleft - mass_last << G4endl;
      }
 
      if (eleft <= mass_last) break;
 
      modules[i] = pmod;
    }
 
    if (i < multiplicity-1) continue;   // Failed to generate full kinematics
 
    G4double plast = eleft * eleft - masses.back()*masses.back();
    if (GetVerboseLevel() > 2) G4cout << " plast ** 2 " << plast << G4endl;
    
    if (plast <= small) continue;   // Not enough momentum left over
 
    plast = std::sqrt(plast);       // Final momentum is what's left over
    modules.back() = plast;
    
    if (multiplicity > 3 || satisfyTriangle(modules)) break;    // Successful
  } // while (itry < itry_max)
 
  if (itry >= itry_max) {       // Too many attempts
    if (GetVerboseLevel() > 2)
      G4cerr << " Unable to generate momenta for multiplicity "
         << multiplicity << G4endl;
 
    modules.clear();        // Something went wrong, throw away partial
  }
}

Referenced by GenerateMultiBody().

◆ FillUsingKopylov()

void G4CascadeFinalStateAlgorithm::FillUsingKopylov	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protected

Definition at line 463 of file G4CascadeFinalStateAlgorithm.cc.

                                                 {
  if (GetVerboseLevel()>2)
    G4cout << " >>> " << GetName() << "::FillUsingKopylov" << G4endl;
 
  finalState.clear();
 
  std::size_t N = masses.size();
  finalState.resize(N);
 
  G4double mtot = std::accumulate(masses.begin(), masses.end(), 0.0);
  G4double mu = mtot;
  G4double Mass = initialMass;
  G4double T = Mass-mtot;
  G4double recoilMass = 0.0;
  G4ThreeVector momV, boostV;       // Buffers to reduce memory churn
  G4LorentzVector recoil(0.0,0.0,0.0,Mass);
 
  for (std::size_t k=N-1; k>0; --k) {
    mu -= masses[k];
    T *= (k>1) ? BetaKopylov((G4int)k) : 0.;
    
    recoilMass = mu + T;
 
    boostV = recoil.boostVector();  // Previous system's rest frame
    
    // Create momentum with a random direction isotropically distributed
    // FIXME:  Should theta distribution should use Bertini fit function?
    momV.setRThetaPhi(TwoBodyMomentum(Mass,masses[k],recoilMass),
              UniformTheta(), UniformPhi());
    
    finalState[k].setVectM(momV,masses[k]);
    recoil.setVectM(-momV,recoilMass);
 
    finalState[k].boost(boostV);
    recoil.boost(boostV);
    Mass = recoilMass;
  }
  
  finalState[0] = recoil;
}

Referenced by GenerateMultiBody().

◆ GenerateCosTheta()

G4double G4CascadeFinalStateAlgorithm::GenerateCosTheta	(	G4int	ptype,
		G4double	pmod ) const

protected

Definition at line 412 of file G4CascadeFinalStateAlgorithm.cc.

                                                   {
  if (GetVerboseLevel() > 2) {
    G4cout << " >>> " << GetName() << "::GenerateCosTheta " << ptype
       << " " << pmod << G4endl;
  }
 
  if (multiplicity == 3) {      // Use distribution for three-body
    return angDist->GetCosTheta(bullet_ekin, ptype);
  }
 
  // Throw multi-body distribution
  G4double p0 = ptype<3 ? 0.36 : 0.25;  // Nucleon vs. everything else
  G4double alf = 1.0 / p0 / (p0 - (pmod+p0)*G4Exp(-pmod / p0));
 
  G4double sinth = 2.0;
 
  G4int itry1 = -1;     /* Loop checking 08.06.2015 MHK */
  while (std::fabs(sinth) > maxCosTheta && ++itry1 < itry_max) {
    G4double s1 = pmod * inuclRndm();
    G4double s2 = alf * oneOverE * p0 * inuclRndm();
    G4double salf = s1 * alf * G4Exp(-s1 / p0);
    if (GetVerboseLevel() > 3) {
      G4cout << " s1 * alf * G4Exp(-s1 / p0) " << salf
         << " s2 " << s2 << G4endl;
    }
    
    if (salf > s2) sinth = s1 / pmod;
  }
  
  if (GetVerboseLevel() > 3)
    G4cout << " itry1 " << itry1 << " sinth " << sinth << G4endl;
  
  if (itry1 == itry_max) {
    if (GetVerboseLevel() > 2)
      G4cout << " high energy angles generation: itry1 " << itry1 << G4endl;
    
    sinth = 0.5 * inuclRndm();
  }
 
  // Convert generated sin(theta) to cos(theta) with random sign
  G4double costh = std::sqrt(1.0 - sinth * sinth);
  if (inuclRndm() > 0.5) costh = -costh;
 
  return costh;
}

Referenced by FillDirManyBody(), and FillDirThreeBody().

◆ GenerateMultiBody()

void G4CascadeFinalStateAlgorithm::GenerateMultiBody	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protectedvirtual

Implements G4VHadDecayAlgorithm.

Definition at line 204 of file G4CascadeFinalStateAlgorithm.cc.

                                                  {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::GenerateMultiBody" << G4endl;
 
  if (G4CascadeParameters::usePhaseSpace()) {
    FillUsingKopylov(initialMass, masses, finalState);
    return;
  }
 
  finalState.clear();       // Initialization and sanity checks
  if (multiplicity < 3) return;
  if (!momDist) return;
 
  G4int itry = -1;      /* Loop checking 08.06.2015 MHK */
  while ((G4int)finalState.size() != multiplicity && ++itry < itry_max) {
    FillMagnitudes(initialMass, masses);
    FillDirections(initialMass, masses, finalState);
  }
}

◆ GenerateTwoBody()

void G4CascadeFinalStateAlgorithm::GenerateTwoBody	(	G4double	initialMass,
		const std::vector< G4double > &	masses,
		std::vector< G4LorentzVector > &	finalState )

protectedvirtual

Implements G4VHadDecayAlgorithm.

Definition at line 163 of file G4CascadeFinalStateAlgorithm.cc.

                                                {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::GenerateTwoBody" << G4endl;
 
  finalState.clear();       // Initialization and sanity checks
 
  if (multiplicity != 2) return;
 
  // Generate momentum vector in CMS for back-to-back particles
  G4double pscm = TwoBodyMomentum(initialMass, masses[0], masses[1]);
 
  G4double costh = angDist ? angDist->GetCosTheta(bullet_ekin, pscm)
                           : (2.*G4UniformRand() - 1.);
 
  mom.setRThetaPhi(pscm, std::acos(costh), UniformPhi());
 
  if (GetVerboseLevel()>3) {        // Copied from old G4EPCollider
    G4cout << " Particle kinds = " << kinds[0] << " , " << kinds[1]
       << "\n pmod " << pscm
       << "\n before rotation px " << mom.x() << " py " << mom.y()
       << " pz " << mom.z() << G4endl;
  }
 
  finalState.resize(2);             // Allows filling by index
 
  finalState[0].setVectM(mom, masses[0]);
  finalState[0] = toSCM.rotate(finalState[0]);
 
  if (GetVerboseLevel()>3) {        // Copied from old G4EPCollider
    G4cout << " after rotation px " << finalState[0].x() << " py "
       << finalState[0].y() << " pz " << finalState[0].z() << G4endl;
  }
 
  finalState[1].setVectM(-finalState[0].vect(), masses[1]);
}

◆ satisfyTriangle()

G4bool G4CascadeFinalStateAlgorithm::satisfyTriangle ( const std::vector< G4double > & pmod ) const

protected

Definition at line 295 of file G4CascadeFinalStateAlgorithm.cc.

                                                       {
  if (GetVerboseLevel() > 3) 
    G4cout << " >>> " << GetName() << "::satisfyTriangle" << G4endl;
 
  return ( (pmod.size() != 3) ||
       !(pmod[0] < std::fabs(pmod[1] - pmod[2]) ||
         pmod[0] > pmod[1] + pmod[2] ||
         pmod[1] < std::fabs(pmod[0] - pmod[2]) ||
         pmod[1] > pmod[0] + pmod[2] ||
         pmod[2] < std::fabs(pmod[0] - pmod[1]) ||
         pmod[2] > pmod[1] + pmod[0])
       );
}

Referenced by FillMagnitudes().

◆ SaveKinematics()

void G4CascadeFinalStateAlgorithm::SaveKinematics	(	G4InuclElementaryParticle *	bullet,
		G4InuclElementaryParticle *	target )

protected

Definition at line 114 of file G4CascadeFinalStateAlgorithm.cc.

                                              {
  if (GetVerboseLevel()>1) 
    G4cout << " >>> " << GetName() << "::SaveKinematics" << G4endl;
 
  if (target->nucleon()) {  // Which particle originated in nucleus?
    toSCM.setBullet(bullet);
    toSCM.setTarget(target);
  } else {
    toSCM.setBullet(target);
    toSCM.setTarget(bullet);
  }
 
  toSCM.toTheCenterOfMass();
 
  bullet_ekin = toSCM.getKinEnergyInTheTRS();
}

Referenced by Configure().

◆ SetVerboseLevel()

void G4CascadeFinalStateAlgorithm::SetVerboseLevel ( G4int verbose )

virtual

Reimplemented from G4VHadDecayAlgorithm.

Definition at line 81 of file G4CascadeFinalStateAlgorithm.cc.

                                                                {
  G4VHadDecayAlgorithm::SetVerboseLevel(verbose);
  G4MultiBodyMomentumDist::setVerboseLevel(verbose);
  G4TwoBodyAngularDist::setVerboseLevel(verbose);
  toSCM.setVerbose(verbose);
}

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

Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ G4CascadeFinalStateAlgorithm()

◆ ~G4CascadeFinalStateAlgorithm()

Member Function Documentation

◆ BetaKopylov()

◆ ChooseGenerators()

◆ Configure()

◆ FillDirections()

◆ FillDirManyBody()

◆ FillDirThreeBody()

◆ FillMagnitudes()

◆ FillUsingKopylov()

◆ GenerateCosTheta()

◆ GenerateMultiBody()

◆ GenerateTwoBody()

◆ satisfyTriangle()

◆ SaveKinematics()

◆ SetVerboseLevel()