-#include "TCanvas.h"\r
-#include "TFile.h"\r
-#include "TGraphErrors.h"\r
-#include "TGraphAsymmErrors.h"\r
-#include "TGraph.h"\r
-#include "TF1.h"\r
-#include "TH1D.h"\r
-#include "TLegend.h"\r
-#include "TMath.h"\r
-#include "TString.h"\r
-#include "TStyle.h"\r
-\r
-#include "AliPID.h"\r
-\r
-#include <iostream>\r
-#include <iomanip>\r
-\r
-#include "SystematicErrorUtils.h"\r
-\r
-const Int_t numSpecies = 5;\r
-\r
-//________________________________________________________\r
-TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t /*nSigma*/,\r
- const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr,\r
- Bool_t ignoreSigmaErrors)\r
-{\r
- // For every bin:\r
- // Since the method with the root finding already takes into account the statistical error,\r
- // there is no need to use nSigma > 0.\r
- // If the statistical error is ignored, nevertheless don't use nSigma > 0 because this might\r
- // give zero systematic error for high pT, which is usually not accepted by people, although\r
- // the natural point of view "no systematic visible for given statistical error" is reasonable to me.\r
- \r
- Double_t ymax = 0;\r
- Double_t ymin = 0;\r
- \r
- \r
- // Just for drawing\r
- for (Int_t j = 0; j < numHistos; j++) {\r
- hSystematics[j] = new TH1F(*histos[j]);\r
- hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));\r
- hSystematics[j]->Reset(); \r
- hSystematics[j]->GetXaxis()->SetRange(0, -1);\r
- \r
- for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {\r
- hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));\r
- hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) - \r
- TMath::Power(histos[j]->GetBinError(bin), 2))));\r
- \r
- if (hSystematics[j]->GetBinError(bin) == 0)\r
- hSystematics[j]->SetBinError(bin, 1e-10);\r
- Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);\r
- if (temp > ymax)\r
- ymax = temp;\r
- \r
- temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);\r
- if (temp < ymin)\r
- ymin = temp;\r
- }\r
- }\r
- \r
- TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);\r
- canv->SetGridy(1);\r
- \r
- hSystematics[reference]->Draw("e p");\r
- hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);\r
- for (Int_t j = 0; j < numHistos; j++) {\r
- if (j == reference)\r
- continue;\r
- \r
- hSystematics[j]->SetMarkerStyle(20 + j);\r
- hSystematics[j]->Draw("e p same");\r
- }\r
- \r
- TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932); \r
- legend->SetBorderSize(0);\r
- legend->SetFillColor(0);\r
- \r
- for (Int_t j = 0; j < numHistos; j++) {\r
- legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");\r
- }\r
- legend->Draw();\r
- \r
- \r
- const Int_t nBins = histos[reference]->GetNbinsX();\r
- Double_t x[nBins];\r
- Double_t y[nBins];\r
- Double_t xerr[nBins];\r
- Double_t yerrl[nBins];\r
- Double_t yerrh[nBins];\r
- \r
- Double_t meansForFit[numHistos];\r
- Double_t sigmasForFit[numHistos];\r
- \r
- for (Int_t bin = 0; bin < nBins; bin++) {\r
- x[bin] = histos[reference]->GetBinCenter(bin + 1);\r
- xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;\r
- y[bin] = histos[reference]->GetBinContent(bin + 1);\r
- \r
- for (Int_t j = 0; j < numHistos; j++) {\r
- meansForFit[j] = histos[j]->GetBinContent(bin + 1);\r
- sigmasForFit[j] = histos[j]->GetBinError(bin + 1);\r
- }\r
- \r
- yerrl[bin] = yerrh[bin] = findSystematicError(numHistos, meansForFit, sigmasForFit, ignoreSigmaErrors);\r
- }\r
- \r
- TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);\r
- *gr = gTemp;\r
- (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));\r
- (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());\r
- //(*gr)->SetFillColor(kGray);\r
- (*gr)->SetFillStyle(0);//3004 + reference);\r
- \r
- return canv;\r
-}\r
-\r
-\r
-/*OLD\r
-//________________________________________________________\r
-TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t nSigma,\r
- const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr)\r
-{\r
- Double_t ymax = 0;\r
- Double_t ymin = 0;\r
- \r
- for (Int_t j = 0; j < numHistos; j++) {\r
- hSystematics[j] = new TH1F(*histos[j]);\r
- hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));\r
- hSystematics[j]->Reset(); \r
- hSystematics[j]->GetXaxis()->SetRange(0, -1);\r
- \r
- for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {\r
- hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));\r
- hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) - \r
- TMath::Power(histos[j]->GetBinError(bin), 2))));\r
- \r
- if (hSystematics[j]->GetBinError(bin) == 0)\r
- hSystematics[j]->SetBinError(bin, 1e-10);\r
- Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);\r
- if (temp > ymax)\r
- ymax = temp;\r
- \r
- temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);\r
- if (temp < ymin)\r
- ymin = temp;\r
- }\r
- }\r
- \r
- TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);\r
- canv->SetGridy(1);\r
- \r
- hSystematics[reference]->Draw("e p");\r
- hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);\r
- for (Int_t j = 0; j < numHistos; j++) {\r
- if (j == reference)\r
- continue;\r
- \r
- hSystematics[j]->SetMarkerStyle(20 + j);\r
- hSystematics[j]->Draw("e p same");\r
- }\r
- \r
- TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932); \r
- legend->SetBorderSize(0);\r
- legend->SetFillColor(0);\r
- \r
- for (Int_t j = 0; j < numHistos; j++) {\r
- legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");\r
- }\r
- legend->Draw();\r
- \r
- \r
- const Int_t nBins = histos[reference]->GetNbinsX();\r
- Double_t x[nBins];\r
- Double_t y[nBins];\r
- Double_t xerr[nBins];\r
- Double_t yerrl[nBins];\r
- Double_t yerrh[nBins];\r
- \r
- for (Int_t bin = 0; bin < nBins; bin++) {\r
- x[bin] = histos[reference]->GetBinCenter(bin + 1);\r
- xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;\r
- y[bin] = histos[reference]->GetBinContent(bin + 1);\r
- \r
- // Take all points that are more than nSigma sigma away from 0.\r
- // If there are at least 2 such points, take the difference between\r
- // the extreme values (i.e. maximum and minimum) as a measure of\r
- // the systematics\r
- Int_t count = 0;\r
- Double_t deltaMin = 0;\r
- Double_t deltaMax = 0;\r
- \r
- for (Int_t j = 0; j < numHistos; j++) {\r
- if (hSystematics[j]->GetBinError(bin + 1) == 0) // Per definition always true for reference histo\r
- continue;\r
- \r
- Double_t delta = hSystematics[j]->GetBinContent(bin + 1);\r
- if (TMath::Abs(delta / hSystematics[j]->GetBinError(bin + 1)) > nSigma) {\r
- //if (count == 0) {\r
- // deltaMin = delta;\r
- // deltaMax = delta;\r
- //}\r
- //else {\r
- if (delta < deltaMin)\r
- deltaMin = delta;\r
- if (delta > deltaMax)\r
- deltaMax = delta;\r
- //}\r
- count++;\r
- }\r
- }\r
- \r
- //if (deltaMax > 0.) \r
- // yerrh[bin] = deltaMax;\r
- //else\r
- // yerrh[bin] = 0.;\r
- // \r
- //if (deltaMin < 0.)\r
- // yerrl[bin] = -deltaMin;\r
- //else\r
- // yerrl[bin] = 0.;\r
- \r
- if (count < 1) // Reference histo is not counted. One can only do systematics if there is at least one other histogram\r
- yerrl[bin] = yerrh[bin] = 0.;\r
- else\r
- yerrl[bin] = yerrh[bin] = (deltaMax - deltaMin) / TMath::Sqrt(2);\r
- \r
- }\r
- \r
- TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);\r
- *gr = gTemp;\r
- (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));\r
- (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());\r
- //(*gr)->SetFillColor(kGray);\r
- (*gr)->SetFillStyle(0);//3004 + reference);\r
- \r
- return canv;\r
-}*/\r
-\r
-\r
-//________________________________________________________\r
-TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TH1F*** hist, Int_t reference,\r
- TGraphAsymmErrors** gr)\r
-{\r
- TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);\r
- canv->SetGridx(1);\r
- canv->SetGridy(1);\r
- canv->SetLogx(1);\r
- for (Int_t i = 0; i < numSpecies; i++) {\r
- hist[i][reference]->GetYaxis()->SetRangeUser(0.0, 1.0);\r
- hist[i][reference]->GetYaxis()->SetTitle(canvTitle.Data());\r
- hist[i][reference]->GetXaxis()->SetRangeUser(pLow, pHigh);\r
- //hist[i][reference]->SetFillStyle(3004 + i);\r
- //hist[i][reference]->SetFillColor(kGray);\r
- hist[i][reference]->SetFillStyle(0);\r
- hist[i][reference]->SetFillColor(hist[i][reference]->GetMarkerColor());\r
- hist[i][reference]->SetLineColor(hist[i][reference]->GetMarkerColor());\r
- }\r
- hist[2][reference]->SetMarkerStyle(20);\r
- hist[2][reference]->Draw("e p");\r
- hist[0][reference]->SetMarkerStyle(21);\r
- hist[0][reference]->Draw("e p same");\r
- hist[1][reference]->SetMarkerStyle(22);\r
- hist[1][reference]->Draw("e p same");\r
- hist[3][reference]->SetMarkerStyle(29);\r
- hist[3][reference]->Draw("e p same");\r
- hist[4][reference]->SetMarkerStyle(30);\r
- hist[4][reference]->Draw("e p same");\r
- \r
- gr[0]->Draw("2 same");\r
- gr[1]->Draw("2 same");\r
- gr[2]->Draw("2 same");\r
- gr[3]->Draw("2 same");\r
- gr[4]->Draw("2 same");\r
- \r
- TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932); \r
- legend->SetBorderSize(0);\r
- legend->SetFillColor(0);\r
- legend->AddEntry(hist[2][reference], "#pi", "flp");\r
- legend->AddEntry(hist[0][reference], "e", "flp");\r
- legend->AddEntry(hist[1][reference], "K", "flp");\r
- legend->AddEntry(hist[3][reference], "p", "flp");\r
- legend->AddEntry(hist[4][reference], "#mu", "flp");\r
- legend->Draw();\r
- \r
- return canv;\r
-}\r
-\r
-\r
-//________________________________________________________\r
-TH1F* loadHisto(const TString histName, TFile* f)\r
-{\r
- if (!f) {\r
- std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;\r
- return 0x0;\r
- }\r
- \r
- TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));\r
- if (!hTemp) {\r
- std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;\r
- return 0x0;\r
- } \r
- \r
- return hTemp;\r
-}\r
-\r
-\r
-//________________________________________________________\r
-Int_t SystematicErrorEstimation(const TString path, const TString outFileTitle, const TString* fileNames, const TString* histTitles,\r
- const Int_t numFiles, const Double_t nSigma, const Bool_t ignoreSigmaErrors) \r
-{ \r
- if (!fileNames || numFiles < 1)\r
- return -1;\r
- \r
- TFile* f[numFiles];\r
- TH1F** hFractions[numSpecies];\r
- for (Int_t i = 0; i < numSpecies; i++) \r
- hFractions[i] = new TH1F*[numFiles];\r
- \r
- const Int_t reference = 0;\r
- TH1F* hYields[numSpecies]; // Only the reference yields\r
- \r
- const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };\r
- \r
- const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };\r
- \r
- for (Int_t iFile = 0; iFile < numFiles; iFile++) {\r
- f[iFile] = TFile::Open(fileNames[iFile].Data());\r
- if (!f[iFile]) {\r
- std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;\r
- return -1;\r
- }\r
- \r
- // Extract the data histograms\r
- for (Int_t i = 0; i < numSpecies; i++) {\r
- hFractions[i][iFile] = loadHisto(histNames[i], f[iFile]);\r
- if (!hFractions[i][iFile])\r
- return -1;\r
- \r
- if (iFile == reference) {\r
- hYields[i] = loadHisto(histNamesYields[i], f[iFile]);\r
- if (!hYields[i])\r
- return -1;\r
- }\r
- }\r
- }\r
- \r
- \r
- TGraphAsymmErrors* grSysErrors[numSpecies] = {0x0,};\r
- TGraphAsymmErrors* grSysErrorsYields[numSpecies] = {0x0,};\r
- \r
- TH1F* hSystematicsPions[numFiles];\r
- TCanvas* cSystematicsPions = calculateSystematics("cSystematicsPions", "Systematics Pions", hFractions[2], numFiles,\r
- AliPID::kPion, nSigma,\r
- histTitles, reference, hSystematicsPions, &grSysErrors[2], ignoreSigmaErrors);\r
-\r
- TH1F* hSystematicsElectrons[numFiles];\r
- TCanvas* cSystematicsElectrons = calculateSystematics("cSystematicsElectrons", "Systematics Electrons", hFractions[0], numFiles,\r
- AliPID::kElectron,\r
- nSigma, histTitles, reference, hSystematicsElectrons,\r
- &grSysErrors[0], ignoreSigmaErrors);\r
- \r
- TH1F* hSystematicsKaons[numFiles];\r
- TCanvas* cSystematicsKaons = calculateSystematics("cSystematicsKaons", "Systematics Kaons", hFractions[1], numFiles, AliPID::kKaon, nSigma,\r
- histTitles, reference, hSystematicsKaons, &grSysErrors[1], ignoreSigmaErrors);\r
- \r
- TH1F* hSystematicsProtons[numFiles]; \r
- TCanvas* cSystematicsProtons = calculateSystematics("cSystematicsProtons", "Systematics Protons", hFractions[3], numFiles,\r
- AliPID::kProton, nSigma,\r
- histTitles, reference, hSystematicsProtons, &grSysErrors[3], ignoreSigmaErrors);\r
- \r
- TH1F* hSystematicsMuons[numFiles];\r
- TCanvas* cSystematicsMuons = calculateSystematics("cSystematicsMuons", "Systematics Muons", hFractions[4], numFiles,\r
- AliPID::kMuon, nSigma,\r
- histTitles, reference, hSystematicsMuons, &grSysErrors[4], ignoreSigmaErrors);\r
- \r
- Double_t pLow = 0.15;\r
- Double_t pHigh = 50.;\r
- TCanvas* cFractionsWithSystematicError = DrawFractionHistos("cFractionsWithSystematicError", "Particle fractions", pLow, pHigh, hFractions, reference, \r
- grSysErrors);\r
- \r
- \r
- //TODO At the moment, the error of the fractions and the yield is just a constant factor (number of tracks in that bin)\r
- // (-> But this can change in future (I have to think about it a little bit more carefully)).\r
- // Thus, the relative errors are the same for fractions and yields and I can just use this fact to\r
- // transform the errors from the fractions to those of the yields.\r
- // However, this causes trouble in case of fraction = 0. Therefore, sum up the yields to the total yield and use this for scaling\r
- for (Int_t i = 0; i < numSpecies; i++) {\r
- grSysErrorsYields[i] = new TGraphAsymmErrors(*grSysErrors[i]);\r
- TString name = grSysErrors[i]->GetName();\r
- name.ReplaceAll("systematicError_", "systematicErrorYields_");\r
- grSysErrorsYields[i]->SetName(name.Data());\r
- \r
- for (Int_t ind = 0; ind < grSysErrorsYields[i]->GetN(); ind++) {\r
- Double_t totalYield = 0;\r
- for (Int_t j = 0; j < numSpecies; j++)\r
- totalYield += hYields[j]->GetBinContent(ind + 1);\r
- \r
- const Double_t yield = hYields[i]->GetBinContent(ind + 1);\r
- const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);\r
- const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);\r
- \r
- grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);\r
- grSysErrorsYields[i]->SetPointEYhigh(ind, totalYield * sysErrorHigh);\r
- grSysErrorsYields[i]->SetPointEYlow(ind, totalYield * sysErrorLow);\r
- \r
- /*\r
- Double_t totalYield = 0;\r
- for (Int_t j = 0; j < numSpecies; j++)\r
- totalYield += hYields[j]->GetBinContent(ind + 1);\r
- \r
- const Double_t yield = hYields[i]->GetBinContent(ind + 1);\r
- const Double_t fraction = hFractions[i][reference]->GetBinContent(ind + 1);\r
- const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);\r
- const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);\r
- \r
- if (fraction <= 0.) {\r
- printf("Error: Fraction = 0 for species %d. Cannot transform error....\n", i);\r
- return -1;\r
- }\r
- const Double_t sysErrorLowRel = sysErrorLow / fraction;\r
- const Double_t sysErrorHighRel = sysErrorHigh / fraction;\r
- \r
- grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);\r
- grSysErrorsYields[i]->SetPointEYhigh(ind, sysErrorHighRel * yield);\r
- grSysErrorsYields[i]->SetPointEYlow(ind, sysErrorLowRel * yield);\r
- */\r
- }\r
- }\r
- \r
- \r
- // Output file\r
- TFile* fSave = 0x0;\r
- TDatime daTime;\r
- TString saveFileName;\r
- \r
- saveFileName = Form("outputSystematics_%s_nSigma%.1f__%04d_%02d_%02d.root", outFileTitle.Data(), nSigma, daTime.GetYear(),\r
- daTime.GetMonth(), daTime.GetDay());\r
- \r
- fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");\r
- if (!fSave) {\r
- std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;\r
- return -1;\r
- }\r
- \r
- // Save final results\r
- fSave->cd();\r
- \r
- for (Int_t i = 0; i < numFiles; i++) {\r
- if (hSystematicsElectrons[i])\r
- hSystematicsElectrons[i]->Write();\r
- \r
- if (hSystematicsPions[i])\r
- hSystematicsPions[i]->Write();\r
- \r
- if (hSystematicsKaons[i])\r
- hSystematicsKaons[i]->Write();\r
- \r
- if (hSystematicsProtons[i])\r
- hSystematicsProtons[i]->Write();\r
- if (hSystematicsMuons[i])\r
- hSystematicsMuons[i]->Write();\r
- }\r
- \r
- if (cSystematicsElectrons)\r
- cSystematicsElectrons->Write();\r
- \r
- if (cSystematicsPions)\r
- cSystematicsPions->Write();\r
- \r
- if (cSystematicsKaons)\r
- cSystematicsKaons->Write();\r
- \r
- if (cSystematicsProtons)\r
- cSystematicsProtons->Write();\r
- \r
- if (cSystematicsMuons)\r
- cSystematicsMuons->Write();\r
- \r
- if (cFractionsWithSystematicError)\r
- cFractionsWithSystematicError->Write();\r
- \r
- for (Int_t i = 0; i < numSpecies; i++) {\r
- if (hFractions[i][reference])\r
- hFractions[i][reference]->Write();\r
- \r
- if (grSysErrors[i])\r
- grSysErrors[i]->Write();\r
- \r
- if (hYields[i])\r
- hYields[i]->Write();\r
- \r
- if (grSysErrorsYields[i])\r
- grSysErrorsYields[i]->Write();\r
- }\r
- \r
- // Save list of file names in output file\r
- TString listOfFileNames = "";\r
- for (Int_t i = 0; i < numFiles; i++) {\r
- listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));\r
- }\r
- \r
- TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),\r
- Form("Used files for systematics: %s\n", listOfFileNames.Data()));\r
- settings->Write();\r
- \r
- fSave->Close();\r
- \r
- delete cSystematicsElectrons;\r
- delete cSystematicsPions;\r
- delete cSystematicsKaons;\r
- delete cSystematicsMuons;\r
- delete cSystematicsProtons;\r
- delete cFractionsWithSystematicError;\r
- \r
- return 0;\r
-}\r
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TGraphErrors.h"
+#include "TGraphAsymmErrors.h"
+#include "TGraph.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TLegend.h"
+#include "TMath.h"
+#include "TString.h"
+#include "TStyle.h"
+
+#include "AliPID.h"
+
+#include <iostream>
+#include <iomanip>
+
+#include "SystematicErrorUtils.h"
+
+const Int_t numSpecies = 5;
+
+//________________________________________________________
+TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t /*nSigma*/,
+ const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr,
+ Bool_t ignoreSigmaErrors)
+{
+ // For every bin:
+ // Since the method with the root finding already takes into account the statistical error,
+ // there is no need to use nSigma > 0.
+ // If the statistical error is ignored, nevertheless don't use nSigma > 0 because this might
+ // give zero systematic error for high pT, which is usually not accepted by people, although
+ // the natural point of view "no systematic visible for given statistical error" is reasonable to me.
+
+ Double_t ymax = 0;
+ Double_t ymin = 0;
+
+
+ // Just for drawing
+ for (Int_t j = 0; j < numHistos; j++) {
+ hSystematics[j] = new TH1F(*histos[j]);
+ hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));
+ hSystematics[j]->Reset();
+ hSystematics[j]->GetXaxis()->SetRange(0, -1);
+
+ for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {
+ hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));
+ hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) -
+ TMath::Power(histos[j]->GetBinError(bin), 2))));
+
+ if (hSystematics[j]->GetBinError(bin) == 0)
+ hSystematics[j]->SetBinError(bin, 1e-10);
+ Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);
+ if (temp > ymax)
+ ymax = temp;
+
+ temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);
+ if (temp < ymin)
+ ymin = temp;
+ }
+ }
+
+ TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
+ canv->SetGridy(1);
+
+ hSystematics[reference]->Draw("e p");
+ hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);
+ for (Int_t j = 0; j < numHistos; j++) {
+ if (j == reference)
+ continue;
+
+ hSystematics[j]->SetMarkerStyle(20 + j);
+ hSystematics[j]->Draw("e p same");
+ }
+
+ TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
+ legend->SetBorderSize(0);
+ legend->SetFillColor(0);
+
+ for (Int_t j = 0; j < numHistos; j++) {
+ legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");
+ }
+ legend->Draw();
+
+
+ const Int_t nBins = histos[reference]->GetNbinsX();
+ Double_t x[nBins];
+ Double_t y[nBins];
+ Double_t xerr[nBins];
+ Double_t yerrl[nBins];
+ Double_t yerrh[nBins];
+
+ Double_t meansForFit[numHistos];
+ Double_t sigmasForFit[numHistos];
+
+ for (Int_t bin = 0; bin < nBins; bin++) {
+ x[bin] = histos[reference]->GetBinCenter(bin + 1);
+ xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;
+ y[bin] = histos[reference]->GetBinContent(bin + 1);
+
+ for (Int_t j = 0; j < numHistos; j++) {
+ meansForFit[j] = histos[j]->GetBinContent(bin + 1);
+ sigmasForFit[j] = histos[j]->GetBinError(bin + 1);
+ }
+
+ yerrl[bin] = yerrh[bin] = findSystematicError(numHistos, meansForFit, sigmasForFit, ignoreSigmaErrors);
+ }
+
+ TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);
+ *gr = gTemp;
+ (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));
+ (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());
+ //(*gr)->SetFillColor(kGray);
+ (*gr)->SetFillStyle(0);//3004 + reference);
+
+ return canv;
+}
+
+
+/*OLD
+//________________________________________________________
+TCanvas* calculateSystematics(TString canvName, TString canvTitle, TH1F** histos, Int_t numHistos, Int_t speciesID, Double_t nSigma,
+ const TString* systematicsHistosName, Int_t reference, TH1F** hSystematics, TGraphAsymmErrors** gr)
+{
+ Double_t ymax = 0;
+ Double_t ymin = 0;
+
+ for (Int_t j = 0; j < numHistos; j++) {
+ hSystematics[j] = new TH1F(*histos[j]);
+ hSystematics[j]->SetName(Form("%s_%s", systematicsHistosName[j].Data(), AliPID::ParticleName(speciesID)));
+ hSystematics[j]->Reset();
+ hSystematics[j]->GetXaxis()->SetRange(0, -1);
+
+ for (Int_t bin = 1; bin <= histos[j]->GetNbinsX(); bin++) {
+ hSystematics[j]->SetBinContent(bin, histos[reference]->GetBinContent(bin) - histos[j]->GetBinContent(bin));
+ hSystematics[j]->SetBinError(bin, TMath::Sqrt(TMath::Abs(TMath::Power(histos[reference]->GetBinError(bin), 2) -
+ TMath::Power(histos[j]->GetBinError(bin), 2))));
+
+ if (hSystematics[j]->GetBinError(bin) == 0)
+ hSystematics[j]->SetBinError(bin, 1e-10);
+ Double_t temp = hSystematics[j]->GetBinContent(bin) + hSystematics[j]->GetBinError(bin);
+ if (temp > ymax)
+ ymax = temp;
+
+ temp = hSystematics[j]->GetBinContent(bin) - hSystematics[j]->GetBinError(bin);
+ if (temp < ymin)
+ ymin = temp;
+ }
+ }
+
+ TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
+ canv->SetGridy(1);
+
+ hSystematics[reference]->Draw("e p");
+ hSystematics[reference]->GetYaxis()->SetRangeUser(ymin, ymax);
+ for (Int_t j = 0; j < numHistos; j++) {
+ if (j == reference)
+ continue;
+
+ hSystematics[j]->SetMarkerStyle(20 + j);
+ hSystematics[j]->Draw("e p same");
+ }
+
+ TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
+ legend->SetBorderSize(0);
+ legend->SetFillColor(0);
+
+ for (Int_t j = 0; j < numHistos; j++) {
+ legend->AddEntry(hSystematics[j], Form("%s", systematicsHistosName[j].Data()), "p");
+ }
+ legend->Draw();
+
+
+ const Int_t nBins = histos[reference]->GetNbinsX();
+ Double_t x[nBins];
+ Double_t y[nBins];
+ Double_t xerr[nBins];
+ Double_t yerrl[nBins];
+ Double_t yerrh[nBins];
+
+ for (Int_t bin = 0; bin < nBins; bin++) {
+ x[bin] = histos[reference]->GetBinCenter(bin + 1);
+ xerr[bin] = histos[reference]->GetBinWidth(bin + 1) / 2.;
+ y[bin] = histos[reference]->GetBinContent(bin + 1);
+
+ // Take all points that are more than nSigma sigma away from 0.
+ // If there are at least 2 such points, take the difference between
+ // the extreme values (i.e. maximum and minimum) as a measure of
+ // the systematics
+ Int_t count = 0;
+ Double_t deltaMin = 0;
+ Double_t deltaMax = 0;
+
+ for (Int_t j = 0; j < numHistos; j++) {
+ if (hSystematics[j]->GetBinError(bin + 1) == 0) // Per definition always true for reference histo
+ continue;
+
+ Double_t delta = hSystematics[j]->GetBinContent(bin + 1);
+ if (TMath::Abs(delta / hSystematics[j]->GetBinError(bin + 1)) > nSigma) {
+ //if (count == 0) {
+ // deltaMin = delta;
+ // deltaMax = delta;
+ //}
+ //else {
+ if (delta < deltaMin)
+ deltaMin = delta;
+ if (delta > deltaMax)
+ deltaMax = delta;
+ //}
+ count++;
+ }
+ }
+
+ //if (deltaMax > 0.)
+ // yerrh[bin] = deltaMax;
+ //else
+ // yerrh[bin] = 0.;
+ //
+ //if (deltaMin < 0.)
+ // yerrl[bin] = -deltaMin;
+ //else
+ // yerrl[bin] = 0.;
+
+ if (count < 1) // Reference histo is not counted. One can only do systematics if there is at least one other histogram
+ yerrl[bin] = yerrh[bin] = 0.;
+ else
+ yerrl[bin] = yerrh[bin] = (deltaMax - deltaMin) / TMath::Sqrt(2);
+
+ }
+
+ TGraphAsymmErrors* gTemp = new TGraphAsymmErrors(nBins, x, y, xerr, xerr, yerrl, yerrh);
+ *gr = gTemp;
+ (*gr)->SetName(Form("systematicError_%s", AliPID::ParticleName(speciesID)));
+ (*gr)->SetLineColor(hSystematics[0]->GetMarkerColor());
+ //(*gr)->SetFillColor(kGray);
+ (*gr)->SetFillStyle(0);//3004 + reference);
+
+ return canv;
+}*/
+
+
+//________________________________________________________
+TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TH1F*** hist, Int_t reference,
+ TGraphAsymmErrors** gr)
+{
+ TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
+ canv->SetGridx(1);
+ canv->SetGridy(1);
+ canv->SetLogx(1);
+ for (Int_t i = 0; i < numSpecies; i++) {
+ hist[i][reference]->GetYaxis()->SetRangeUser(0.0, 1.0);
+ hist[i][reference]->GetYaxis()->SetTitle(canvTitle.Data());
+ hist[i][reference]->GetXaxis()->SetRangeUser(pLow, pHigh);
+ //hist[i][reference]->SetFillStyle(3004 + i);
+ //hist[i][reference]->SetFillColor(kGray);
+ hist[i][reference]->SetFillStyle(0);
+ hist[i][reference]->SetFillColor(hist[i][reference]->GetMarkerColor());
+ hist[i][reference]->SetLineColor(hist[i][reference]->GetMarkerColor());
+ }
+ hist[2][reference]->SetMarkerStyle(20);
+ hist[2][reference]->Draw("e p");
+ hist[0][reference]->SetMarkerStyle(21);
+ hist[0][reference]->Draw("e p same");
+ hist[1][reference]->SetMarkerStyle(22);
+ hist[1][reference]->Draw("e p same");
+ hist[3][reference]->SetMarkerStyle(29);
+ hist[3][reference]->Draw("e p same");
+ hist[4][reference]->SetMarkerStyle(30);
+ hist[4][reference]->Draw("e p same");
+
+ gr[0]->Draw("2 same");
+ gr[1]->Draw("2 same");
+ gr[2]->Draw("2 same");
+ gr[3]->Draw("2 same");
+ gr[4]->Draw("2 same");
+
+ TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
+ legend->SetBorderSize(0);
+ legend->SetFillColor(0);
+ legend->AddEntry(hist[2][reference], "#pi", "flp");
+ legend->AddEntry(hist[0][reference], "e", "flp");
+ legend->AddEntry(hist[1][reference], "K", "flp");
+ legend->AddEntry(hist[3][reference], "p", "flp");
+ legend->AddEntry(hist[4][reference], "#mu", "flp");
+ legend->Draw();
+
+ return canv;
+}
+
+
+//________________________________________________________
+TH1F* loadHisto(const TString histName, TFile* f)
+{
+ if (!f) {
+ std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;
+ return 0x0;
+ }
+
+ TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));
+ if (!hTemp) {
+ std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;
+ return 0x0;
+ }
+
+ return hTemp;
+}
+
+
+//________________________________________________________
+Int_t SystematicErrorEstimation(const TString path, const TString outFileTitle, const TString* fileNames, const TString* histTitles,
+ const Int_t numFiles, const Double_t nSigma, const Bool_t ignoreSigmaErrors)
+{
+ if (!fileNames || numFiles < 1)
+ return -1;
+
+ TFile* f[numFiles];
+ TH1F** hFractions[numSpecies];
+ for (Int_t i = 0; i < numSpecies; i++)
+ hFractions[i] = new TH1F*[numFiles];
+
+ const Int_t reference = 0;
+ TH1F* hYields[numSpecies]; // Only the reference yields
+
+ const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };
+
+ const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };
+
+ for (Int_t iFile = 0; iFile < numFiles; iFile++) {
+ f[iFile] = TFile::Open(fileNames[iFile].Data());
+ if (!f[iFile]) {
+ std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;
+ return -1;
+ }
+
+ // Extract the data histograms
+ for (Int_t i = 0; i < numSpecies; i++) {
+ hFractions[i][iFile] = loadHisto(histNames[i], f[iFile]);
+ if (!hFractions[i][iFile])
+ return -1;
+
+ if (iFile == reference) {
+ hYields[i] = loadHisto(histNamesYields[i], f[iFile]);
+ if (!hYields[i])
+ return -1;
+ }
+ }
+ }
+
+
+ TGraphAsymmErrors* grSysErrors[numSpecies] = {0x0,};
+ TGraphAsymmErrors* grSysErrorsYields[numSpecies] = {0x0,};
+
+ TH1F* hSystematicsPions[numFiles];
+ TCanvas* cSystematicsPions = calculateSystematics("cSystematicsPions", "Systematics Pions", hFractions[2], numFiles,
+ AliPID::kPion, nSigma,
+ histTitles, reference, hSystematicsPions, &grSysErrors[2], ignoreSigmaErrors);
+
+ TH1F* hSystematicsElectrons[numFiles];
+ TCanvas* cSystematicsElectrons = calculateSystematics("cSystematicsElectrons", "Systematics Electrons", hFractions[0], numFiles,
+ AliPID::kElectron,
+ nSigma, histTitles, reference, hSystematicsElectrons,
+ &grSysErrors[0], ignoreSigmaErrors);
+
+ TH1F* hSystematicsKaons[numFiles];
+ TCanvas* cSystematicsKaons = calculateSystematics("cSystematicsKaons", "Systematics Kaons", hFractions[1], numFiles, AliPID::kKaon, nSigma,
+ histTitles, reference, hSystematicsKaons, &grSysErrors[1], ignoreSigmaErrors);
+
+ TH1F* hSystematicsProtons[numFiles];
+ TCanvas* cSystematicsProtons = calculateSystematics("cSystematicsProtons", "Systematics Protons", hFractions[3], numFiles,
+ AliPID::kProton, nSigma,
+ histTitles, reference, hSystematicsProtons, &grSysErrors[3], ignoreSigmaErrors);
+
+ TH1F* hSystematicsMuons[numFiles];
+ TCanvas* cSystematicsMuons = calculateSystematics("cSystematicsMuons", "Systematics Muons", hFractions[4], numFiles,
+ AliPID::kMuon, nSigma,
+ histTitles, reference, hSystematicsMuons, &grSysErrors[4], ignoreSigmaErrors);
+
+ Double_t pLow = 0.15;
+ Double_t pHigh = 50.;
+ TCanvas* cFractionsWithSystematicError = DrawFractionHistos("cFractionsWithSystematicError", "Particle fractions", pLow, pHigh, hFractions, reference,
+ grSysErrors);
+
+
+ //TODO At the moment, the error of the fractions and the yield is just a constant factor (number of tracks in that bin)
+ // (-> But this can change in future (I have to think about it a little bit more carefully)).
+ // Thus, the relative errors are the same for fractions and yields and I can just use this fact to
+ // transform the errors from the fractions to those of the yields.
+ // However, this causes trouble in case of fraction = 0. Therefore, sum up the yields to the total yield and use this for scaling
+ for (Int_t i = 0; i < numSpecies; i++) {
+ grSysErrorsYields[i] = new TGraphAsymmErrors(*grSysErrors[i]);
+ TString name = grSysErrors[i]->GetName();
+ name.ReplaceAll("systematicError_", "systematicErrorYields_");
+ grSysErrorsYields[i]->SetName(name.Data());
+
+ for (Int_t ind = 0; ind < grSysErrorsYields[i]->GetN(); ind++) {
+ Double_t totalYield = 0;
+ for (Int_t j = 0; j < numSpecies; j++)
+ totalYield += hYields[j]->GetBinContent(ind + 1);
+
+ const Double_t yield = hYields[i]->GetBinContent(ind + 1);
+ const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);
+ const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);
+
+ grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);
+ grSysErrorsYields[i]->SetPointEYhigh(ind, totalYield * sysErrorHigh);
+ grSysErrorsYields[i]->SetPointEYlow(ind, totalYield * sysErrorLow);
+
+ /*
+ Double_t totalYield = 0;
+ for (Int_t j = 0; j < numSpecies; j++)
+ totalYield += hYields[j]->GetBinContent(ind + 1);
+
+ const Double_t yield = hYields[i]->GetBinContent(ind + 1);
+ const Double_t fraction = hFractions[i][reference]->GetBinContent(ind + 1);
+ const Double_t sysErrorLow = grSysErrors[i]->GetErrorYlow(ind);
+ const Double_t sysErrorHigh = grSysErrors[i]->GetErrorYhigh(ind);
+
+ if (fraction <= 0.) {
+ printf("Error: Fraction = 0 for species %d. Cannot transform error....\n", i);
+ return -1;
+ }
+ const Double_t sysErrorLowRel = sysErrorLow / fraction;
+ const Double_t sysErrorHighRel = sysErrorHigh / fraction;
+
+ grSysErrorsYields[i]->SetPoint(ind, grSysErrorsYields[i]->GetX()[ind], yield);
+ grSysErrorsYields[i]->SetPointEYhigh(ind, sysErrorHighRel * yield);
+ grSysErrorsYields[i]->SetPointEYlow(ind, sysErrorLowRel * yield);
+ */
+ }
+ }
+
+
+ // Output file
+ TFile* fSave = 0x0;
+ TDatime daTime;
+ TString saveFileName;
+
+ saveFileName = Form("outputSystematics_%s_nSigma%.1f__%04d_%02d_%02d.root", outFileTitle.Data(), nSigma, daTime.GetYear(),
+ daTime.GetMonth(), daTime.GetDay());
+
+ fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");
+ if (!fSave) {
+ std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;
+ return -1;
+ }
+
+ // Save final results
+ fSave->cd();
+
+ for (Int_t i = 0; i < numFiles; i++) {
+ if (hSystematicsElectrons[i])
+ hSystematicsElectrons[i]->Write();
+
+ if (hSystematicsPions[i])
+ hSystematicsPions[i]->Write();
+
+ if (hSystematicsKaons[i])
+ hSystematicsKaons[i]->Write();
+
+ if (hSystematicsProtons[i])
+ hSystematicsProtons[i]->Write();
+ if (hSystematicsMuons[i])
+ hSystematicsMuons[i]->Write();
+ }
+
+ if (cSystematicsElectrons)
+ cSystematicsElectrons->Write();
+
+ if (cSystematicsPions)
+ cSystematicsPions->Write();
+
+ if (cSystematicsKaons)
+ cSystematicsKaons->Write();
+
+ if (cSystematicsProtons)
+ cSystematicsProtons->Write();
+
+ if (cSystematicsMuons)
+ cSystematicsMuons->Write();
+
+ if (cFractionsWithSystematicError)
+ cFractionsWithSystematicError->Write();
+
+ for (Int_t i = 0; i < numSpecies; i++) {
+ if (hFractions[i][reference])
+ hFractions[i][reference]->Write();
+
+ if (grSysErrors[i])
+ grSysErrors[i]->Write();
+
+ if (hYields[i])
+ hYields[i]->Write();
+
+ if (grSysErrorsYields[i])
+ grSysErrorsYields[i]->Write();
+ }
+
+ // Save list of file names in output file
+ TString listOfFileNames = "";
+ for (Int_t i = 0; i < numFiles; i++) {
+ listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));
+ }
+
+ TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),
+ Form("Used files for systematics: %s\n", listOfFileNames.Data()));
+ settings->Write();
+
+ fSave->Close();
+
+ delete cSystematicsElectrons;
+ delete cSystematicsPions;
+ delete cSystematicsKaons;
+ delete cSystematicsMuons;
+ delete cSystematicsProtons;
+ delete cFractionsWithSystematicError;
+
+ return 0;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff