54 nFinal(0), totalMass(0.), massExcess(0.), weightMax(0.), nTrials(0) {;}
61 const std::vector<G4double>& masses,
62 std::vector<G4LorentzVector>& finalState) {
83 nFinal = masses.size();
84 msum.resize(nFinal, 0.);
85 msq.resize(nFinal, 0.);
87 std::partial_sum(masses.begin(), masses.end(), msum.begin());
88 std::transform(masses.begin(), masses.end(), masses.begin(), msq.begin(),
89 std::multiplies<G4double>());
90 totalMass = msum.back();
91 massExcess = initialMass - totalMass;
96 G4cout <<
" totalMass " << totalMass <<
" massExcess " << massExcess
109 rndm.resize(nFinal-2,0.);
110 std::generate(rndm.begin(), rndm.end(), uniformRand);
120 const std::vector<G4double>& masses) {
126 meff.push_back(masses[0]);
127 for (
size_t i=1; i<nFinal-1; i++) {
128 meff.push_back(rndm[i-1]*massExcess + msum[i]);
131 meff.push_back(initialMass);
132 pd.push_back(
TwoBodyMomentum(meff[nFinal-1], meff[nFinal-2], masses[nFinal-1]));
149 for (
size_t i=1; i<nFinal; i++) {
150 weightMax *=
TwoBodyMomentum(massExcess+msum[i], msum[i-1], masses[i]);
162 return (std::accumulate(pd.begin(), pd.end(), 1./weightMax,
163 std::multiplies<G4double>()));
178 std::vector<G4LorentzVector>& finalState) {
181 finalState.resize(nFinal);
183 for (
size_t i=0; i<nFinal; i++) {
186 G4cout <<
" finalState[" << i <<
"] " << finalState[i] <<
G4endl;
194 const std::vector<G4double>& masses,
195 std::vector<G4LorentzVector>& finalState) {
200 finalState[i].setVectM(
G4ThreeVector(0.,pd[i],0.),masses[i]);
204 finalState[i].setVectM(
G4ThreeVector(0.,-pd[i-1],0.),masses[i]);
209 G4cout <<
" initialized Py " << -pd[i-1] <<
" phi " << phi
210 <<
" theta " << theta <<
G4endl;
215 esys = std::sqrt(pd[i]*pd[i]+meff[i]*meff[i]);
217 gamma = esys / meff[i];
220 G4cout <<
" esys " << esys <<
" beta " << beta <<
" gamma " << gamma
224 for (
size_t j=0; j<=i; j++) {
225 finalState[j].rotateZ(theta).rotateY(phi);
226 finalState[j].setY(gamma*(finalState[j].
y() + beta*finalState[j].e()));
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
void PrintVector(const std::vector< G4double > &v, const G4String &name, std::ostream &os) const
const G4String & GetName() const
void Initialize(G4double initialMass, const std::vector< G4double > &masses)
G4bool AcceptEvent() const
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const