PWGPP/MUON/lite/PlotEfficiency.C

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * All rights reserved.                                                   *
   4  *                                                                        *
   5  * Primary Authors:                                                       *
   6  *   Artur Szostak <artursz@iafrica.com>                                  *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /// \ingroup macros
  18 /// \file PlotEfficiency.C
  19 /// \brief Generates efficiency plots from the histograms generated by AliAnalysisTaskLinkToMC.
  20 /// \author Artur Szostak <artursz@iafrica.com>
  21 ///
  22 /// This macro is used to generate efficiency plots, fake track ratio plots
  23 /// and also calculate the total integrated efficiency and fake track ratio.
  24 /// The histograms generated by the AliAnalysisTaskLinkToMC analysis task is
  25 /// used as input. The macro can be run with root as follows:
  26 ///
  27 ///  $ root PlotEfficiency.C\(\"hists.root\"\)
  28 ///
  29 /// where hists.root is the name of the output file containing the histograms
  30 /// generated by the analysis task. Note that the '\' character before the '('
  31 /// and '"' characters is required. Alternatively run the macros as follows from
  32 /// the root command prompt:
  33 ///
  34 ///  $ root
  35 ///  .L PlotEfficiency.C
  36 ///  PlotEfficiency("hists.root")
  37
  38
  39 #if !defined(__CINT__) || defined(__MAKECINT__)
  40 #include "Riostream.h"
  41 #include "TH1.h"
  42 #include "TH2.h"
  43 #include "TFile.h"
  44 #include "TList.h"
  45 #include "TGraphAsymmErrors.h"
  46 #include "TCanvas.h"
  47 #endif
  48
  49
  50 void EfficiencyIntegral(const char* prefix, TH1* pass, TH1* total, double& eff, double& error, bool includeOverflow = false)
  51 {
  52         /// This routine calculates the integrated efficiency for a 1D histogram.
  53         /// [in] \param prefix  The string to print before the calculated efficiency result.
  54         /// [in] \param pass  The histogram of reconstructed tracks.
  55         /// [in] \param total  The histogram of findable (reconstructable) tracks.
  56         /// [out] \param eff  Filled by the total integrated efficiency calculated.
  57         /// [out] \param error  Filled by the uncertainty in the efficiency value.
  58         /// [in] \param includeOverflow  If set then the overflow bins are also used
  59         ///     in the calculation.
  60
  61         TH1D passhist("temppasshist", "", 1, 0, 1);
  62         TH1D totalhist("temptotalhist", "", 1, 0, 1);
  63         if (includeOverflow)
  64         {
  65                 passhist.SetBinContent(1, pass->Integral(0, pass->GetNbinsX()+1));
  66                 totalhist.SetBinContent(1, total->Integral(0, total->GetNbinsX()+1));
  67         }
  68         else
  69         {
  70                 passhist.SetBinContent(1, pass->Integral(1, pass->GetNbinsX()));
  71                 totalhist.SetBinContent(1, total->Integral(1, total->GetNbinsX()));
  72         }
  73         TGraphAsymmErrors* effgraph = new TGraphAsymmErrors(&passhist, &totalhist);
  74         eff = effgraph->GetY()[0];
  75         double effhigh = effgraph->GetEYhigh()[0];
  76         double efflow = effgraph->GetEYlow()[0];
  77         error = effhigh > efflow ? effhigh : efflow;
  78         cout << prefix << eff << " + " << effhigh << " - " << efflow << endl;
  79 }
  80
  81
  82 void PlotEfficiency(const char* histfilename = "hists.root", bool includeOverflow = true)
  83 {
  84         /// Opens the file containing the histograms generated by the AliAnalysisTaskLinkToMC
  85         /// analysis task and tries to find the first TList in the file.
  86         /// From the TList we try to find TH2 classes with the following names:
  87         /// "findableTracksHist", "foundTracksHistMC" and "fakeTracksHist"
  88         /// From these we plot the efficiency plots for the pT and rapidity dimension
  89         /// by projecting the 2D histograms. Also, the fake track ratio histograms are
  90         /// created. Finally the total integrated efficiency and fake track ratio is
  91         /// calculated and printed.
  92         /// \param histfilename  The name of the root file containing the histograms
  93         ///    generated by the AliAnalysisTaskLinkToMC analysis task.
  94         /// \param includeOverflow  Indicates if the overflow bins should be used for
  95         ///    for the total integral efficiency and fake track ratio calculations.
  96         ///    The default is not to use the overflow bins.
  97
  98         TFile file(histfilename);
  99
 100         TList* histolist = NULL;
 101         TIter next(file.GetListOfKeys());
 102         for (TObject* key = next(); key != NULL; key = next())
 103         {
 104                 TList* list = dynamic_cast<TList*>( file.Get(key->GetName()) );
 105                 if (list != NULL)
 106                 {
 107                         histolist = list;
 108                         break;
 109                 }
 110         }
 111
 112         if (histolist == NULL)
 113         {
 114                 cout << "ERROR: Could not find the histogram list in file '" << histfilename << "'" << endl;
 115                 return;
 116         }
 117
 118         TH2* findableTracksHist = dynamic_cast<TH2*>( histolist->FindObject("findableTracksHist") );
 119         TH2* foundTracksHistMC = dynamic_cast<TH2*>( histolist->FindObject("foundTracksHistMC") );
 120         TH2* fakeTracksHist = dynamic_cast<TH2*>( histolist->FindObject("fakeTracksHist") );
 121
 122         if (findableTracksHist == NULL)
 123         {
 124                 cout << "ERROR: Could not find the histogram findableTracksHist in TList '" << histolist->GetName() << "'" << endl;
 125                 return;
 126         }
 127         if (foundTracksHistMC == NULL)
 128         {
 129                 cout << "ERROR: Could not find the histogram foundTracksHistMC in TList '" << histolist->GetName() << "'" << endl;
 130                 return;
 131         }
 132         if (fakeTracksHist == NULL)
 133         {
 134                 cout << "ERROR: Could not find the histogram fakeTracksHist in TList '" << histolist->GetName() << "'" << endl;
 135                 return;
 136         }
 137
 138         TH1* projeffx = foundTracksHistMC->ProjectionX();
 139         TH1* projfindx = findableTracksHist->ProjectionX();
 140         TGraphAsymmErrors* effpt = new TGraphAsymmErrors(projeffx, projfindx);
 141         effpt->Draw("ap");
 142         effpt->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
 143         effpt->GetHistogram()->SetXTitle(projeffx->GetXaxis()->GetTitle());
 144         effpt->GetHistogram()->SetYTitle("#epsilon");
 145         effpt->GetYaxis()->SetRangeUser(0, 1.1);
 146
 147         new TCanvas;
 148         TH1* projfakex = fakeTracksHist->ProjectionX();
 149         TGraphAsymmErrors* fakept = new TGraphAsymmErrors(projfakex, projfindx);
 150         fakept->Draw("ap");
 151         fakept->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
 152         fakept->GetHistogram()->SetXTitle(projfakex->GetXaxis()->GetTitle());
 153         fakept->GetHistogram()->SetYTitle("fake tracks / findable tracks");
 154
 155         new TCanvas;
 156         TH1* projeffy = foundTracksHistMC->ProjectionY();
 157         TH1* projfindy = findableTracksHist->ProjectionY();
 158         TGraphAsymmErrors* effy = new TGraphAsymmErrors(projeffy, projfindy);
 159         effy->Draw("ap");
 160         effy->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
 161         effy->GetHistogram()->SetXTitle(projeffy->GetXaxis()->GetTitle());
 162         effy->GetHistogram()->SetYTitle("#epsilon");
 163         effy->GetYaxis()->SetRangeUser(0, 1.1);
 164
 165         new TCanvas;
 166         TH1* projfakey = fakeTracksHist->ProjectionY();
 167         TGraphAsymmErrors* fakey = new TGraphAsymmErrors(projfakey, projfindy);
 168         fakey->Draw("ap");
 169         fakey->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
 170         fakey->GetHistogram()->SetXTitle(projfakey->GetXaxis()->GetTitle());
 171         fakey->GetHistogram()->SetYTitle("fake tracks / findable tracks");
 172
 173         if (includeOverflow)
 174         {
 175                 cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
 176                 cout << "found tracks    = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
 177                 cout << "fake tracks     = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
 178         }
 179         else
 180         {
 181                 cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
 182                 cout << "found tracks    = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
 183                 cout << "fake tracks     = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
 184         }
 185         double eff, fake, deff, dfake;
 186         EfficiencyIntegral("Total efficiency = ", projeffx, projfindx, eff, deff, includeOverflow);
 187         EfficiencyIntegral("Total fake ratio = ", projfakex, projfindx, fake, dfake, includeOverflow);
 188         cout << "Sum = " << eff + fake << " +/- " << sqrt(deff*deff + dfake*dfake) << endl;
 189 }
 190