Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 G4StringGetName () 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
 

Detailed Description

Definition at line 38 of file G4HadPhaseSpaceGenbod.hh.

Constructor & Destructor Documentation

G4HadPhaseSpaceGenbod::G4HadPhaseSpaceGenbod ( G4int  verbose = 0)

Definition at line 52 of file G4HadPhaseSpaceGenbod.cc.

53  : G4VHadPhaseSpaceAlgorithm("G4HadPhaseSpaceGenbod",verbose),
54  nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}
G4VHadPhaseSpaceAlgorithm(const G4String &algName, G4int verbose=0)
virtual G4HadPhaseSpaceGenbod::~G4HadPhaseSpaceGenbod ( )
inlinevirtual

Definition at line 41 of file G4HadPhaseSpaceGenbod.hh.

41 {;}

Member Function Documentation

G4bool G4HadPhaseSpaceGenbod::AcceptEvent ( ) const
protected

Definition at line 171 of file G4HadPhaseSpaceGenbod.cc.

171  {
172  if (GetVerboseLevel()>1)
173  G4cout << GetName() << "::AcceptEvent? " << nTrials << G4endl;
174 
175  return (G4UniformRand() <= ComputeWeight());
176 }
G4int GetVerboseLevel() const
G4double ComputeWeight() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 198 of file G4HadPhaseSpaceGenbod.cc.

200  {
201  if (GetVerboseLevel()>2)
202  G4cout << GetName() << "::AccumulateFinalState " << i << G4endl;
203 
204  if (i==0) { // First final state particle left alone
205  finalState[i].setVectM(G4ThreeVector(0.,pd[i],0.),masses[i]);
206  return;
207  }
208 
209  finalState[i].setVectM(G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
210  G4double phi = G4UniformRand() * twopi;
211  G4double theta = std::acos(2.*G4UniformRand() - 1.);
212 
213  if (GetVerboseLevel() > 2) {
214  G4cout << " initialized Py " << -pd[i-1] << " phi " << phi
215  << " theta " << theta << G4endl;
216  }
217 
218  G4double esys=0.,beta=0.,gamma=1.;
219  if (i < nFinal-1) { // Do not boost final particle
220  esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
221  beta = pd[i] / esys;
222  gamma = esys / meff[i];
223 
224  if (GetVerboseLevel()>2)
225  G4cout << " esys " << esys << " beta " << beta << " gamma " << gamma
226  << G4endl;
227  }
228 
229  for (size_t j=0; j<=i; j++) { // Accumulate rotations
230  finalState[j].rotateZ(theta).rotateY(phi);
231  finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
232  if (GetVerboseLevel()>2) G4cout << " j " << j << " " << finalState[j] << G4endl;
233  }
234 }
CLHEP::Hep3Vector G4ThreeVector
G4int GetVerboseLevel() const
static constexpr double twopi
Definition: G4SIunits.hh:76
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4HadPhaseSpaceGenbod::ComputeWeight ( ) const
protected

Definition at line 164 of file G4HadPhaseSpaceGenbod.cc.

164  {
165  if (GetVerboseLevel()>1) G4cout << GetName() << "::ComputeWeight" << G4endl;
166 
167  return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
168  std::multiplies<G4double>()));
169 }
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 149 of file G4HadPhaseSpaceGenbod.cc.

149  {
150  if (GetVerboseLevel()>1)
151  G4cout << GetName() << "::ComputeWeightScale" << G4endl;
152 
153  weightMax = 1.;
154  for (size_t i=1; i<nFinal; i++) {
155  weightMax *= TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
156  }
157 
158  if (GetVerboseLevel()>2) G4cout << " weightMax = " << weightMax << G4endl;
159 }
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 124 of file G4HadPhaseSpaceGenbod.cc.

125  {
126  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillEnergySteps" << G4endl;
127 
128  meff.clear();
129  pd.clear();
130 
131  meff.push_back(masses[0]);
132  for (size_t i=1; i<nFinal-1; i++) {
133  meff.push_back(rndm[i-1]*massExcess + msum[i]);
134  pd.push_back(TwoBodyMomentum(meff[i], meff[i-1], masses[i]));
135  }
136  meff.push_back(initialMass);
137  pd.push_back(TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[nFinal-1]));
138 
139  if (GetVerboseLevel()>2) {
140  PrintVector(meff,"meff",G4cout);
141  PrintVector(pd,"pd",G4cout);
142  }
143 }
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
const G4String & GetName() const
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

void G4HadPhaseSpaceGenbod::FillRandomBuffer ( )
protected

Definition at line 111 of file G4HadPhaseSpaceGenbod.cc.

111  {
112  if (GetVerboseLevel()>1) G4cout << GetName() << "::FillRandomBuffer" << G4endl;
113 
114  rndm.resize(nFinal-2,0.); // Final states generated in sorted order
115  std::generate(rndm.begin(), rndm.end(), uniformRand);
116  std::sort(rndm.begin(), rndm.end());
117  if (GetVerboseLevel()>2) PrintVector(rndm, "rndm", G4cout);
118 }
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
const G4String & GetName() const
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 182 of file G4HadPhaseSpaceGenbod.cc.

183  {
184  if (GetVerboseLevel()>1) G4cout << GetName() << "::GenerateMomenta" << G4endl;
185 
186  finalState.resize(nFinal); // Preallocate vectors for convenience below
187 
188  for (size_t i=0; i<nFinal; i++) {
189  AccumulateFinalState(i, masses, finalState);
190  if (GetVerboseLevel()>2)
191  G4cout << " finalState[" << i << "] " << finalState[i] << G4endl;
192  }
193 }
void AccumulateFinalState(size_t i, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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.

62  {
63  if (GetVerboseLevel()) G4cout << GetName() << "::GenerateMultiBody" << G4endl;
64 
65  finalState.clear();
66 
67  Initialize(initialMass, masses);
68 
69  const G4int maxNumberOfLoops = 10000;
70  nTrials = 0;
71  do { // Apply accept/reject to get distribution
72  ++nTrials;
74  FillEnergySteps(initialMass, masses);
75  } while ( (!AcceptEvent()) && nTrials < maxNumberOfLoops ); /* Loop checking, 02.11.2015, A.Ribon */
76  if ( nTrials >= maxNumberOfLoops ) {
78  ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
79  G4Exception( " G4HadPhaseSpaceGenbod::GenerateMultiBody ", "HAD_GENBOD_001", FatalException, ed );
80  }
81  GenerateMomenta(masses, finalState);
82 }
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void FillEnergySteps(G4double initialMass, const std::vector< G4double > &masses)
void GenerateMomenta(const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4String & GetName() const
void Initialize(G4double initialMass, const std::vector< G4double > &masses)
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

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

Definition at line 85 of file G4HadPhaseSpaceGenbod.cc.

85  {
86  if (GetVerboseLevel()>1) G4cout << GetName() << "::Initialize" << G4endl;
87 
88  nFinal = masses.size();
89  msum.resize(nFinal, 0.); // Initialize buffers for filling
90  msq.resize(nFinal, 0.);
91 
92  std::partial_sum(masses.begin(), masses.end(), msum.begin());
93  std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
94  std::multiplies<G4double>());
95  totalMass = msum.back();
96  massExcess = initialMass - totalMass;
97 
98  if (GetVerboseLevel()>2) {
99  PrintVector(msum, "msum", G4cout);
100  PrintVector(msq, "msq", G4cout);
101  G4cout << " totalMass " << totalMass << " massExcess " << massExcess
102  << G4endl;
103  }
104 
105  ComputeWeightScale(masses);
106 }
void ComputeWeightScale(const std::vector< G4double > &masses)
G4int GetVerboseLevel() const
G4GLOB_DLL std::ostream G4cout
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
const G4String & GetName() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:


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