The 4-body phase-space can be rather large and chance to obtain events in the corners of m12-m34 space is rather small. There are specific use cases where one wants to generate only events in the restricted range in m12-m34 plane and as one might want to select only small part or because small probability in the corners, just selecting events from full phase-space decays can be very inefficient. This model solves that by generating only events in wanted part of m12-m34 plane and sets maximum probability appropriately for the given window. Full phase-space can be generated with this model, but might not be best option for that.
In the process of testing, results of this decay model is compared to those from generic PHSP model including cases of restricted window.
Algorithm looks OK overall, although it would be good if we could try to avoid using double pointers in various functions.
This function can be constant.
I don't think we need to use double pointers here. Just define two phspFactor() functions, where the 2nd one includes q, p1 and p3 as extra arguments?
I'm not sure I understand. The reason why I use pointers here is that I'm actually returning also individual momenta so I do not have to calculate them twice. An alternative would be to avoid this function completely as it is used just in two places (when searching for maximum probability and at one point in decay).
This is not class member, just standalone function, so const does not apply here.
After a bit of thinking, I decided against removing this function completely as in the search for maximum probability it is used in a few places and that would get a bit messier by repeating several lines every time.
I'm afraid I don't have time right now to really fully review the physics. But hopefully your tests should have covered that. If we think it's necessary, perhaps we can ask the analysts who were requesting this to do some further validations?
I've added a follow-up on the issue of the various double* arguments - in short I'm also not a fan but let's leave them as they are for now.
I've also added one other small question on initialisation of member variables.
I don't see any other problems.
I see you've given the last 3 of these initial values. Can we give the others sensible initial values as well?
I'm also not a fan of this sort of construct. There are some alternatives, such as returning an array of 4 doubles or a tuple of 4 doubles.