Non stability of fitting in roofit

VDiBella · July 10, 2026, 2:55pm

Hi,

I’m trying to fit double dimension data in roofit. I’m trying to calculate the raw yield.

I’m studying principally two particles, the D+ and D-. My main issue is that the data for the D+ and D- are really similar, but when I’m using the code for the D+ I got a raw yield kinda good, (an efficiency of around 1.5) but when I use the same on D- that is really the same I got a raw yield really bad (an efficiency of around 3) and I don’t understand because the code should be the same for the two.

Can someone help me and have an idea ?

Here is my code :

#include "TF2.h"
#include "TH2.h"
#include "TMath.h"

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"


TH1 *makeTH1(TRandom &trnd);
TTree *makeTTree(TRandom &trnd);

TString DataFile = Form("Root_Files/AnalysisResults_symmetric.root");
TString Table = Form("hf-task-correlation-dplus-dplus-reduced");

//Choose the data we will look into : 1 for Dplus, 2 : for Dminus, 3 : for Dplusminus
int Choose_type = 1;
bool Fitting_on_matched = false;

TString Data_type = Form("");
if (Choose_type == 1) {Data_type = Form("DplusPair");}
else if (Choose_type == 2) {Data_type = Form("DminusPair");}
else if (Choose_type == 3) {Data_type = Form("DplusminusPair");}
else {
  cout<<"Please choose a correct choose_type, it should be 1, 2 or 3"<<endl;
  exit(0);
}
TString Data_name = Form("hMass");
Data_name = Data_name + Data_type;
if (Fitting_on_matched){Data_name = Data_name + "sMatched";}
TString Save_name = Form("Root_Files/Results");
Save_name = Save_name + Data_type;
Save_name = Save_name + Form(".root");


void Fitting() {
  // O p e n   t h e   d a t a   f i l e
  // -------------------------------------------------------------------
	TFile *myFile = new TFile(DataFile);
	TDirectoryFile *Folder = (TDirectoryFile*)myFile->Get(Table);
	TH2F *Graph = (TH2F*)Folder->Get(Data_name.Data());
  using namespace RooFit;


  // P r e p a r a t i o n   o f   t h e   f i t t i n g
  // -------------------------------------------------------------------

  //Creation of the variables
  RooRealVar invmassx("invmassx", "inv_mass_first_candidate (Gev/c²)", 1.80, 1.95);
  RooRealVar invmassy("invmassy", "inv_mass_second_candidate (Gev/c²)", 1.80, 1.95);

  RooRealVar mean1("mean1", "mean", 1.87113, 1.85, 1.88);
  RooRealVar sigma1("sigma1", "sigma", 0.00903759, 0.007, 0.015);
  RooRealVar sigma2("sigma2", "sigma", 0.00903917, 0.007, 0.015);
  RooRealVar bkgfirst1("bkgfirst1", "background first parameter for back sig", 295.624, -100, 1000);
  RooRealVar bkgsecond1("bkgsecond1", "background second parameter for back sig", 37.6065, -1, 300);
  RooRealVar bkgthird1("bkgthird1", "background third parameter for back sig", -69.8315, -300, 2);
  RooRealVar bkgfirst2("bkgfirst2", "background first parameter for back back", 280.982, -100, 1000);
  RooRealVar bkgsecond2("bkgsecond2", "background second parameter for back back", 24.063, -1, 300);
  RooRealVar bkgthird2("bkgthird2", "background third parameter for back back", -63.3258, -300, 2);

  // mean1.setConstant(kTRUE);
  // sigma1.setConstant(kTRUE);
  // sigma2.setConstant(kTRUE);
  // bkgfirst1.setConstant(kTRUE);
  // bkgsecond1.setConstant(kTRUE);
  // bkgthird1.setConstant(kTRUE);
  // bkgfirst2.setConstant(kTRUE);
  // bkgsecond2.setConstant(kTRUE);
  // bkgthird2.setConstant(kTRUE);


  //Normalisation variable for the additional pdf
  RooRealVar nsig("nsig", "#signal events", 858., 0., 30000);
  RooRealVar nsigbkg("nsigbkg", "#signal events", 83598., 100., 1000000);
  RooRealVar nbkgsig("nbkgsig", "#signal events", 83598., 100., 1000000);
  RooRealVar nbkg("nbkg", "#signal events", 449101., 0., 1000000);

  //Creation of the frame
  RooDataHist dh("dh", "dh", {invmassx,invmassy}, Graph);
  RooPlot *frame = new RooPlot(invmassx,invmassy,1.80,1.95,1.80,1.95);
  dh.plotOn(frame,Project({invmassx,invmassy}), Name("dh"));

  //Creation of the function
  RooGaussian gaussX("gaussX", "gauss", invmassx, mean1, sigma1);
  RooPolynomial backgroundX("backgroundX", "Polynomial", invmassx, {bkgfirst1,bkgsecond1,bkgthird1});

  RooGaussian gaussY("gaussY", "gauss", invmassy, mean1, sigma1);
  RooPolynomial backgroundY("backgroundY", "Polynomial", invmassy, {bkgfirst1,bkgsecond1,bkgthird1});

  RooPolynomial backgroundbackX("backgroundbackX", "Polynomial", invmassx, {bkgfirst2,bkgsecond2,bkgthird2});
  RooPolynomial backgroundbackY("backgroundbackY", "Polynomial", invmassy, {bkgfirst2,bkgsecond2,bkgthird2});

  // C o n s t r u c t   u n c o r r e l a t e d   p r o d u c t   p d f
  // -------------------------------------------------------------------

  //Multiply gaussx and gaussy into a two-dimensional p.d.f. gaussxy
  RooProdPdf  gaussxy("gaussxy","gaussX*gaussY",RooArgList(gaussX,gaussY));
  RooProdPdf  gaussxbkgy("gaussxbkgy","gaussX*backgroundY",RooArgList(gaussX,backgroundY));
  RooProdPdf  bkgxgaussy("bkgxgaussy","backgroundX*gaussY",RooArgList(backgroundX,gaussY));
  RooProdPdf  bkgxy("bkgxy","backgroundX*backgroundY",RooArgList(backgroundbackX,backgroundbackY));

  //Addition of all the functions
  RooAddPdf total("total", "g+a", {gaussxy,gaussxbkgy,bkgxgaussy,bkgxy}, {nsig, nsigbkg, nbkgsig, nbkg});
  RooFitResult *fit_result_data = total.fitTo(dh, Save());
  total.plotOn(frame,Project({invmassx,invmassy}), Name("model_curve"));

  TCanvas *c = new TCanvas("rf304_uncorrprod","rf304_uncorrprod",800, 400);
  
  Graph->GetXaxis()->SetRangeUser(1.8, 1.95);
  Graph->GetYaxis()->SetRangeUser(1.8, 1.95);

  
  // I n t e g r a t i o n   o f   t h e   d o u b l e   g a u s s i a n
  // -------------------------------------------------------------------
  invmassx.setRange("range",mean1.getVal()-sigma1.getVal()*3,mean1.getVal()+sigma1.getVal()*3);
  invmassy.setRange("range",mean1.getVal()-sigma1.getVal()*3,mean1.getVal()+sigma1.getVal()*3);

  RooAbsReal *Integral = gaussxy.createIntegral(RooArgSet(invmassx,invmassy),NormSet(RooArgSet(invmassx,invmassy)),Range("range"));

  // G e t   t h e   n u m b e r   o f   m a t c h e d
  // -------------------------------------------------------------------

  int Matched_value = 0;
  if (!Fitting_on_matched) {
  TString Matched_name = Data_name + "sMatched";
	TH2F *Matched_graph = (TH2F*)Folder->Get(Matched_name.Data());
  Matched_value = Matched_graph->GetEntries();
  }
  else {Matched_value = -1;}

  // C r e a t i o n    o f    p r o j e c t i o n    t o    a c c e s s    t h e    c h i s q u a r e
  // -------------------------------------------------------------------

  RooPlot *xframe = invmassx.frame(Title("Projection on x"));
  dh.plotOn(xframe);
  total.plotOn(xframe);

  RooPlot *yframe = invmassy.frame(Title("Projection on y"));
  dh.plotOn(yframe);
  total.plotOn(yframe);

  // P r e s e n t a t i o n   o f   r e s u l t s
  // -------------------------------------------------------------------

  TH2* pdf_2d = (TH2*) total.createHistogram("invmassx,invmassy",Graph->FindBin(1.95,1.80)-Graph->FindBin(1.80,1.80),Graph->FindBin(1.95,1.80)-Graph->FindBin(1.80,1.80)) ;
  pdf_2d->SetLineColor(2);
  Graph->Draw("lego2");
  pdf_2d->Draw("surf same");
  c->Print(Save_name);

  cout<<endl<<endl<<endl<<endl<<"E n d   r e s u l t s :"<<endl<<endl<<endl<<endl;
  cout<<"Number of events :"<<pdf_2d->Integral(pdf_2d->GetXaxis()->FindBin(mean1.getVal()-sigma1.getVal()*3),pdf_2d->GetXaxis()->FindBin(mean1.getVal()+sigma1.getVal()*3),pdf_2d->GetYaxis()->FindBin(mean1.getVal()-sigma1.getVal()*3),pdf_2d->GetYaxis()->FindBin(mean1.getVal()-sigma1.getVal()*3))<<endl;
  cout<<"Sigma :"<<sigma1.getVal()<<"+-"<<sigma1.getError()<<endl;
  cout<<"Mean :"<<mean1.getVal()<<"+-"<<mean1.getError()<<endl;
  cout<<"Normalisation factor :"<<nsig.getVal()<<"+-"<<nsig.getError()<<endl;
  cout<<"MonteCarlo truth : "<<Matched_value<<endl;
  cout<<"Raw yield :"<<nsig.getVal()*Integral->getVal()<<"+-"<<nsig.getError()*Integral->getVal()+nsig.getVal()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})+nsig.getError()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})<<endl;
  cout<<"Efficiency : "<<nsig.getVal()*Integral->getVal()/Matched_value<<"+-"<<(nsig.getError()*Integral->getVal()+nsig.getVal()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})+nsig.getError()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy}))/Matched_value<<endl;
  cout<<"chisquare for x projection :"<<xframe->chiSquare()<<endl;
  cout<<"chisquare for y projection :"<<yframe->chiSquare()<<endl;

  // cout<<endl<<endl<<endl<<endl<<"E n d   r e s u l t s :"<<endl<<endl<<endl<<endl;
  // cout<<"Sigma1 :"<<sigma1.getVal()<<endl;
  // cout<<"Mean1 :"<<mean1.getVal()<<endl;
  // cout<<"Sigma2 :"<<sigma2.getVal()<<endl;
  // cout<<"bkgfirst1 :"<<bkgfirst1.getVal()<<endl;
  // cout<<"bkgsecond1 :"<<bkgsecond1.getVal()<<endl;
  // cout<<"bkgthird1 :"<<bkgthird1.getVal()<<endl;
  // cout<<"bkgfirst2 :"<<bkgfirst2.getVal()<<endl;
  // cout<<"bkgsecond2 :"<<bkgsecond2.getVal()<<endl;
  // cout<<"bkgthird2 :"<<bkgthird2.getVal()<<endl;
  // cout<<"alpha :"<<alpha.getVal()<<endl;
  // cout<<"n :"<<n.getVal()<<endl;
}

As a new user I can’t link my data, tell me if you need them.

Thanks a lot in advance.

bellenot · July 10, 2026, 3:34pm

Let’s see if @jonas can help

jonas · July 12, 2026, 7:15pm

Hi @VDiBella,

Thanks for the post! Nothing in your code treats D⁺ and D⁻ differently, so the difference must come from how the fit converges on the two datasets. Two suggestions:

Use RooChebychev or RooBernstein instead of RooPolynomial for the backgrounds. With your coefficient ranges, the polynomial can become negative inside the fit range during minimization. The PDF is then invalid, and the minimizer has to navigate around these invalid regions of parameter space, which can make it end up in a bad minimum for one dataset and not the other, even if the data are very similar. Chebychev or Bernstein parametrizations are much better behaved in this respect (Bernstein polynomials with positive coefficients are positive by construction).
Please post the full fit result for both the D⁺ and D⁻ fits, i.e. the output of:
```
fit_result_data->Print("v");
```
That prints the fit status, the covariance matrix quality, and all parameter values with their uncertainties. Things to look out for: a parameter that ended up at one of its limits (for example sigma1, which is only allowed in the narrow window [0.007, 0.015] tuned on D⁺), a covariance quality below 3, or an uncertainty on nsig that is comparable to the value itself. Any of these would explain a bad yield.

With that output we can continue the diagnosis.

I hope this helps and we can eventually make the fit converge!
Jonas

VDiBella · July 13, 2026, 4:07pm

Hello, thanks a lot for your answer. I’m not able to test the first proposition today I will do it as possible.

Here are the result you proposed for the second proposition :

For D+ :

RooFitResult: minimized FCN value: -9.75305e+06, estimated distance to minimum: 0.0583859
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=-1 HESSE=3

Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.

         bkgfirst1    2.9562e+02    2.9792e+02 +/-  3.11e+01  <none>
         bkgfirst2    2.8098e+02    2.8110e+02 +/-  1.36e+01  <none>
        bkgsecond1    3.7606e+01    3.6730e+01 +/-  1.40e+01  <none>
        bkgsecond2    2.4063e+01    2.4074e+01 +/-  6.27e+00  <none>
         bkgthird1   -6.9832e+01   -6.9988e+01 +/-  6.21e+00  <none>
         bkgthird2   -6.3326e+01   -6.3333e+01 +/-  2.81e+00  <none>
             mean1    1.8711e+00    1.8711e+00 +/-  4.52e-05  <none>
              nbkg    4.4910e+05    4.4245e+05 +/-  9.81e+02  <none>
           nbkgsig    8.3598e+04    8.0801e+04 +/-  5.70e+02  <none>
              nsig    8.5800e+02    1.2496e+03 +/-  2.94e+02  <none>
           nsigbkg    8.3598e+04    8.0801e+04 +/-  5.70e+02  <none>
            sigma1    9.0376e-03    9.0416e-03 +/-  5.00e-05  <none>

For D- :

RooFitResult: minimized FCN value: -9.95474e+06, estimated distance to minimum: 0.0611924
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=-1 HESSE=3

Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.

         bkgfirst1    2.9562e+02    3.0086e+02 +/-  3.53e+01  <none>
         bkgfirst2    2.8098e+02    2.7995e+02 +/-  1.30e+01  <none>
        bkgsecond1    3.7606e+01    3.9063e+01 +/-  1.49e+01  <none>
        bkgsecond2    2.4063e+01    2.3598e+01 +/-  6.06e+00  <none>
         bkgthird1   -6.9832e+01   -6.9113e+01 +/-  6.86e+00  <none>
         bkgthird2   -6.3326e+01   -6.3537e+01 +/-  2.69e+00  <none>
             mean1    1.8711e+00    1.8712e+00 +/-  4.80e-05  <none>
              nbkg    4.4910e+05    4.4996e+05 +/-  1.01e+03  <none>
           nbkgsig    8.3598e+04    8.2350e+04 +/-  5.92e+02  <none>
              nsig    8.5800e+02    2.5462e+03 +/-  3.16e+02  <none>
           nsigbkg    8.3598e+04    8.2350e+04 +/-  5.92e+02  <none>
            sigma1    9.0376e-03    9.6062e-03 +/-  5.23e-05  <none>

The results seems really similar, the only difference seems to be the nsig, and its precisly the point that matter the most for me.

Again, thanks a lot.

jonas · July 16, 2026, 2:24pm

Hi @VDiBella,

Thanks, this output is actually very informative! The key thing is in the first two lines, not in the parameter table:

estimated distance to minimum: 0.0583859
Status : MINIMIZE=-1 HESSE=3

The MINIMIZE=-1 status and the estimated distance to minimum (EDM) of ~0.06 tell us that neither fit actually converged. For a converged fit the EDM should be below roughly 1e-3. The “Full, accurate covariance matrix” message is not a contradiction: it only means HESSE could compute a good covariance matrix at the point where MIGRAD gave up, not that this point is the minimum. So both results are essentially “wherever the minimizer stalled”, and it’s then not surprising that two very similar datasets end up in different places. Neither nsig value should be trusted at this point.

There is a likely explanation for why MIGRAD stalls: your minimized FCN value is about -9.75e6. With double precision, a number of that magnitude can only be resolved to an absolute precision of ~0.01 to 0.1, which is exactly the size of your EDM. In other words, the minimizer can no longer see the shape of the likelihood through the numerical noise. RooFit has a built-in remedy for this: likelihood offsetting, which subtracts the initial likelihood value so the minimizer works with numbers close to zero:

RooFitResult *fit_result_data = total.fitTo(dh, Save(), Offset(true), PrintLevel(1), Strategy(2));

The Offset(true) is the important part (in recent ROOT versions you can also write Offset("initial")). I added Strategy(2) for a more careful minimization and PrintLevel(1) so you can see what MIGRAD is doing.

Two more suggestions to make the minimizer’s life easier:

In both fit results, nsigbkg and nbkgsig come out exactly identical, including their uncertainties. That makes sense: your 2D histogram is symmetric under exchange of the two axes, so these two yields are maximally correlated and the fit has an (almost) flat direction. Since they must be equal by symmetry anyway, use a single parameter for both:
```
RooRealVar ncross("ncross", "#cross-term events", 83598., 100., 1000000);
RooAddPdf total("total", "g+a", {gaussxy, gaussxbkgy, bkgxgaussy, bkgxy}, {nsig, ncross, ncross, nbkg});
```
The switch from RooPolynomial to RooChebychev from my previous post is still worth doing: invalid (negative) PDF regions during minimization are another common reason for MIGRAD to fail.

Once you have fits with EDM < 1e-3 and status 0, a good stability cross-check is to fit the D- data starting from the D+ result and vice versa. If both converge to the same respective minima, you can be reasonably confident in the yields.

Two small unrelated things I noticed in the code: sigma2 is declared but never used, and in the pdf_2d->Integral(...) call that prints “Number of events”, the last argument uses mean1.getVal()-sigma1.getVal()*3 for the upper edge of the y range (that should be +). And finally, once the fits converge, I’d take a critical look at why the “efficiency” is above 1 for both charges. A raw yield 50% larger than the MC-matched count even in the “good” D+ fit suggests the yield extraction or the matched-count normalization deserves a second look, independently of the fit stability.

I hope this gets the fits to converge!

Jonas

VDiBella · July 16, 2026, 3:02pm

Hello,

Thanks a lot for all the information !

I was trying on implementing the chebychev. But I had some issue for now. I will continue to work on it.

I have apply all your advice (thanks a lot for noticing the error such as the - for the sigma instead of the +).

When i apply your proposition for a correction this does not change anything at all.
The Strategy(2) change the result a bit (~multiply by 2 the EDM for example), but the offset(true) does not change anything.

I have run the code and if I put this line :

RooFitResult *fit_result_data = total.fitTo(dh, Save(), Offset(true));

I have the exact same results for everything as using this line :

RooFitResult *fit_result_data = total.fitTo(dh, Save());

Its important I think to note that if I say that the parameter for the background can’t be negative, the result as worst, but the Minimize became egal to 1.

I can note too that if I add the Strategy(2) for D- the covariance matrix can’t even be calculated.

Do you think I should just try to use the chebychev ? Or there was something I did not understood in your answer ?

Thanks a lot again.

VDiBella · July 22, 2026, 2:23pm

Hello,

I’ve tried a lot of things. I have a fit for some pairs that is quite stable, when trying different starting point my results are in the same error range for example.

But it still imply that this is not converging. Even with the Offset being on true or “initial”.

Do you have an idea of an other solution ?

Thanks a lot for all your help

Edit of my post : It seems that the chebychev function was the solution now I have :

RooFitResult: minimized FCN value: -8.17543e+06, estimated distance to minimum: 0.00109702
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0

Floating Parameter  InitialValue    FinalValue +/-  Error     GblCorr.

         bkgfirst1   -1.2440e-01   -1.2438e-01 +/-  1.29e-02  <none>
         bkgfirst2   -1.1391e-01   -1.1392e-01 +/-  2.45e-03  <none>
        bkgsecond1   -6.6532e-02   -6.6460e-02 +/-  1.34e-02  <none>
        bkgsecond2   -2.9873e-02   -2.9876e-02 +/-  2.77e-03  <none>
         bkgthird1    2.7553e-02    2.7672e-02 +/-  1.12e-02  <none>
         bkgthird2    2.5822e-02    2.5801e-02 +/-  2.20e-03  <none>
             mean1    1.8711e+00    1.8714e+00 +/-  6.80e-05  <none>
              nbkg    4.2055e+05    4.2058e+05 +/-  1.11e+03  <none>
              nsig    3.0000e+04    4.7228e+03 +/-  3.33e+02  <none>
           nsigbkg    4.4117e+04    4.4099e+04 +/-  5.13e+02  <none>
            sigma1    9.0376e-03    8.7610e-03 +/-  8.24e-05  <none>

As a result.

And for the D+ D- issue, it seems I just don’t have enough data, so in the end it should not be able to converge at all. Here is the result for an other set of data with more data.

Thanks a lot again for all your help.
But I’m still curious because the offset was not changing any results.