Reported by Carla Gobel:

This is my understanding: (suppose we have the decay D-> R 2, with R -> 1 3)
There is indeed a problem when the mass of the resonance m_R is below the 13 threshold (m_1+m_3) since the triangle function gets negative (so is good to have that protection), but there is not such a problem when m_R is above the upper limit (mD – m_2); the “decay momenta” q(m_R) can still be calculated (lambda is positive).
It is true that above the upper limit you cannot calculate what it would be the “decay momentum” of the bachelor 2, calculated at m_R, in the D frame because lambda(m_D^2, m_R^2,m_2^2) is negative. But there is no *need* to calculate this quantity (you would need this, for instance, to calculate the normalized barrier factor in the vertex D -> R m_3, but you don’t need to normalize it since it just multiplies the amplitude). Explicitly:

A = F_D(J) * F_R(J) * BW (s13, m_0, gamma_R) * f_theta

The form factors F_D and F_R are not normalized as factors in the amplitude (they don’t need to, since any normalization there can be absorbed, so just don’t divide them by their value at q(m_R)). In the mass dependent width gamma_R the decay momentum calculated at m_R is needed (both in the normalized F_R(m13)/F_R(m_R) and in the factor [q/q(m_R)]^2J+1) but the quantity q(m_R) is is well defined for any value m_R > (m_1+m_3) (does not matter if m_R > (m_D – m_2))

From inspection of the code it does indeed seem unnecessary to explicitly normalise the two Blatt-Weisskopf factors that multiply the amplitude.
Should modify the code to remove it and check that behaviour is unchanged for both generation and fitting using the various examples.