Dear experts,
I am performing a fit to the sidebands of a linear function, I’m using RooPolynomial.
The problem is that I’m getting some non-sense values when computing the error on the integral over the “signal” region. The error is huge, orders of magnitude larger than the integral itself.
What is even more strange is that this problem doesn’t appear when the slope of the linear function is negative.
I report here the output of the macro
positive_slope_fit.pdf (20.7 KB)
Extrapolation in signal region with RooFit (positive slope), integral = 335.416 ± 4258.17
negative_slope_fit.pdf (20.5 KB)
Extrapolation in signal region with RooFit (negative slope), integral =289.03 ± 14.2727
Do you have any idea of what is happening here?
I’m using ROOT 5.34.18
Thanks for your help,
Andrea
Below is the code
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooPolynomial.h"
#include "RooPlot.h"
#include "RooNLLVar.h"
#include "RooAbsData.h"
#include "RooAbsPdf.h"
#include "RooAbsArg.h"
#include "RooAbsReal.h"
#include "RooDataSet.h"
#include "RooFitResult.h"
#include "RooAddPdf.h"
#include "RooArgList.h"
#include "RooMinuit.h"
using namespace RooFit;
RooDataSet* PolyData(double coeff){
RooRealVar mass("mass", "mass", 0., 1);
RooRealVar c0("c0","c0", coeff);
RooPolynomial poly("poly","poly",mass,RooArgList(c0));
RooDataSet* data = poly.generate(mass, 3e+03);
return data;
}
////////////////////////////////////////////////////////////////////////////////
// this function do the fit //
////////////////////////////////////////////////////////////////////////////////
RooPlot* plot_fit(double coeff){
double min_signal = 0.5;
double max_signal = 0.6;
double low_sb_min = 0.3;
double low_sb_max = min_signal;
double high_sb_min = max_signal;
double high_sb_max = 0.8;
TString str_low_sb, str_high_sb, str_sig;
TCut cut_low_sb, cut_high_sb, cut_sig;
// Declare observables
RooRealVar mass("mass","mass", 0., 1) ;
// Define signal, sideband range
mass.setRange("signal", min_signal, max_signal) ;
mass.setRange("low_sb", low_sb_min, low_sb_max) ;
mass.setRange("high_sb", high_sb_min, high_sb_max) ;
mass.setRange("fullRange", low_sb_min, high_sb_max) ;
// Define coefficient
RooRealVar c0("c0","c0", 0.01, -1, 1e6);
RooPolynomial poly("poly","poly",mass,RooArgList(c0));
RooRealVar norm("norm","number of signal events",2000,1.,1000000) ;
RooExtendPdf extpoly("extpoly","extended signal p.d.f",poly,norm,"fullRange") ;
RooDataSet* data = PolyData(coeff);
// Perform unbinned extended ML fit to data
RooFitResult * rf = extpoly.fitTo(*data,Extended(kTRUE), Range("low_sb,high_sb"), Save()) ;
// LM: you need to use the parameters of the fitted function
RooArgList pars(* extpoly.getParameters(RooArgSet(mass) ) );
//LM: make fitted function (need to add normalized term to pdf)
RooArgSet prodSet(extpoly); //prodSet.add(norm);
RooProduct unNormPdf("fitted Function", "fitted Function", prodSet);
TF1 * f2 = unNormPdf.asTF(RooArgList(mass), pars);
float norm1 = ((RooRealVar*) pars.find("norm"))->getVal();
float dnorm1 = ((RooRealVar*) pars.find("norm"))->getError();
float coef0 = ((RooRealVar*) pars.find("c0"))->getVal();
float dcoef0 = ((RooRealVar*) pars.find("c0"))->getError();
Double_t integ2_full = f2->Integral(low_sb_min, high_sb_max);
Double_t integ2 = norm1*f2->Integral(min_signal, max_signal, 0)/integ2_full;
Double_t dinteg2 = norm1*f2->IntegralError(min_signal, max_signal, 0, rf->covarianceMatrix().GetMatrixArray())/integ2_full; //
RooPlot * xframe = mass.frame(Title("Roofit fit"));
data->plotOn(xframe, Binning(50));
extpoly.plotOn(xframe,Range("low_sb,high_sb"));
extpoly.plotOn(xframe,Range("signal"),LineColor(kRed),LineStyle(kDashed));
extpoly.paramOn(xframe,Layout(0.55));
if(coeff > 0)
cout << endl << "---------------- Positive slope ----------------" << endl << endl;
else
cout << endl << "---------------- Negative slope ----------------" << endl << endl;
cout << " norm = " << norm1 << " +- " << dnorm1 << endl;
cout << " c0 = " << coef0 << " +- " << dcoef0 << endl;
cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << endl;
cout << " Total # of entries in dataset = " << data->sumEntries() << endl;
str_low_sb = TString(Form("mass>%f && mass<%f",low_sb_min, low_sb_max));
cut_low_sb = str_low_sb;
cout << " Actual # of entries in low SB [" << low_sb_min << ", " << low_sb_max << "] = " << data->sumEntries(cut_low_sb) << endl;
str_high_sb = TString(Form("mass>%f && mass<%f",high_sb_min, high_sb_max));
cut_high_sb = str_high_sb;
cout << " Actual # of entries in high SB [" << high_sb_min << ", " << high_sb_max << "] = " << data->sumEntries(cut_high_sb) << endl;
str_sig = TString(Form("mass>%f && mass<%f",min_signal, max_signal));
cut_sig = str_sig;
cout << " Actual # of entries between [" << min_signal << ", " << max_signal << "] = " << data->sumEntries(cut_sig) << endl;
cout << " Extrapolation between [" << min_signal << ", " << max_signal << "] with RooFit, integral = " << integ2 << " +- " << dinteg2 << endl;
cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << endl << endl;
return xframe;
}
////////////////////////////////////////////////////////////////////////////////
// main //
////////////////////////////////////////////////////////////////////////////////
void fit(){
RooPlot* xframe_plus = plot_fit(60);
TCanvas* c1 = new TCanvas("positive slope fit","positive slope fit",900,600) ;
xframe_plus->Draw();
c1->Update();
RooPlot* xframe_minus = plot_fit(-0.8);
TCanvas* c2 = new TCanvas("negative slope fit","negative slope fit",900,600) ;
xframe_minus->Draw();
c2->Update();
}