Logo ROOT   6.21/01
Reference Guide
h1analysisProxy.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_tree
3 /// Example of analysis class for the H1 data using code generated by MakeProxy.
4 ///
5 /// This file uses 4 large data sets from the H1 collaboration at DESY Hamburg.
6 /// One can access these data sets (277 MBytes) from the standard Root web site
7 /// at: `ftp:///root.cern.ch/root/h1analysis`
8 /// The Physics plots below generated by this example cannot be produced when
9 /// using smaller data sets.
10 ///
11 /// There are several ways to analyze data stored in a Root Tree
12 /// - Using TTree::Draw: This is very convenient and efficient for small tasks.
13 /// A TTree::Draw call produces one histogram at the time. The histogram
14 /// is automatically generated. The selection expression may be specified
15 /// in the command line.
16 ///
17 /// - Using the TTreeViewer: This is a graphical interface to TTree::Draw
18 /// with the same functionality.
19 ///
20 /// - Using the code generated by TTree::MakeClass: In this case, the user
21 /// creates an instance of the analysis class. They have the control over
22 /// the event loop and he can generate an unlimited number of histograms.
23 ///
24 /// - Using the code generated by TTree::MakeSelector. Like for the code
25 /// generated by TTree::MakeClass, the user can do complex analysis.
26 /// However, they cannot control the event loop. The event loop is controlled
27 /// by TTree::Process called by the user. This solution is illustrated
28 /// by the current code. The advantage of this method is that it can be run
29 /// in a parallel environment using PROOF (the Parallel Root Facility).
30 ///
31 /// A chain of 4 files (originally converted from PAW ntuples) is used
32 /// to illustrate the various ways to loop on Root data sets.
33 /// Each data set contains a Root Tree named "h42"
34 ///
35 /// h1analysProxy.C can be used either via TTree::Draw:
36 /// ~~~{.cpp}
37 /// h42->Draw("h1analysisProxy.C");
38 /// ~~~
39 /// or it can be used directly with TTree::MakeProxy, for example to generate a
40 /// shared library. TTree::MakeProxy will generate a TSelector skeleton that
41 /// include h1analysProxy.C:
42 /// ~~~{.cpp}
43 /// h42->MakeProxy("h1sel","h1analysisProxy.C");
44 /// ~~~
45 /// This produces one file: h1sel.h which does a #include "h1analysProxy.C"
46 /// The h1sel class is derived from the Root class TSelector and can then
47 /// be used as:
48 /// ~~~{.cpp}
49 /// h42->Process("h1sel.h+");
50 /// ~~~
51 ///
52 /// The following members functions are called by the TTree::Process functions.
53 /// - **h1analysProxy_Begin()**: Called every time a loop on the tree starts.
54 /// a convenient place to create your histograms.
55 ///
56 /// - **h1analysProxy_Notify()**: This function is called at the first entry of a new Tree
57 /// in a chain.
58 /// - **h1analysProxy_Process()**: called to analyze each entry.
59 ///
60 /// - **h1analysProxy_Terminate()**: called at the end of a loop on a TTree.
61 /// a convenient place to draw/fit your histograms.
62 ///
63 /// To use this file, try the following session
64 /// ~~~{.cpp}
65 /// Root > gROOT->Time(); ///will show RT & CPU time per command
66 /// ~~~
67 /// ### Case A: Create a TChain with the 4 H1 data files
68 /// The chain can be created by executed the short macro h1chain.C below:
69 /// ~~~{.cpp}
70 /// {
71 /// TChain chain("h42");
72 /// chain.Add("$H1/dstarmb.root"); /// 21330730 bytes 21920 events
73 /// chain.Add("$H1/dstarp1a.root"); /// 71464503 bytes 73243 events
74 /// chain.Add("$H1/dstarp1b.root"); /// 83827959 bytes 85597 events
75 /// chain.Add("$H1/dstarp2.root"); /// 100675234 bytes 103053 events
76 /// ///where $H1 is a system symbol pointing to the H1 data directory.
77 /// }
78 /// ~~~
79 ///
80 /// ### Case B: Loop on all events
81 /// ~~~{.cpp}
82 /// Root > chain.Draw("h1analysisProxy.C")
83 /// ~~~
84 ///
85 /// ### Case C: Same as B, but in addition fill the event list with selected entries.
86 /// The event list is saved to a file "elist.root" by the Terminate function.
87 /// To see the list of selected events, you can do elist->Print("all").
88 /// The selection function has selected 7525 events out of the 283813 events
89 /// in the chain of files. (2.65 per cent)
90 /// ~~~{.cpp}
91 /// Root > chain.Draw("h1analysisProxy.C","","fillList")
92 /// ~~~
93 /// ### Case D: Process only entries in the event list
94 /// The event list is read from the file in elist.root generated by step C
95 /// ~~~{.cpp}
96 /// Root > chain.Draw("h1analysisProxy.C","","useList")
97 /// ~~~
98 ///
99 /// The commands executed with the 3 different methods B,C and D
100 /// produce two canvases shown below:
101 /// begin_html <a href="gif/h1analysis_dstar.gif" >the Dstar plot</a> end_html
102 /// begin_html <a href="gif/h1analysis_tau.gif" >the Tau D0 plot</a> end_html
103 ///
104 /// \macro_code
105 ///
106 /// \author Philippe Canal from original h1analysis.C by Rene Brun
107 
108 TEntryList *elist;
109 Bool_t useList, fillList;
110 TH1F *hdmd;
111 TH2F *h2;
112 
113 
114 void h1analysisProxy_Begin(TTree *tree)
115 {
116 // function called before starting the event loop
117 // -it performs some cleanup
118 // -it creates histograms
119 // -it sets some initialisation for the event list
120 
121  //print the option specified in the Process function.
122  TString option = GetOption();
123  printf("Starting (begin) h1analysis with process option: %s\n",option.Data());
124 
125  //process cases with event list
126  fillList = kFALSE;
127  useList = kFALSE;
128  if (fChain) fChain->SetEntryList(0);
129  delete gDirectory->GetList()->FindObject("elist");
130 
131  // case when one creates/fills the event list
132  if (option.Contains("fillList")) {
133  fillList = kTRUE;
134  elist = new TEntryList("elist","H1 selection from Cut");
135  // Add to the input list for processing in PROOF, if needed
136  if (fInput) {
137  fInput->Add(new TNamed("fillList",""));
138  fInput->Add(elist);
139  }
140  } else elist = 0;
141 
142  // case when one uses the event list generated in a previous call
143  if (option.Contains("useList")) {
144  useList = kTRUE;
145  if (fInput) {
146  tree->SetEntryList(elist);
147  TFile f("elist.root");
148  elist = (TEntryList*)f.Get("elist");
149  if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
150  } else {
151  // Option "useList" not supported in PROOF directly
152  Warning("Begin", "option 'useList' not supported in PROOF - ignoring");
153  Warning("Begin", "the entry list must be set on the chain *before* calling Process");
154  }
155  }
156 }
157 
158 
159 void h1analysisProxy_SlaveBegin(TTree *tree)
160 {
161 // function called before starting the event loop
162 // -it performs some cleanup
163 // -it creates histograms
164 // -it sets some initialisation for the entry list
165 
166  //initialize the Tree branch addresses
167  Init(tree);
168 
169  //print the option specified in the Process function.
170  TString option = GetOption();
171  printf("Starting (slave) h1analysis with process option: %s\n",option.Data());
172 
173  //create histograms
174  hdmd = new TH1F("hdmd","dm_d",40,0.13,0.17);
175  h2 = new TH2F("h2","ptD0 vs dm_d",30,0.135,0.165,30,-3,6);
176 
177  fOutput->Add(hdmd);
178  fOutput->Add(h2);
179 
180  //process cases with entry list
181  fillList = kFALSE;
182  useList = kFALSE;
183 
184  // case when one creates/fills the entry list
185  if (option.Contains("fillList")) {
186  fillList = kTRUE;
187  // Get the list
188  if (fInput) {
189  if ((elist = (TEntryList *) fInput->FindObject("elist")))
190  // Need to clone to avoid problems when destroying the selector
191  elist = (TEntryList *) elist->Clone();
192  }
193  if (elist)
194  fOutput->Add(elist);
195  else
196  fillList = kFALSE;
197  } else elist = 0;
198 
199  // case when one uses the entry list generated in a previous call
200  if (option.Contains("useList")) {
201  useList = kTRUE;
202  TFile f("elist.root");
203  elist = (TEntryList*)f.Get("elist");
204  if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
205  if (tree) tree->SetEntryList(elist);
206  else {
207  // Option "useList" not supported in PROOF directly
208  Warning("Begin", "option 'useList' not supported in PROOF - ignoring");
209  Warning("Begin", "the entry list must be set on the chain *before* calling Process");
210  }
211  }
212 
213 }
214 
215 Double_t h1analysisProxy() {
216  return 0;
217 }
218 
219 
220 Bool_t h1analysisProxy_Process(Long64_t entry)
221 {
222 // entry is the entry number in the current Tree
223 // Selection function to select D* and D0.
224 
225  //in case one entry list is given in input, the selection has already been done.
226  if (!useList) {
227 
228  float f1 = md0_d;
229  float f2 = md0_d-1.8646;
230  bool test = TMath::Abs(md0_d-1.8646) >= 0.04;
231  if (gDebug>0) fprintf(stderr,"entry #%lld f1=%f f2=%f test=%d\n",
232  fChain->GetReadEntry(),f1,f2,test);
233 
234  if (TMath::Abs(md0_d-1.8646) >= 0.04) return kFALSE;
235  if (ptds_d <= 2.5) return kFALSE;
236  if (TMath::Abs(etads_d) >= 1.5) return kFALSE;
237 
238  int cik = ik-1; //original ik used f77 convention starting at 1
239  int cipi = ipi-1; //original ipi used f77 convention starting at 1
240 
241  f1 = nhitrp[cik];
242  f2 = nhitrp[cipi];
243  test = nhitrp[cik]*nhitrp[cipi] <= 1;
244  if (gDebug>0) fprintf(stderr,"entry #%lld f1=%f f2=%f test=%d\n",
245  fChain->GetReadEntry(),f1,f2,test);
246 
247  if (nhitrp[cik]*nhitrp[cipi] <= 1) return kFALSE;
248  if (rend[cik] -rstart[cik] <= 22) return kFALSE;
249  if (rend[cipi]-rstart[cipi] <= 22) return kFALSE;
250  if (nlhk[cik] <= 0.1) return kFALSE;
251  if (nlhpi[cipi] <= 0.1) return kFALSE;
252  // fix because read-only
253  if (nlhpi[ipis-1] <= 0.1) return kFALSE;
254  if (njets < 1) return kFALSE;
255 
256  }
257  // if option fillList, fill the event list
258  if (fillList) elist->Enter(entry);
259 
260  //fill some histograms
261  hdmd->Fill(dm_d);
262  h2->Fill(dm_d,rpd0_t/0.029979*1.8646/ptd0_d);
263 
264  return kTRUE;
265 }
266 
267 
268 
269 void h1analysisProxy_SlaveTerminate()
270 {
271  // nothing to be done
272  printf("Terminate (slave) h1analysis\n");
273 }
274 
275 
276 void h1analysisProxy_Terminate()
277 {
278  printf("Terminate (final) h1analysis\n");
279 
280  // function called at the end of the event loop
281 
282  hdmd = dynamic_cast<TH1F*>(fOutput->FindObject("hdmd"));
283  h2 = dynamic_cast<TH2F*>(fOutput->FindObject("h2"));
284 
285  if (hdmd == 0 || h2 == 0) {
286  Error("Terminate", "hdmd = %p , h2 = %p", hdmd, h2);
287  return;
288  }
289 
290  //create the canvas for the h1analysis fit
291  gStyle->SetOptFit();
292  TCanvas *c1 = new TCanvas("c1","h1analysis analysis",10,10,800,600);
293  c1->SetBottomMargin(0.15);
294  hdmd->GetXaxis()->SetTitle("m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]");
295  hdmd->GetXaxis()->SetTitleOffset(1.4);
296 
297  //fit histogram hdmd with function f5 using the log-likelihood option
298  TF1 *f5 = new TF1("f5",fdm5,0.139,0.17,5);
299  f5->SetParameters(1000000, .25, 2000, .1454, .001);
300  hdmd->Fit("f5","lr");
301 
302  //create the canvas for tau d0
303  gStyle->SetOptFit(0);
304  gStyle->SetOptStat(1100);
305  TCanvas *c2 = new TCanvas("c2","tauD0",100,100,800,600);
306  c2->SetGrid();
307  c2->SetBottomMargin(0.15);
308 
309  // Project slices of 2-d histogram h2 along X , then fit each slice
310  // with function f2 and make a histogram for each fit parameter
311  // Note that the generated histograms are added to the list of objects
312  // in the current directory.
313  TF1 *f2 = new TF1("f2",fdm2,0.139,0.17,2);
314  f2->SetParameters(10000, 10);
315  h2->FitSlicesX(f2,0,-1,1,"qln");
316  TH1D *h2_1 = (TH1D*)gDirectory->Get("h2_1");
317  h2_1->GetXaxis()->SetTitle("#tau[ps]");
318  h2_1->SetMarkerStyle(21);
319  h2_1->Draw();
320  c2->Update();
321  TLine *line = new TLine(0,0,0,c2->GetUymax());
322  line->Draw();
323 
324  // Have the number of entries on the first histogram (to cross check when running
325  // with entry lists)
326  TPaveStats *psdmd = (TPaveStats *)hdmd->GetListOfFunctions()->FindObject("stats");
327  psdmd->SetOptStat(1110);
328  c1->Modified();
329 
330  //save the entry list to a Root file if one was produced
331  if (fillList) {
332  elist = dynamic_cast<TEntryList*>(fOutput->FindObject("elist"));
333  if (elist) {
334  TFile efile("elist.root","recreate");
335  elist->Write();
336  } else {
337  Error("Terminate", "entry list requested but not found in output");
338  }
339  }
340 }
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:785
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
long long Long64_t
Definition: RtypesCore.h:69
TLine * line
return c1
Definition: legend1.C:41
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
Double_t fdm2(Double_t *xx, Double_t *par)
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:48
Definition: test.py:1
Basic string class.
Definition: TString.h:131
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
#define f(i)
Definition: RSha256.hxx:104
bool Bool_t
Definition: RtypesCore.h:59
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
The histogram statistics painter class.
Definition: TPaveStats.h:18
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
virtual void SetDirectory(TDirectory *dir)
Add reference to directory dir. dir can be 0.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
void Init(TClassEdit::TInterpreterLookupHelper *helper)
Definition: TClassEdit.cxx:155
void Error(const char *location, const char *msgfmt,...)
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition: TObject.cxx:321
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
virtual void FitSlicesX(TF1 *f1=0, Int_t firstybin=0, Int_t lastybin=-1, Int_t cut=0, Option_t *option="QNR", TObjArray *arr=0)
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:858
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
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:1404
A simple line.
Definition: TLine.h:23
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
void Warning(const char *location, const char *msgfmt,...)
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
const Bool_t kFALSE
Definition: RtypesCore.h:88
The Canvas class.
Definition: TCanvas.h:31
return c2
Definition: legend2.C:14
double Double_t
Definition: RtypesCore.h:55
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
virtual Bool_t Enter(Long64_t entry, TTree *tree=0)
Add entry #entry to the list.
Definition: TEntryList.cxx:558
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:74
Double_t fdm5(Double_t *xx, Double_t *par)
1-Dim function class
Definition: TF1.h:211
TF1 * f1
Definition: legend1.C:11
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:1452
R__EXTERN Int_t gDebug
Definition: Rtypes.h:91
Definition: tree.py:1
A TTree represents a columnar dataset.
Definition: TTree.h:72
#define gDirectory
Definition: TDirectory.h:223
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
A List of entry numbers in a TTree or TChain.
Definition: TEntryList.h:25
const Bool_t kTRUE
Definition: RtypesCore.h:87
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
void SetOptStat(Int_t stat=1)
Set the stat option.
Definition: TPaveStats.cxx:302
const char * Data() const
Definition: TString.h:364