duttaANI
commented
2 years ago Invariant mass for N - body decay is calculated here. Also two-body decay here. The approach seems different from the one mentioned in stackoverflow
The problem and the solution is documented here: https://stackoverflow.com/questions/70758079/numerically-stable-calculation-of-invariant-mass-in-particle-physics/70775687#70775687
I discovered a numerically stable calculation of the invariant mass for a two-body decay. It should be trivial to generalize this to N-body decay.
The problem and the solution is documented here: https://stackoverflow.com/questions/70758079/numerically-stable-calculation-of-invariant-mass-in-particle-physics/70775687#70775687
Briefly, when one computes (1D and two-body decay for simplicity)
m_squared = (sqrt(m1 ** 2 + p1 ** 2) + sqrt(m2 ** 2 + p2 ** 2)) ** 2 - (p1 + p2) ** 2one gets catastrophic cancellation from the subtraction, when p_i >> m_i. I found an equivalent formula which is perfectly numerically stable and works accurately even in single precision.ROOT seems to use the equivalent of this naive formula: https://root.cern/root/html522/src/ROOT__Math__PxPyPzE4D_double_.h.html#D1mppD The fundamental issue of catastrophic cancellation is partially mitigated by computing this in double precision. This increases the range of values where the problem appears, but the problem is even there in double precision. This is evident from the fact that there is a check for negative values of M2() in the implementation of the method M(). The numerically stable formula is guarantees to never produce a negative value for M2().
I propose to implement the numerically stable calculation in ROOT.