Strange behaviour for multidimensional PDF with conditional observables

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

I guess @moneta can help.

Hi @achernov,

sorry for the late reply!

Yes, there is indeed something wrong in the RooFit implementation then, more specifically with fitTo.

I tried to generate toy datasets from both versions of your model, and they are identical. It’s only the fitting that goes wrong in case 2.

The problem is that when you do a conditional fit, you generally need to specify what are the conditional observables:

total.fitTo(*_arr, ConditionalObservables(DMass));

Only then does RooFit know what do normalize the pdf over.

In case 2, you where just lucky because you hit a code path where the pdfs were evaluated in such a order that the caching of normalization sets somehow worked out to get it right in the end, but that was just luck :slight_smile:

Actually I also like your “case 2” formulation more, it’s more concise.

Thanks!
Jonas

Dear @jonas

Thank you for a reply. “Case 2” was the original one, and that’s one we had problems with. Case 1 was found out to not reproduce the same error, and that’s when I started wondering why these two mathemaitcally equivalent ways of constructing composite PDFs end with very different results.

I’m afraid I don’t quite understand your solution, even if it seems to work. The pdf is declared conditional in the mass difference between D0 and D* (aka dM), since resolution in that variable is a function of mass.

Or does ConditionalObservables() argument in fitTo function work like Conditional(pdf,observables, depsAreCond=true) (latter argument defaulting to 0, and it’s not set to 1 in the example)?

Hi @achernov, thanks for following up!

Hmmm I need to study again in detail what the ConditionalObservables() does…

Anyway, I wanted to let you know that I came across a similar problem, which might fundamentally be the same that you have! I realized that when I create a RooAddPdf of products, I get inconsistent results:

using namespace RooFit;

// Create observables
RooRealVar x("x", "x", -5, 5);
RooRealVar y("y", "y", -5, 5);

// Gaussian signal in x and y
RooGaussian gx("gx", "gx", x, RooConst(0), RooConst(1));
RooGaussian gy("gy", "gy", y, RooConst(0), RooConst(1));

// Polynomial background in x and y
RooPolynomial px("px", "px", x, RooArgSet(-0.1, 0.004));
RooPolynomial py("py", "py", y, RooArgSet(-0.1, 0.004));

// Normalization sets to try for model evaluation
RooArgSet nsetx{x};
RooArgSet nsety{y};
RooArgSet nsetxy{x, y};

{
   // Formulate model as a product of sums
   RooAddPdf modelx("modelx", "modelx", {gx, px}, RooConst(0.1));
   RooAddPdf modely("modely", "modely", {gy, py}, RooConst(0.1));
   RooProdPdf model("model", "model", {modelx, modely});

   std::cout << "Product of sums:" << std::endl;

   std::cout << "p(x)        = " << modelx.getVal(nsetx) << std::endl;
   std::cout << "p(y)        = " << modely.getVal(nsety) << std::endl;
   std::cout << "p(x) * p(y) = " << model.getVal(nsetxy) << std::endl;
   std::cout << std::endl;
}

{
   // Formulate model as a sum of products
   RooProdPdf sig("sig", "sig", {gx, gy});
   RooProdPdf bkg("bkg", "bkg", {px, py});
   RooAddPdf model("model", "model", {sig, bkg}, RooConst(0.1));

   std::cout << "Sum of products:" << std::endl;

   std::cout << "p(x)        = " << model.getVal(nsetx) << std::endl;
   std::cout << "p(y)        = " << model.getVal(nsety) << std::endl;
   std::cout << "p(x) * p(y) = " << model.getVal(nsetxy)
                  << " -> wrong value, should be p(x) * p(y)!" << std::endl;
   std::cout << std::endl;
}

This snippet gives the following ouput:

Processing simple_repro.C...
Product of sums:
p(x)        = 0.126991
p(y)        = 0.126991
p(x) * p(y) = 0.0161267

Sum of products:
p(x)        = 0.126991
p(y)        = 0.126991
p(x) * p(y) = 0.0243442 -> wrong value, should be p(x) * p(y)!

The model that is formulated as a sum of products is wrong, as it doesn’t factorize correctly. I think this is a bug in RooFit that I’ll further investigate. I hope your problem is related and will also be fixed then! In any case, I will keep an eye on your code, use it to investigate the problem, and keep you updated :+1:

After reading through documentation for declaration of Conditional Observables for product of PDFs and for the fit function, I think they (with default parameters) mean the opposite - normalize over all observables except ones in the set/normalize only over observables in the set. So, if one naively declares a Conditional(observable) and then uses ConditionalObservables(observable), it kills normalization and of course results in an error.

Passing an optional argument to a PDF product solves this problem, however, in the current code I’m still looking at different behaviour between sum of products and product of sums, even after using @jonas suggestions… will post a minimal working example later once I’m sure there’s nothing else in the current code that could cause similar problems.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.