PWGJE/macros/PID/AddUpSystematicErrors.C

   1 #include "TH1D.h"\r
   2 #include "TCanvas.h"\r
   3 #include "TStyle.h"\r
   4 #include "TString.h"\r
   5 #include "TLegend.h"\r
   6 #include "TFile.h"\r
   7 #include "TGraphErrors.h"\r
   8 #include "TGraphAsymmErrors.h"\r
   9 #include "TGraph.h"\r
  10 #include "TMath.h"\r
  11 \r
  12 #include "AliPID.h"\r
  13 \r
  14 #include "THnSparseDefinitions.h"\r
  15 \r
  16 #include <iostream>\r
  17 #include <iomanip>\r
  18 \r
  19 const Int_t numSpecies = 5;\r
  20 \r
  21 //________________________________________________________\r
  22 TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TGraphAsymmErrors** gr,  TH1F** histRef)\r
  23 {\r
  24   TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);\r
  25   canv->SetGridx(1);\r
  26   canv->SetGridy(1);\r
  27   canv->SetLogx(1);\r
  28   \r
  29   for (Int_t i = 0; i < numSpecies; i++) {\r
  30     histRef[i]->GetYaxis()->SetRangeUser(0.0, 1.0);\r
  31     histRef[i]->GetYaxis()->SetTitle(canvTitle.Data());\r
  32     histRef[i]->GetXaxis()->SetRangeUser(pLow, pHigh);\r
  33     //histRef[i]->SetFillStyle(3004 + i);\r
  34     //histRef[i]->SetFillColor(kGray);\r
  35     histRef[i]->SetFillStyle(0);\r
  36     histRef[i]->SetFillColor(histRef[i]->GetMarkerColor());\r
  37     histRef[i]->SetLineColor(histRef[i]->GetMarkerColor());\r
  38   }\r
  39   histRef[2]->SetMarkerStyle(20);\r
  40   histRef[2]->Draw("e p");\r
  41   histRef[0]->SetMarkerStyle(21);\r
  42   histRef[0]->Draw("e p same");\r
  43   histRef[1]->SetMarkerStyle(22);\r
  44   histRef[1]->Draw("e p same");\r
  45   histRef[3]->SetMarkerStyle(29);\r
  46   histRef[3]->Draw("e p same");\r
  47   histRef[4]->SetMarkerStyle(30);\r
  48   histRef[4]->Draw("e p same");\r
  49   \r
  50   gr[0]->GetHistogram()->GetXaxis()->SetRangeUser(pLow, pHigh);\r
  51   gr[0]->GetHistogram()->GetYaxis()->SetRangeUser(0., 1.0);\r
  52   gr[0]->Draw("2 same");\r
  53   gr[1]->Draw("2 same");\r
  54   gr[2]->Draw("2 same");\r
  55   gr[3]->Draw("2 same");\r
  56   gr[4]->Draw("2 same");\r
  57   \r
  58   TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);    \r
  59   legend->SetBorderSize(0);\r
  60   legend->SetFillColor(0);\r
  61   legend->AddEntry(histRef[2], "#pi", "flp");\r
  62   legend->AddEntry(histRef[0], "e", "flp");\r
  63   legend->AddEntry(histRef[1], "K", "flp");\r
  64   legend->AddEntry(histRef[3], "p", "flp");\r
  65   legend->AddEntry(histRef[4], "#mu", "flp");\r
  66   legend->Draw();\r
  67   \r
  68   ClearTitleFromHistoInCanvas(canv);\r
  69   \r
  70   return canv;\r
  71 }\r
  72 \r
  73 \r
  74 //________________________________________________________\r
  75 TGraphAsymmErrors* loadGraph(const TString graphName, TFile* f)\r
  76 {\r
  77   if (!f) {\r
  78     std::cout << "No file. Cannot load graph \"" << graphName.Data() << "\n!" << std::endl;\r
  79     return 0x0;\r
  80   }\r
  81   \r
  82   TGraphAsymmErrors* grTemp = dynamic_cast<TGraphAsymmErrors*>(f->Get(graphName.Data()));\r
  83   if (!grTemp) {\r
  84     std::cout << "Failed to load histo \"" << graphName.Data() << "\"!" << std::endl;\r
  85     return 0x0;\r
  86   } \r
  87   \r
  88   return grTemp;\r
  89 }\r
  90 \r
  91 \r
  92 //________________________________________________________\r
  93 TH1F* loadHisto(const TString histName, TFile* f)\r
  94 {\r
  95   if (!f) {\r
  96     std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;\r
  97     return 0x0;\r
  98   }\r
  99   \r
 100   TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));\r
 101   if (!hTemp) {\r
 102     std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;\r
 103     return 0x0;\r
 104   } \r
 105   \r
 106   return hTemp;\r
 107 }\r
 108 \r
 109 \r
 110 //________________________________________________________\r
 111 Int_t AddUpSystematicErrors(const TString path, const TString outFileTitle, const TString* fileNames, const Int_t numFiles,\r
 112                             const TString fileNameReference) \r
 113 { \r
 114   if (!fileNames || numFiles < 1)\r
 115     return -1;\r
 116   \r
 117   TFile* f[numFiles];\r
 118   TGraphAsymmErrors** grSysError[numSpecies];\r
 119   TString grNames[numSpecies] = {"systematicError_electron", "systematicError_kaon", "systematicError_pion", "systematicError_proton",\r
 120                                  "systematicError_muon" };\r
 121   \r
 122   const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };\r
 123   \r
 124   \r
 125   TGraphAsymmErrors** grSysErrorYields[numSpecies];\r
 126   TString grNamesYields[numSpecies] = {"systematicErrorYields_electron", "systematicErrorYields_kaon", "systematicErrorYields_pion",\r
 127                                        "systematicErrorYields_proton", "systematicErrorYields_muon" };\r
 128                                  \r
 129   const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };\r
 130   \r
 131   for (Int_t i = 0; i < numSpecies; i++) {\r
 132     grSysError[i] = new TGraphAsymmErrors*[numFiles];\r
 133     grSysErrorYields[i] = new TGraphAsymmErrors*[numFiles];\r
 134   }\r
 135     \r
 136   for (Int_t iFile = 0; iFile < numFiles; iFile++) {\r
 137     f[iFile] = TFile::Open(fileNames[iFile].Data());\r
 138     if (!f[iFile])  {\r
 139       std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;\r
 140       return -1;\r
 141     }\r
 142         \r
 143     // Extract the data graphs\r
 144     for (Int_t i = 0; i < numSpecies; i++) {\r
 145       grSysError[i][iFile] = loadGraph(grNames[i], f[iFile]);\r
 146       if (!grSysError[i][iFile])\r
 147         return -1;\r
 148       \r
 149       grSysError[i][iFile]->SetName(Form("%s_fileID%d", grSysError[i][iFile]->GetName(), iFile));\r
 150       \r
 151       \r
 152       grSysErrorYields[i][iFile] = loadGraph(grNamesYields[i], f[iFile]);\r
 153       if (!grSysErrorYields[i][iFile])\r
 154         return -1;\r
 155       \r
 156       grSysErrorYields[i][iFile]->SetName(Form("%s_fileID%d", grSysErrorYields[i][iFile]->GetName(), iFile));\r
 157     }\r
 158   }\r
 159   \r
 160   // Load data points with statistical errors\r
 161   TFile* fReferenceData = TFile::Open(fileNameReference.Data());\r
 162   if (!fReferenceData) {\r
 163     std::cout << "Failed to open file \"" << fileNameReference.Data() << "\"!" << std::endl;\r
 164     return -1;\r
 165   }\r
 166     \r
 167   TH1F* hReferenceFractions[numSpecies];\r
 168   TH1F* hReferenceYields[numSpecies];\r
 169   \r
 170   for (Int_t i = 0; i < numSpecies; i++) {\r
 171     hReferenceFractions[i] = loadHisto(histNames[i], fReferenceData);\r
 172     if (!hReferenceFractions[i])\r
 173       return -1;\r
 174     \r
 175     hReferenceYields[i] = loadHisto(histNamesYields[i], fReferenceData);\r
 176     if (!hReferenceYields[i])\r
 177       return -1;\r
 178   }\r
 179   \r
 180   const Int_t reference = 0;\r
 181   \r
 182   // The x,y coordinates should be those of the reference graph. Then, all corresponding sys. errors are added in quadrature.\r
 183   TGraphAsymmErrors* grTotSysError[numSpecies];\r
 184   TGraphAsymmErrors* grTotSysErrorYields[numSpecies];\r
 185   Double_t sysErrorTotalSquared = 0;\r
 186   Double_t temp = 0;\r
 187   \r
 188   for (Int_t i = 0; i < numSpecies; i++) {\r
 189     grTotSysError[i] = new TGraphAsymmErrors(*grSysError[i][reference]);\r
 190     grTotSysError[i]->SetName(grNames[i]);\r
 191     \r
 192     for (Int_t iPoint = 0; iPoint < grTotSysError[i]->GetN(); iPoint++) {\r
 193       sysErrorTotalSquared = 0;\r
 194       for (Int_t iFile = 0; iFile < numFiles; iFile++) {\r
 195         // Already averages high and low value -> Since they are by now the same, this is ok.\r
 196         temp = grSysError[i][iFile]->GetErrorY(iPoint); \r
 197         if (temp > 0) {\r
 198           sysErrorTotalSquared += temp * temp;\r
 199         }\r
 200       }\r
 201       \r
 202       grTotSysError[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));\r
 203       grTotSysError[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));\r
 204     }\r
 205     \r
 206     // Same for the yields\r
 207     grTotSysErrorYields[i] = new TGraphAsymmErrors(*grSysErrorYields[i][reference]);\r
 208     grTotSysErrorYields[i]->SetName(grNamesYields[i]);\r
 209     \r
 210     for (Int_t iPoint = 0; iPoint < grTotSysErrorYields[i]->GetN(); iPoint++) {\r
 211       sysErrorTotalSquared = 0;\r
 212       for (Int_t iFile = 0; iFile < numFiles; iFile++) {\r
 213         // Already averages high and low value -> Since they are by now the same, this is ok.\r
 214         temp = grSysErrorYields[i][iFile]->GetErrorY(iPoint); \r
 215         if (temp > 0) {\r
 216           sysErrorTotalSquared += temp * temp;\r
 217         }\r
 218       }\r
 219       \r
 220       grTotSysErrorYields[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));\r
 221       grTotSysErrorYields[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));\r
 222     }\r
 223   }\r
 224   \r
 225   const Double_t pLow = 0.15;\r
 226   const Double_t pHigh = 50.;\r
 227   TCanvas* cFractionsWithTotalSystematicError = DrawFractionHistos("cFractionsWithTotalSystematicError", "Particle fractions", pLow, pHigh,\r
 228                                                                    grTotSysError, hReferenceFractions);\r
 229   \r
 230     \r
 231   // Output file\r
 232   TFile* fSave = 0x0;\r
 233   TDatime daTime;\r
 234   TString saveFileName;\r
 235     \r
 236   saveFileName = Form("outputSystematicsTotal_%s__%04d_%02d_%02d.root", outFileTitle.Data(), daTime.GetYear(),\r
 237                       daTime.GetMonth(), daTime.GetDay());\r
 238     \r
 239   fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");\r
 240   if (!fSave) {\r
 241     std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;\r
 242     return -1;\r
 243   }\r
 244   \r
 245   // Save final results\r
 246   fSave->cd();\r
 247   \r
 248   if (cFractionsWithTotalSystematicError)\r
 249       cFractionsWithTotalSystematicError->Write();\r
 250     \r
 251     for (Int_t i = 0; i < numSpecies; i++) {\r
 252       if (grTotSysError[i])\r
 253         grTotSysError[i]->Write();\r
 254       if (hReferenceFractions[i])\r
 255         hReferenceFractions[i]->Write();\r
 256       \r
 257       if (grTotSysErrorYields[i])\r
 258         grTotSysErrorYields[i]->Write();\r
 259       if (hReferenceYields[i])\r
 260         hReferenceYields[i]->Write();\r
 261   }\r
 262   \r
 263   // Save list of file names in output file\r
 264   TString listOfFileNames = "";\r
 265   for (Int_t i = 0; i < numFiles; i++) {\r
 266     listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));\r
 267   }\r
 268   \r
 269   TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),\r
 270                                 Form("Used files for systematics: %s\n", listOfFileNames.Data()));\r
 271   settings->Write();\r
 272   \r
 273   delete cFractionsWithTotalSystematicError;\r
 274     \r
 275   fSave->Close();\r
 276      \r
 277   return 0;\r
 278 }\r