When we do sequential decays, there is part of the phase-space factor missing compared to treating it as decay directly to final state with amplitude corresponding to resonance from sequential decay. Particularly well it is visible with wide resonance close to kinematic boundary. I understand issue itself and in principle understand what exactly is missing, but when looking to it I did not immediately see how to fix it in generic way (possibly by not having enough time to dig through couple of test cases).

Original communication from Tomasz Skwarnicki:

I personally don’t have direct experience with 3-body PHSP generation in EvtGen, but it is likely that omega mass is generated from non-relativistic Breit-Wigner according to point 6. on page 65 of https://evtgen.hepforge.org/doc/EvtGenGuide.pdf , and then 3-body B->J/psi omega K decay is generated.

This will be fine for most of phase-space. Unfortunately, the X(3872) is below J/psi omega threshold if you take the omega pole mass. That means X(3872) is near the J/psi omega K Dalitz plot boundary, where EvtGen algorithm will not do a good job. The correct way would have been to generate flat omega mass (in relevant range), generate B->J/psi omega K 3-body phase-space and then reweight generated distribution in omega mass, with the resonant omega shape as weight. If you are using this MC in the narrow region near X(3872) mass, then probably reweighting MC with J/psi momentum in J/psi omega rest frame should give you a good approximation.

If you are using this MC to do full amplitude analysis of all J/psi omega and K omega structures, then you could do toy simulations of B->J/psi omega K the way I described (phase-space first, omega mass weights second) and the way EvtGen does it (omega mass first, phase-space second).

You will probably need to reweight MC by the ratio of these two (correct/EvtGen-way) in 3D space of m(J/psi omega), cos(helicity_angle), m(omega) where helicity_angle is the angle between J/psi and K direction in J/psi omega rest frame. The differences in m(omega) shape should show up near lower and upper boundaries of m(J/psi omega) and cos(helicity_angle). You can use such toys to study, how much of the distortion of omega shape you get between the two approaches. You could also exclude the regions near Dalitz plot boundaries from the amplitude fit, but this will exclude X(3872) itself. What you need to do depends on how exactly you are analyzing data. If you like, we can take this discussion offline.