TGenPhaseSpace Issue with d+d->3He+n

[physRR]

Hello,

I am trying to write a simple Monte Carlo routine to get the neutron energy distribution from d-d fusion. I was hoping to use TGenPhaseSpace to do so, and wanted to make sure I understood it so I started with the simple case of two deuterons at rest, which should yield a 2.45 MeV neutron. I have the follwing wirtten in python:

  # Two deutrium at rest
  d1 = ROOT.TLorentzVector(0.0,0.0,0.0,1.8756)
  d2 = ROOT.TLorentzVector(0.0,0.0,0.0,1.8756)
  total = d1+d2
 
  # end product masses - neutron and 3He
  masses = np.array([0.93956,2.809])

  event = ROOT.TGenPhaseSpace()
  event.SetDecay(total,2,masses)
  
  for it in range(0,10):
    event.Generate();
    pNeut = event.GetDecay(0);
    pHe = event.GetDecay(1);

    #kinetic energy is total minus rest mass
    nKE = pNeut.Energy() - pNeut.Mag()
    print "Neutron KE is "+str(nKE)

However, I get 1.97 MeV for the neutron energy instead of the expected 2.45MeV. Am I misinterpreting how to use this class? Thanks so much for any help!

_ROOT Version: 6.16/00
_Platform: Ubuntu
Compiler: Not Provided

---

[physRR]

Ok nevermind, I think I figured it out. I just needed to add more digits to the masses. Since they are in GeV they are pretty sensitive. Using 2.8084 for the 3He mass gets me much closer to the expected 2.45MeV for the neutron energy.

[etejedor]

Thanks for reporting the solution!