]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/macros/PID/SystematicErrorEstimation.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / macros / PID / SystematicErrorEstimation.C
CommitLineData
c683985a 1#include "TCanvas.h"
2#include "TFile.h"
3#include "TGraphErrors.h"
4#include "TGraphAsymmErrors.h"
5#include "TGraph.h"
6#include "TF1.h"
7#include "TH1D.h"
8#include "TLegend.h"
9#include "TMath.h"
10#include "TString.h"
11#include "TStyle.h"
12
13#include "AliPID.h"
14
15#include <iostream>
16#include <iomanip>
17
18#include "SystematicErrorUtils.h"
19
20const Int_t numSpecies = 5;
21
22//________________________________________________________
23TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t /*nSigma*/,
24 const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr,
25 Bool_t ignoreSigmaErrors)
26{
27 // For every bin:
28 // Since the method with the root finding already takes into account the statistical error,
29 // there is no need to use nSigma > 0.
30 // If the statistical error is ignored, nevertheless don't use nSigma > 0 because this might
31 // give zero systematic error for high pT, which is usually not accepted by people, although
32 // the natural point of view "no systematic visible for given statistical error" is reasonable to me.
33
34 Double_t ymax = 0;
35 Double_t ymin = 0;
36
37
38 // Just for drawing
39 for (Int_t j = 0; j < numHistos; j++) {
40 hSystematics[j] = new TH1F(*histos[j]);
41 hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));
42 hSystematics[j]->Reset();
43 hSystematics[j]->GetXaxis()->SetRange(0, -1);
44
45 for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {
46 hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));
47 hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) -
48 TMath::Power(histos[j]->GetBinError(bin), 2))));
49
50 if (hSystematics[j]->GetBinError(bin) == 0)
51 hSystematics[j]->SetBinError(bin, 1e-10);
52 Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);
53 if (temp > ymax)
54 ymax = temp;
55
56 temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);
57 if (temp < ymin)
58 ymin = temp;
59 }
60 }
61
62 TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
63 canv->SetGridy(1);
64
65 hSystematics[reference]->Draw("e p");
66 hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);
67 for (Int_t j = 0; j < numHistos; j++) {
68 if (j == reference)
69 continue;
70
71 hSystematics[j]->SetMarkerStyle(20 + j);
72 hSystematics[j]->Draw("e p same");
73 }
74
75 TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
76 legend->SetBorderSize(0);
77 legend->SetFillColor(0);
78
79 for (Int_t j = 0; j < numHistos; j++) {
80 legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");
81 }
82 legend->Draw();
83
84
85 const Int_t nBins = histos[reference]->GetNbinsX();
86 Double_t x[nBins];
87 Double_t y[nBins];
88 Double_t xerr[nBins];
89 Double_t yerrl[nBins];
90 Double_t yerrh[nBins];
91
92 Double_t meansForFit[numHistos];
93 Double_t sigmasForFit[numHistos];
94
95 for (Int_t bin = 0; bin < nBins; bin++) {
96 x[bin] = histos[reference]->GetBinCenter(bin + 1);
97 xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;
98 y[bin] = histos[reference]->GetBinContent(bin + 1);
99
100 for (Int_t j = 0; j < numHistos; j++) {
101 meansForFit[j] = histos[j]->GetBinContent(bin + 1);
102 sigmasForFit[j] = histos[j]->GetBinError(bin + 1);
103 }
104
105 yerrl[bin] = yerrh[bin] = findSystematicError(numHistos, meansForFit, sigmasForFit, ignoreSigmaErrors);
106 }
107
108 TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);
109 *gr = gTemp;
110 (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));
111 (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());
112 //(*gr)->SetFillColor(kGray);
113 (*gr)->SetFillStyle(0);//3004 + reference);
114
115 return canv;
116}
117
118
119/*OLD
120//________________________________________________________
121TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t nSigma,
122 const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr)
123{
124 Double_t ymax = 0;
125 Double_t ymin = 0;
126
127 for (Int_t j = 0; j < numHistos; j++) {
128 hSystematics[j] = new TH1F(*histos[j]);
129 hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));
130 hSystematics[j]->Reset();
131 hSystematics[j]->GetXaxis()->SetRange(0, -1);
132
133 for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {
134 hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));
135 hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) -
136 TMath::Power(histos[j]->GetBinError(bin), 2))));
137
138 if (hSystematics[j]->GetBinError(bin) == 0)
139 hSystematics[j]->SetBinError(bin, 1e-10);
140 Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);
141 if (temp > ymax)
142 ymax = temp;
143
144 temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);
145 if (temp < ymin)
146 ymin = temp;
147 }
148 }
149
150 TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
151 canv->SetGridy(1);
152
153 hSystematics[reference]->Draw("e p");
154 hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);
155 for (Int_t j = 0; j < numHistos; j++) {
156 if (j == reference)
157 continue;
158
159 hSystematics[j]->SetMarkerStyle(20 + j);
160 hSystematics[j]->Draw("e p same");
161 }
162
163 TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
164 legend->SetBorderSize(0);
165 legend->SetFillColor(0);
166
167 for (Int_t j = 0; j < numHistos; j++) {
168 legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");
169 }
170 legend->Draw();
171
172
173 const Int_t nBins = histos[reference]->GetNbinsX();
174 Double_t x[nBins];
175 Double_t y[nBins];
176 Double_t xerr[nBins];
177 Double_t yerrl[nBins];
178 Double_t yerrh[nBins];
179
180 for (Int_t bin = 0; bin < nBins; bin++) {
181 x[bin] = histos[reference]->GetBinCenter(bin + 1);
182 xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;
183 y[bin] = histos[reference]->GetBinContent(bin + 1);
184
185 // Take all points that are more than nSigma sigma away from 0.
186 // If there are at least 2 such points, take the difference between
187 // the extreme values (i.e. maximum and minimum) as a measure of
188 // the systematics
189 Int_t count = 0;
190 Double_t deltaMin = 0;
191 Double_t deltaMax = 0;
192
193 for (Int_t j = 0; j < numHistos; j++) {
194 if (hSystematics[j]->GetBinError(bin + 1) == 0) // Per definition always true for reference histo
195 continue;
196
197 Double_t delta = hSystematics[j]->GetBinContent(bin + 1);
198 if (TMath::Abs(delta / hSystematics[j]->GetBinError(bin + 1)) > nSigma) {
199 //if (count == 0) {
200 // deltaMin = delta;
201 // deltaMax = delta;
202 //}
203 //else {
204 if (delta < deltaMin)
205 deltaMin = delta;
206 if (delta > deltaMax)
207 deltaMax = delta;
208 //}
209 count++;
210 }
211 }
212
213 //if (deltaMax > 0.)
214 // yerrh[bin] = deltaMax;
215 //else
216 // yerrh[bin] = 0.;
217 //
218 //if (deltaMin < 0.)
219 // yerrl[bin] = -deltaMin;
220 //else
221 // yerrl[bin] = 0.;
222
223 if (count < 1) // Reference histo is not counted. One can only do systematics if there is at least one other histogram
224 yerrl[bin] = yerrh[bin] = 0.;
225 else
226 yerrl[bin] = yerrh[bin] = (deltaMax - deltaMin) / TMath::Sqrt(2);
227
228 }
229
230 TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);
231 *gr = gTemp;
232 (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));
233 (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());
234 //(*gr)->SetFillColor(kGray);
235 (*gr)->SetFillStyle(0);//3004 + reference);
236
237 return canv;
238}*/
239
240
241//________________________________________________________
242TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TH1F*** hist, Int_t reference,
243 TGraphAsymmErrors** gr)
244{
245 TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
246 canv->SetGridx(1);
247 canv->SetGridy(1);
248 canv->SetLogx(1);
249 for (Int_t i = 0; i < numSpecies; i++) {
250 hist[i][reference]->GetYaxis()->SetRangeUser(0.0, 1.0);
251 hist[i][reference]->GetYaxis()->SetTitle(canvTitle.Data());
252 hist[i][reference]->GetXaxis()->SetRangeUser(pLow, pHigh);
253 //hist[i][reference]->SetFillStyle(3004 + i);
254 //hist[i][reference]->SetFillColor(kGray);
255 hist[i][reference]->SetFillStyle(0);
256 hist[i][reference]->SetFillColor(hist[i][reference]->GetMarkerColor());
257 hist[i][reference]->SetLineColor(hist[i][reference]->GetMarkerColor());
258 }
259 hist[2][reference]->SetMarkerStyle(20);
260 hist[2][reference]->Draw("e p");
261 hist[0][reference]->SetMarkerStyle(21);
262 hist[0][reference]->Draw("e p same");
263 hist[1][reference]->SetMarkerStyle(22);
264 hist[1][reference]->Draw("e p same");
265 hist[3][reference]->SetMarkerStyle(29);
266 hist[3][reference]->Draw("e p same");
267 hist[4][reference]->SetMarkerStyle(30);
268 hist[4][reference]->Draw("e p same");
269
270 gr[0]->Draw("2 same");
271 gr[1]->Draw("2 same");
272 gr[2]->Draw("2 same");
273 gr[3]->Draw("2 same");
274 gr[4]->Draw("2 same");
275
276 TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
277 legend->SetBorderSize(0);
278 legend->SetFillColor(0);
279 legend->AddEntry(hist[2][reference], "#pi", "flp");
280 legend->AddEntry(hist[0][reference], "e", "flp");
281 legend->AddEntry(hist[1][reference], "K", "flp");
282 legend->AddEntry(hist[3][reference], "p", "flp");
283 legend->AddEntry(hist[4][reference], "#mu", "flp");
284 legend->Draw();
285
286 return canv;
287}
288
289
290//________________________________________________________
291TH1F* loadHisto(const TString histName, TFile* f)
292{
293 if (!f) {
294 std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;
295 return 0x0;
296 }
297
298 TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));
299 if (!hTemp) {
300 std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;
301 return 0x0;
302 }
303
304 return hTemp;
305}
306
307
308//________________________________________________________
309Int_t SystematicErrorEstimation(const TString path, const TString outFileTitle, const TString* fileNames, const TString* histTitles,
310 const Int_t numFiles, const Double_t nSigma, const Bool_t ignoreSigmaErrors)
311{
312 if (!fileNames || numFiles < 1)
313 return -1;
314
315 TFile* f[numFiles];
316 TH1F** hFractions[numSpecies];
317 for (Int_t i = 0; i < numSpecies; i++)
318 hFractions[i] = new TH1F*[numFiles];
319
320 const Int_t reference = 0;
321 TH1F* hYields[numSpecies]; // Only the reference yields
322
323 const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };
324
325 const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };
326
327 for (Int_t iFile = 0; iFile < numFiles; iFile++) {
328 f[iFile] = TFile::Open(fileNames[iFile].Data());
329 if (!f[iFile]) {
330 std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;
331 return -1;
332 }
333
334 // Extract the data histograms
335 for (Int_t i = 0; i < numSpecies; i++) {
336 hFractions[i][iFile] = loadHisto(histNames[i], f[iFile]);
337 if (!hFractions[i][iFile])
338 return -1;
339
340 if (iFile == reference) {
341 hYields[i] = loadHisto(histNamesYields[i], f[iFile]);
342 if (!hYields[i])
343 return -1;
344 }
345 }
346 }
347
348
349 TGraphAsymmErrors* grSysErrors[numSpecies] = {0x0,};
350 TGraphAsymmErrors* grSysErrorsYields[numSpecies] = {0x0,};
351
352 TH1F* hSystematicsPions[numFiles];
353 TCanvas* cSystematicsPions = calculateSystematics("cSystematicsPions", "Systematics Pions", hFractions[2], numFiles,
354 AliPID::kPion, nSigma,
355 histTitles, reference, hSystematicsPions, &grSysErrors[2], ignoreSigmaErrors);
356
357 TH1F* hSystematicsElectrons[numFiles];
358 TCanvas* cSystematicsElectrons = calculateSystematics("cSystematicsElectrons", "Systematics Electrons", hFractions[0], numFiles,
359 AliPID::kElectron,
360 nSigma, histTitles, reference, hSystematicsElectrons,
361 &grSysErrors[0], ignoreSigmaErrors);
362
363 TH1F* hSystematicsKaons[numFiles];
364 TCanvas* cSystematicsKaons = calculateSystematics("cSystematicsKaons", "Systematics Kaons", hFractions[1], numFiles, AliPID::kKaon, nSigma,
365 histTitles, reference, hSystematicsKaons, &grSysErrors[1], ignoreSigmaErrors);
366
367 TH1F* hSystematicsProtons[numFiles];
368 TCanvas* cSystematicsProtons = calculateSystematics("cSystematicsProtons", "Systematics Protons", hFractions[3], numFiles,
369 AliPID::kProton, nSigma,
370 histTitles, reference, hSystematicsProtons, &grSysErrors[3], ignoreSigmaErrors);
371
372 TH1F* hSystematicsMuons[numFiles];
373 TCanvas* cSystematicsMuons = calculateSystematics("cSystematicsMuons", "Systematics Muons", hFractions[4], numFiles,
374 AliPID::kMuon, nSigma,
375 histTitles, reference, hSystematicsMuons, &grSysErrors[4], ignoreSigmaErrors);
376
377 Double_t pLow = 0.15;
378 Double_t pHigh = 50.;
379 TCanvas* cFractionsWithSystematicError = DrawFractionHistos("cFractionsWithSystematicError", "Particle fractions", pLow, pHigh, hFractions, reference,
380 grSysErrors);
381
382
383 //TODO At the moment, the error of the fractions and the yield is just a constant factor (number of tracks in that bin)
384 // (-> But this can change in future (I have to think about it a little bit more carefully)).
385 // Thus, the relative errors are the same for fractions and yields and I can just use this fact to
386 // transform the errors from the fractions to those of the yields.
387 // However, this causes trouble in case of fraction = 0. Therefore, sum up the yields to the total yield and use this for scaling
388 for (Int_t i = 0; i < numSpecies; i++) {
389 grSysErrorsYields[i] = new TGraphAsymmErrors(*grSysErrors[i]);
390 TString name = grSysErrors[i]->GetName();
391 name.ReplaceAll("systematicError_", "systematicErrorYields_");
392 grSysErrorsYields[i]->SetName(name.Data());
393
394 for (Int_t ind = 0; ind < grSysErrorsYields[i]->GetN(); ind++) {
395 Double_t totalYield = 0;
396 for (Int_t j = 0; j < numSpecies; j++)
397 totalYield += hYields[j]->GetBinContent(ind + 1);
398
399 const Double_t yield = hYields[i]->GetBinContent(ind + 1);
400 const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);
401 const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);
402
403 grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);
404 grSysErrorsYields[i]->SetPointEYhigh(ind, totalYield * sysErrorHigh);
405 grSysErrorsYields[i]->SetPointEYlow(ind, totalYield * sysErrorLow);
406
407 /*
408 Double_t totalYield = 0;
409 for (Int_t j = 0; j < numSpecies; j++)
410 totalYield += hYields[j]->GetBinContent(ind + 1);
411
412 const Double_t yield = hYields[i]->GetBinContent(ind + 1);
413 const Double_t fraction = hFractions[i][reference]->GetBinContent(ind + 1);
414 const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);
415 const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);
416
417 if (fraction <= 0.) {
418 printf("Error: Fraction = 0 for species %d. Cannot transform error....\n", i);
419 return -1;
420 }
421 const Double_t sysErrorLowRel = sysErrorLow / fraction;
422 const Double_t sysErrorHighRel = sysErrorHigh / fraction;
423
424 grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);
425 grSysErrorsYields[i]->SetPointEYhigh(ind, sysErrorHighRel * yield);
426 grSysErrorsYields[i]->SetPointEYlow(ind, sysErrorLowRel * yield);
427 */
428 }
429 }
430
431
432 // Output file
433 TFile* fSave = 0x0;
434 TDatime daTime;
435 TString saveFileName;
436
437 saveFileName = Form("outputSystematics_%s_nSigma%.1f__%04d_%02d_%02d.root", outFileTitle.Data(), nSigma, daTime.GetYear(),
438 daTime.GetMonth(), daTime.GetDay());
439
440 fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");
441 if (!fSave) {
442 std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;
443 return -1;
444 }
445
446 // Save final results
447 fSave->cd();
448
449 for (Int_t i = 0; i < numFiles; i++) {
450 if (hSystematicsElectrons[i])
451 hSystematicsElectrons[i]->Write();
452
453 if (hSystematicsPions[i])
454 hSystematicsPions[i]->Write();
455
456 if (hSystematicsKaons[i])
457 hSystematicsKaons[i]->Write();
458
459 if (hSystematicsProtons[i])
460 hSystematicsProtons[i]->Write();
461 if (hSystematicsMuons[i])
462 hSystematicsMuons[i]->Write();
463 }
464
465 if (cSystematicsElectrons)
466 cSystematicsElectrons->Write();
467
468 if (cSystematicsPions)
469 cSystematicsPions->Write();
470
471 if (cSystematicsKaons)
472 cSystematicsKaons->Write();
473
474 if (cSystematicsProtons)
475 cSystematicsProtons->Write();
476
477 if (cSystematicsMuons)
478 cSystematicsMuons->Write();
479
480 if (cFractionsWithSystematicError)
481 cFractionsWithSystematicError->Write();
482
483 for (Int_t i = 0; i < numSpecies; i++) {
484 if (hFractions[i][reference])
485 hFractions[i][reference]->Write();
486
487 if (grSysErrors[i])
488 grSysErrors[i]->Write();
489
490 if (hYields[i])
491 hYields[i]->Write();
492
493 if (grSysErrorsYields[i])
494 grSysErrorsYields[i]->Write();
495 }
496
497 // Save list of file names in output file
498 TString listOfFileNames = "";
499 for (Int_t i = 0; i < numFiles; i++) {
500 listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));
501 }
502
503 TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),
504 Form("Used files for systematics: %s\n", listOfFileNames.Data()));
505 settings->Write();
506
507 fSave->Close();
508
509 delete cSystematicsElectrons;
510 delete cSystematicsPions;
511 delete cSystematicsKaons;
512 delete cSystematicsMuons;
513 delete cSystematicsProtons;
514 delete cFractionsWithSystematicError;
515
516 return 0;
517}