PWGHF/hfe/AliHFEInclusiveSpectrumQA.cxx

   1
   2 /**************************************************************************
   3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4 *                                                                        *
   5 * Author: The ALICE Off-line Project.                                    *
   6 * Contributors are mentioned in the code where appropriate.              *
   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 // Class for spectrum correction
  18 // Subtraction of hadronic background, Unfolding of the data and
  19 // Renormalization done here
  20 // The following containers have to be set:
  21 //  - Correction framework container for real data
  22 //  - Correction framework container for MC (Efficiency Map)
  23 //  - Correction framework container for background coming from data
  24 //  - Correction framework container for background coming from MC
  25 //
  26 //  Author:
  27 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
  28 //
  29
  30 #include <TArrayD.h>
  31 #include <TH1.h>
  32 #include <TList.h>
  33 #include <TObjArray.h>
  34 #include <TObject.h>
  35 #include <TROOT.h>
  36 #include <TCanvas.h>
  37 #include <TLegend.h>
  38 #include <TStyle.h>
  39 #include <TMath.h>
  40 #include <TAxis.h>
  41 #include <TGraphErrors.h>
  42 #include <TFile.h>
  43 #include <TPad.h>
  44 #include <TH2D.h>
  45 #include <TF1.h>
  46
  47 #include "AliPID.h"
  48 #include "AliCFContainer.h"
  49 #include "AliCFDataGrid.h"
  50 #include "AliCFEffGrid.h"
  51 #include "AliCFGridSparse.h"
  52 #include "AliCFUnfolding.h"
  53 #include "AliLog.h"
  54
  55 #include "AliHFEInclusiveSpectrumQA.h"
  56 #include "AliHFECorrectSpectrumBase.h"
  57 #include "AliHFEcuts.h"
  58 #include "AliHFEcontainer.h"
  59 #include "AliHFEtools.h"
  60
  61 ClassImp(AliHFEInclusiveSpectrumQA)
  62
  63 const Char_t *AliHFEInclusiveSpectrumQA::fgkNameCanvas[AliHFEInclusiveSpectrumQA::kNTypeEfficiency] = {
  64   "V0Efficiency",
  65   "MCEfficiency",
  66   "ParametrizedEfficiency"
  67 };
  68
  69 const Char_t *AliHFEInclusiveSpectrumQA::fgkNameCanvasND[AliHFEInclusiveSpectrumQA::kNTypeEfficiency] = {
  70   "V0EfficiencyND",
  71   "MCEfficiencyND",
  72   "ParametrizedEfficiencyND"
  73 };
  74
  75 //____________________________________________________________
  76 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA():
  77   TNamed(),
  78   fPtMax(10.0),
  79   fListOfResult(),
  80   fWriteToFile(kTRUE)
  81 {
  82   //
  83   // Default constructor
  84   //
  85
  86   fListOfResult = new TObjArray(kNResults);
  87   fListOfResult->SetName("ListOfResults");
  88
  89
  90 }
  91 //____________________________________________________________
  92 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA(const char *name):
  93   TNamed(name, ""),
  94   fPtMax(10.0),
  95   fListOfResult(),
  96   fWriteToFile(kTRUE)
  97 {
  98   //
  99   // Default constructor
 100   //
 101
 102   fListOfResult = new TObjArray(kNResults);
 103   fListOfResult->SetName("ListOfResults");
 104
 105
 106 }
 107
 108 //____________________________________________________________
 109 AliHFEInclusiveSpectrumQA::~AliHFEInclusiveSpectrumQA(){
 110   //
 111   // Destructor
 112   //
 113   if(fListOfResult) delete fListOfResult;
 114
 115 }
 116 //____________________________________________________________
 117 void AliHFEInclusiveSpectrumQA::AddResultAt(TObject *obj,Int_t index)
 118 {
 119   //
 120   // Init what we need for the correction:
 121   //
 122
 123   if(fListOfResult) fListOfResult->AddAt(obj,index);
 124
 125 }
 126 //____________________________________________________________
 127 TObject *AliHFEInclusiveSpectrumQA::GetResult(Int_t index)
 128 {
 129   //
 130   // Get result
 131   //
 132
 133   if(fListOfResult) return fListOfResult->UncheckedAt(index);
 134   else return 0x0;
 135
 136 }
 137 //____________________________________________________________
 138 void AliHFEInclusiveSpectrumQA::DrawProjections() const
 139 {
 140   //
 141   // get spectrum for beauty 2nd method
 142   //
 143   //
 144   AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
 145   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
 146   if(!data || !correlation) return;
 147
 148   Int_t ndimcont = data->GetNVar();
 149   Int_t ndimcor = correlation->GetNdimensions();
 150   Int_t charge = 3;
 151   Int_t centrality = 5;
 152   Int_t eta = 1;
 153
 154   TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
 155   Int_t n = 0;
 156   if(charge < ndimcont) n++;
 157   if(centrality < ndimcont) n++;
 158   if(eta < ndimcont) n++;
 159   canvas->Divide(2,n);
 160   Int_t counter = 1;
 161
 162   if(charge < ndimcont) {
 163
 164     canvas->cd(counter);
 165     TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
 166     checkcharge->Draw();
 167     counter++;
 168     canvas->cd(counter);
 169     TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
 170     projectioncharge->Draw("colz");
 171     counter++;
 172
 173   }
 174
 175   if(centrality < ndimcont) {
 176     canvas->cd(counter);
 177     TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
 178     checkcentrality->Draw();
 179     counter++;
 180     canvas->cd(counter);
 181     TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
 182     projectioncentrality->Draw("colz");
 183     counter++;
 184   }
 185
 186   if(eta < ndimcont) {
 187     canvas->cd(counter);
 188     TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
 189     checketa->Draw();
 190     counter++;
 191     canvas->cd(counter);
 192     TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
 193     projectioneta->Draw("colz");
 194   }
 195
 196  if(fWriteToFile) canvas->SaveAs("Projections.png");
 197
 198
 199
 200 }
 201 //____________________________________________________________
 202 void AliHFEInclusiveSpectrumQA::DrawSubtractContamination() const
 203 {
 204   //
 205   // get spectrum for beauty 2nd method
 206   //
 207   //
 208   TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
 209   TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
 210   if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
 211
 212   SetStyle();
 213
 214   TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
 215   cbackgroundsubtraction->Divide(2,1);
 216   cbackgroundsubtraction->cd(1);
 217   gPad->SetLogy();
 218   gPad->SetTicks();
 219   measuredTH1Daftersubstraction->SetStats(0);
 220   measuredTH1Daftersubstraction->SetTitle("");
 221   measuredTH1Daftersubstraction->GetYaxis()->SetTitleOffset(1.5);
 222   measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 223   measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 224   measuredTH1Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 225   measuredTH1Daftersubstraction->SetMarkerStyle(25);
 226   measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
 227   measuredTH1Daftersubstraction->SetLineColor(kBlack);
 228   measuredTH1Dbeforesubstraction->SetStats(0);
 229   measuredTH1Dbeforesubstraction->SetTitle("");
 230   measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 231   measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 232   measuredTH1Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 233   measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
 234   measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
 235   measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
 236   measuredTH1Daftersubstraction->Draw();
 237   measuredTH1Dbeforesubstraction->Draw("same");
 238   TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
 239   legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
 240   legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
 241   legsubstraction->SetFillStyle(0);
 242   legsubstraction->SetLineStyle(0);
 243   legsubstraction->SetLineColor(0);
 244   legsubstraction->Draw("same");
 245   cbackgroundsubtraction->cd(2);
 246   gPad->SetLogy();
 247   gPad->SetTicks();
 248   TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
 249   ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
 250   ratiomeasuredcontamination->SetTitle("");
 251   ratiomeasuredcontamination->GetYaxis()->SetTitleOffset(1.5);
 252   ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
 253   ratiomeasuredcontamination->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 254   ratiomeasuredcontamination->GetYaxis()->SetRangeUser(0.8,1.2);
 255   ratiomeasuredcontamination->GetXaxis()->SetRangeUser(0.0,fPtMax);
 256   ratiomeasuredcontamination->Sumw2();
 257   ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
 258   ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
 259   ratiomeasuredcontamination->SetStats(0);
 260   ratiomeasuredcontamination->SetMarkerStyle(26);
 261   ratiomeasuredcontamination->SetMarkerColor(kBlack);
 262   ratiomeasuredcontamination->SetLineColor(kBlack);
 263   for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
 264     ratiomeasuredcontamination->SetBinError(k+1,0.0);
 265   }
 266   ratiomeasuredcontamination->Draw("P");
 267   if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
 268
 269 }
 270
 271 //____________________________________________________________
 272 void AliHFEInclusiveSpectrumQA::DrawSubtractContaminationND() const
 273 {
 274   //
 275   // subtract the hadron contamination
 276   //
 277   //
 278   AliCFDataGrid *afterE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kAfterSCND);
 279   AliCFDataGrid *beforeE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kBeforeSCND);
 280   AliCFDataGrid *contamination = (AliCFDataGrid *) fListOfResult->UncheckedAt(kHadronContaminationND);
 281
 282   if(!afterE || !beforeE || !contamination) return;
 283
 284   SetStyle();
 285
 286   TCanvas * cbackgroundsubtractionND = new TCanvas("backgroundsubtractionND","backgroundsubtractionND",1000,700);
 287   cbackgroundsubtractionND->Divide(3,1);
 288   cbackgroundsubtractionND->cd(1);
 289   gPad->SetLogz();
 290   gPad->SetTicks();
 291   TH2D *measuredTH2Dbeforesubstraction = (TH2D *) beforeE->Project(0,1);
 292   measuredTH2Dbeforesubstraction->SetStats(0);
 293   measuredTH2Dbeforesubstraction->SetTitle("Before contamination");
 294   measuredTH2Dbeforesubstraction->GetZaxis()->SetTitleOffset(1.5);
 295   measuredTH2Dbeforesubstraction->GetZaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 296   measuredTH2Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 297   measuredTH2Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 298   measuredTH2Dbeforesubstraction->Draw("lego");
 299   cbackgroundsubtractionND->cd(2);
 300   gPad->SetLogz();
 301   gPad->SetTicks();
 302   TH2D *measuredTH2Daftersubstraction = (TH2D *) afterE->Project(0,1);
 303   measuredTH2Daftersubstraction->SetStats(0);
 304   measuredTH2Daftersubstraction->SetTitle("After contamination");
 305   measuredTH2Daftersubstraction->GetZaxis()->SetTitleOffset(1.5);
 306   measuredTH2Daftersubstraction->GetZaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 307   measuredTH2Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 308   measuredTH2Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 309   measuredTH2Daftersubstraction->Draw("lego");
 310   cbackgroundsubtractionND->cd(3);
 311   gPad->SetLogz();
 312   gPad->SetTicks();
 313   TH2D *measuredsubstraction = (TH2D *) contamination->Project(0,1);
 314   measuredsubstraction->SetStats(0);
 315   measuredsubstraction->SetTitle("Contamination");
 316   measuredsubstraction->GetZaxis()->SetTitleOffset(1.5);
 317   measuredsubstraction->GetZaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 318   measuredsubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 319   measuredsubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 320   measuredsubstraction->SetMarkerStyle(25);
 321   measuredsubstraction->SetMarkerColor(kBlack);
 322   measuredsubstraction->SetLineColor(kBlack);
 323   measuredsubstraction->Draw("colz");
 324   if(fWriteToFile) cbackgroundsubtractionND->SaveAs("BackgroundSubtractedND.png");
 325
 326 }
 327
 328 //____________________________________________________________
 329 void AliHFEInclusiveSpectrumQA::DrawSubtractPhotonicBackground() const
 330 {
 331   //
 332   // get spectrum
 333   //
 334
 335   TH1D *measuredTH1Dafterphotonicsubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSPB);
 336   TH1D *measuredTH1Dbeforephotonicsubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSPB);
 337   if(!measuredTH1Dafterphotonicsubstraction || !measuredTH1Dbeforephotonicsubstraction) return;
 338
 339   SetStyle();
 340
 341   TCanvas * cphotonic = new TCanvas("Photonic Subtraction","Photonic Subtraction",1000,700);
 342   cphotonic->Divide(2,1);
 343   cphotonic->cd(1);
 344   gPad->SetLogy();
 345   gPad->SetTicks();
 346   measuredTH1Dafterphotonicsubstraction->SetStats(0);
 347   measuredTH1Dafterphotonicsubstraction->SetTitle("");
 348   measuredTH1Dafterphotonicsubstraction->GetYaxis()->SetTitleOffset(1.5);
 349   measuredTH1Dafterphotonicsubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 350   measuredTH1Dafterphotonicsubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 351   measuredTH1Dafterphotonicsubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 352   measuredTH1Dafterphotonicsubstraction->SetMarkerStyle(25);
 353   measuredTH1Dafterphotonicsubstraction->SetMarkerColor(kBlack);
 354   measuredTH1Dafterphotonicsubstraction->SetLineColor(kBlack);
 355   measuredTH1Dbeforephotonicsubstraction->SetStats(0);
 356   measuredTH1Dbeforephotonicsubstraction->SetTitle("");
 357   measuredTH1Dbeforephotonicsubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 358   measuredTH1Dbeforephotonicsubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 359   measuredTH1Dbeforephotonicsubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 360   measuredTH1Dbeforephotonicsubstraction->SetMarkerStyle(24);
 361   measuredTH1Dbeforephotonicsubstraction->SetMarkerColor(kBlue);
 362   measuredTH1Dbeforephotonicsubstraction->SetLineColor(kBlue);
 363   measuredTH1Dafterphotonicsubstraction->Draw();
 364   measuredTH1Dbeforephotonicsubstraction->Draw("same");
 365   TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
 366   legsubstraction->AddEntry(measuredTH1Dbeforephotonicsubstraction,"With photonic background","p");
 367   legsubstraction->AddEntry(measuredTH1Dafterphotonicsubstraction,"Without photonic background","p");
 368   legsubstraction->SetFillStyle(0);
 369   legsubstraction->SetLineStyle(0);
 370   legsubstraction->SetLineColor(0);
 371   legsubstraction->Draw("same");
 372   cphotonic->cd(2);
 373   gPad->SetLogy();
 374   gPad->SetTicks();
 375   TH1D* ratiomeasuredphotonic = (TH1D*)measuredTH1Dbeforephotonicsubstraction->Clone();
 376   ratiomeasuredphotonic->SetName("ratiomeasuredphotonic");
 377   ratiomeasuredphotonic->SetTitle("");
 378   ratiomeasuredphotonic->GetYaxis()->SetTitleOffset(1.5);
 379   ratiomeasuredphotonic->GetYaxis()->SetTitle("(with photonic background - without photonic background) / with photonic background");
 380   ratiomeasuredphotonic->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 381   ratiomeasuredphotonic->GetYaxis()->SetRangeUser(0.8,1.2);
 382   ratiomeasuredphotonic->GetXaxis()->SetRangeUser(0.0,fPtMax);
 383   ratiomeasuredphotonic->Sumw2();
 384   ratiomeasuredphotonic->Add(measuredTH1Dafterphotonicsubstraction,-1.0);
 385   ratiomeasuredphotonic->Divide(measuredTH1Dbeforephotonicsubstraction);
 386   ratiomeasuredphotonic->SetStats(0);
 387   ratiomeasuredphotonic->SetMarkerStyle(26);
 388   ratiomeasuredphotonic->SetMarkerColor(kBlack);
 389   ratiomeasuredphotonic->SetLineColor(kBlack);
 390   for(Int_t k=0; k < ratiomeasuredphotonic->GetNbinsX(); k++){
 391     ratiomeasuredphotonic->SetBinError(k+1,0.0);
 392   }
 393   ratiomeasuredphotonic->Draw("P");
 394   if(fWriteToFile) cphotonic->SaveAs("PhotonicSubtracted.png");
 395
 396 }
 397
 398 //____________________________________________________________
 399 void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
 400 {
 401   //
 402   // Correct the spectrum for efficiency and unfolding
 403   // with both method and compare
 404   //
 405
 406   TH1D *afterE = 0x0;
 407   TH1D *beforeE = 0x0;
 408   TH1D *efficiencyDproj = 0x0;
 409   TF1 *efficiencyparametrized = 0x0;
 410
 411   if(typeeff== kV0) {
 412     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterV0);
 413     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeV0);
 414     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kV0Efficiency);
 415   }
 416   if(typeeff== kMC) {
 417     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
 418     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
 419     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
 420   }
 421  if(typeeff== kParametrized) {
 422     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
 423     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
 424     efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
 425   }
 426
 427  if(!afterE || !beforeE) return;
 428
 429  if((typeeff==kV0 || typeeff==kMC) && (!efficiencyDproj)) return;
 430  if(typeeff==kParametrized && (!efficiencyparametrized)) return;
 431
 432   SetStyle();
 433
 434   TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],1000,700);
 435   cEfficiency->Divide(2,1);
 436   cEfficiency->cd(1);
 437   gPad->SetLogy();
 438   gPad->SetTicks();
 439   afterE->SetStats(0);
 440   afterE->SetTitle("");
 441   afterE->GetYaxis()->SetTitleOffset(1.5);
 442   afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 443   afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 444   afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
 445   afterE->SetMarkerStyle(25);
 446   afterE->SetMarkerColor(kBlack);
 447   afterE->SetLineColor(kBlack);
 448   beforeE->SetStats(0);
 449   beforeE->SetTitle("");
 450   beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 451   beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 452   beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
 453   beforeE->SetMarkerStyle(24);
 454   beforeE->SetMarkerColor(kBlue);
 455   beforeE->SetLineColor(kBlue);
 456   afterE->Draw();
 457   beforeE->Draw("same");
 458   TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
 459   legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
 460   legefficiency->AddEntry(afterE,"After Efficiency correction","p");
 461   legefficiency->SetFillStyle(0);
 462   legefficiency->SetLineStyle(0);
 463   legefficiency->SetLineColor(0);
 464   legefficiency->Draw("same");
 465   cEfficiency->cd(2);
 466   gPad->SetTicks();
 467   if((typeeff==kV0 || typeeff==kMC)) {
 468     if(efficiencyDproj) {
 469       efficiencyDproj->SetTitle("");
 470       efficiencyDproj->SetStats(0);
 471       efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
 472       efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
 473       efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
 474       efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 475       efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
 476       efficiencyDproj->SetMarkerStyle(25);
 477       efficiencyDproj->Draw();
 478     }
 479   }
 480   if(typeeff==kParametrized) {
 481     if(efficiencyparametrized) {
 482       efficiencyparametrized->GetYaxis()->SetTitleOffset(1.5);
 483       efficiencyparametrized->GetYaxis()->SetRangeUser(0.0,1.0);
 484       efficiencyparametrized->GetYaxis()->SetTitle("Efficiency");
 485       efficiencyparametrized->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 486       efficiencyparametrized->GetXaxis()->SetRangeUser(0.0,fPtMax);
 487       efficiencyparametrized->Draw();
 488     }
 489   }
 490
 491   if(fWriteToFile) {
 492     if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0.png");
 493     if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
 494     if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
 495   }
 496
 497 }
 498 //____________________________________________________________
 499 void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiencyND(Int_t typeeff) const
 500 {
 501   //
 502   // Correct the spectrum for efficiency and unfolding
 503   // with both method and compare
 504   //
 505
 506   AliCFDataGrid *afterE = 0x0;
 507   AliCFDataGrid *beforeE = 0x0;
 508   AliCFEffGrid *efficiencyND = 0x0;
 509   TF1 *efficiencyparametrized = 0x0;
 510
 511   if(typeeff== kV0) {
 512     afterE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kAfterV0ND);
 513     beforeE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kBeforeV0ND);
 514     efficiencyND = (AliCFEffGrid *) fListOfResult->UncheckedAt(kV0EfficiencyND);
 515   }
 516   if(typeeff== kMC) {
 517     afterE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kAfterMCEND);
 518     beforeE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kBeforeMCEND);
 519     efficiencyND = (AliCFEffGrid *) fListOfResult->UncheckedAt(kMCEfficiencyND);
 520   }
 521  if(typeeff== kParametrized) {
 522     afterE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kAfterPEND);
 523     beforeE = (AliCFDataGrid *) fListOfResult->UncheckedAt(kBeforePEND);
 524     efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiencyND);
 525   }
 526
 527  if(!afterE || !beforeE) return;
 528
 529  if((typeeff==kV0 || typeeff==kMC) && (!efficiencyND)) return;
 530  if(typeeff==kParametrized && (!efficiencyparametrized)) return;
 531
 532   SetStyle();
 533
 534   TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvasND[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvasND[typeeff],1000,700);
 535   cEfficiency->Divide(3,1);
 536   cEfficiency->cd(1);
 537   gPad->SetLogz();
 538   gPad->SetTicks();
 539   TH2D *b2D = (TH2D *) beforeE->Project(0,1);
 540   b2D->SetStats(0);
 541   b2D->SetTitle("Before efficiency correction");
 542   b2D->GetZaxis()->SetTitleOffset(1.5);
 543   b2D->GetZaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 544   b2D->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 545   b2D->GetXaxis()->SetRangeUser(0.0,fPtMax);
 546   b2D->Draw("lego");
 547   cEfficiency->cd(2);
 548   gPad->SetLogz();
 549   gPad->SetTicks();
 550   TH2D *a2D = (TH2D *) afterE->Project(0,1);
 551   a2D->SetStats(0);
 552   a2D->SetTitle("After efficiency correction");
 553   a2D->GetZaxis()->SetTitleOffset(1.5);
 554   a2D->GetZaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 555   a2D->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 556   a2D->GetXaxis()->SetRangeUser(0.0,fPtMax);
 557   a2D->Draw("lego");
 558   cEfficiency->cd(3);
 559   gPad->SetTicks();
 560   if((typeeff==kV0 || typeeff==kMC)) {
 561     if(efficiencyND) {
 562       THnSparseF *gride = (THnSparseF *) efficiencyND->GetGrid();
 563       TH2D *e2D = (TH2D *) gride->Projection(1,0);
 564       e2D->SetStats(0);
 565       e2D->SetTitle("");
 566       e2D->SetStats(0);
 567       e2D->GetZaxis()->SetTitleOffset(1.5);
 568       e2D->GetZaxis()->SetRangeUser(0.0,1.0);
 569       e2D->GetZaxis()->SetTitle("Efficiency");
 570       e2D->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 571       e2D->GetXaxis()->SetRangeUser(0.0,fPtMax);
 572       e2D->Draw("lego");
 573     }
 574   }
 575   if(typeeff==kParametrized) {
 576     if(efficiencyparametrized) {
 577       efficiencyparametrized->GetYaxis()->SetTitleOffset(1.5);
 578       efficiencyparametrized->GetYaxis()->SetRangeUser(0.0,1.0);
 579       efficiencyparametrized->GetYaxis()->SetTitle("Efficiency");
 580       efficiencyparametrized->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 581       efficiencyparametrized->GetXaxis()->SetRangeUser(0.0,fPtMax);
 582       efficiencyparametrized->Draw();
 583     }
 584   }
 585
 586   if(fWriteToFile) {
 587     if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0ND.png");
 588     if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMCND.png");
 589     if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrizedND.png");
 590   }
 591
 592 }
 593
 594 //____________________________________________________________
 595 void AliHFEInclusiveSpectrumQA::DrawUnfolding() const
 596 {
 597   //
 598   // Draw unfolding
 599   //
 600   TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
 601   TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
 602   TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
 603   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
 604
 605   if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
 606
 607   Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
 608
 609   SetStyle();
 610
 611   TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
 612   cunfolding->Divide(2,2);
 613   cunfolding->cd(1);
 614   gPad->SetLogy();
 615   gPad->SetTicks();
 616   residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 617   residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 618   residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
 619   residualspectrumD->SetStats(0);
 620   residualspectrumD->SetTitle("");
 621   residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
 622   residualspectrumD->SetMarkerStyle(26);
 623   residualspectrumD->SetMarkerColor(kBlue);
 624   residualspectrumD->SetLineColor(kBlue);
 625   residualspectrumD->Sumw2();
 626   residualspectrumD->Draw("P");
 627   measuredspectrumD->SetStats(0);
 628   measuredspectrumD->SetTitle("");
 629   measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
 630   measuredspectrumD->SetMarkerStyle(25);
 631   measuredspectrumD->SetMarkerColor(kBlack);
 632   measuredspectrumD->SetLineColor(kBlack);
 633   measuredspectrumD->Draw("same");
 634   TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
 635   legres->AddEntry(residualspectrumD,"Residual","p");
 636   legres->AddEntry(measuredspectrumD,"Measured","p");
 637   legres->SetFillStyle(0);
 638   legres->SetLineStyle(0);
 639   legres->SetLineColor(0);
 640   legres->Draw("same");
 641   cunfolding->cd(2);
 642   gPad->SetTicks();
 643   TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
 644   ratioresidual->SetName("ratioresidual");
 645   ratioresidual->SetTitle("");
 646   ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
 647   ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
 648   ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 649   ratioresidual->Divide(measuredspectrumD);
 650   ratioresidual->SetStats(0);
 651   ratioresidual->Draw();
 652   cunfolding->cd(3);
 653   gPad->SetTicks();
 654   efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
 655   efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
 656   efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
 657   efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
 658   efficiencyDproj->SetStats(0);
 659   efficiencyDproj->SetTitle("");
 660   efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
 661   efficiencyDproj->SetMarkerStyle(26);
 662   efficiencyDproj->SetMarkerColor(kBlue);
 663   efficiencyDproj->SetLineColor(kBlue);
 664   efficiencyDproj->Sumw2();
 665   efficiencyDproj->Draw("P");
 666   cunfolding->cd(4);
 667   TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
 668   projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
 669   projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
 670   projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
 671   projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
 672   projectioncorr->SetStats(0);
 673   projectioncorr->SetTitle("");
 674   projectioncorr->Draw("colz");
 675
 676   if(fWriteToFile){
 677     cunfolding->SaveAs("Unfolding.png");
 678   }
 679
 680 }
 681 //____________________________________________________________
 682 void AliHFEInclusiveSpectrumQA::DrawResult()
 683 {
 684   //
 685   // Draw Results
 686   //
 687   TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
 688   TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt( kFinalResultDirectEfficiency);
 689   if(!correctedspectrumD || !alltogetherspectrumD) return;
 690
 691
 692   correctedspectrumD->SetTitle("");
 693   correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
 694   correctedspectrumD->GetYaxis()->SetRangeUser(0.000001,100.0);
 695   correctedspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
 696   correctedspectrumD->SetMarkerStyle(26);
 697   correctedspectrumD->SetMarkerColor(kBlue);
 698   correctedspectrumD->SetLineColor(kBlue);
 699   alltogetherspectrumD->SetTitle("");
 700   alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
 701   alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 702   alltogetherspectrumD->SetMarkerStyle(25);
 703   alltogetherspectrumD->SetMarkerColor(kBlack);
 704   alltogetherspectrumD->SetLineColor(kBlack);
 705   TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
 706   legcorrected->AddEntry(correctedspectrumD,"Unfolded","p");
 707   legcorrected->AddEntry(alltogetherspectrumD,"Direct corrected","p");
 708   legcorrected->SetFillStyle(0);
 709   legcorrected->SetLineStyle(0);
 710   legcorrected->SetLineColor(0);
 711
 712
 713   TH1D* ratiocorrected = DivideSpectra(correctedspectrumD,alltogetherspectrumD);
 714   if(ratiocorrected) {
 715     ratiocorrected->SetName("ratiocorrected");
 716     ratiocorrected->SetTitle("");
 717     ratiocorrected->GetYaxis()->SetTitleOffset(1.5);
 718     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 719     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 720     ratiocorrected->GetXaxis()->SetRangeUser(0.0,fPtMax);
 721     ratiocorrected->GetYaxis()->SetRangeUser(0.4,1.4);
 722     ratiocorrected->SetStats(0);
 723   }
 724
 725
 726   TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
 727   if(ratiocorrected) ccorrected->Divide(2,1);
 728   SetStyle();
 729   ccorrected->cd(1);
 730   gPad->SetLogy();
 731   gPad->SetTicks();
 732   correctedspectrumD->Draw("AP");
 733   alltogetherspectrumD->Draw("P");
 734   legcorrected->Draw("same");
 735   if(ratiocorrected) {
 736     ccorrected->cd(2);
 737     gPad->SetTicks();
 738     ratiocorrected->Draw();
 739   }
 740   if(fWriteToFile)ccorrected->SaveAs("CorrectedResults.png");
 741
 742 }
 743 //____________________________________________________________
 744 TH1D *AliHFEInclusiveSpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb)
 745 {
 746   //
 747   // Divide Spectra
 748   //
 749
 750   TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
 751   if(!afterE) return 0x0;
 752
 753   TH1D *histoB = (TH1D*) afterE->Clone();
 754   histoB->Sumw2();
 755   histoB->SetName("ratio");
 756   TH1D *histoa = (TH1D*) afterE->Clone();
 757   histoa->Sumw2();
 758   histoa->SetName("a");
 759   TH1D *histob = (TH1D*) afterE->Clone();
 760   histob->Sumw2();
 761   histob->SetName("b");
 762
 763   double xa,ya,xb,yb,eya,eyb;
 764   Int_t npointsa = ga->GetN();
 765   Int_t npointsb = gb->GetN();
 766   if(npointsa != npointsb) {
 767     printf("Problem the two spectra have not the same number of points\n");
 768     return 0x0;
 769   }
 770   for(Int_t k = 0; k < npointsa; k++){
 771     ga->GetPoint(k,xa,ya);
 772     gb->GetPoint(k,xb,yb);
 773     //
 774     Double_t centerhisto = histoa->GetBinCenter(k+1);
 775     //
 776     if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
 777       printf("Problem not the same x axis\n");
 778       return 0x0;
 779     }
 780     histoa->SetBinContent(k+1,ya);
 781     histob->SetBinContent(k+1,yb);
 782     //
 783     eya = ga->GetErrorY(k);
 784     eyb = gb->GetErrorY(k);
 785     //
 786     histoa->SetBinError(k+1,eya);
 787     histob->SetBinError(k+1,eyb);
 788
 789   }
 790
 791   histoB->Sumw2();
 792   histoB->Divide(histoa,histob,1.0,1.0,"B");
 793
 794   return histoB;
 795
 796 }
 797 //__________________________________________
 798 void AliHFEInclusiveSpectrumQA::SetStyle() const
 799 {
 800   //
 801   // Set style
 802   //
 803
 804   gStyle->SetPalette(1);
 805   gStyle->SetOptStat(1111);
 806   gStyle->SetPadBorderMode(0);
 807   gStyle->SetCanvasColor(10);
 808   gStyle->SetPadLeftMargin(0.13);
 809   gStyle->SetPadRightMargin(0.13);
 810
 811 }