Basic functionality: importing data from ROOT TTrees and THx histograms.
Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:14 PM.
%%cpp -d
#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"
using namespace RooFit;
TH1 *makeTH1();
TTree *makeTTree();
Create ROOT TH1 filled with a Gaussian distribution
%%cpp -d
TH1 *makeTH1()
{
TH1D *hh = new TH1D("hh", "hh", 25, -10, 10);
for (int i = 0; i < 100; i++) {
hh->Fill(gRandom->Gaus(0, 3));
}
return hh;
}
Create ROOT TTree filled with a Gaussian distribution in x and a uniform distribution in y
%%cpp -d
TTree *makeTTree()
{
TTree *tree = new TTree("tree", "tree");
double *px = new double;
double *py = new double;
tree->Branch("x", px, "x/D");
tree->Branch("y", py, "y/D");
for (int i = 0; i < 100; i++) {
*px = gRandom->Gaus(0, 3);
*py = gRandom->Uniform() * 30 - 15;
tree->Fill();
}
return tree;
}
Create a ROOT TH1 histogram
TH1 *hh = makeTH1();
Declare observable x
RooRealVar x("x", "x", -10, 10);
Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
RooDataHist dh("dh", "dh", x, Import(*hh));
Make plot of binned dataset showing Poisson error bars (RooFit default)
RooPlot *frame = x.frame(Title("Imported TH1 with Poisson error bars"));
dh.plotOn(frame);
Fit a Gaussian pdf to the data
RooRealVar mean("mean", "mean", 0, -10, 10);
RooRealVar sigma("sigma", "sigma", 3, 0.1, 10);
RooGaussian gauss("gauss", "gauss", x, mean, sigma);
gauss.fitTo(dh, PrintLevel(-1));
gauss.plotOn(frame);
[#1] INFO:Fitting -- RooAbsPdf::fitTo(gauss_over_gauss_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_gauss_over_gauss_Int[x]_dh) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
If histogram has custom error (i.e. its contents is does not originate from a Poisson process but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead (same error bars as shown by ROOT)
RooPlot *frame2 = x.frame(Title("Imported TH1 with internal errors"));
dh.plotOn(frame2, DataError(RooAbsData::SumW2));
gauss.plotOn(frame2);
Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used in a maximum likelihood fit
A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly fitted with a chi^2 fit (see rf602_chi2fit.C)
TTree *tree = makeTTree();
Define 2nd observable y
RooRealVar y("y", "y", -10, 10);
Construct unbinned dataset importing tree branches x and y matching between branches and RooRealVars is done by name of the branch/RRV
Note that ONLY entries for which x,y have values within their allowed ranges as defined in RooRealVar x and y are imported. Since the y values in the import tree are in the range [-15,15] and RRV y defines a range [-10,10] this means that the RooDataSet below will have less entries than the TTree 'tree'
RooDataSet ds("ds", "ds", RooArgSet(x, y), Import(*tree));
[#1] INFO:DataHandling -- RooTreeDataStore::loadValues(ds) Skipping event #0 because y cannot accommodate the value 14.424 [#1] INFO:DataHandling -- RooTreeDataStore::loadValues(ds) Skipping event #3 because y cannot accommodate the value -12.0022 [#1] INFO:DataHandling -- RooTreeDataStore::loadValues(ds) Skipping event #5 because y cannot accommodate the value 13.8261 [#1] INFO:DataHandling -- RooTreeDataStore::loadValues(ds) Skipping event #6 because y cannot accommodate the value -14.9925 [#1] INFO:DataHandling -- RooTreeDataStore::loadValues(ds) Skipping ... [#0] WARNING:DataHandling -- RooTreeDataStore::loadValues(ds) Ignored 36 out-of-range events
{
// Write data to output stream
std::ofstream outstream("rf102_testData.txt");
// Optionally, adjust the stream here (e.g. std::setprecision)
ds.write(outstream);
outstream.close();
}
Read data from input stream. The variables of the dataset need to be supplied to the RooDataSet::read() function.
std::cout << "\n-----------------------\nReading data from ASCII\n";
RooDataSet *dataReadBack =
RooDataSet::read("rf102_testData.txt",
RooArgList(x, y), // variables to be read. If the file has more fields, these are ignored.
"D"); // Prints if a RooFit message stream listens for debug messages. Use Q for quiet.
dataReadBack->Print("V");
std::cout << "\nOriginal data, line 20:\n";
ds.get(20)->Print("V");
std::cout << "\nRead-back data, line 20:\n";
dataReadBack->get(20)->Print("V");
----------------------- Reading data from ASCII [#1] INFO:DataHandling -- RooDataSet::read: reading file rf102_testData.txt [#1] INFO:DataHandling -- RooDataSet::read: read 64 events (ignored 0 out of range events) DataStore dataset (rf102_testData.txt) Contains 64 entries Observables: 1) x = 0.0174204 L(-10 - 10) "x" 2) y = 9.46654 L(-10 - 10) "y" 3) blindState = Normal(idx = 0) "Blinding State" Original data, line 20: 1) RooRealVar:: x = -0.79919 2) RooRealVar:: y = 0.0106407 Read-back data, line 20: 1) RooRealVar:: x = -0.79919 2) RooRealVar:: y = 0.0106407 3) RooCategory:: blindState = Normal(idx = 0)
Print number of events in dataset
ds.Print();
RooDataSet::ds[x,y] = 64 entries
Print unbinned dataset with default frame binning (100 bins)
RooPlot *frame3 = y.frame(Title("Unbinned data shown in default frame binning"));
ds.plotOn(frame3);
Print unbinned dataset with custom binning choice (20 bins)
RooPlot *frame4 = y.frame(Title("Unbinned data shown with custom binning"));
ds.plotOn(frame4, Binning(20));
RooPlot *frame5 = y.frame(Title("Unbinned data read back from ASCII file"));
ds.plotOn(frame5, Binning(20));
dataReadBack->plotOn(frame5, Binning(20), MarkerColor(kRed), MarkerStyle(5));
Draw all frames on a canvas
TCanvas *c = new TCanvas("rf102_dataimport", "rf102_dataimport", 1000, 800);
c->Divide(3, 2);
c->cd(1);
gPad->SetLeftMargin(0.15);
frame->GetYaxis()->SetTitleOffset(1.4);
frame->Draw();
c->cd(2);
gPad->SetLeftMargin(0.15);
frame2->GetYaxis()->SetTitleOffset(1.4);
frame2->Draw();
c->cd(4);
gPad->SetLeftMargin(0.15);
frame3->GetYaxis()->SetTitleOffset(1.4);
frame3->Draw();
c->cd(5);
gPad->SetLeftMargin(0.15);
frame4->GetYaxis()->SetTitleOffset(1.4);
frame4->Draw();
c->cd(6);
gPad->SetLeftMargin(0.15);
frame4->GetYaxis()->SetTitleOffset(1.4);
frame5->Draw();
Draw all canvases
%jsroot on
gROOT->GetListOfCanvases()->Draw()