In CERN Root analysis software, there is a class called TGenPhaseSpace. This class can generate all possible solution of a nuclear reaction ( elastic, inelastic, 2-body, 3-body, decay, etc. ) base on isotropic distribution and balancing four-momenta.
The unit of mass, energy, and momentum are GeV/c^2, GeV, and GeV/c, respectively. The typical usage is
TLorentzVector target; TLorentzVector beam; Double_t masses[3] = { m1, m2, m3}; // what particles after reaction TGenPhaseSpace event; event.SetDecay(target + beam, 3, masses); for( Int_t n = 0; n < numEvent; n++){ Double_t xsec = event.Generate(); // cross-section TLorentzVector *p1 = event.GetDecay(0); // get the 4-momentum of m1 TLorentzVector *p2 = event.GetDecay(1); TLorentzVector *p3 = event.GetDecay(2); // fill histogram; }
In the proton-proton elastic reaction at 300 MeV, the calculation is very nice.
However, in the 23F(p,2p)22O inverse knockout reaction at 300A MeV , the Fermi-momentum of the bound proton, which is given by
where 
The Fermi-energy of a bound nucleon has maximum 40 MeV, or the momentum can be at most 400 MeV/c. But the TGenPhaseSpace only calculate all possible solution that balancing the 4-momenta, the Fermi-momentum can be more than 1 GeV/c.
This gives an un-realistic result. For example, the proton KE can be as large as 1200 MeV, while only 300A MeV of 23F KE. In realistic situation, the maximum KE of the scattered proton is at most 300 MeV at zero degree.
Now, if we restricted the Fermi-momentum to be 100 MeV/c +- 10 MeV/c, we have,
This is more realistic result.
If we can control the distribution of Fermi-momentum, and also understand the estimation of the cross-section, we can use it to simulate many reactions.
Nuclear Physics 101 