59 using namespace G4InuclSpecialFunctions;
64 const G4double G4CascadeFinalStateAlgorithm::maxCosTheta = 0.9999;
65 const G4double G4CascadeFinalStateAlgorithm::oneOverE = 0.3678794;
66 const G4double G4CascadeFinalStateAlgorithm::small = 1.e-10;
67 const G4int G4CascadeFinalStateAlgorithm::itry_max = 10;
74 momDist(0), angDist(0), multiplicity(0), bullet_ekin(0.) {;}
91 const std::vector<G4int>& particle_kinds) {
96 multiplicity = particle_kinds.size();
98 G4int fs = (multiplicity==2) ? particle_kinds[0]*particle_kinds[1] : 0;
106 kinds = particle_kinds;
136 <<
" is " << is <<
" fs " << fs <<
G4endl;
142 if (fs > 0 && multiplicity == 2) {
143 G4int kw = (fs==is) ? 1 : 2;
145 }
else if (multiplicity == 3) {
153 << (angDist?angDist->
GetName():
"") << G4endl;
162 std::vector<G4LorentzVector>& finalState) {
168 if (multiplicity != 2)
return;
179 G4cout <<
" Particle kinds = " << kinds[0] <<
" , " << kinds[1]
180 <<
"\n pmod " << pscm
181 <<
"\n before rotation px " << mom.
x() <<
" py " << mom.
y()
182 <<
" pz " << mom.
z() <<
G4endl;
185 finalState.resize(2);
187 finalState[0].setVectM(mom, masses[0]);
188 finalState[0] = toSCM.
rotate(finalState[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) {
213 if (multiplicity < 3)
return;
214 if (!momDist)
return;
217 while ((
G4int)finalState.size() != multiplicity && ++itry < itry_max) {
230 if (!momDist)
return;
232 modules.resize(multiplicity,0.);
238 G4cout <<
" knd_last " << kinds.back() <<
" mass_last "
243 while (++itry < itry_max) {
251 for (i=0; i < multiplicity-1; i++) {
252 pmod = momDist->
GetMomentum(kinds[i], bullet_ekin);
254 if (pmod < small)
break;
255 eleft -= std::sqrt(pmod*pmod + masses[i]*masses[i]);
258 G4cout <<
" kp " << kinds[i] <<
" pmod " << pmod
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);
276 modules.back() = plast;
281 if (itry >= itry_max) {
283 G4cerr <<
" Unable to generate momenta for multiplicity "
284 << multiplicity <<
G4endl;
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) {
316 if ((
G4int)modules.size() != multiplicity)
return;
319 if (multiplicity == 3)
327 std::vector<G4LorentzVector>& finalState) {
331 finalState.resize(3);
335 finalState[2] = toSCM.
rotate(finalState[2]);
338 costh = -0.5 * (modules[2]*modules[2] + modules[0]*modules[0] -
339 modules[1]*modules[1]) / modules[2] / modules[0];
341 if (std::fabs(costh) >= maxCosTheta) {
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) {
367 finalState.resize(multiplicity);
369 for (
G4int i=0; i<multiplicity-2; i++) {
372 finalState[i] = toSCM.
rotate(finalState[i]);
377 std::accumulate(finalState.begin(), finalState.end()-2,
G4LorentzVector());
380 costh = -0.5 * (pmod*pmod +
381 modules[multiplicity-2]*modules[multiplicity-2] -
382 modules[multiplicity-1]*modules[multiplicity-1])
383 / pmod / modules[multiplicity-2];
387 if (std::fabs(costh) >= maxCosTheta) {
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
416 if (multiplicity == 3) {
421 G4double p0 = ptype<3 ? 0.36 : 0.25;
422 G4double alf = 1.0 / p0 / (p0 - (pmod+p0)*std::exp(-pmod / p0));
427 while (std::fabs(sinth) > maxCosTheta && ++itry1 < itry_max) {
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;
442 if (itry1 == itry_max) {
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) {
493 finalState[k].setVectM(momV,masses[k]);
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 G4VTwoBodyAngDst * GetDist(G4int is, G4int fs, G4int kw)
Hep3Vector boostVector() const
virtual void GenerateMultiBody(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
G4double UniformTheta() const
G4double powN(G4double x, G4int n) const
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)
void SaveKinematics(G4InuclElementaryParticle *bullet, G4InuclElementaryParticle *target)
void setVectM(const Hep3Vector &spatial, double mass)
void FillMagnitudes(G4double initialMass, const std::vector< G4double > &masses)
G4GLOB_DLL std::ostream G4cout
void setVerbose(G4int vb=0)
HepLorentzVector & boost(double, double, double)
G4bool satisfyTriangle(const std::vector< G4double > &pmod) const
G4double getKinEnergyInTheTRS() const
void setRThetaPhi(double r, double theta, double phi)
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 G4String & GetName() const
void set(double x, double y, double z, double t)
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
G4double TwoBodyMomentum(G4double M0, G4double M1, G4double M2) const
void FillDirections(G4double initialMass, const std::vector< G4double > &masses, std::vector< G4LorentzVector > &finalState)
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