75 momDist(0), angDist(0), multiplicity(0), bullet_ekin(0.) {;}
92 const std::vector<G4int>& particle_kinds) {
107 kinds = particle_kinds;
137 <<
" is " << is <<
" fs " << fs <<
G4endl;
144 G4int kw = (fs==is) ? 1 : 2;
163 std::vector<G4LorentzVector>& finalState) {
181 <<
"\n pmod " << pscm
182 <<
"\n before rotation px " <<
mom.x() <<
" py " <<
mom.y()
186 finalState.resize(2);
188 finalState[0].setVectM(
mom, masses[0]);
192 G4cout <<
" after rotation px " << finalState[0].x() <<
" py "
193 << finalState[0].y() <<
" pz " << finalState[0].z() <<
G4endl;
196 finalState[1].setVectM(-finalState[0].vect(), masses[1]);
204 std::vector<G4LorentzVector>& finalState) {
239 G4cout <<
" knd_last " <<
kinds.back() <<
" mass_last "
255 if (pmod <
small)
break;
256 eleft -= std::sqrt(pmod*pmod + masses[i]*masses[i]);
260 <<
" mass2 " << masses[i]*masses[i] <<
" eleft " << eleft
261 <<
"\n x1 " << eleft - mass_last <<
G4endl;
264 if (eleft <= mass_last)
break;
269 if (i < multiplicity-1)
continue;
271 G4double plast = eleft * eleft - masses.back()*masses.back();
274 if (plast <=
small)
continue;
276 plast = std::sqrt(plast);
284 G4cerr <<
" Unable to generate momenta for multiplicity "
298 return ( (pmod.size() != 3) ||
299 !(pmod[0] < std::fabs(pmod[1] - pmod[2]) ||
300 pmod[0] > pmod[1] + pmod[2] ||
301 pmod[1] < std::fabs(pmod[0] - pmod[2]) ||
302 pmod[1] > pmod[0] + pmod[2] ||
303 pmod[2] < std::fabs(pmod[0] - pmod[1]) ||
304 pmod[2] > pmod[1] + pmod[0])
312 std::vector<G4LorentzVector>& finalState) {
328 std::vector<G4LorentzVector>& finalState) {
332 finalState.resize(3);
352 finalState[0] =
toSCM.
rotate(finalState[2], finalState[0]);
355 finalState[1].set(0.,0.,0.,initialMass);
356 finalState[1] -= finalState[0] + finalState[2];
361 std::vector<G4LorentzVector>& finalState) {
378 std::accumulate(finalState.begin(), finalState.end()-2,
G4LorentzVector());
381 costh = -0.5 * (pmod*pmod +
384 / pmod /
modules[multiplicity-2];
397 finalState[multiplicity-2] =
399 masses[multiplicity-2]);
400 finalState[multiplicity-2] =
toSCM.
rotate(psum, finalState[multiplicity-2]);
403 finalState[multiplicity-1].set(0.,0.,0.,initialMass);
404 finalState[multiplicity-1] -= psum + finalState[multiplicity-2];
413 G4cout <<
" >>> " <<
GetName() <<
"::GenerateCosTheta " << ptype
422 G4double p0 = ptype<3 ? 0.36 : 0.25;
423 G4double alf = 1.0 / p0 / (p0 - (pmod+p0)*std::exp(-pmod / p0));
431 G4double salf = s1 * alf * std::exp(-s1 / p0);
433 G4cout <<
" s1 * alf * std::exp(-s1 / p0) " << salf
434 <<
" s2 " << s2 <<
G4endl;
437 if (salf > s2) sinth = s1 / pmod;
441 G4cout <<
" itry1 " << itry1 <<
" sinth " << sinth <<
G4endl;
445 G4cout <<
" high energy angles generation: itry1 " << itry1 <<
G4endl;
451 G4double costh = std::sqrt(1.0 - sinth * sinth);
463 const std::vector<G4double>& masses,
464 std::vector<G4LorentzVector>& finalState) {
470 size_t N = masses.size();
471 finalState.resize(N);
473 G4double mtot = std::accumulate(masses.begin(), masses.end(), 0.0);
481 for (
size_t k=N-1; k>0; --k) {
487 boostV = recoil.boostVector();
494 finalState[k].setVectM(momV,masses[k]);
495 recoil.setVectM(-momV,recoilMass);
497 finalState[k].boost(boostV);
498 recoil.boost(boostV);
502 finalState[0] = recoil;
510 G4double Fmax = std::sqrt(g4pow->
powN(xN/(xN+1.),N)/(xN+1.));
515 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