PWGPP/EvTrkSelection/macros/MakeSensitivityPlots.C

   1 #include "TH1D.h"
   2 #include "TCanvas.h"
   3 #include "TFile.h"
   4 #include "THn.h"
   5 #include "TList.h"
   6 #include "TStyle.h"
   7 #include "TLegend.h"
   8 #include "TLatex.h"
   9
  10 //loadlibs();
  11 //gSystem->Load("libANALYSIS");
  12 //gSystem->Load("libANALYSISalice");
  13 //#include "AliESDtrackCuts.h"
  14
  15 void MakeSensitivityPlots2();
  16 //
  17 TCanvas * GetSensitivityPlot(TString cutname = "Ncl",
  18                              Int_t projectionAxis = 1,
  19                              Int_t particleType = 5,
  20                              Double_t lowCut = 80.,
  21                              Double_t highCut = 1e4,
  22                              TString inFileNameData = "data_LHC10d.root",
  23                              TString inFileNameMC = "mc_LHC10d.root");
  24 //
  25 TH1D   * GetAcceptedFraction(TString cutname = "Ncl",
  26                              Int_t projectionAxis = 1,
  27                              Int_t particleType = 5,
  28                              Double_t lowCut = 80.,
  29                              Double_t highCut = 1e4,
  30                              TString inFileName = "data_LHC10d.root");
  31 //
  32 void PrintCutGallery(TString inFileNameData = "data_LHC10d.root",
  33                      TString inFileNameMC   = "mc_LHC10d.root",
  34                      TString outFileName    = "plots/10d.pdf");
  35
  36 //
  37 TH1D *GetITSTPCMatchingEff(TString inFileNameData = "data_LHC10d.root",
  38                               Int_t particleType = 5,
  39                               Int_t projectionAxis = 1);
  40
  41 //
  42 TH1D *GetITSTPCMatchingHisto(TString inFileNameData = "data_LHC10d.root",
  43                              Int_t particleType = 5,
  44                              Bool_t isMatched = kTRUE,
  45                              Int_t projectionAxis = 1);
  46
  47
  48 //______________________________________________________________________________
  49 void MakeSensitivityPlots2() {
  50   //
  51   // make all the senstivity plots
  52   //
  53   //PrintCutGallery("output/LHC10b_data.root", "output/LHC10b_MC.root", "plots/10b.pdf");
  54   PrintCutGallery("data_LHC10d.root", "mc_LHC10d.root", "plots/10d.pdf");
  55   //PrintCutGallery("output/LHC10e_data.root", "output/LHC10e_MC.root", "plots/10e.pdf");
  56   //PrintCutGallery("output/LHC10h_data.root", "output/LHC10h_MC.root", "plots/10h.pdf");
  57
  58 }
  59
  60 //______________________________________________________________________________
  61 void PrintCutGallery(TString inFileNameData, TString inFileNameMC, TString outFileName) {
  62   //
  63   // print a gallery of cuts
  64   //
  65   TCanvas * canvMaster = new TCanvas("canvMaster","canvMaster");
  66   canvMaster->Divide(3,2);
  67   //TCanvas * canv = GetSensitivityPlot();
  68   //
  69   TCanvas * canvSens[6];
  70   for(Int_t i=0; i < 6; i++) {
  71     canvSens[i]= GetSensitivityPlot("Ncl", 1, 5, 50 + i*10, 1e4, inFileNameData.Data(),  inFileNameMC.Data());
  72     //canvSens[i]= GetSensitivityPlot("Chi2Tpc", 1, 5, -1e4, 2.5 + i*0.5, inFileNameData.Data(),  inFileNameMC.Data());
  73     canvSens[i]->SetName(Form("canvSens_%i",i));
  74     canvMaster->cd(i+1);
  75     canvSens[i]->DrawClonePad();
  76
  77   }
  78   canvMaster->Print(outFileName.Data());
  79
  80 }
  81
  82
  83 //______________________________________________________________________________
  84 TCanvas * GetSensitivityPlot(TString cutname,
  85                              Int_t projectionAxis,
  86                              Int_t particleType,
  87                              Double_t lowCut,
  88                              Double_t highCut,
  89                              TString inFileNameData,
  90                              TString inFileNameMC) {
  91   //
  92   // make a single plot
  93   //
  94   TH1D * nclAcceptedData = GetAcceptedFraction(cutname.Data(), projectionAxis, particleType, lowCut, highCut, inFileNameData);
  95   nclAcceptedData->SetNameTitle(Form("%s_DATA",nclAcceptedData->GetName()), Form("%s_DATA",nclAcceptedData->GetName()));
  96   //
  97   //
  98   nclAcceptedData->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
  99   nclAcceptedData->GetYaxis()->SetTitle("accepted fraction");
 100   nclAcceptedData->GetYaxis()->SetLabelSize(0.04);
 101   nclAcceptedData->GetYaxis()->SetTitleSize(0.06);
 102   nclAcceptedData->GetYaxis()->SetTitleOffset(1);
 103   nclAcceptedData->SetMaximum(1.07);
 104   nclAcceptedData->SetMinimum(0.625);
 105   nclAcceptedData->SetLineColor(kRed -3);
 106   //
 107   TH1D * nclAcceptedMc   = GetAcceptedFraction(cutname.Data(), projectionAxis, particleType, lowCut, highCut, inFileNameMC);
 108   nclAcceptedMc->SetNameTitle(Form("%s_MC",nclAcceptedMc->GetName()), Form("%s_MC",nclAcceptedMc->GetName()));
 109   nclAcceptedMc->SetLineColor(kBlue -3);
 110   //
 111   // create the complicated splitted canvas
 112   //
 113   gStyle->SetOptStat(0);
 114   gStyle->SetOptTitle(0);
 115   TCanvas * canvCut = new TCanvas("canvCut",Form("sensitivity to %s cut",cutname.Data()),600,600);
 116   canvCut->Divide(1,2);
 117   //
 118   TPad* canvas_up = (TPad*) canvCut->GetListOfPrimitives()->FindObject(Form("canvCut_1"));
 119   TPad* canvas_dw = (TPad*) canvCut->GetListOfPrimitives()->FindObject(Form("canvCut_2"));
 120   //
 121   canvas_up->SetLogx();
 122   canvas_dw->SetLogx();
 123   //
 124   // define the size
 125   double up_height     = 0.75; // please tune so that the upper figures size will meet your requirement
 126   double dw_correction = 1.30; // please tune so that the smaller canvas size will work in your environment
 127   //
 128   //double font_size_dw  = 0.1; // please tune the font size parameter for bottom figure
 129   double dw_height    = (1. - up_height) * dw_correction;
 130   // set pad size
 131   canvas_up->SetPad(0., 1 - up_height, 1., 1.);
 132   canvas_dw->SetPad(0., 0., 1., dw_height);
 133   canvas_up->SetFrameFillColor(0);
 134   canvas_up->SetFillColor(0);
 135   canvas_dw->SetFillColor(0);
 136   canvas_dw->SetFrameFillColor(0);
 137   //
 138   // set top margin 0 for bottom figure
 139   canvas_dw->SetTopMargin(0);
 140   canvas_dw->SetLeftMargin(0.12);
 141   canvas_dw->SetBottomMargin(0.3);
 142   //
 143   //set margins for top figure
 144   canvas_up->SetLeftMargin(0.12);
 145   //
 146   canvas_up->cd();
 147   nclAcceptedData->GetYaxis()->SetLabelFont(62);
 148   nclAcceptedData->SetObjectStat(0);
 149   nclAcceptedMc->SetObjectStat(0);
 150   //
 151   nclAcceptedData->DrawCopy();
 152   nclAcceptedMc->DrawCopy("SAME");
 153   //
 154   TLegend *leg = new TLegend(0.65,0.2,.85,0.4);
 155   leg->AddEntry(nclAcceptedData,"Data","f");
 156   leg->AddEntry(nclAcceptedMc,"MC","f");
 157   leg->SetBorderSize(0);
 158   gStyle->SetFillColor(0);
 159   leg->Draw();
 160   Float_t lowCutTitle = lowCut < -999. ? 0 : lowCut;
 161   TLatex *   tex = new TLatex(0.9654768,1.032485,Form("%4.2f < %s < %4.2f", lowCutTitle, cutname.Data(), highCut));
 162   tex->Draw();
 163   //
 164   // Draw the ratio
 165   //
 166   TH1D * nclAcceptedMcDataRatio = (TH1D*)nclAcceptedData->Clone();
 167   //
 168   canvCut->cd(2)->SetLogx();
 169   // gPad->SetTicky(2);
 170   gPad->SetFillStyle(0);
 171   //
 172   nclAcceptedMcDataRatio->Divide(nclAcceptedMc);
 173   nclAcceptedMcDataRatio->GetYaxis()->SetTitle("ratio");
 174   nclAcceptedMcDataRatio->GetXaxis()->SetLabelSize(0.1);
 175   nclAcceptedMcDataRatio->GetXaxis()->SetTitleSize(0.11);
 176   nclAcceptedMcDataRatio->GetYaxis()->SetLabelSize(0.07);
 177   nclAcceptedMcDataRatio->GetYaxis()->SetTitleSize(0.11);
 178   nclAcceptedMcDataRatio->GetXaxis()->SetTitleOffset(1.23);
 179   nclAcceptedMcDataRatio->GetYaxis()->SetTitleOffset(0.5);
 180   nclAcceptedMcDataRatio->GetYaxis()->CenterTitle();
 181   nclAcceptedMcDataRatio->SetMaximum(1.18);
 182   nclAcceptedMcDataRatio->SetMinimum(0.8);
 183   nclAcceptedMcDataRatio->GetYaxis()->SetLabelFont(62);
 184   nclAcceptedMcDataRatio->SetLineColor(kGreen+2);
 185   nclAcceptedMcDataRatio->Draw();
 186   canvas_dw->SetFillStyle(1001);
 187   //
 188   //
 189   //
 190   return canvCut;
 191
 192 }
 193
 194
 195 //______________________________________________________________________________
 196 TH1D * GetAcceptedFraction(TString cutname,
 197                            Int_t projectionAxis,
 198                            Int_t particleType,
 199                            Double_t lowCut,
 200                            Double_t highCut,
 201                            TString inFileName) {
 202   //
 203   // accepted fraction of tracks for ncl cut vs. pT
 204   //
 205   TFile * inFileData = TFile::Open(inFileName);
 206   TList * l = (TList * ) inFileData->Get("akalweit_TrackingUncert");
 207   THnF * histNcl = (THnF *) l->FindObject(Form("hist%s",cutname.Data()));
 208   //
 209   // select particleType
 210   //
 211   histNcl->GetAxis(4)->SetRangeUser(particleType, particleType);
 212   //
 213   // restrict eta-range
 214   //
 215   histNcl->GetAxis(2)->SetRangeUser(-0.75, 0.75);
 216   //
 217   // determine sensitivities
 218   //
 219   TH1D * hAll = histNcl->Projection(projectionAxis);
 220   hAll->SetNameTitle("hAll","hAll");
 221   //
 222   const Int_t kVeryBig = 10000;
 223   histNcl->GetAxis(0)->SetRangeUser(lowCut, highCut);
 224   TH1D * hAccepted = histNcl->Projection(1);
 225   hAccepted->SetNameTitle(Form("hAccepted_%s_%f_%f",cutname.Data(),lowCut,highCut),Form("hAccepted_%s_%f_%f",cutname.Data(),lowCut,highCut));
 226   //
 227   hAccepted->Divide(hAll);
 228   //
 229   // some cosmetics
 230   //
 231   hAccepted->SetLineWidth(3);
 232   hAccepted->SetDirectory(0);
 233   //
 234   delete l;
 235   inFileData->Close();
 236   //
 237   return hAccepted;
 238
 239
 240 }
 241
 242 //______________________________________________________________________________
 243 TH1D *GetITSTPCMatchingEff(TString inFileNameData,
 244                               Int_t particleType,
 245                               Int_t projectionAxis){
 246
 247     //
 248     // make a single plot
 249     //
 250     TH1D * hMatching = GetITSTPCMatchingHisto(inFileNameData, particleType, kTRUE, projectionAxis);
 251     hMatching->SetNameTitle(Form("PID: %d\t proj: %d",particleType,projectionAxis), Form("PID: %d\t proj: %d",particleType,projectionAxis));
 252     hMatching->Sumw2();
 253     //
 254     TH1D *hNoMatching = GetITSTPCMatchingHisto(inFileNameData, particleType, kFALSE, projectionAxis);
 255     hNoMatching->SetNameTitle(Form("PID: %d\t proj: %d",particleType,projectionAxis), Form("PID: %d\t proj: %d",particleType,projectionAxis));
 256     hNoMatching->Sumw2();
 257     //
 258     hMatching->Divide(hMatching,hNoMatching,1,1,"B");
 259
 260     //TCanvas * canvEff = new TCanvas("canvEff","matching efficiency",600,600);
 261     //hMatching->Draw();
 262
 263     return hMatching;
 264
 265
 266 }
 267
 268
 269
 270 //______________________________________________________________________________
 271 TH1D *GetITSTPCMatchingHisto(TString inFileNameData,
 272                              Int_t particleType,
 273                              Bool_t isMatched,
 274                              Int_t projectionAxis){
 275
 276   //
 277   // ITS-TPC matching histograms as a funct. of pT, eta, phi, for each species
 278   //
 279   TFile * inFileData = TFile::Open(inFileNameData);
 280   TList * l = (TList * ) inFileData->Get("akalweit_TrackingUncert");
 281   THnF * histITSTPC = (THnF *) l->FindObject("histTpcItsMatch");
 282
 283   //
 284   // select particleType
 285   //
 286   histITSTPC->GetAxis(4)->SetRangeUser(particleType, particleType);
 287
 288   //
 289   // extract isMatched = kFALSE or kTRUE
 290   //
 291   if(isMatched)
 292     histITSTPC->GetAxis(0)->SetRangeUser(1,1);
 293   else
 294     histITSTPC->GetAxis(0)->SetRangeUser(0,0);
 295
 296   TH1D * hProj  = histITSTPC->Projection(projectionAxis);
 297   hProj->SetDirectory(0);
 298   hProj->SetNameTitle(Form("h%d",projectionAxis),Form("h%d",projectionAxis));
 299
 300   //
 301   delete l;
 302   inFileData->Close();
 303   //
 304
 305
 306   return hProj;
 307
 308
 309 }