]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/macros/PID/SystematicErrorEstimation.C
Initial Tasks and macros for PIDed Fragmentation Function (Benjamin Hess)
[u/mrichter/AliRoot.git] / PWGJE / macros / PID / SystematicErrorEstimation.C
diff --git a/PWGJE/macros/PID/SystematicErrorEstimation.C b/PWGJE/macros/PID/SystematicErrorEstimation.C
new file mode 100644 (file)
index 0000000..0585f56
--- /dev/null
@@ -0,0 +1,517 @@
+#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