Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
h1analysis.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_tree
3/// \notebook -header -nodraw
4/// Example of analysis class for the H1 data.
5///
6/// This file uses 4 large data sets from the H1 collaboration at DESY Hamburg.
7/// One can access these data sets (277 MBytes) from the standard Root web site
8/// at: [https://root.cern/files/h1/](https://root.cern/files/h1/)
9/// The Physics plots below generated by this example cannot be produced when
10/// using smaller data sets.
11///
12/// There are several ways to analyze data stored in a Root Tree
13/// - Using TTree::Draw: This is very convenient and efficient for small tasks.
14/// A TTree::Draw call produces one histogram at the time. The histogram
15/// is automatically generated. The selection expression may be specified
16/// in the command line.
17///
18/// - Using the TTreeViewer: This is a graphical interface to TTree::Draw
19/// with the same functionality.
20///
21/// - Using the code generated by TTree::MakeClass: In this case, the user
22/// creates an instance of the analysis class. They have the control over
23/// the event loop and he can generate an unlimited number of histograms.
24///
25/// - Using the code generated by TTree::MakeSelector. Like for the code
26/// generated by TTree::MakeClass, the user can do complex analysis.
27/// However, they cannot control the event loop. The event loop is controlled
28/// by TTree::Process called by the user. This solution is illustrated
29/// by the current code. The advantage of this method is that it can be run
30/// in a parallel environment using PROOF (the Parallel Root Facility).
31///
32/// A chain of 4 files (originally converted from PAW ntuples) is used
33/// to illustrate the various ways to loop on Root data sets.
34/// Each data set contains a Root Tree named "h42"
35/// The class definition in h1analysis.h has been generated automatically
36/// by the Root utility TTree::MakeSelector using one of the files with the
37/// following statement:
38///
39/// ~~~{.cpp}
40/// h42->MakeSelector("h1analysis");
41/// ~~~
42///
43/// This produces two files: h1analysis.h and h1analysis.C (skeleton of this file)
44/// The h1analysis class is derived from the Root class TSelector.
45///
46/// The following members functions are called by the TTree::Process functions.
47/// - **Begin()**: Called every time a loop on the tree starts.
48/// A convenient place to create your histograms.
49/// - **Notify()**: This function is called at the first entry of a new Tree
50/// in a chain.
51/// - **Process()**: Called to analyze each entry.
52///
53/// - **Terminate()**: Called at the end of a loop on a TTree.
54/// A convenient place to draw/fit your histograms.
55///
56/// To use this file, try the following sessions
57///
58/// ~~~{.cpp}
59/// Root > gROOT->Time(); /// will show RT & CPU time per command
60/// ~~~
61///
62/// ### Case A: Create a TChain with the 4 H1 data files
63///
64/// The chain can be created by executed the short macro h1chain.C below:
65///
66/// ~~~{.cpp}
67/// {
68/// TChain chain("h42");
69/// chain.Add("$H1/dstarmb.root"); /// 21330730 bytes 21920 events
70/// chain.Add("$H1/dstarp1a.root"); /// 71464503 bytes 73243 events
71/// chain.Add("$H1/dstarp1b.root"); /// 83827959 bytes 85597 events
72/// chain.Add("$H1/dstarp2.root"); /// 100675234 bytes 103053 events
73/// /// where $H1 is a system symbol pointing to the H1 data directory.
74/// }
75/// ~~~
76///
77/// ### Case B: Loop on all events
78///
79/// ~~~{.cpp}
80/// Root > chain.Process("h1analysis.C")
81/// ~~~
82///
83/// ### Case C: Same as B, but in addition fill the entry list with selected entries.
84///
85/// The entry list is saved to a file "elist.root" by the Terminate function.
86/// To see the list of selected events, you can do `elist->Print("all")`.
87/// The selection function has selected 7525 events out of the 283813 events
88/// in the chain of files. (2.65 per cent)
89///
90/// ~~~{.cpp}
91/// Root > chain.Process("h1analysis.C","fillList")
92/// ~~~
93///
94/// ### Case D: Process only entries in the entry list
95///
96/// The entry list is read from the file in elist.root generated by step C
97///
98/// ~~~{.cpp}
99/// Root > chain.Process("h1analysis.C","useList")
100/// ~~~
101///
102/// ### Case E: The above steps have been executed via the interpreter.
103/// You can repeat the steps B, C and D using the script compiler
104/// by replacing "h1analysis.C" by "h1analysis.C+" or "h1analysis.C++"
105/// in a new session (see F).
106///
107/// ### Case F: Create the chain as in A, then execute
108///
109/// ~~~{.cpp}
110/// Root > chain.Process("h1analysis.C+","useList")
111/// ~~~
112///
113/// The same analysis can be run on PROOF. For a quick try start a PROOF-Lite
114/// session
115///
116/// ~~~{.cpp}
117/// Root > TProof *p = TProof::Open("")
118/// ~~~
119///
120/// create (if not already done) the chain by executing the 'h1chain.C' macro
121/// mentioned above, and then tell ROOT to use PROOF to process the chain:
122///
123/// ~~~{.cpp}
124/// Root > chain.SetProof()
125/// ~~~
126///
127/// You can then repeat step B above. Step C can also be executed in PROOF. However,
128/// step D cannot be executed in PROOF as in the local session (i.e. just passing
129/// option 'useList'): to use the entry list you have to
130///
131/// ### Case G: Load first in the session the list form the file
132///
133/// ~~~{.cpp}
134/// Root > TFile f("elist.root")
135/// Root > TEntryList *elist = (TEntryList *) f.Get("elist")
136/// ~~~
137///
138/// set it on the chain:
139///
140/// ~~~{.cpp}
141/// Root > chain.SetEntryList(elist)
142/// ~~~
143///
144/// call Process as in step B. Of course this works also for local processing.
145///
146/// \macro_code
147///
148/// \author Rene Brun
149
150#include "h1analysis.h"
151#include "TH2.h"
152#include "TF1.h"
153#include "TStyle.h"
154#include "TBranch.h"
155#include "TCanvas.h"
156#include "TPaveStats.h"
157#include "TLine.h"
158#include "TMath.h"
159
160const Double_t dxbin = (0.17-0.13)/40; // Bin-width
161const Double_t sigma = 0.0012;
162
163
165{
166 Double_t x = xx[0];
167 if (x <= 0.13957) return 0;
168 Double_t xp3 = (x-par[3])*(x-par[3]);
169 Double_t res = dxbin*(par[0]*TMath::Power(x-0.13957, par[1])
170 + par[2] / 2.5066/par[4]*TMath::Exp(-xp3/2/par[4]/par[4]));
171 return res;
172}
173
174
176{
177 Double_t x = xx[0];
178 if (x <= 0.13957) return 0;
179 Double_t xp3 = (x-0.1454)*(x-0.1454);
180 Double_t res = dxbin*(par[0]*TMath::Power(x-0.13957, 0.25)
181 + par[1] / 2.5066/sigma*TMath::Exp(-xp3/2/sigma/sigma));
182 return res;
183}
184
185
186void h1analysis::Begin(TTree * /*tree*/)
187{
188// function called before starting the event loop
189// -it performs some cleanup
190// -it creates histograms
191// -it sets some initialisation for the entry list
192
193 // This is needed when re-processing the object
194 Reset();
195
196 //print the option specified in the Process function.
198 Info("Begin", "starting h1analysis with process option: %s", option.Data());
199
200 //process cases with entry list
201 if (fChain) fChain->SetEntryList(nullptr);
202 delete gDirectory->GetList()->FindObject("elist");
203
204 // case when one creates/fills the entry list
205 if (option.Contains("fillList")) {
206 fillList = kTRUE;
207 elist = new TEntryList("elist", "H1 selection from Cut");
208 // Add to the input list for processing in PROOF, if needed
209 if (fInput) {
210 fInput->Add(new TNamed("fillList",""));
211 // We send a clone to avoid double deletes when importing the result
212 fInput->Add(elist);
213 // This is needed to avoid warnings from output-to-members mapping
214 elist = nullptr;
215 }
216 Info("Begin", "creating an entry-list");
217 }
218 // case when one uses the entry list generated in a previous call
219 if (option.Contains("useList")) {
220 useList = kTRUE;
221 if (fInput) {
222 // In PROOF option "useList" is processed in SlaveBegin and we do not need
223 // to do anything here
224 } else {
225 TFile f("elist.root");
226 elist = (TEntryList*)f.Get("elist");
227 if (elist) elist->SetDirectory(nullptr); //otherwise the file destructor will delete elist
228 }
229 }
230}
231
232
234{
235// function called before starting the event loop
236// -it performs some cleanup
237// -it creates histograms
238// -it sets some initialisation for the entry list
239
240 //initialize the Tree branch addresses
241 Init(tree);
242
243 //print the option specified in the Process function.
245 Info("SlaveBegin",
246 "starting h1analysis with process option: %s (tree: %p)", option.Data(), tree);
247
248 //create histograms
249 hdmd = new TH1F("hdmd","dm_d",40,0.13,0.17);
250 h2 = new TH2F("h2","ptD0 vs dm_d",30,0.135,0.165,30,-3,6);
251
252 fOutput->Add(hdmd);
253 fOutput->Add(h2);
254
255 // Entry list stuff (re-parse option because on PROOF only SlaveBegin is called)
256 if (option.Contains("fillList")) {
257 fillList = kTRUE;
258 // Get the list
259 if (fInput) {
260 if ((elist = (TEntryList *) fInput->FindObject("elist")))
261 // Need to clone to avoid problems when destroying the selector
262 elist = (TEntryList *) elist->Clone();
263 if (elist)
264 fOutput->Add(elist);
265 else
267 }
268 }
269 if (fillList) Info("SlaveBegin", "creating an entry-list");
270 if (option.Contains("useList")) useList = kTRUE;
271}
272
273
275{
276// entry is the entry number in the current Tree
277// Selection function to select D* and D0.
278
279 fProcessed++;
280 //in case one entry list is given in input, the selection has already been done.
281 if (!useList) {
282 // Read only the necessary branches to select entries.
283 // return as soon as a bad entry is detected
284 // to read complete event, call fChain->GetTree()->GetEntry(entry)
285 b_md0_d->GetEntry(entry); if (TMath::Abs(md0_d-1.8646) >= 0.04) return kFALSE;
286 b_ptds_d->GetEntry(entry); if (ptds_d <= 2.5) return kFALSE;
287 b_etads_d->GetEntry(entry); if (TMath::Abs(etads_d) >= 1.5) return kFALSE;
288 b_ik->GetEntry(entry); ik--; //original ik used f77 convention starting at 1
289 b_ipi->GetEntry(entry); ipi--;
290 b_ntracks->GetEntry(entry);
291 b_nhitrp->GetEntry(entry);
292 if (nhitrp[ik]*nhitrp[ipi] <= 1) return kFALSE;
293 b_rend->GetEntry(entry);
294 b_rstart->GetEntry(entry);
295 if (rend[ik] -rstart[ik] <= 22) return kFALSE;
296 if (rend[ipi]-rstart[ipi] <= 22) return kFALSE;
297 b_nlhk->GetEntry(entry); if (nlhk[ik] <= 0.1) return kFALSE;
298 b_nlhpi->GetEntry(entry); if (nlhpi[ipi] <= 0.1) return kFALSE;
299 b_ipis->GetEntry(entry); ipis--; if (nlhpi[ipis] <= 0.1) return kFALSE;
300 b_njets->GetEntry(entry); if (njets < 1) return kFALSE;
301 }
302 // if option fillList, fill the entry list
303 if (fillList) elist->Enter(entry);
304
305 // to read complete event, call fChain->GetTree()->GetEntry(entry)
306 // read branches not processed in ProcessCut
307 b_dm_d->GetEntry(entry); //read branch holding dm_d
308 b_rpd0_t->GetEntry(entry); //read branch holding rpd0_t
309 b_ptd0_d->GetEntry(entry); //read branch holding ptd0_d
310
311 //fill some histograms
312 hdmd->Fill(dm_d);
313 h2->Fill(dm_d,rpd0_t/0.029979*1.8646/ptd0_d);
314
315 // Count the number of selected events
316 fStatus++;
317
318 return kTRUE;
319}
320
321
322
324{
325 // nothing to be done
326}
327
328
330{
331// function called at the end of the event loop
332
333 hdmd = dynamic_cast<TH1F*>(fOutput->FindObject("hdmd"));
334 h2 = dynamic_cast<TH2F*>(fOutput->FindObject("h2"));
335
336 if (hdmd == nullptr || h2 == nullptr) {
337 Error("Terminate", "hdmd = %p , h2 = %p", hdmd, h2);
338 return;
339 }
340
341 //create the canvas for the h1analysis fit
342 gStyle->SetOptFit();
343 TCanvas *c1 = new TCanvas("c1","h1analysis analysis",10,10,800,600);
344 c1->SetBottomMargin(0.15);
345 hdmd->GetXaxis()->SetTitle("m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]");
347
348 //fit histogram hdmd with function f5 using the log-likelihood option
349 if (gROOT->GetListOfFunctions()->FindObject("f5"))
350 delete gROOT->GetFunction("f5");
351 TF1 *f5 = new TF1("f5",fdm5,0.139,0.17,5);
352 f5->SetParameters(1000000, .25, 2000, .1454, .001);
353 hdmd->Fit("f5","lr");
354
355 //create the canvas for tau d0
356 gStyle->SetOptFit(0);
357 gStyle->SetOptStat(1100);
358 TCanvas *c2 = new TCanvas("c2","tauD0",100,100,800,600);
359 c2->SetGrid();
360 c2->SetBottomMargin(0.15);
361
362 // Project slices of 2-d histogram h2 along X , then fit each slice
363 // with function f2 and make a histogram for each fit parameter
364 // Note that the generated histograms are added to the list of objects
365 // in the current directory.
366 if (gROOT->GetListOfFunctions()->FindObject("f2"))
367 delete gROOT->GetFunction("f2");
368 TF1 *f2 = new TF1("f2",fdm2,0.139,0.17,2);
369 f2->SetParameters(10000, 10);
370 h2->FitSlicesX(f2,0,-1,1,"qln");
371 TH1D *h2_1 = (TH1D*)gDirectory->Get("h2_1");
372 h2_1->GetXaxis()->SetTitle("#tau[ps]");
373 h2_1->SetMarkerStyle(21);
374 h2_1->Draw();
375 c2->Update();
376 TLine *line = new TLine(0,0,0,c2->GetUymax());
377 line->Draw();
378
379 // Have the number of entries on the first histogram (to cross check when running
380 // with entry lists)
381 TPaveStats *psdmd = (TPaveStats *)hdmd->GetListOfFunctions()->FindObject("stats");
382 psdmd->SetOptStat(1110);
383 c1->Modified();
384
385 //save the entry list to a Root file if one was produced
386 if (fillList) {
387 if (!elist)
388 elist = dynamic_cast<TEntryList*>(fOutput->FindObject("elist"));
389 if (elist) {
390 Printf("Entry list 'elist' created:");
391 elist->Print();
392 TFile efile("elist.root","recreate");
393 elist->Write();
394 } else {
395 Error("Terminate", "entry list requested but not found in output");
396 }
397 }
398 // Notify the amount of processed events
399 if (!fInput) Info("Terminate", "processed %lld events", fProcessed);
400}
#define f(i)
Definition RSha256.hxx:104
bool Bool_t
Definition RtypesCore.h:63
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
double Double_t
Definition RtypesCore.h:59
long long Long64_t
Definition RtypesCore.h:69
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
#define gDirectory
Definition TDirectory.h:384
Option_t Option_t option
#define gROOT
Definition TROOT.h:406
void Printf(const char *fmt,...)
Formats a string in a circular formatting buffer and prints the string.
Definition TString.cxx:2503
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all leaves of entry and return total number of bytes read.
Definition TBranch.cxx:1706
The Canvas class.
Definition TCanvas.h:23
A List of entry numbers in a TTree or TChain.
Definition TEntryList.h:26
virtual bool Enter(Long64_t entry, TTree *tree=nullptr)
Add entry #entry to the list.
virtual void SetDirectory(TDirectory *dir)
Add reference to directory dir. dir can be 0.
void Print(const Option_t *option="") const override
Print this list.
1-Dim function class
Definition TF1.h:233
virtual void SetParameters(const Double_t *params)
Definition TF1.h:677
A ROOT file is an on-disk file, usually with extension .root, that stores objects in a file-system-li...
Definition TFile.h:53
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:671
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:623
TAxis * GetXaxis()
Definition TH1.h:325
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:3906
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3346
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3068
TList * GetListOfFunctions() const
Definition TH1.h:245
2-D histogram with a float per channel (see TH1 documentation)
Definition TH2.h:308
Int_t Fill(Double_t) override
Invalid Fill method.
Definition TH2.cxx:393
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:1013
TObject * FindObject(const char *name) const override
Find object using its name.
Use the TLine constructor to create a simple line.
Definition TLine.h:22
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:576
void Add(TObject *obj) override
Definition TList.h:81
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
TObject * Clone(const char *newname="") const override
Make a clone of an object using the Streamer facility.
Definition TNamed.cxx:74
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:898
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1005
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:292
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:979
The histogram statistics painter class.
Definition TPaveStats.h:18
void SetOptStat(Int_t stat=1)
Set the stat option.
TList * fInput
List of objects available during processing.
Definition TSelector.h:41
TSelectorList * fOutput
! List of objects created during processing
Definition TSelector.h:42
Long64_t fStatus
Selector status.
Definition TSelector.h:37
const char * GetOption() const override
Definition TSelector.h:57
Basic string class.
Definition TString.h:139
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1640
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:1593
A TTree represents a columnar dataset.
Definition TTree.h:79
virtual void SetEntryList(TEntryList *list, Option_t *opt="")
Set an EntryList.
Definition TTree.cxx:9036
TBranch * b_md0_d
Definition h1analysis.h:265
TBranch * b_dm_d
Definition h1analysis.h:254
Float_t rpd0_t
Definition h1analysis.h:131
void SlaveTerminate() override
TBranch * b_nlhpi
Definition h1analysis.h:320
TBranch * b_ipis
Definition h1analysis.h:261
TBranch * b_ipi
Definition h1analysis.h:260
Float_t nlhpi[200]
Definition h1analysis.h:166
Bool_t useList
Definition h1analysis.h:30
Int_t nhitrp[200]
Definition h1analysis.h:157
TBranch * b_ik
Definition h1analysis.h:259
TBranch * b_rstart
Definition h1analysis.h:315
TBranch * b_rend
Definition h1analysis.h:316
TBranch * b_rpd0_t
Definition h1analysis.h:285
Bool_t fillList
Definition h1analysis.h:31
Float_t md0_d
Definition h1analysis.h:111
void Init(TTree *tree) override
Definition h1analysis.h:390
TBranch * b_etads_d
Definition h1analysis.h:253
TTree * fChain
Definition h1analysis.h:35
TBranch * b_nhitrp
Definition h1analysis.h:311
Long64_t fProcessed
Definition h1analysis.h:33
TH1F * hdmd
Definition h1analysis.h:27
Float_t dm_d
Definition h1analysis.h:100
void Reset()
Definition h1analysis.h:376
Bool_t Process(Long64_t entry) override
The Process() function is called for each entry in the tree (or possibly keyed object in the case of ...
Float_t rend[200]
Definition h1analysis.h:162
TBranch * b_nlhk
Definition h1analysis.h:319
TBranch * b_ptds_d
Definition h1analysis.h:252
TBranch * b_ntracks
Definition h1analysis.h:303
Float_t rstart[200]
Definition h1analysis.h:161
TBranch * b_njets
Definition h1analysis.h:328
void Begin(TTree *tree) override
Float_t ptd0_d
Definition h1analysis.h:109
Int_t njets
Definition h1analysis.h:174
Float_t nlhk[200]
Definition h1analysis.h:165
TEntryList * elist
Definition h1analysis.h:32
void SlaveBegin(TTree *tree) override
Float_t etads_d
Definition h1analysis.h:99
TH2F * h2
Definition h1analysis.h:28
TBranch * b_ptd0_d
Definition h1analysis.h:263
Float_t ptds_d
Definition h1analysis.h:98
void Terminate() override
TLine * line
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
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Definition TMath.h:713
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:725
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123