G4HadPhaseSpaceGenbod Class Reference

#include <G4HadPhaseSpaceGenbod.hh>

Inheritance diagram for G4HadPhaseSpaceGenbod:

Collaboration diagram for G4HadPhaseSpaceGenbod:

Public Member Functions
	G4HadPhaseSpaceGenbod (G4int verbose=0)
virtual	~G4HadPhaseSpaceGenbod ()
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)
void	Initialize (G4double initialMass, const std::vector< G4double > &masses)
void	FillRandomBuffer ()
void	ComputeWeightScale (const std::vector< G4double > &masses)
void	FillEnergySteps (G4double initialMass, const std::vector< G4double > &masses)
void	GenerateMomenta (const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
void	AccumulateFinalState (size_t i, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4bool	AcceptEvent () const
G4double	ComputeWeight () const
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

Private Attributes
size_t	nFinal
G4double	totalMass
G4double	massExcess
G4double	weightMax
G4int	nTrials
std::vector< G4double >	msum
std::vector< G4double >	msq
std::vector< G4double >	rndm
std::vector< G4double >	meff
std::vector< G4double >	pd

Detailed Description

Definition at line 38 of file G4HadPhaseSpaceGenbod.hh.

Constructor & Destructor Documentation

G4HadPhaseSpaceGenbod()

G4HadPhaseSpaceGenbod::G4HadPhaseSpaceGenbod

(

G4int

verbose = 0

)

Definition at line 52 of file G4HadPhaseSpaceGenbod.cc.

  : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceGenbod",verbose),
    nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}

Here is the call graph for this function:

~G4HadPhaseSpaceGenbod()

virtual G4HadPhaseSpaceGenbod::~G4HadPhaseSpaceGenbod ( )

inlinevirtual

Definition at line 41 of file G4HadPhaseSpaceGenbod.hh.

{;}

Here is the call graph for this function:

Member Function Documentation

AcceptEvent()

G4bool G4HadPhaseSpaceGenbod::AcceptEvent ( ) const

protected

Definition at line 171 of file G4HadPhaseSpaceGenbod.cc.

                                                {
  if (GetVerboseLevel()>1) 
    G4cout << GetName() << "::AcceptEvent? " << nTrials << G4endl;

  return (G4UniformRand() <= ComputeWeight());
}

Here is the call graph for this function:

Here is the caller graph for this function:

AccumulateFinalState()

void G4HadPhaseSpaceGenbod::AccumulateFinalState ( size_t  i, const std::vector< G4double > &  masses, std::vector< G4LorentzVector > &  finalState  )

protected

Definition at line 198 of file G4HadPhaseSpaceGenbod.cc.

                                                     {
  if (GetVerboseLevel()>2)
    G4cout << GetName() << "::AccumulateFinalState " << i << G4endl;

  if (i==0) {           // First final state particle left alone
    finalState[i].setVectM(G4ThreeVector(0.,pd[i],0.),masses[i]);
    return;
  }
  
  finalState[i].setVectM(G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
  G4double phi = G4UniformRand() * twopi;
  G4double theta = std::acos(2.*G4UniformRand() - 1.);

  if (GetVerboseLevel() > 2) {
    G4cout << " initialized Py " << -pd[i-1] << " phi " << phi
       << " theta " << theta << G4endl;
  }

  G4double esys=0.,beta=0.,gamma=1.;
  if (i < nFinal-1) {           // Do not boost final particle
    esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
    beta = pd[i] / esys;
    gamma = esys / meff[i];

    if (GetVerboseLevel()>2)
      G4cout << " esys " << esys << " beta " << beta << " gamma " << gamma
         << G4endl;
  }

  for (size_t j=0; j<=i; j++) {     // Accumulate rotations
    finalState[j].rotateZ(theta).rotateY(phi);
    finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
    if (GetVerboseLevel()>2) G4cout << " j " << j << " " << finalState[j] << G4endl;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeWeight()

G4double G4HadPhaseSpaceGenbod::ComputeWeight ( ) const

protected

Definition at line 164 of file G4HadPhaseSpaceGenbod.cc.

                                                    {
  if (GetVerboseLevel()>1) G4cout << GetName() << "::ComputeWeight" << G4endl;

  return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
              std::multiplies<G4double>()));
}

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeWeightScale()

void G4HadPhaseSpaceGenbod::ComputeWeightScale ( const std::vector< G4double > &  masses )

protected

Definition at line 149 of file G4HadPhaseSpaceGenbod.cc.

                                                                           {
  if (GetVerboseLevel()>1) 
    G4cout << GetName() << "::ComputeWeightScale" << G4endl;

  weightMax = 1.;
  for (size_t i=1; i<nFinal; i++) {
    weightMax *= TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
  }

  if (GetVerboseLevel()>2) G4cout << " weightMax = " << weightMax << G4endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

FillEnergySteps()

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

protected

Definition at line 124 of file G4HadPhaseSpaceGenbod.cc.

                                                          {
  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillEnergySteps" << G4endl;

  meff.clear();
  pd.clear();

  meff.push_back(masses[0]);
  for (size_t i=1; i<nFinal-1; i++) {
    meff.push_back(rndm[i-1]*massExcess + msum[i]);
    pd.push_back(TwoBodyMomentum(meff[i], meff[i-1], masses[i]));
  }
  meff.push_back(initialMass);
  pd.push_back(TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[nFinal-1]));

  if (GetVerboseLevel()>2) {
    PrintVector(meff,"meff",G4cout);
    PrintVector(pd,"pd",G4cout);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

FillRandomBuffer()

void G4HadPhaseSpaceGenbod::FillRandomBuffer ( )

protected

Definition at line 111 of file G4HadPhaseSpaceGenbod.cc.

                                             {
  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillRandomBuffer" << G4endl;

  rndm.resize(nFinal-2,0.); // Final states generated in sorted order
  std::generate(rndm.begin(), rndm.end(), uniformRand);
  std::sort(rndm.begin(), rndm.end());
  if (GetVerboseLevel()>2) PrintVector(rndm, "rndm", G4cout);
}

Here is the call graph for this function:

Here is the caller graph for this function:

GenerateMomenta()

void G4HadPhaseSpaceGenbod::GenerateMomenta ( const std::vector< G4double > &  masses, std::vector< G4LorentzVector > &  finalState  )

protected

Definition at line 182 of file G4HadPhaseSpaceGenbod.cc.

                                                {
  if (GetVerboseLevel()>1) G4cout << GetName() << "::GenerateMomenta" << G4endl;

  finalState.resize(nFinal);    // Preallocate vectors for convenience below

  for (size_t i=0; i<nFinal; i++) {
    AccumulateFinalState(i, masses, finalState);
    if (GetVerboseLevel()>2)
      G4cout << " finalState[" << i << "] " << finalState[i] << G4endl;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

GenerateMultiBody()

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

protectedvirtual

Implements G4VHadDecayAlgorithm.

Definition at line 60 of file G4HadPhaseSpaceGenbod.cc.

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

  finalState.clear();

  Initialize(initialMass, masses);

  const G4int maxNumberOfLoops = 10000;
  nTrials = 0;
  do {              // Apply accept/reject to get distribution
    ++nTrials;
    FillRandomBuffer();
    FillEnergySteps(initialMass, masses);
  } while ( (!AcceptEvent()) && nTrials < maxNumberOfLoops );  /* Loop checking, 02.11.2015, A.Ribon */
  if ( nTrials >= maxNumberOfLoops ) {
    G4ExceptionDescription ed;
    ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
    G4Exception( " G4HadPhaseSpaceGenbod::GenerateMultiBody ", "HAD_GENBOD_001", FatalException, ed );
  }
  GenerateMomenta(masses, finalState);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Initialize()

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

protected

Definition at line 85 of file G4HadPhaseSpaceGenbod.cc.

                                                                    {
  if (GetVerboseLevel()>1) G4cout << GetName() << "::Initialize" << G4endl;

  nFinal = masses.size();
  msum.resize(nFinal, 0.);      // Initialize buffers for filling
  msq.resize(nFinal, 0.);

  std::partial_sum(masses.begin(), masses.end(), msum.begin());
  std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
         std::multiplies<G4double>());
  totalMass  = msum.back();
  massExcess = initialMass - totalMass;

  if (GetVerboseLevel()>2) {
    PrintVector(msum, "msum", G4cout);
    PrintVector(msq, "msq", G4cout);
    G4cout << " totalMass " << totalMass << " massExcess " << massExcess
       << G4endl;
  }

  ComputeWeightScale(masses);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

massExcess

G4double G4HadPhaseSpaceGenbod::massExcess

private

Definition at line 72 of file G4HadPhaseSpaceGenbod.hh.

meff

std::vector<G4double> G4HadPhaseSpaceGenbod::meff

private

Definition at line 79 of file G4HadPhaseSpaceGenbod.hh.

msq

std::vector<G4double> G4HadPhaseSpaceGenbod::msq

private

Definition at line 77 of file G4HadPhaseSpaceGenbod.hh.

msum

std::vector<G4double> G4HadPhaseSpaceGenbod::msum

private

Definition at line 76 of file G4HadPhaseSpaceGenbod.hh.

nFinal

size_t G4HadPhaseSpaceGenbod::nFinal

private

Definition at line 70 of file G4HadPhaseSpaceGenbod.hh.

nTrials

G4int G4HadPhaseSpaceGenbod::nTrials

private

Definition at line 74 of file G4HadPhaseSpaceGenbod.hh.

pd

std::vector<G4double> G4HadPhaseSpaceGenbod::pd

private

Definition at line 80 of file G4HadPhaseSpaceGenbod.hh.

rndm

std::vector<G4double> G4HadPhaseSpaceGenbod::rndm

private

Definition at line 78 of file G4HadPhaseSpaceGenbod.hh.

totalMass

G4double G4HadPhaseSpaceGenbod::totalMass

private

Definition at line 71 of file G4HadPhaseSpaceGenbod.hh.

weightMax

G4double G4HadPhaseSpaceGenbod::weightMax

private

Definition at line 73 of file G4HadPhaseSpaceGenbod.hh.

---

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

• G4HadPhaseSpaceGenbod.hh
• G4HadPhaseSpaceGenbod.cc