PWGJE/macros/PID/SystematicErrorEstimation.C

   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
  20 const Int_t numSpecies = 5;
  21
  22 //________________________________________________________
  23 TCanvas* 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 //________________________________________________________
 121 TCanvas* 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 //________________________________________________________
 242 TCanvas* 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 //________________________________________________________
 291 TH1F* 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 //________________________________________________________
 309 Int_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 }