Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
df101_h1Analysis.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_dataframe
3/// \notebook
4/// Show how to express ROOT's standard H1 analysis with RDataFrame.
5///
6/// \macro_code
7/// \macro_image
8///
9/// \date December 2016
10/// \authors Axel Naumann, Danilo Piparo (CERN)
11
12auto Select = [](ROOT::RDataFrame &dataFrame) {
13 using namespace ROOT;
14
15 auto ret = dataFrame.Filter("TMath::Abs(md0_d - 1.8646) < 0.04")
16 .Filter("ptds_d > 2.5")
17 .Filter("TMath::Abs(etads_d) < 1.5")
18 .Filter([](int ik, int ipi, RVecI& nhitrp) { return nhitrp[ik - 1] * nhitrp[ipi - 1] > 1; },
19 {"ik", "ipi", "nhitrp"})
20 .Filter([](int ik, RVecF& rstart, RVecF& rend) { return rend[ik - 1] - rstart[ik - 1] > 22; },
21 {"ik", "rstart", "rend"})
22 .Filter([](int ipi, RVecF& rstart, RVecF& rend) { return rend[ipi - 1] - rstart[ipi - 1] > 22; },
23 {"ipi", "rstart", "rend"})
24 .Filter([](int ik, RVecF& nlhk) { return nlhk[ik - 1] > 0.1; }, {"ik", "nlhk"})
25 .Filter([](int ipi, RVecF& nlhpi) { return nlhpi[ipi - 1] > 0.1; }, {"ipi", "nlhpi"})
26 .Filter([](int ipis, RVecF& nlhpi) { return nlhpi[ipis - 1] > 0.1; }, {"ipis", "nlhpi"})
27 .Filter("njets >= 1");
28
29 return ret;
30};
31
32const Double_t dxbin = (0.17 - 0.13) / 40; // Bin-width
33
35{
36 Double_t x = xx[0];
37 if (x <= 0.13957)
38 return 0;
39 Double_t xp3 = (x - par[3]) * (x - par[3]);
40 Double_t res =
41 dxbin * (par[0] * pow(x - 0.13957, par[1]) + par[2] / 2.5066 / par[4] * exp(-xp3 / 2 / par[4] / par[4]));
42 return res;
43}
44
46{
47 static const Double_t sigma = 0.0012;
48 Double_t x = xx[0];
49 if (x <= 0.13957)
50 return 0;
51 Double_t xp3 = (x - 0.1454) * (x - 0.1454);
52 Double_t res = dxbin * (par[0] * pow(x - 0.13957, 0.25) + par[1] / 2.5066 / sigma * exp(-xp3 / 2 / sigma / sigma));
53 return res;
54}
55
56void FitAndPlotHdmd(TH1 &hdmd)
57{
58 // create the canvas for the h1analysis fit
60 auto c1 = new TCanvas("c1", "h1analysis analysis", 10, 10, 800, 600);
61 hdmd.GetXaxis()->SetTitleOffset(1.4);
62
63 // fit histogram hdmd with function f5 using the loglikelihood option
64 auto f5 = new TF1("f5", fdm5, 0.139, 0.17, 5);
65 f5->SetParameters(1000000, .25, 2000, .1454, .001);
66 hdmd.Fit("f5", "lr");
67
68 hdmd.DrawClone();
69}
70
71void FitAndPlotH2(TH2 &h2)
72{
73 // create the canvas for tau d0
74 auto c2 = new TCanvas("c2", "tauD0", 100, 100, 800, 600);
75
76 c2->SetGrid();
77 c2->SetBottomMargin(0.15);
78
79 // Project slices of 2-d histogram h2 along X , then fit each slice
80 // with function f2 and make a histogram for each fit parameter
81 // Note that the generated histograms are added to the list of objects
82 // in the current directory.
83 auto f2 = new TF1("f2", fdm2, 0.139, 0.17, 2);
84 f2->SetParameters(10000, 10);
85 h2.FitSlicesX(f2, 0, -1, 1, "qln");
86
87 // See TH2::FitSlicesX documentation
88 auto h2_1 = (TH1D *)gDirectory->Get("h2_1");
89 h2_1->GetXaxis()->SetTitle("#tau [ps]");
90 h2_1->SetMarkerStyle(21);
91 h2_1->DrawClone();
92 c2->Update();
93
94 auto line = new TLine(0, 0, 0, c2->GetUymax());
95 line->Draw();
96}
97
98void df101_h1Analysis()
99{
100 TChain chain("h42");
101 chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarmb.root");
102 chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp1a.root");
103 chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp1b.root");
104 chain.Add("root://eospublic.cern.ch//eos/root-eos/h1/dstarp2.root");
105
107
108 ROOT::RDataFrame dataFrame(chain);
109 auto selected = Select(dataFrame);
110 // Note: The title syntax is "<Title>;<Label x axis>;<Label y axis>"
111 auto hdmdARP = selected.Histo1D({"hdmd", "Dm_d;m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]", 40, 0.13, 0.17}, "dm_d");
112 auto selectedAddedBranch = selected.Define("h2_y", "rpd0_t / 0.029979f * 1.8646f / ptd0_d");
113 auto h2ARP = selectedAddedBranch.Histo2D({"h2", "ptD0 vs Dm_d", 30, 0.135, 0.165, 30, -3, 6}, "dm_d", "h2_y");
114
115 FitAndPlotHdmd(*hdmdARP);
116 FitAndPlotH2(*h2ARP);
117}
double Double_t
Definition RtypesCore.h:59
#define gDirectory
Definition TDirectory.h:384
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
ROOT's RDataFrame offers a modern, high-level interface for analysis of data stored in TTree ,...
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
The Canvas class.
Definition TCanvas.h:23
A chain is a collection of files containing TTree objects.
Definition TChain.h:33
1-Dim function class
Definition TF1.h:214
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
TAxis * GetXaxis()
Definition TH1.h:322
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition TH1.cxx:3901
Service class for 2-D histogram classes.
Definition TH2.h:30
virtual void FitSlicesX(TF1 *f1=nullptr, Int_t firstybin=0, Int_t lastybin=-1, Int_t cut=0, Option_t *option="QNR", TObjArray *arr=nullptr)
Project slices along X in case of a 2-D histogram, then fit each slice with function f1 and make a hi...
Definition TH2.cxx:970
Use the TLine constructor to create a simple line.
Definition TLine.h:22
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:299
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition TStyle.cxx:1589
TLine * line
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition RVec.hxx:1809
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition RVec.hxx:1800
RVec< T > Filter(const RVec< T > &v, F &&f)
Create a new collection with the elements passing the filter expressed by the predicate.
Definition RVec.hxx:2145
Double_t fdm5(Double_t *xx, Double_t *par)
const Double_t sigma
const Double_t dxbin
Double_t fdm2(Double_t *xx, Double_t *par)
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
return c2
Definition legend2.C:14
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
void EnableImplicitMT(UInt_t numthreads=0)
Enable ROOT's implicit multi-threading for all objects and methods that provide an internal paralleli...
Definition TROOT.cxx:539