Consider initialising float fields to zero, not NAN. e.g. in smeared output, particles outside acceptance retain default values, therefore in e.g. jacquet-blondel calculations, those not detected screw up the calculation because they have E, p = NAN, not zero!
Revisit implementation, giving option for using particle energy or momentum when computing "energy", and think how to handle mass for smeared particles.
This could do with some reworking, to make it compatible with an arbitrary function via a functor (which I think is the best way to go). Move to its own files.
Consider the implementation here. If the particle is identified (either correctly or incorrectly), it could return the PDG mass of that particle. If the particle is not identified, should it return sqrt(E^2 - p^2)? What about if p or E are not known? Return E or p? Or zero (e.g. that would work for photons, E known, p not).