Unable to plot pull distribution of Nsignal when fitting model1(signal+baground) to background only MC toy data

ataxeidi · September 29, 2024, 7:12pm

Dear experts,
In my macro i have model_0 ,which respesents background and its a composite model of both fractional and composite models and model_1 which represent signal+background .I perform 3 fits with 3 RooMCstudies.I n the 2nd case where fit model1 to generated data from model_0,pull distribution is not showing up,just an empty canvas and i think it has to do with output messege :

RooAbsMinimizerFcn: Minimized function has error status.
Returning maximum FCN so far (-inf) to force MIGRAD to back out of this region. Error log follows.
Parameter values: 	pullMean=-2	pullSigma=1
RooNLLVar::nll_pullGauss_fitParData_model_1_model_0[ parameters=(pullMean,pullSigma) ]
     function value is NAN @ parameters=(pullMean = -2 +/- 0,pullSigma = 1 +/- 0)
RooGaussian::pullGauss[ x=Nexp_sigpull mean=pullMean sigma=pullSigma ]
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4196.94, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4861.19, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4353.27, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4397.93, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4228.23, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4741.67, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4035.26, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4347.29, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4352.36, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4402.53, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4405.51, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4598.3, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0
    ... (remaining 2038 messages suppressed)

and so on.
All other pulls seem to be ok and fitting works well,even thought i get many error messenges. Do you have any hint about what causes this beheaviour and how to avoid it?
I attach my code as well as output log
Thank you in advance,
Artemis.
fit2.C (37.6 KB)

Danilo · September 30, 2024, 5:14am

Hi,

Thanks for sharing also the macro. Would you be able to reduce the code to something minimal (10-20 lines)? Right now it’s about 1000 lines which makes it hard to read.

Cheers,
Danilo

ataxeidi · September 30, 2024, 12:34pm

Nsig=h_m20->Integral();
Ntt1=h_tt1->Integral();
Ntt4=h_tt4->Integral();
Ntt5=h_tt5->Integral();
Ntt=Ntt1+Ntt4+Ntt5;

Nw=h_w->Integral();;
Nqcd=h_qcd->Integral();
Nw_qcd=Nw+Nqcd;
Nznn=h_znn->Integral();
Nbkg=Ntt1+Ntt4+Ntt5+Nznn+Nw+Nqcd;
Nsb=Nsig+Nbkg;

RooRealVar output_BDT("output_BDT","BDT score",-0.9,0.6);

RooRealVar Nexp_sig("Nexp_sig","Expected number of signal events",Nsig,-3*Nsig,3*Nsig);
  RooRealVar Nexp_sig_th("Nexp_sig_th","Expected number of signal events",Nsig);

  
  RooRealVar Nexp_tt("Nexp_tt", "Expected number of TT bkg", Ntt, -2*Ntt, 2*Ntt);
   RooRealVar Nexp_tt_th("Nexp_tt_th", "Expected number of TT  bkg", Ntt);
    RooRealVar Nexp_w_qcd("Nexp_w_qcd", "Expected number of W+QCD bkg", Nw_qcd, -2Nw_qcd, 2*Nw_qcd);
    RooRealVar Nexp_w_qcd_th("Nexp_w_qcd_th", "Expected number of  W+QCD bkg", Nw_qcd);
    RooRealVar Nexp_znn("Nexp_znn", "Expected number of Zvv bkg", Nznn,-2*Nznn,2*Nznn);
     RooRealVar Nexp_znn_th("Nexp_znn", "Expected number of Zvv bkg", Nznn);
float f_tt1=Ntt1/(Ntt1+Ntt4+Ntt5);
    float f_tt4=Ntt4/(Ntt1+Ntt4+Ntt5);
    float f_tt5=Ntt5/(Ntt1+Ntt4+Ntt5);
    float f_w=Nw/(Nw+Nqcd);
// bkg fractions
    RooRealVar fr_tt1("fraction_tt1", "fraction_tt1",f_tt1);
    RooRealVar fr_tt4("fraction_tt4", "fraction_tt4",f_tt4);
    RooRealVar fr_tt5("fraction_tt5", "fraction_tt5",f_tt5);
    RooRealVar fr_w("fraction_w", "fraction_w",f_w);
fr_tt1.setConstant(kTRUE);
    fr_tt4.setConstant(kTRUE);
    fr_tt5.setConstant(kTRUE);
    fr_w.setConstant(kTRUE);

RooDataHist sig20("sig20","sig20",output_BDT,h_m20_r);
RooDataHist tt1("tt","tt1",output_BDT,h_tt1_r);
    RooDataHist tt4("tt","tt4",output_BDT,h_tt4_r);
    RooDataHist tt5("tt","tt5",output_BDT,h_tt5_r);
   RooDataHist w("w","w",output_BDT,h_w_r);
    RooDataHist qcd("qcd","qcd",output_BDT,h_qcd_r);
  RooDataHist znn("znn","znn",output_BDT,h_znn_r);

RooHistPdf znn_pdf("znn_pdf","znn_pdf",output_BDT,znn);
  RooHistPdf w_pdf("w_pdf","w_pdf",output_BDT,w);
  RooHistPdf qcd_pdf("qcd_pdf","qcd_pdf",output_BDT,qcd);
  RooHistPdf tt1_pdf("tt1_pdf","tt1_pdf",output_BDT,tt1);
   RooHistPdf tt4_pdf("tt4_pdf","tt4_pdf",output_BDT,tt4);
    RooHistPdf tt5_pdf("tt5_pdf","tt5_pdf",output_BDT,tt5);
RooHistPdf sig20("sig20","sig20",output_BDT,sig20);

  //define  bkg frac
   RooAddPdf w_qcd_pdf("w_qcd_pdf", "w_qcd_pdf", RooArgList(w_pdf, qcd_pdf), RooArgList(fr_w));
   RooAddPdf tt_pdf("tt_pdf", "tt_pdf", RooArgList(tt1_pdf,tt4_pdf,tt5_pdf), RooArgList(fr_tt1,fr_tt4));
   //define total bkg
   RooAddPdf model_0("model_0", "Background", RooArgList( tt_pdf,w_qcd_pdf, znn_pdf), RooArgList(  Nexp_tt,Nexp_w_qcd,Nexp_znn));
   //define s+b
   RooAddPdf model_1("model_1", "Signal + Background", RooArgList(sig_pdf,tt_pdf,w_qcd_pdf, znn_pdf), RooArgList(Nexp_sig,Nexp_tt,Nexp_w_qcd,Nexp_znn));

   RooDataSet *data_B = model_0.generate(output_BDT,Nbkg);
   RooDataSet *data_SB = model_1.generate(output_BDT,Nsb);

RooMCStudy * mcstudy_0_1 = new RooMCStudy(model_0,output_BDT,Binned(kTRUE), FitModel(model_1),EvalErrorWall(false), FitOptions(Save(kTRUE),PrintEvalErrors(0)),Extended());
   mcstudy_0_1->generateAndFit(Ntoys);
  
 TCanvas *c_tt_01 = new TCanvas("c_tt_01", " model 1 fit on B-only data", 1600, 600);
 c_tt_01 ->Divide(3,1);
c_tt_01 ->cd(1);
 RooPlot* framett_01_1 = mcstudy_0_1->plotParam(Nexp_tt, Bins(50));
 framett_01_1->Draw();
 c_tt_01 ->cd(2);
 RooPlot* framett_01_2= mcstudy_0_1->plotError(Nexp_tt, Bins(50));
 framett_01_2->Draw();
c_tt_01 ->cd(3);
 RooPlot* framett_01_3 = mcstudy_0_1->plotPull(Nexp_tt, Bins(50),FitGauss(kTRUE)) ;
 framett_01_3->Draw();
 
  TCanvas *c_wqcd_01 = new TCanvas("c_wqcd_01", " model 1 fit on B-only data", 1600, 600);
 c_wqcd_01 ->Divide(3,1);
c_wqcd_01 ->cd(1);
 RooPlot* framewqcd_01_1 = mcstudy_0_1->plotParam(Nexp_w_qcd, Bins(50));
 framewqcd_01_1->Draw();
 c_wqcd_01 ->cd(2);
 RooPlot* framewqcd_01_2= mcstudy_0_1->plotError(Nexp_w_qcd, Bins(50));
 framewqcd_01_2->Draw();
c_wqcd_01 ->cd(3);
 RooPlot* framewqcd_01_3 = mcstudy_0_1->plotPull(Nexp_w_qcd, Bins(50),FitGauss(kTRUE)) ;
 framewqcd_01_3->Draw();
 

  TCanvas *c_znn_01 = new TCanvas("c_znn_01"," model 1 fit on B-only data", 1600, 600);
 c_znn_01 ->Divide(3,1);
c_znn_01 ->cd(1);
 RooPlot* frameznn_01_1 = mcstudy_0_1->plotParam(Nexp_znn, Bins(50));
 frameznn_01_1->Draw();
 c_znn_01 ->cd(2);
 RooPlot* frameznn_01_2= mcstudy_0_1->plotError(Nexp_znn, Bins(50));
 frameznn_01_2->Draw();
c_znn_01 ->cd(3);
 RooPlot* frameznn_01_3 = mcstudy_0_1->plotPull(Nexp_znn, Bins(50),FitGauss(kTRUE)) ;
 frameznn_01_3->Draw();
  
     TCanvas *c_1s = new TCanvas("c_1s", "Nsig: model 1 fit on B-only data ", 1600, 600);
   c_1s->Divide(3,1);
   c_1s->cd(1);
   RooPlot* frameM11_1 = mcstudy_0_1->plotParam(Nexp_sig,Bins(50));
   // mcstudy_0_1.fitParDataSet().statOn(frameM11_1) ;
    frameM11_1->Draw();
    c_1s->cd(2);
    RooPlot* frameM11_2= mcstudy_0_1->plotError(Nexp_sig,Bins(50));
    frameM11_2->Draw();
    c_1s->cd(3);
    RooPlot* frameM11_3 = mcstudy_0_1->plotPull(Nexp_sig,Bins(50)) ;
    frameM11_3->Draw();
     TCanvas *c_1sl = new TCanvas("c_1sl", "NLL: model 1 fit on B-only data ", 1600, 600);
      RooPlot* frameM01_4 = mcstudy_0_1->plotNLL(Bins(100)) ;

ataxeidi · September 30, 2024, 12:36pm

I sent some relevant lines.The pull distribution of Nexp_sig is empty all other plots are fine

ataxeidi · September 30, 2024, 12:44pm

nsig01.pdf (26.8 KB)
nll01.pdf (24.0 KB)

nbkg2_01.pdf (30.0 KB)
nbkg3_01.pdf (30.8 KB)

StephanH · October 1, 2024, 7:32am

Hello @ataxeidi,

when I look at the error messages, it seems that the pull distribution is maybe just plotted in the wrong range:

     getLogVal() top-level p.d.f evaluates to zero @ x=Nexp_sigpull=-4196.94, mean=pullMean=-2 +/- 0, sigma=pullSigma=1 +/- 0

It looks like it’s trying to fit a Gaussian(x, mean=-2, sigma=1) at x=-4196, so that’s very much zero within machine precision. It’s hard to compute the logarithm of zero, and that’s the error you see.

Could it be those lines?

RooPlot* frameznn_00_3 = mcstudy_0_0->plotPull(Nexp_znn, Bins(50),FitGauss(kTRUE)) ;

It seems you are trying to fit a gaussian into the pull distribution. Maybe first try to create the pull distribution without that fit, check if the range of the variable makes sense, and then constrain the range a bit, and repeat the fit.

That’s my best guess without being able to run it myself. Let us know if this helps.