root-project / root

The official repository for ROOT: analyzing, storing and visualizing big data, scientifically
https://root.cern
Other
2.68k stars 1.27k forks source link

[RF] RooDataHist from a variable binning TH1D: different shape, different bin content #16028

Open santocch opened 3 months ago

santocch commented 3 months ago

Check duplicate issues.

Description

This problem has been already mentioned here [0] but no solution has been proposed. When you compare a RooDataHist obtained from a TH1 with fixed size binning everything looks ok but when you do the same from a variable size binning the comparison fails. See code below to reproduce the comparison.

[0] https://root-forum.cern.ch/t/roodatahist-and-th1d-why-different-shape-bin-content/56971

Reproducer

just run root mytest.cpp

using namespace RooFit;

void mytest() {

  const int nBin=100;
  double xbins[nBin+1];
  const float left_Val = 150;
  const float rightVal = 1150;  
  const float logxmin = TMath::Log10(left_Val);
  const float logxmax = TMath::Log10(rightVal);
  const double dxLog = (logxmax-logxmin)/nBin;
  for (int i=0;i<=nBin;i++) xbins[i] = TMath::Power(10., logxmin + i * dxLog);

  TF1 *f1 = new TF1("f1","expo(0)",left_Val,rightVal);
  f1->SetParameters(1e2,-1e-2);

  TH1D *h1 = new TH1D("h1","",nBin,xbins);
  h1->FillRandom("f1",10000000);
  RooRealVar mInv("mInv","m [GeV]",150,1150);
  RooDataHist hist_test1("hist_test1","hist_test1",mInv,Import(*h1,false));
  TCanvas *cExpo1 = new TCanvas();
  gPad->SetLogx();
  gPad->SetLogy();
  auto plot_test1 = mInv.frame();
  hist_test1.plotOn(plot_test1,DataError(RooAbsData::SumW2));
  plot_test1->Draw();
  h1->Draw("same");

  TH1D *h2 = new TH1D("h2","",nBin,left_Val,rightVal);
  h2->FillRandom("f1",10000000);
  RooDataHist hist_test2("hist_test2","hist_test2",mInv,Import(*h2,false));
  TCanvas *cExpo2 = new TCanvas();
  gPad->SetLogx();
  gPad->SetLogy();
  auto plot_test2 = mInv.frame();
  hist_test2.plotOn(plot_test2,DataError(RooAbsData::SumW2));
  plot_test2->Draw();
  h2->Draw("same");

}

ROOT version

ROOT 6.30/06

Installation method

homebrew

Operating system

OSX Sonoma 14.5

Additional context

No response

santocch commented 2 months ago

Hi few comments after digging a bit in the code...

See the attached modified version of the previous code to see the PDF begaviour mentioned above.

Best

Attilio

[0] https://github.com/root-project/root/blob/b8b0a8150325be271b45038dd81b751cb8a7a41d/roofit/roofitcore/src/RooHist.cxx#L434


int nEntries = 1000000;

void myTestRooDataHist() {

  const int nBin=100;
  double xbins[nBin+1];
  const float left_Val = 150;
  const float rightVal = 1150;  
  const float logxmin = TMath::Log10(left_Val);
  const float logxmax = TMath::Log10(rightVal);
  const double dxLog = (logxmax-logxmin)/nBin;
  for (int i=0;i<=nBin;i++) xbins[i] = TMath::Power(10., logxmin + i * dxLog);

  TF1 *f1 = new TF1("f1","expo(0)",left_Val,rightVal);
  f1->SetParameters(1e2,-1e-2);

  RooRealVar mInv("mInv","m [GeV]",150,1150);
  RooRealVar expo("expo", "expo", -1.0, -0.00001);
  RooExponential fitExpoFun("background", "background", mInv, expo);

  TH1D *h1 = new TH1D("h1","",nBin,left_Val,rightVal);
  h1->FillRandom("f1",nEntries);
  RooDataHist hist_test1("hist_test1","hist_test1",mInv,h1,1);
  TCanvas *cExpo1 = new TCanvas("cExpo1","",600,600);
  gPad->SetLogx();
  gPad->SetLogy();
  auto plot_test1 = mInv.frame();
  hist_test1.plotOn(plot_test1);
  plot_test1->Draw();
  h1->SetLineColor(2);
  h1->Draw("same");

  TH1D *h2 = new TH1D("h2","",nBin,xbins);
  h2->FillRandom("f1",nEntries);
  RooDataHist hist_test2("hist_test2","hist_test2",mInv,h2,1);
  TCanvas *cExpo2 = new TCanvas("cExpo2","",650,0,600,600);
  gPad->SetLogx();
  gPad->SetLogy();
  auto plot_test2 = mInv.frame();
  hist_test2.plotOn(plot_test2);
  plot_test2->Draw();
  h2->SetLineColor(2);
  h2->Draw("same");

  TCanvas *cExpo4 = new TCanvas("cExpo4","",1300,0,600,600);
  gPad->SetLogx();
  gPad->SetLogy();  
  RooHistPdf h1_PDF("h1_PDF","h1_PDF",mInv,hist_test1,0);
  auto plot_test1_pdf = mInv.frame();
  h1_PDF.plotOn(plot_test1_pdf);
  plot_test1_pdf->Draw();
  TH1D *h1clone = (TH1D*)h1->Clone("h1clone");
  h1clone->Scale(1./h1clone->Integral());
  h1clone->SetMarkerSize(0.5);
  h1clone->SetMarkerStyle(20);
  h1clone->SetMarkerColor(2);
  h1clone->SetLineColor(2);
  h1clone->Draw("same EP");
  auto result1 = fitExpoFun.fitTo(hist_test1,PrintLevel(0),RooFit::Save(true));
  result1->Print();

  TCanvas *cExpo5 = new TCanvas("cExpo5","",0,650,600,600);
  gPad->SetLogx();
  gPad->SetLogy();  
  RooHistPdf h2_PDF("h2_PDF","h2_PDF",mInv,hist_test2,0);
  auto plot_test2_pdf = mInv.frame();
  h2_PDF.plotOn(plot_test2_pdf);
  plot_test2_pdf->Draw();
  TH1D *h2clone = (TH1D*)h2->Clone("h2clone");
  h2clone->Scale(1./h2clone->Integral());
  h2clone->SetMarkerSize(0.5);
  h2clone->SetMarkerStyle(20);
  h2clone->SetMarkerColor(2);
  h2clone->SetLineColor(2);
  h2clone->Draw("same EP");
  auto result2 = fitExpoFun.fitTo(hist_test2,PrintLevel(0),RooFit::Save(true));
  result2->Print();

}