Problem with RooGenericPdf fit

Roofit and RooStats
1

Hi everybody,
I’m pretty new of RooFit. I’m trying to fit a dataset with a RooGenericPdf but the parameters of the function remain fixed at the initial value and their errors are too small (about 8 magnitude order smaller…). I give you the code:

Float_t alpha=1.2;
Float_t beta=3.4;
Float_t min=0.900;

RooRealVar aalpha("aalpha","aalpha", alpha,0.40,1.80);
RooRealVar bbeta("bbeta","bbeta", beta, 1.00,11.00);
RooRealVar mmin("mmin","mmin",min,0.700,1.200);

RooGenericPdf Bk("Bk","Generic PDF 2","((Mass-mmin)**aalpha)*((2.100-Mass)**bbeta)", RooArgSet(*Mass,aalpha,bbeta,mmin) );

RooRealVar Bk_nbkg("nbkg","#background events",3500,0.,6000);

RooAddPdf low_sum("low_sum","low_sum",RooArgList( psi2S_sigModel, square, Bk  ), RooArgList( psi2S_nsig, square_nsig, Bk_nbkg ) );     // psi2S_sigModel and square are previously defined

low_sum.fitTo(*dataw_phi,Range(1.070,2.070),PrintLevel(-1));
low_sum.plotOn(frame2);
low_sum.plotOn(frame2,Components(Bk),LineStyle(kDashed));
low_sum.paramOn(frame2);
frame2->Draw();

Is there any error in the code or I badly define my function?
Thank you in advance.

Emiliano

2

Hi Emiliano,

your reproducer is not complete. What is “Mass” for example?

Cheers,
Danilo

3

Thank you for your replay, Mass is a RooRealVar and in particular is the B mass (my data). What do you mean my reproducer is not complete?

Thanks again,
Emiliano

4

Hi Emiliano,

I mean that it would be very useful if you could post a stand alone macro/program that reproduces your issue that other people could try out. Your is incomplete since the absence some variables, for example Mass, prevents others to run your example and help you debugging.

Cheers,
Danilo

5

Hi again,
it took few time but finally I produced a stand alone macro: I had to cut many lines of the original program containing the pre-analysis I made before the fit. The code is:

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooStats/SPlot.h"
#include "RooAbsData.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "TCanvas.h"
#include "RooAbsPdf.h"
#include "RooFit.h"
#include "RooFitResult.h"
#include "RooWorkspace.h"
#include "RooConstVar.h"
#include "TMath.h"
#include "TTree.h"
#include "TH1F.h"
#include "TRandom.h"
#include "RooPolynomial.h"
#include "RooDataHist.h"
#include "RooPlot.h"
#include "TH1D.h"

using namespace RooFit ;
using namespace RooStats ;
#include <iostream>
#include <vector>
using namespace std;
using namespace ROOT::Math;
#pragma link C++ class vector<double>+;

void AddModel(RooWorkspace*);
void AddData(RooWorkspace*);
void DoSPlot(RooWorkspace*);
void MakePlots(RooWorkspace*);

void test(){

  TFile *f = new TFile("sPlot_to_fit.root","read"); 
  RooWorkspace *ws = (RooWorkspace *) f->Get("nWS") ;
  f->Close();
  RooRealVar* Mass = ws->var("Mass");
  RooDataSet * data = (RooDataSet*) ws->data("dataWithSWeights");
  RooDataSet * dataw_phi = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"nsig_sw") ;

  TCanvas* cdata1 = new TCanvas("sPlot1","sPlot1", 600, 450);
  RooDataSet * dataw_phi = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"nsig_sw") ;
  RooPlot* frame2 = Mass->frame() ; 
  dataw_phi->plotOn(frame2, Binning(150), DataError(RooAbsData::SumW2) ) ; 
    
  frame2->SetTitle("LHCb preliminary");
  frame2->Draw() ;
  
  RooRealVar  psi2S_mean("psi2S_mean","mean of the signal",1.660,1.600,1.800);
  RooRealVar  psi2S_width("psi2S_width","width/sigma of the signal",0.0512,0.020,0.070);
  RooGaussian psi2S_sigModel("psi2S_sigModel","signal model",*Mass,psi2S_mean,psi2S_width);

  RooRealVar par0("par0","par0", 0.035 , 0. , 0.150);
  RooRealVar par1("par1","par1", 0.0106 , 0. , 0.027 );
  RooRealVar par2("par2","par2", 1.995 , 1.900 ,2.005 );
  RooRealVar par3("par3","par3", 1.686 , 1.500 ,1.800 );

  par3.setConstant(kTRUE);
  // par2.setConstant(kTRUE);
  par1.setConstant(kTRUE);
  par0.setConstant(kTRUE);

  RooGenericPdf square("square","Generic PDF","0.5*par0*(tanh((Mass-par3)/par1)-tanh((Mass-par2)/par1))", RooArgSet(*Mass,par0,par1,par2,par3));

  RooRealVar psi2S_nsig("psi2S_nsig","#signal events",2000,100,10000);
  RooRealVar square_nsig("square_nsig","square_nsig",1000,0.,300000);
 
  Double_t alpha=1.2;
  Double_t beta=3.4;
  Double_t min=0.900;

  RooRealVar aalpha("aalpha","aalpha", alpha,0.40,1.80);
  RooRealVar bbeta("bbeta","bbeta", beta, 0.0,2.00);
  RooRealVar mmin("mmin","mmin",min,0.700,1.200);

  RooGenericPdf Bk("Bk","Generic PDF 2","((  (Mass-mmin))**aalpha)*(( (2.100-Mass) )**bbeta)", RooArgSet(*Mass,aalpha,bbeta,mmin) );
  RooRealVar Bk_nbkg("nbkg","#background events",3400,3100,3800);

  RooAddPdf low_sum("low_sum","low_sum",RooArgList( psi2S_sigModel, square, Bk  ), RooArgList( psi2S_nsig, square_nsig, Bk_nbkg ) );

  low_sum.fitTo(*dataw_phi,Range(1.070,2.070),PrintLevel(-1));

  low_sum.plotOn(frame2);
  Double_t chi2 = frame2->chiSquare();

  low_sum.paramOn(frame2); 
  low_sum.plotOn(frame2,Components(psi2S_sigModel),LineColor(kRed));
  low_sum.plotOn(frame2,Components(square),LineColor(kCyan));
  low_sum.plotOn(frame2,Components(Bk),LineStyle(kDashed));

  frame2->Draw();
}

Thank you again,
Emiliano
sPlot_to_fit.root (519 KB)

6

Hi Emiliano,

are you sure that everything is well defined in your model?
I get this error:

[#0] WARNING:Minization -- RooFitGlue: Minimized function has error status.
Returning maximum FCN so far (-46865) to force MIGRAD to back out of this region. Error log follows
Parameter values: aalpha=1.22103, bbeta=3.4, mmin=0.9, nbkg=3400, par2=1.995, psi2S_mean=1.66, psi2S_nsig=2000, psi2S_width=0.0512, square_nsig=1000
RooGenericPdf::Bk[ actualVars=(Mass,aalpha,bbeta,mmin) formula="((  (Mass-mmin))**aalpha)*(( (2.100-Mass) )**bbeta)" ]
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.5,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.125,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.3125,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.21875,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.40625,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.17188,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.26562,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.35938,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.45312,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.10156,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.14844,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ actualVars=(Mass = 2.19531,aalpha = 1.22103,bbeta = 3.4,mmin = 0.9)
    ... (remaining 8 messages suppressed)

I had to manipulate 2 or 3 lines in your code, please find it here:

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooStats/SPlot.h"
#include "RooAbsData.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooGenericPdf.h"
#include "RooProdPdf.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "TCanvas.h"
#include "RooAbsPdf.h"
#include "RooFit.h"
#include "RooFitResult.h"
#include "RooWorkspace.h"
#include "RooConstVar.h"
#include "TMath.h"
#include "TTree.h"
#include "TH1F.h"
#include "TRandom.h"
#include "RooPolynomial.h"
#include "RooDataHist.h"
#include "RooPlot.h"
#include "TH1D.h"
#include "TFile.h"

using namespace RooFit ;
using namespace RooStats ;
#include <iostream>
#include <vector>
using namespace std;
// #include "MathCore.h"
// using namespace ROOT::Math;
// #pragma link C++ class vector<double>+;

// void AddModel(RooWorkspace*);
// void AddData(RooWorkspace*);
// void DoSPlot(RooWorkspace*);
// void MakePlots(RooWorkspace*);

void test(){

  TFile *f = new TFile("sPlot_to_fit.root","read"); 
  RooWorkspace *ws = (RooWorkspace *) f->Get("nWS") ;
  f->Close();
  RooRealVar* Mass = ws->var("Mass");
  RooDataSet * data = (RooDataSet*) ws->data("dataWithSWeights");
  RooDataSet * dataw_phi = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"nsig_sw") ;

  TCanvas* cdata1 = new TCanvas("sPlot1","sPlot1", 600, 450);
  RooPlot* frame2 = Mass->frame() ; 
  dataw_phi->plotOn(frame2, Binning(150), DataError(RooAbsData::SumW2) ) ; 
    
  frame2->SetTitle("LHCb preliminary");
  frame2->Draw() ;
  
  RooRealVar  psi2S_mean("psi2S_mean","mean of the signal",1.660,1.600,1.800);
  RooRealVar  psi2S_width("psi2S_width","width/sigma of the signal",0.0512,0.020,0.070);
  RooGaussian psi2S_sigModel("psi2S_sigModel","signal model",*Mass,psi2S_mean,psi2S_width);

  RooRealVar par0("par0","par0", 0.035 , 0. , 0.150);
  RooRealVar par1("par1","par1", 0.0106 , 0. , 0.027 );
  RooRealVar par2("par2","par2", 1.995 , 1.900 ,2.005 );
  RooRealVar par3("par3","par3", 1.686 , 1.500 ,1.800 );

  par3.setConstant(kTRUE);
  // par2.setConstant(kTRUE);
  par1.setConstant(kTRUE);
  par0.setConstant(kTRUE);

  RooGenericPdf square("square","Generic PDF","0.5*par0*(tanh((Mass-par3)/par1)-tanh((Mass-par2)/par1))", RooArgSet(*Mass,par0,par1,par2,par3));

  RooRealVar psi2S_nsig("psi2S_nsig","#signal events",2000,100,10000);
  RooRealVar square_nsig("square_nsig","square_nsig",1000,0.,300000);
 
  Double_t alpha=1.2;
  Double_t beta=3.4;
  Double_t min=0.900;

  RooRealVar aalpha("aalpha","aalpha", alpha,0.40,1.80);
  RooRealVar bbeta("bbeta","bbeta", beta, 0.0,2.00);
  RooRealVar mmin("mmin","mmin",min,0.700,1.200);

  RooGenericPdf Bk("Bk","Generic PDF 2","((  (Mass-mmin))**aalpha)*(( (2.100-Mass) )**bbeta)", RooArgSet(*Mass,aalpha,bbeta,mmin) );
  RooRealVar Bk_nbkg("nbkg","#background events",3400,3100,3800);

  RooAddPdf low_sum("low_sum","low_sum",RooArgList( psi2S_sigModel, square, Bk  ), RooArgList( psi2S_nsig, square_nsig, Bk_nbkg ) );

  low_sum.fitTo(*dataw_phi,Range(1.070,2.070),PrintLevel(-1),SumW2Error(kTRUE));

  low_sum.plotOn(frame2);
  Double_t chi2 = frame2->chiSquare();

  low_sum.paramOn(frame2); 
  low_sum.plotOn(frame2,Components(psi2S_sigModel),LineColor(kRed));
  low_sum.plotOn(frame2,Components(square),LineColor(kCyan));
  low_sum.plotOn(frame2,Components(Bk),LineStyle(kDashed));

  frame2->Draw();
}
7

Thanks for your replay. I fix the error but unfortunately the results are the same. I give you the modified code:

 #ifndef __CINT__
    #include "RooGlobalFunc.h"
    #endif
    #include "RooRealVar.h"
    #include "RooStats/SPlot.h"
    #include "RooAbsData.h"
    #include "RooDataSet.h"
    #include "RooRealVar.h"
    #include "RooGaussian.h"
    #include "RooExponential.h"
    #include "RooChebychev.h"
    #include "RooAddPdf.h"
    #include "RooGenericPdf.h"
    #include "RooProdPdf.h"
    #include "RooAddition.h"
    #include "RooProduct.h"
    #include "TCanvas.h"
    #include "RooAbsPdf.h"
    #include "RooFit.h"
    #include "RooFitResult.h"
    #include "RooWorkspace.h"
    #include "RooConstVar.h"
    #include "TMath.h"
    #include "TTree.h"
    #include "TH1F.h"
    #include "TRandom.h"
    #include "RooPolynomial.h"
    #include "RooDataHist.h"
    #include "RooPlot.h"
    #include "TH1D.h"
    #include "TFile.h"

    using namespace RooFit ;
    using namespace RooStats ;
    #include <iostream>
    #include <vector>
    using namespace std;
    // #include "MathCore.h"
    // using namespace ROOT::Math;
    // #pragma link C++ class vector<double>+;

    // void AddModel(RooWorkspace*);
    // void AddData(RooWorkspace*);
    // void DoSPlot(RooWorkspace*);
    // void MakePlots(RooWorkspace*);

    void test(){

      TFile *f = new TFile("sPlot_to_fit.root","read");
      RooWorkspace *ws = (RooWorkspace *) f->Get("nWS") ;
      f->Close();
      RooRealVar* Mass = ws->var("Mass");
      RooDataSet * data = (RooDataSet*) ws->data("dataWithSWeights");
      RooDataSet * dataw_phi = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"nsig_sw") ;

      TCanvas* cdata1 = new TCanvas("sPlot1","sPlot1", 600, 450);
      RooPlot* frame2 = Mass->frame() ;
      dataw_phi->plotOn(frame2, Binning(150), DataError(RooAbsData::SumW2) ) ;
       
      frame2->SetTitle("LHCb preliminary");
      frame2->Draw() ;
     
      RooRealVar  psi2S_mean("psi2S_mean","mean of the signal",1.660,1.600,1.800);
      RooRealVar  psi2S_width("psi2S_width","width/sigma of the signal",0.0512,0.020,0.070);
      RooGaussian psi2S_sigModel("psi2S_sigModel","signal model",*Mass,psi2S_mean,psi2S_width);

      RooRealVar par0("par0","par0", 0.035 , 0. , 0.150);
      RooRealVar par1("par1","par1", 0.0106 , 0. , 0.027 );
      RooRealVar par2("par2","par2", 1.995 , 1.900 ,2.005 );
      RooRealVar par3("par3","par3", 1.686 , 1.500 ,1.800 );

      par3.setConstant(kTRUE);
      // par2.setConstant(kTRUE);
      par1.setConstant(kTRUE);
      par0.setConstant(kTRUE);

      RooGenericPdf square("square","Generic PDF","0.5*par0*(tanh((Mass-par3)/par1)-tanh((Mass-par2)/par1))", RooArgSet(*Mass,par0,par1,par2,par3));

      RooRealVar psi2S_nsig("psi2S_nsig","#signal events",2000,100,10000);
      RooRealVar square_nsig("square_nsig","square_nsig",1000,0.,300000);
     
      Double_t alpha=0.6;  //////modified value
      Double_t beta=1.7;   //////modified value
      Double_t min=0.900;  //////modified value

      RooRealVar aalpha("aalpha","aalpha", alpha,0.40,1.80);
      RooRealVar bbeta("bbeta","bbeta", beta, 0.0,2.00);
      RooRealVar mmin("mmin","mmin",min,0.700,1.200);

      RooGenericPdf Bk("Bk","Generic PDF 2","(( (Mass-mmin)^2)**aalpha)*(( (2.100-Mass)^2 )**bbeta)", RooArgSet(*Mass,aalpha,bbeta,mmin) );   //////modified PDF
      RooRealVar Bk_nbkg("nbkg","#background events",3400,3100,3800);

      RooAddPdf low_sum("low_sum","low_sum",RooArgList( psi2S_sigModel, square, Bk  ), RooArgList( psi2S_nsig, square_nsig, Bk_nbkg ) );

      low_sum.fitTo(*dataw_phi,Range(1.070,2.070),PrintLevel(-1),SumW2Error(kTRUE));

      low_sum.plotOn(frame2);
      Double_t chi2 = frame2->chiSquare();

      low_sum.paramOn(frame2);
      low_sum.plotOn(frame2,Components(psi2S_sigModel),LineColor(kRed));
      low_sum.plotOn(frame2,Components(square),LineColor(kCyan));
      low_sum.plotOn(frame2,Components(Bk),LineStyle(kDashed));

      frame2->Draw();
    }