Dear ROOT experts;

I’m at my wits’ end after discovering a seemingly major bug affecting our fits, and all my collegues are likewise perplexed.

Short story - we’re doing an analysis of radiative charm decays, and performing three-dimensional fits with the following observables: Mass of a D0 meson, difference in mass between D* and D0 and a helicity angle between one of D0 daughters and the original particle. While the latter factorises nearly, M(D0) and DM = M(D*)-M(D0) are obviously correlated.

We use MC to perform fits to the DM in bins of M(D0) [our dominant peaking background has a radiative tail in low mass that cancels out to the first order in DM, so it’s easier to parametrise DM rather than M(D0) ] to determine effective resolution in DM as a function of mass.

So, the final PDF is then:

PDF(M0) X PDF(DM|M0) X PDF(helicity angle)

The simplest MC-based models that were giving a passable description of Run1 data were just triple Gaussians. Normally, you’d expect that addition and miltiplication are commutative:

[G1(DM|M0)+G2(DM|M0)+G3(DM|M0)] X F(M0) = G1(DM|M0)XF(M0) + G2(DM|M0)XF(M0) + G3(DM|M0)XF(M0).

However, in practice we observe that explicitly defining RooProdPdf for each sub-component of DM triple Gaussian seems to result in a convergent fit, while constructing triple Gaussians first and then multiplying them results in a resolution parameters always hitting the limit and the overall fit diverging.

I’m attaching a mininal working example using a sample of Run1 MC after pre-selection completed, just uncomment one or the other fitter if you want to see the difference.

Minimal_WE.C (7.8 KB)

Would be nice to hear from experts if there’s some subtlety in statistics we don’t quite understand or if there’s indeed a bug in RooFit implementation.

Best wishes,

Aleksei