G4HadPhaseSpaceNBodyAsai Class Reference

#include <G4HadPhaseSpaceNBodyAsai.hh>

Inheritance diagram for G4HadPhaseSpaceNBodyAsai:

Public Member Functions
	G4HadPhaseSpaceNBodyAsai (G4int verbose=0)
virtual	~G4HadPhaseSpaceNBodyAsai ()
Public Member Functions inherited from G4VHadPhaseSpaceAlgorithm
	G4VHadPhaseSpaceAlgorithm (const G4String &algName, G4int verbose=0)
virtual	~G4VHadPhaseSpaceAlgorithm ()
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)
virtual void	SetVerboseLevel (G4int verbose)
G4int	GetVerboseLevel () const
const G4String &	GetName () const

Protected Member Functions
virtual void	GenerateMultiBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
Protected Member Functions inherited from G4VHadPhaseSpaceAlgorithm
virtual void	GenerateTwoBody (G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4ThreeVector	UniformVector (G4double mag=1.) 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 38 of file G4HadPhaseSpaceNBodyAsai.hh.

Constructor & Destructor Documentation

G4HadPhaseSpaceNBodyAsai::G4HadPhaseSpaceNBodyAsai ( G4int  verbose = 0 )

inline

Definition at line 40 of file G4HadPhaseSpaceNBodyAsai.hh.

    : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceNBodyAsai",verbose) {;}

virtual G4HadPhaseSpaceNBodyAsai::~G4HadPhaseSpaceNBodyAsai ( )

inlinevirtual

Definition at line 42 of file G4HadPhaseSpaceNBodyAsai.hh.

{;}

Member Function Documentation

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

protectedvirtual

Implements G4VHadDecayAlgorithm.

Definition at line 52 of file G4HadPhaseSpaceNBodyAsai.cc.

References G4cout, G4endl, G4UniformRand, G4VHadDecayAlgorithm::GetName(), G4VHadDecayAlgorithm::GetVerboseLevel(), python.hepunit::GeV, G4VHadDecayAlgorithm::PrintVector(), G4InuclParticleNames::sm, sort(), G4VHadDecayAlgorithm::TwoBodyMomentum(), and G4VHadPhaseSpaceAlgorithm::UniformVector().

                                                  {
  if (GetVerboseLevel()) G4cout << GetName() << "::GenerateMultiBody" << G4endl;

  finalState.clear();

  //daughters' mass
  G4int numberOfDaughters = masses.size();
  G4double sumofmasses =
    std::accumulate(masses.begin(), masses.end(), 0.);
  
  //Calculate daughter momentum
  std::vector<G4double> daughtermomentum(numberOfDaughters);
  std::vector<G4double> sm(numberOfDaughters);
  G4double tmas;
  G4double weight = 1.0;
  G4int numberOfTry = 0; 
  G4int i;

  std::vector<G4double> rd(numberOfDaughters);
  do {
    //Generate random number in descending order
    rd[0] = 1.0;
    std::generate(rd.begin()+1, rd.end(), uniformRand);
    std::sort(rd.begin(), rd.end(), std::greater<G4double>());

    if (GetVerboseLevel()>1) PrintVector(rd,"rd",G4cout);

    //calcurate virtual mass 
    tmas = initialMass -  sumofmasses;
    G4double temp = sumofmasses; 
    for(i =0; i < numberOfDaughters; i++) {
      sm[i] = rd[i]*tmas + temp;
      temp -= masses[i];
      if (GetVerboseLevel()>1) {
        G4cout << i << " random number:" << rd[i]
           << " virtual mass:" << sm[i]/GeV << " GeV/c2" <<G4endl; 
      }
    }

    //Calculate daughter momentum
    weight = 1.0;
    i = numberOfDaughters-1;
    daughtermomentum[i] = TwoBodyMomentum(sm[i-1],masses[i-1],sm[i]);
    if (GetVerboseLevel()>1) {
      G4cout << " daughter " << i << ": momentum "
         << daughtermomentum[i]/GeV << " GeV/c" <<G4endl;
    }
    for(i =numberOfDaughters-2; i>=0; i--) {
      // calculate 
      daughtermomentum[i] = TwoBodyMomentum(sm[i],masses[i],sm[i+1]);
      if(daughtermomentum[i] < 0.0) {
        // !!! illegal momentum !!!
        if (GetVerboseLevel()>0) {
          G4cout << "G4HadPhaseSpaceNBodyAsai::Generate "
         << " can not calculate daughter momentum "
         << "\n initialMass " << initialMass/GeV << " GeV/c2"
         << "\n daughter " << i << ": mass "
         << masses[i]/GeV << " GeV/c2; momentum "
         << daughtermomentum[i]/GeV << " GeV/c" << G4endl;
        }
    return;         // Error detection
      }

      // calculate weight of this events
      weight *= daughtermomentum[i]/sm[i];
      if (GetVerboseLevel()>1) {
    G4cout << " daughter " << i << ": momentum "
           << daughtermomentum[i]/GeV << " GeV/c" <<G4endl;
      }
    }
    if (GetVerboseLevel()>1) {
      G4cout << " weight: " << weight <<G4endl;
    }
    
    // exit if number of Try exceeds 100
    if (numberOfTry++ > 100) {
      if (GetVerboseLevel()>0) {
        G4cout << "G4HadPhaseSpaceNBodyAsai::Generate "
           << " can not determine Decay Kinematics " << G4endl;
      }
      return;           // Error detection
    }
  } while (weight > G4UniformRand());

  if (GetVerboseLevel()>1) {
      G4cout << "Start calulation of daughters momentum vector "<<G4endl;
  }
  
  G4double beta;

  finalState.resize(numberOfDaughters);

  i = numberOfDaughters-2;

  G4ThreeVector direction = UniformVector(daughtermomentum[i]);

  finalState[i].setVectM(direction, masses[i]);
  finalState[i+1].setVectM(-direction, masses[i+1]);

  for (i = numberOfDaughters-3;  i >= 0; i--) {
    direction = UniformVector();

    //create daughter particle
    finalState[i].setVectM(-daughtermomentum[i]*direction, masses[i]);

    // boost already created particles 
    beta = daughtermomentum[i];
    beta /= std::sqrt(beta*beta + sm[i+1]*sm[i+1]);
    for (G4int j = i+1; j<numberOfDaughters; j++) {
      finalState[j].boost(beta*direction);
    }
  }
}

---

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

• G4HadPhaseSpaceNBodyAsai.hh
• G4HadPhaseSpaceNBodyAsai.cc