74 momDist(0), angDist(0), multiplicity(0), bullet_ekin(0.) {;}
91 const std::vector<G4int>& particle_kinds) {
106 kinds = particle_kinds;
136 <<
" is " << is <<
" fs " << fs <<
G4endl;
143 G4int kw = (fs==is) ? 1 : 2;
162 std::vector<G4LorentzVector>& finalState) {
180 <<
"\n pmod " << pscm
181 <<
"\n before rotation px " <<
mom.x() <<
" py " <<
mom.y()
185 finalState.resize(2);
187 finalState[0].setVectM(
mom, masses[0]);
191 G4cout <<
" after rotation px " << finalState[0].x() <<
" py "
192 << finalState[0].y() <<
" pz " << finalState[0].z() <<
G4endl;
195 finalState[1].setVectM(-finalState[0].vect(), masses[1]);
203 std::vector<G4LorentzVector>& finalState) {
238 G4cout <<
" knd_last " <<
kinds.back() <<
" mass_last "
254 if (pmod <
small)
break;
255 eleft -= std::sqrt(pmod*pmod + masses[i]*masses[i]);
259 <<
" mass2 " << masses[i]*masses[i] <<
" eleft " << eleft
260 <<
"\n x1 " << eleft - mass_last <<
G4endl;
263 if (eleft <= mass_last)
break;
268 if (i < multiplicity-1)
continue;
270 G4double plast = eleft * eleft - masses.back()*masses.back();
273 if (plast <=
small)
continue;
275 plast = std::sqrt(plast);
283 G4cerr <<
" Unable to generate momenta for multiplicity "
297 return ( (pmod.size() != 3) ||
298 !(pmod[0] < std::fabs(pmod[1] - pmod[2]) ||
299 pmod[0] > pmod[1] + pmod[2] ||
300 pmod[1] < std::fabs(pmod[0] - pmod[2]) ||
301 pmod[1] > pmod[0] + pmod[2] ||
302 pmod[2] < std::fabs(pmod[0] - pmod[1]) ||
303 pmod[2] > pmod[1] + pmod[0])
311 std::vector<G4LorentzVector>& finalState) {
327 std::vector<G4LorentzVector>& finalState) {
331 finalState.resize(3);
351 finalState[0] =
toSCM.
rotate(finalState[2], finalState[0]);
354 finalState[1].set(0.,0.,0.,initialMass);
355 finalState[1] -= finalState[0] + finalState[2];
360 std::vector<G4LorentzVector>& finalState) {
377 std::accumulate(finalState.begin(), finalState.end()-2,
G4LorentzVector());
380 costh = -0.5 * (pmod*pmod +
383 / pmod /
modules[multiplicity-2];
396 finalState[multiplicity-2] =
398 masses[multiplicity-2]);
399 finalState[multiplicity-2] =
toSCM.
rotate(psum, finalState[multiplicity-2]);
402 finalState[multiplicity-1].set(0.,0.,0.,initialMass);
403 finalState[multiplicity-1] -= psum + finalState[multiplicity-2];
412 G4cout <<
" >>> " <<
GetName() <<
"::GenerateCosTheta " << ptype
421 G4double p0 = ptype<3 ? 0.36 : 0.25;
422 G4double alf = 1.0 / p0 / (p0 - (pmod+p0)*std::exp(-pmod / p0));
430 G4double salf = s1 * alf * std::exp(-s1 / p0);
432 G4cout <<
" s1 * alf * std::exp(-s1 / p0) " << salf
433 <<
" s2 " << s2 <<
G4endl;
436 if (salf > s2) sinth = s1 / pmod;
440 G4cout <<
" itry1 " << itry1 <<
" sinth " << sinth <<
G4endl;
444 G4cout <<
" high energy angles generation: itry1 " << itry1 <<
G4endl;
450 G4double costh = std::sqrt(1.0 - sinth * sinth);
462 const std::vector<G4double>& masses,
463 std::vector<G4LorentzVector>& finalState) {
469 size_t N = masses.size();
470 finalState.resize(N);
472 G4double mtot = std::accumulate(masses.begin(), masses.end(), 0.0);
480 for (
size_t k=N-1; k>0; --k) {
486 boostV = recoil.boostVector();
493 finalState[k].setVectM(momV,masses[k]);
494 recoil.setVectM(-momV,recoilMass);
496 finalState[k].boost(boostV);
497 recoil.boost(boostV);
501 finalState[0] = recoil;
509 G4double Fmax = std::sqrt(g4pow->
powN(xN/(xN+1.),N)/(xN+1.));
514 F = std::sqrt(g4pow->
powN(chi,N)*(1.-chi));
virtual const G4String & GetName() const
static G4Pow * GetInstance()
static void setVerboseLevel(G4int vb=0)
static const G4double small
static const G4VTwoBodyAngDst * GetDist(G4int is, G4int fs, G4int kw)
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4double UniformTheta() const
G4double powN(G4double x, G4int n) const
CLHEP::Hep3Vector G4ThreeVector
void FillDirThreeBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4int GetVerboseLevel() const
G4LorentzVector rotate(const G4LorentzVector &mom) const
void setBullet(const G4InuclParticle *bullet)
std::vector< G4int > kinds
void SaveKinematics(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target)
const G4VTwoBodyAngDst * angDist
void FillMagnitudes(G4double initialMass, const std::vector< G4double > &masses)
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
static const G4int itry_max
G4bool satisfyTriangle(const std::vector< G4double > &pmod) const
std::vector< G4double > modules
G4double getKinEnergyInTheTRS() const
G4double GenerateCosTheta(G4int ptype, G4double pmod) const
virtual const G4String & GetName() const
virtual G4double GetMomentum(G4int ptype, const G4double &ekin) const =0
static void setVerboseLevel(G4int vb=0)
static G4bool usePhaseSpace()
virtual void SetVerboseLevel(G4int verbose)
void Configure(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target, const std::vector< G4int > &particle_kinds)
const G4VMultiBodyMomDst * momDist
const G4String & GetName() const
virtual void GenerateTwoBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
G4double UniformPhi() const
static const G4double oneOverE
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
void FillDirections(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
static const G4double maxCosTheta
void FillUsingKopylov(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4CascadeFinalStateAlgorithm()
virtual void SetVerboseLevel(G4int verbose)
void ChooseGenerators(G4int is, G4int fs)
void FillDirManyBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
virtual ~G4CascadeFinalStateAlgorithm()
virtual G4double GetCosTheta(const G4double &ekin, const G4double &pcm) const =0
static const G4VMultiBodyMomDst * GetDist(G4int is, G4int mult)
G4double BetaKopylov(G4int K) const
G4GLOB_DLL std::ostream G4cerr
void setTarget(const G4InuclParticle *target)
CLHEP::HepLorentzVector G4LorentzVector