Logo ROOT  
Reference Guide
mp_H1_lambdas.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_multicore
3 /// \notebook -nodraw
4 /// Lambdas used to check and fit the result of the H1 analysis.
5 /// Used by mp104_processH1.C, mp105_processEntryList.C and roottest/root/multicore/tProcessExecutorH1Test.C
6 ///
7 /// \macro_code
8 ///
9 /// \author Gerardo Ganis
10 
11 // This function is used to check the result of the H1 analysis
12 auto checkH1 = [](TList *out) {
13 
14  // Make sure the output list is there
15  if (!out) {
16  std::cout << "checkH1 >>> Test failure: output list not found\n";
17  return -1;
18  }
19 
20  // Check the 'hdmd' histo
21  auto hdmd = dynamic_cast<TH1F *>(out->FindObject("hdmd"));
22  if (!hdmd) {
23  std::cout << "checkH1 >>> Test failure: 'hdmd' histo not found\n";
24  return -1;
25  }
26  if ((Int_t)(hdmd->GetEntries()) != 7525) {
27  std::cout << "checkH1 >>> Test failure: 'hdmd' histo: wrong number"
28  " of entries ("
29  << (Int_t)(hdmd->GetEntries()) << ": expected 7525) \n";
30  return -1;
31  }
32  if (TMath::Abs((hdmd->GetMean() - 0.15512023) / 0.15512023) > 0.001) {
33  std::cout << "checkH1 >>> Test failure: 'hdmd' histo: wrong mean (" << hdmd->GetMean()
34  << ": expected 0.15512023) \n";
35  return -1;
36  }
37 
38  auto h2 = dynamic_cast<TH2F *>(out->FindObject("h2"));
39  if (!h2) {
40  std::cout << "checkH1 >>> Test failure: 'h2' histo not found\n";
41  return -1;
42  }
43  if ((Int_t)(h2->GetEntries()) != 7525) {
44  std::cout << "checkH1 >>> Test failure: 'h2' histo: wrong number"
45  " of entries ("
46  << (Int_t)(h2->GetEntries()) << ": expected 7525) \n";
47  return -1;
48  }
49  if (TMath::Abs((h2->GetMean() - 0.15245688) / 0.15245688) > 0.001) {
50  std::cout << "checkH1 >>> Test failure: 'h2' histo: wrong mean (" << h2->GetMean() << ": expected 0.15245688) \n";
51  return -1;
52  }
53 
54  // Done
55  return 0;
56 };
57 
58 // This function is used to fit the result of the analysis with graphics
59 auto doFit = [](TList *out, const char *lfn = 0) -> Int_t {
60 
61  RedirectHandle_t redH;
62  if (lfn)
63  gSystem->RedirectOutput(lfn, "a", &redH);
64 
65  auto hdmd = dynamic_cast<TH1F *>(out->FindObject("hdmd"));
66  auto h2 = dynamic_cast<TH2F *>(out->FindObject("h2"));
67 
68  // function called at the end of the event loop
69  if (hdmd == 0 || h2 == 0) {
70  std::cout << "doFit: hdmd = " << hdmd << " , h2 = " << h2 << "\n";
71  return -1;
72  if (lfn)
73  gSystem->RedirectOutput(0, 0, &redH);
74  }
75 
76  // create the canvas for the h1analysis fit
77  gStyle->SetOptFit();
78  TCanvas *c1 = new TCanvas("c1", "h1analysis analysis", 10, 10, 800, 600);
79  c1->SetBottomMargin(0.15);
80  hdmd->GetXaxis()->SetTitle("m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]");
81  hdmd->GetXaxis()->SetTitleOffset(1.4);
82 
83  // fit histogram hdmd with function f5 using the log-likelihood option
84  if (gROOT->GetListOfFunctions()->FindObject("f5"))
85  delete gROOT->GetFunction("f5");
86 
87  auto fdm5 = [](Double_t *xx, Double_t *par) -> Double_t {
88  const Double_t dxbin = (0.17 - 0.13) / 40; // Bin-width
89  Double_t x = xx[0];
90  if (x <= 0.13957)
91  return 0;
92  Double_t xp3 = (x - par[3]) * (x - par[3]);
93  Double_t res = dxbin * (par[0] * TMath::Power(x - 0.13957, par[1]) +
94  par[2] / 2.5066 / par[4] * TMath::Exp(-xp3 / 2 / par[4] / par[4]));
95  return res;
96  };
97 
98  auto f5 = new TF1("f5", fdm5, 0.139, 0.17, 5);
99  f5->SetParameters(1000000, .25, 2000, .1454, .001);
100  hdmd->Fit("f5", "lr");
101 
102  // Check the result of the fit
103  Double_t ref_f5[4] = {959915.0, 0.351114, 1185.03, 0.145569};
104  for (int i : {0, 1, 2, 3}) {
105  if ((TMath::Abs((f5->GetParameters())[i] - ref_f5[i]) / ref_f5[i]) > 0.001) {
106  std::cout << "\n >>> Test failure: fit to 'f5': parameter '" << f5->GetParName(i) << "' has wrong value ("
107  << (f5->GetParameters())[i] << ": expected" << ref_f5[i] << ") \n";
108  if (lfn)
109  gSystem->RedirectOutput(0, 0, &redH);
110  return -1;
111  }
112  }
113 
114  // create the canvas for tau d0
115  gStyle->SetOptFit(0);
116  gStyle->SetOptStat(1100);
117  auto c2 = new TCanvas("c2", "tauD0", 100, 100, 800, 600);
118  c2->SetGrid();
119  c2->SetBottomMargin(0.15);
120 
121  // Project slices of 2-d histogram h2 along X , then fit each slice
122  // with function f2 and make a histogram for each fit parameter
123  // Note that the generated histograms are added to the list of objects
124  // in the current directory.
125  if (gROOT->GetListOfFunctions()->FindObject("f2"))
126  delete gROOT->GetFunction("f2");
127 
128  auto fdm2 = [](Double_t *xx, Double_t *par) -> Double_t {
129  const auto dxbin = (0.17 - 0.13) / 40; // Bin-width
130  const auto sigma = 0.0012;
131  auto x = xx[0];
132  if (x <= 0.13957)
133  return 0;
134  auto xp3 = (x - 0.1454) * (x - 0.1454);
135  auto res = dxbin * (par[0] * TMath::Power(x - 0.13957, 0.25) +
136  par[1] / 2.5066 / sigma * TMath::Exp(-xp3 / 2 / sigma / sigma));
137  return res;
138  };
139 
140  auto f2 = new TF1("f2", fdm2, 0.139, 0.17, 2);
141  f2->SetParameters(10000, 10);
142 
143  // Restrict to three bins in this example
144  std::cout << "doFit: restricting fit to two bins only in this example...\n";
145 
146  h2->FitSlicesX(f2, 10, 20, 10, "g5 l");
147 
148  // Check the result of the fit
149  Double_t ref_f2[2] = {52432.2, 105.481};
150  for (int i : {0, 1}) {
151  if ((TMath::Abs((f2->GetParameters())[i] - ref_f2[i]) / ref_f2[i]) > 0.001) {
152  std::cout << "\n >>> Test failure: fit to 'f2': parameter '" << f2->GetParName(i) << "' has wrong value ("
153  << (f2->GetParameters())[i] << ": expected" << ref_f2[i] << ") \n";
154  if (lfn)
155  gSystem->RedirectOutput(0, 0, &redH);
156  return -1;
157  }
158  }
159 
160  auto h2_1 = (TH1D *)gDirectory->Get("h2_1");
161  h2_1->GetXaxis()->SetTitle("#tau[ps]");
162  h2_1->SetMarkerStyle(21);
163  h2_1->Draw();
164  c2->Update();
165  auto line = new TLine(0, 0, 0, c2->GetUymax());
166  line->Draw();
167 
168  // Have the number of entries on the first histogram (to cross check when running
169  // with entry lists)
170  auto psdmd = (TPaveStats *)hdmd->GetListOfFunctions()->FindObject("stats");
171  psdmd->SetOptStat(1110);
172  c1->Modified();
173 
174  if (lfn)
175  gSystem->RedirectOutput(0, 0, &redH);
176 
177  return 0;
178 };
179 
180 // This is the function invoked during the processing of the trees.
181 auto doH1 = [](TTreeReader &reader) {
182 
183  // Histograms
184  auto hdmd = new TH1F("hdmd", "Dm_d", 40, 0.13, 0.17);
185  auto h2 = new TH2F("h2", "ptD0 vs Dm_d", 30, 0.135, 0.165, 30, -3, 6);
186 
187  TTreeReaderValue<Float_t> fPtds_d(reader, "ptds_d");
188  TTreeReaderValue<Float_t> fEtads_d(reader, "etads_d");
189  TTreeReaderValue<Float_t> fDm_d(reader, "dm_d");
190  TTreeReaderValue<Int_t> fIk(reader, "ik");
191  TTreeReaderValue<Int_t> fIpi(reader, "ipi");
192  TTreeReaderValue<Int_t> fIpis(reader, "ipis");
193  TTreeReaderValue<Float_t> fPtd0_d(reader, "ptd0_d");
194  TTreeReaderValue<Float_t> fMd0_d(reader, "md0_d");
195  TTreeReaderValue<Float_t> fRpd0_t(reader, "rpd0_t");
196  TTreeReaderArray<Int_t> fNhitrp(reader, "nhitrp");
197  TTreeReaderArray<Float_t> fRstart(reader, "rstart");
198  TTreeReaderArray<Float_t> fRend(reader, "rend");
199  TTreeReaderArray<Float_t> fNlhk(reader, "nlhk");
200  TTreeReaderArray<Float_t> fNlhpi(reader, "nlhpi");
201  TTreeReaderValue<Int_t> fNjets(reader, "njets");
202 
203  while (reader.Next()) {
204 
205  // Return as soon as a bad entry is detected
206  if (TMath::Abs(*fMd0_d - 1.8646) >= 0.04)
207  continue;
208  if (*fPtds_d <= 2.5)
209  continue;
210  if (TMath::Abs(*fEtads_d) >= 1.5)
211  continue;
212  (*fIk)--; // original fIk used f77 convention starting at 1
213  (*fIpi)--;
214 
215  if (fNhitrp.At(*fIk) * fNhitrp.At(*fIpi) <= 1)
216  continue;
217 
218  if (fRend.At(*fIk) - fRstart.At(*fIk) <= 22)
219  continue;
220  if (fRend.At(*fIpi) - fRstart.At(*fIpi) <= 22)
221  continue;
222  if (fNlhk.At(*fIk) <= 0.1)
223  continue;
224  if (fNlhpi.At(*fIpi) <= 0.1)
225  continue;
226  (*fIpis)--;
227  if (fNlhpi.At(*fIpis) <= 0.1)
228  continue;
229  if (*fNjets < 1)
230  continue;
231 
232  // Fill the histograms
233  hdmd->Fill(*fDm_d);
234  h2->Fill(*fDm_d, *fRpd0_t / 0.029979 * 1.8646 / *fPtd0_d);
235  }
236 
237  // Return a list
238  auto l = new TList;
239  l->Add(hdmd);
240  l->Add(h2);
241  l->SetOwner(kFALSE);
242 
243  return l;
244 };
245 
246 // This is the function invoked during the processing of the trees to create a TEntryList
247 auto doH1fillList = [](TTreeReader &reader) {
248 
249  // Entry list
250  auto elist = new TEntryList("elist", "H1 selection from Cut");
251 
252  TTreeReaderValue<Float_t> fPtds_d(reader, "ptds_d");
253  TTreeReaderValue<Float_t> fEtads_d(reader, "etads_d");
254  TTreeReaderValue<Int_t> fIk(reader, "ik");
255  TTreeReaderValue<Int_t> fIpi(reader, "ipi");
256  TTreeReaderValue<Int_t> fIpis(reader, "ipis");
257  TTreeReaderValue<Float_t> fMd0_d(reader, "md0_d");
258  TTreeReaderArray<Int_t> fNhitrp(reader, "nhitrp");
259  TTreeReaderArray<Float_t> fRstart(reader, "rstart");
260  TTreeReaderArray<Float_t> fRend(reader, "rend");
261  TTreeReaderArray<Float_t> fNlhk(reader, "nlhk");
262  TTreeReaderArray<Float_t> fNlhpi(reader, "nlhpi");
263  TTreeReaderValue<Int_t> fNjets(reader, "njets");
264 
265  while (reader.Next()) {
266 
267  // Return as soon as a bad entry is detected
268  if (TMath::Abs(*fMd0_d - 1.8646) >= 0.04)
269  continue;
270  if (*fPtds_d <= 2.5)
271  continue;
272  if (TMath::Abs(*fEtads_d) >= 1.5)
273  continue;
274  (*fIk)--; // original fIk used f77 convention starting at 1
275  (*fIpi)--;
276 
277  if (fNhitrp.At(*fIk) * fNhitrp.At(*fIpi) <= 1)
278  continue;
279 
280  if (fRend.At(*fIk) - fRstart.At(*fIk) <= 22)
281  continue;
282  if (fRend.At(*fIpi) - fRstart.At(*fIpi) <= 22)
283  continue;
284  if (fNlhk.At(*fIk) <= 0.1)
285  continue;
286  if (fNlhpi.At(*fIpi) <= 0.1)
287  continue;
288  (*fIpis)--;
289  if (fNlhpi.At(*fIpis) <= 0.1)
290  continue;
291  if (*fNjets < 1)
292  continue;
293 
294  // Fill the entry list
295  elist->Enter(reader.GetCurrentEntry(), reader.GetTree());
296  }
297 
298  return elist;
299 };
300 
301 // This is the function invoked during the processing of the trees using a TEntryList
302 auto doH1useList = [](TTreeReader &reader) {
303 
304  // Histograms
305  auto hdmd = new TH1F("hdmd", "Dm_d", 40, 0.13, 0.17);
306  auto h2 = new TH2F("h2", "ptD0 vs Dm_d", 30, 0.135, 0.165, 30, -3, 6);
307 
308  TTreeReaderValue<Float_t> fDm_d(reader, "dm_d");
309  TTreeReaderValue<Float_t> fPtd0_d(reader, "ptd0_d");
310  TTreeReaderValue<Float_t> fRpd0_t(reader, "rpd0_t");
311 
312  while (reader.Next()) {
313  // Fill the histograms
314  hdmd->Fill(*fDm_d);
315  h2->Fill(*fDm_d, *fRpd0_t / 0.029979 * 1.8646 / *fPtd0_d);
316  }
317 
318  // Return a list
319  auto l = new TList;
320  l->Add(hdmd);
321  l->Add(h2);
322  l->SetOwner(kFALSE);
323 
324  return l;
325 };
l
auto * l
Definition: textangle.C:4
TLine
A simple line.
Definition: TLine.h:22
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:578
TTreeReaderValue< Float_t >
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
TTreeReaderArray
An interface for reading collections stored in ROOT columnar datasets.
Definition: TTreeReaderArray.h:75
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:727
Int_t
int Int_t
Definition: RtypesCore.h:45
x
Double_t x[n]
Definition: legend1.C:17
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
fdm5
Double_t fdm5(Double_t *xx, Double_t *par)
Definition: h1analysisProxy.h:14
TStyle::SetOptFit
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:1541
TPaveStats
The histogram statistics painter class.
Definition: TPaveStats.h:18
RedirectHandle_t
Definition: TSystem.h:203
TH1::GetMean
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7428
TSystem::RedirectOutput
virtual Int_t RedirectOutput(const char *name, const char *mode="a", RedirectHandle_t *h=nullptr)
Redirect standard output (stdout, stderr) to the specified file.
Definition: TSystem.cxx:1711
gStyle
R__EXTERN TStyle * gStyle
Definition: TStyle.h:412
dxbin
const Double_t dxbin
Definition: h1analysisProxy.h:10
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
TTreeReader
A simple, robust and fast interface to read values from ROOT columnar datasets such as TTree,...
Definition: TTreeReader.h:44
gDirectory
#define gDirectory
Definition: TDirectory.h:290
line
TLine * line
Definition: entrylistblock_figure1.C:235
TStyle::SetOptStat
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:1589
TEntryList
A List of entry numbers in a TTree or TChain.
Definition: TEntryList.h:26
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
TObject::Draw
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:197
Double_t
double Double_t
Definition: RtypesCore.h:59
TCanvas
The Canvas class.
Definition: TCanvas.h:23
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
c2
return c2
Definition: legend2.C:14
TF1
1-Dim function class
Definition: TF1.h:213
fdm2
Double_t fdm2(Double_t *xx, Double_t *par)
Definition: h1analysisProxy.h:25
TList
A doubly linked list.
Definition: TList.h:44
gROOT
#define gROOT
Definition: TROOT.h:406
c1
return c1
Definition: legend1.C:41