54 nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}
61 const std::vector<G4double>& masses,
62 std::vector<G4LorentzVector>& finalState) {
69 const G4int maxNumberOfLoops = 10000;
76 if (
nTrials >= maxNumberOfLoops ) {
78 ed <<
" Failed sampling after maxNumberOfLoops attempts : forced exit" <<
G4endl;
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>());
116 std::sort(
rndm.begin(),
rndm.end());
125 const std::vector<G4double>& masses) {
131 meff.push_back(masses[0]);
132 for (
size_t i=1; i<
nFinal-1; i++) {
136 meff.push_back(initialMass);
154 for (
size_t i=1; i<
nFinal; i++) {
168 std::multiplies<G4double>()));
183 std::vector<G4LorentzVector>& finalState) {
186 finalState.resize(
nFinal);
188 for (
size_t i=0; i<
nFinal; i++) {
191 G4cout <<
" finalState[" << i <<
"] " << finalState[i] <<
G4endl;
199 const std::vector<G4double>& masses,
200 std::vector<G4LorentzVector>& finalState) {
214 G4cout <<
" initialized Py " << -
pd[i-1] <<
" phi " << phi
215 <<
" theta " << theta <<
G4endl;
222 gamma = esys /
meff[i];
225 G4cout <<
" esys " << esys <<
" beta " << beta <<
" gamma " << gamma
229 for (
size_t j=0; j<=i; j++) {
230 finalState[j].rotateZ(theta).rotateY(phi);
231 finalState[j].setY(gamma*(finalState[j].y() + beta*finalState[j].e()));
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
void FillEnergySteps(G4double initialMass, const std::vector< G4double > &masses)
void GenerateMomenta(const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
void ComputeWeightScale(const std::vector< G4double > &masses)
void AccumulateFinalState(size_t i, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
G4double ComputeWeight() const
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4HadPhaseSpaceGenbod(G4int verbose=0)
G4GLOB_DLL std::ostream G4cout
static const double twopi
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
std::vector< G4double > msq
const G4String & GetName() const
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
void Initialize(G4double initialMass, const std::vector< G4double > &masses)
std::vector< G4double > pd
std::vector< G4double > rndm
G4bool AcceptEvent() const
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
std::vector< G4double > msum
std::vector< G4double > meff