Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
12auto 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
59auto doFit = [](TList *out, const char *lfn = 0) -> Int_t {
60
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
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.
181auto 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
247auto 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
302auto 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};
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
#define gDirectory
Definition TDirectory.h:386
#define gROOT
Definition TROOT.h:405
R__EXTERN TStyle * gStyle
Definition TStyle.h:414
R__EXTERN TSystem * gSystem
Definition TSystem.h:560
The Canvas class.
Definition TCanvas.h:23
A List of entry numbers in a TTree or TChain.
Definition TEntryList.h:26
1-Dim function class
Definition TF1.h:213
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:577
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:7445
2-D histogram with a float per channel (see TH1 documentation)}
Definition TH2.h:257
Use the TLine constructor to create a simple line.
Definition TLine.h:22
A doubly linked list.
Definition TList.h:38
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:578
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
The histogram statistics painter class.
Definition TPaveStats.h:18
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
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:1542
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:1716
An interface for reading collections stored in ROOT columnar datasets.
An interface for reading values stored in ROOT columnar datasets.
A simple, robust and fast interface to read values from ROOT columnar datasets such as TTree,...
Definition TTreeReader.h:44
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:707
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:719
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TLine l
Definition textangle.C:4