53 const std::vector<G4double>& masses,
54 std::vector<G4LorentzVector>& finalState) {
60 G4int numberOfDaughters = masses.size();
62 std::accumulate(masses.begin(), masses.end(), 0.);
65 std::vector<G4double> daughtermomentum(numberOfDaughters);
66 std::vector<G4double>
sm(numberOfDaughters);
69 G4int numberOfTry = 0;
72 std::vector<G4double> rd(numberOfDaughters);
76 std::generate(rd.begin()+1, rd.end(), uniformRand);
77 std::sort(rd.begin(), rd.end(), std::greater<G4double>());
82 tmas = initialMass - sumofmasses;
84 for(i =0; i < numberOfDaughters; i++) {
85 sm[i] = rd[i]*tmas + temp;
88 G4cout << i <<
" random number:" << rd[i]
89 <<
" virtual mass:" << sm[i]/
GeV <<
" GeV/c2" <<
G4endl;
95 i = numberOfDaughters-1;
98 G4cout <<
" daughter " << i <<
": momentum "
99 << daughtermomentum[i]/
GeV <<
" GeV/c" <<
G4endl;
101 for(i =numberOfDaughters-2; i>=0; i--) {
104 if(daughtermomentum[i] < 0.0) {
107 G4cout <<
"G4HadPhaseSpaceNBodyAsai::Generate "
108 <<
" can not calculate daughter momentum "
109 <<
"\n initialMass " << initialMass/
GeV <<
" GeV/c2"
110 <<
"\n daughter " << i <<
": mass "
111 << masses[i]/
GeV <<
" GeV/c2; momentum "
112 << daughtermomentum[i]/
GeV <<
" GeV/c" <<
G4endl;
118 weight *= daughtermomentum[i]/sm[i];
120 G4cout <<
" daughter " << i <<
": momentum "
121 << daughtermomentum[i]/
GeV <<
" GeV/c" <<
G4endl;
129 if (numberOfTry++ > 100) {
131 G4cout <<
"G4HadPhaseSpaceNBodyAsai::Generate "
132 <<
" can not determine Decay Kinematics " <<
G4endl;
139 G4cout <<
"Start calulation of daughters momentum vector "<<
G4endl;
144 finalState.resize(numberOfDaughters);
146 i = numberOfDaughters-2;
150 finalState[i].setVectM(direction, masses[i]);
151 finalState[i+1].setVectM(-direction, masses[i+1]);
153 for (i = numberOfDaughters-3; i >= 0; i--) {
157 finalState[i].setVectM(-daughtermomentum[i]*direction, masses[i]);
160 beta = daughtermomentum[i];
161 beta /= std::sqrt(beta*beta + sm[i+1]*sm[i+1]);
162 for (
G4int j = i+1; j<numberOfDaughters; j++) {
163 finalState[j].boost(beta*direction);
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
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4ThreeVector UniformVector(G4double mag=1.) const