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 //____________________________________________________________
  70 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA():
  71   TNamed(),
  72   fPtMax(7.0),
  73   fListOfResult(),
  74   fWriteToFile(kTRUE)
  75 {
  76   //
  77   // Default constructor
  78   //
  79
  80   fListOfResult = new TObjArray(kNResults);
  81   fListOfResult->SetName("ListOfResults");
  82
  83
  84 }
  85 //____________________________________________________________
  86 AliHFEInclusiveSpectrumQA::AliHFEInclusiveSpectrumQA(const char *name):
  87   TNamed(name, ""),
  88   fPtMax(7.0),
  89   fListOfResult(),
  90   fWriteToFile(kTRUE)
  91 {
  92   //
  93   // Default constructor
  94   //
  95
  96   fListOfResult = new TObjArray(kNResults);
  97   fListOfResult->SetName("ListOfResults");
  98
  99
 100 }
 101
 102 //____________________________________________________________
 103 AliHFEInclusiveSpectrumQA::~AliHFEInclusiveSpectrumQA(){
 104   //
 105   // Destructor
 106   //
 107   if(fListOfResult) delete fListOfResult;
 108
 109 }
 110 //____________________________________________________________
 111 void AliHFEInclusiveSpectrumQA::AddResultAt(TObject *obj,Int_t index)
 112 {
 113   //
 114   // Init what we need for the correction:
 115   //
 116
 117   if(fListOfResult) fListOfResult->AddAt(obj,index);
 118
 119 }
 120 //____________________________________________________________
 121 TObject *AliHFEInclusiveSpectrumQA::GetResult(Int_t index)
 122 {
 123   //
 124   // Get result
 125   //
 126
 127   if(fListOfResult) return fListOfResult->UncheckedAt(index);
 128   else return 0x0;
 129
 130 }
 131 //____________________________________________________________
 132 void AliHFEInclusiveSpectrumQA::DrawProjections() const
 133 {
 134   //
 135   // get spectrum for beauty 2nd method
 136   //
 137   //
 138   AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
 139   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
 140   if(!data || !correlation) return;
 141
 142   Int_t ndimcont = data->GetNVar();
 143   Int_t ndimcor = correlation->GetNdimensions();
 144   Int_t charge = 3;
 145   Int_t centrality = 5;
 146   Int_t eta = 1;
 147
 148   TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
 149   Int_t n = 0;
 150   if(charge < ndimcont) n++;
 151   if(centrality < ndimcont) n++;
 152   if(eta < ndimcont) n++;
 153   canvas->Divide(2,n);
 154   Int_t counter = 1;
 155
 156   if(charge < ndimcont) {
 157
 158     canvas->cd(counter);
 159     TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
 160     checkcharge->Draw();
 161     counter++;
 162     canvas->cd(counter);
 163     TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
 164     projectioncharge->Draw("colz");
 165     counter++;
 166
 167   }
 168
 169   if(centrality < ndimcont) {
 170     canvas->cd(counter);
 171     TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
 172     checkcentrality->Draw();
 173     counter++;
 174     canvas->cd(counter);
 175     TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
 176     projectioncentrality->Draw("colz");
 177     counter++;
 178   }
 179
 180   if(eta < ndimcont) {
 181     canvas->cd(counter);
 182     TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
 183     checketa->Draw();
 184     counter++;
 185     canvas->cd(counter);
 186     TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
 187     projectioneta->Draw("colz");
 188   }
 189
 190
 191 }
 192 //____________________________________________________________
 193 void AliHFEInclusiveSpectrumQA::DrawSubtractContamination() const
 194 {
 195   //
 196   // get spectrum for beauty 2nd method
 197   //
 198   //
 199   TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
 200   TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
 201   if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
 202
 203   SetStyle();
 204
 205   TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
 206   cbackgroundsubtraction->Divide(2,1);
 207   cbackgroundsubtraction->cd(1);
 208   gPad->SetLogy();
 209   gPad->SetTicks();
 210   measuredTH1Daftersubstraction->SetStats(0);
 211   measuredTH1Daftersubstraction->SetTitle("");
 212   measuredTH1Daftersubstraction->GetYaxis()->SetTitleOffset(1.5);
 213   measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 214   measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 215   measuredTH1Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 216   measuredTH1Daftersubstraction->SetMarkerStyle(25);
 217   measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
 218   measuredTH1Daftersubstraction->SetLineColor(kBlack);
 219   measuredTH1Dbeforesubstraction->SetStats(0);
 220   measuredTH1Dbeforesubstraction->SetTitle("");
 221   measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 222   measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 223   measuredTH1Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 224   measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
 225   measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
 226   measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
 227   measuredTH1Daftersubstraction->Draw();
 228   measuredTH1Dbeforesubstraction->Draw("same");
 229   TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
 230   legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
 231   legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
 232   legsubstraction->SetFillStyle(0);
 233   legsubstraction->SetLineStyle(0);
 234   legsubstraction->SetLineColor(0);
 235   legsubstraction->Draw("same");
 236   cbackgroundsubtraction->cd(2);
 237   gPad->SetLogy();
 238   gPad->SetTicks();
 239   TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
 240   ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
 241   ratiomeasuredcontamination->SetTitle("");
 242   ratiomeasuredcontamination->GetYaxis()->SetTitleOffset(1.5);
 243   ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
 244   ratiomeasuredcontamination->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 245   ratiomeasuredcontamination->GetYaxis()->SetRangeUser(0.8,1.2);
 246   ratiomeasuredcontamination->GetXaxis()->SetRangeUser(0.0,fPtMax);
 247   ratiomeasuredcontamination->Sumw2();
 248   ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
 249   ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
 250   ratiomeasuredcontamination->SetStats(0);
 251   ratiomeasuredcontamination->SetMarkerStyle(26);
 252   ratiomeasuredcontamination->SetMarkerColor(kBlack);
 253   ratiomeasuredcontamination->SetLineColor(kBlack);
 254   for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
 255     ratiomeasuredcontamination->SetBinError(k+1,0.0);
 256   }
 257   ratiomeasuredcontamination->Draw("P");
 258   if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
 259
 260 }
 261
 262 //____________________________________________________________
 263 void AliHFEInclusiveSpectrumQA::DrawSubtractPhotonicBackground() const
 264 {
 265   //
 266   // get spectrum
 267   //
 268
 269   TH1D *measuredTH1Dafterphotonicsubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSPB);
 270   TH1D *measuredTH1Dbeforephotonicsubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSPB);
 271   if(!measuredTH1Dafterphotonicsubstraction || !measuredTH1Dbeforephotonicsubstraction) return;
 272
 273   SetStyle();
 274
 275   TCanvas * cphotonic = new TCanvas("Photonic Subtraction","Photonic Subtraction",1000,700);
 276   cphotonic->Divide(2,1);
 277   cphotonic->cd(1);
 278   gPad->SetLogy();
 279   gPad->SetTicks();
 280   measuredTH1Dafterphotonicsubstraction->SetStats(0);
 281   measuredTH1Dafterphotonicsubstraction->SetTitle("");
 282   measuredTH1Dafterphotonicsubstraction->GetYaxis()->SetTitleOffset(1.5);
 283   measuredTH1Dafterphotonicsubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 284   measuredTH1Dafterphotonicsubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 285   measuredTH1Dafterphotonicsubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 286   measuredTH1Dafterphotonicsubstraction->SetMarkerStyle(25);
 287   measuredTH1Dafterphotonicsubstraction->SetMarkerColor(kBlack);
 288   measuredTH1Dafterphotonicsubstraction->SetLineColor(kBlack);
 289   measuredTH1Dbeforephotonicsubstraction->SetStats(0);
 290   measuredTH1Dbeforephotonicsubstraction->SetTitle("");
 291   measuredTH1Dbeforephotonicsubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 292   measuredTH1Dbeforephotonicsubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 293   measuredTH1Dbeforephotonicsubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
 294   measuredTH1Dbeforephotonicsubstraction->SetMarkerStyle(24);
 295   measuredTH1Dbeforephotonicsubstraction->SetMarkerColor(kBlue);
 296   measuredTH1Dbeforephotonicsubstraction->SetLineColor(kBlue);
 297   measuredTH1Dafterphotonicsubstraction->Draw();
 298   measuredTH1Dbeforephotonicsubstraction->Draw("same");
 299   TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
 300   legsubstraction->AddEntry(measuredTH1Dbeforephotonicsubstraction,"With photonic background","p");
 301   legsubstraction->AddEntry(measuredTH1Dafterphotonicsubstraction,"Without photonic background","p");
 302   legsubstraction->SetFillStyle(0);
 303   legsubstraction->SetLineStyle(0);
 304   legsubstraction->SetLineColor(0);
 305   legsubstraction->Draw("same");
 306   cphotonic->cd(2);
 307   gPad->SetLogy();
 308   gPad->SetTicks();
 309   TH1D* ratiomeasuredphotonic = (TH1D*)measuredTH1Dbeforephotonicsubstraction->Clone();
 310   ratiomeasuredphotonic->SetName("ratiomeasuredphotonic");
 311   ratiomeasuredphotonic->SetTitle("");
 312   ratiomeasuredphotonic->GetYaxis()->SetTitleOffset(1.5);
 313   ratiomeasuredphotonic->GetYaxis()->SetTitle("(with photonic background - without photonic background) / with photonic background");
 314   ratiomeasuredphotonic->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 315   ratiomeasuredphotonic->GetYaxis()->SetRangeUser(0.8,1.2);
 316   ratiomeasuredphotonic->GetXaxis()->SetRangeUser(0.0,fPtMax);
 317   ratiomeasuredphotonic->Sumw2();
 318   ratiomeasuredphotonic->Add(measuredTH1Dafterphotonicsubstraction,-1.0);
 319   ratiomeasuredphotonic->Divide(measuredTH1Dbeforephotonicsubstraction);
 320   ratiomeasuredphotonic->SetStats(0);
 321   ratiomeasuredphotonic->SetMarkerStyle(26);
 322   ratiomeasuredphotonic->SetMarkerColor(kBlack);
 323   ratiomeasuredphotonic->SetLineColor(kBlack);
 324   for(Int_t k=0; k < ratiomeasuredphotonic->GetNbinsX(); k++){
 325     ratiomeasuredphotonic->SetBinError(k+1,0.0);
 326   }
 327   ratiomeasuredphotonic->Draw("P");
 328   if(fWriteToFile) cphotonic->SaveAs("PhotonicSubtracted.png");
 329
 330 }
 331
 332 //____________________________________________________________
 333 void AliHFEInclusiveSpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
 334 {
 335   //
 336   // Correct the spectrum for efficiency and unfolding
 337   // with both method and compare
 338   //
 339
 340   TH1D *afterE = 0x0;
 341   TH1D *beforeE = 0x0;
 342   TH1D *efficiencyDproj = 0x0;
 343   TF1 *efficiencyparametrized = 0x0;
 344
 345   if(typeeff== kV0) {
 346     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterV0);
 347     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeV0);
 348     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kV0Efficiency);
 349   }
 350   if(typeeff== kMC) {
 351     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
 352     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
 353     efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
 354   }
 355  if(typeeff== kParametrized) {
 356     afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
 357     beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
 358     efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
 359   }
 360
 361  if(!afterE || !beforeE) return;
 362
 363  if((typeeff==kV0 || typeeff==kMC) && !efficiencyDproj) return;
 364  if(typeeff==kParametrized && !efficiencyparametrized) return;
 365
 366   SetStyle();
 367
 368   TCanvas * cEfficiency = new TCanvas(AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],AliHFEInclusiveSpectrumQA::fgkNameCanvas[typeeff],1000,700);
 369   cEfficiency->Divide(2,1);
 370   cEfficiency->cd(1);
 371   gPad->SetLogy();
 372   gPad->SetTicks();
 373   afterE->SetStats(0);
 374   afterE->SetTitle("");
 375   afterE->GetYaxis()->SetTitleOffset(1.5);
 376   afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 377   afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 378   afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
 379   afterE->SetMarkerStyle(25);
 380   afterE->SetMarkerColor(kBlack);
 381   afterE->SetLineColor(kBlack);
 382   beforeE->SetStats(0);
 383   beforeE->SetTitle("");
 384   beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 385   beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 386   beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
 387   beforeE->SetMarkerStyle(24);
 388   beforeE->SetMarkerColor(kBlue);
 389   beforeE->SetLineColor(kBlue);
 390   afterE->Draw();
 391   beforeE->Draw("same");
 392   TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
 393   legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
 394   legefficiency->AddEntry(afterE,"After Efficiency correction","p");
 395   legefficiency->SetFillStyle(0);
 396   legefficiency->SetLineStyle(0);
 397   legefficiency->SetLineColor(0);
 398   legefficiency->Draw("same");
 399   cEfficiency->cd(2);
 400   gPad->SetTicks();
 401   if((typeeff==kV0 || typeeff==kMC)) {
 402     if(efficiencyDproj) {
 403       efficiencyDproj->SetTitle("");
 404       efficiencyDproj->SetStats(0);
 405       efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
 406       efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
 407       efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
 408       efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 409       efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
 410       efficiencyDproj->SetMarkerStyle(25);
 411       efficiencyDproj->Draw();
 412     }
 413   }
 414   if(typeeff==kParametrized) {
 415     if(efficiencyparametrized) efficiencyparametrized->Draw();
 416   }
 417
 418   if(fWriteToFile) {
 419     if(typeeff==kV0) cEfficiency->SaveAs("EfficiencyV0.png");
 420     if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
 421     if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
 422   }
 423
 424 }
 425
 426 //____________________________________________________________
 427 void AliHFEInclusiveSpectrumQA::DrawUnfolding() const
 428 {
 429   //
 430   // Draw unfolding
 431   //
 432   TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
 433   TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
 434   TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
 435   THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
 436
 437   if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
 438
 439   Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
 440
 441   SetStyle();
 442
 443   TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
 444   cunfolding->Divide(2,2);
 445   cunfolding->cd(1);
 446   gPad->SetLogy();
 447   gPad->SetTicks();
 448   residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
 449   residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
 450   residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
 451   residualspectrumD->SetStats(0);
 452   residualspectrumD->SetTitle("");
 453   residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
 454   residualspectrumD->SetMarkerStyle(26);
 455   residualspectrumD->SetMarkerColor(kBlue);
 456   residualspectrumD->SetLineColor(kBlue);
 457   residualspectrumD->Sumw2();
 458   residualspectrumD->Draw("P");
 459   measuredspectrumD->SetStats(0);
 460   measuredspectrumD->SetTitle("");
 461   measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
 462   measuredspectrumD->SetMarkerStyle(25);
 463   measuredspectrumD->SetMarkerColor(kBlack);
 464   measuredspectrumD->SetLineColor(kBlack);
 465   measuredspectrumD->Draw("same");
 466   TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
 467   legres->AddEntry(residualspectrumD,"Residual","p");
 468   legres->AddEntry(measuredspectrumD,"Measured","p");
 469   legres->SetFillStyle(0);
 470   legres->SetLineStyle(0);
 471   legres->SetLineColor(0);
 472   legres->Draw("same");
 473   cunfolding->cd(2);
 474   gPad->SetTicks();
 475   TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
 476   ratioresidual->SetName("ratioresidual");
 477   ratioresidual->SetTitle("");
 478   ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
 479   ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
 480   ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 481   ratioresidual->Divide(measuredspectrumD);
 482   ratioresidual->SetStats(0);
 483   ratioresidual->Draw();
 484   cunfolding->cd(3);
 485   gPad->SetTicks();
 486   efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
 487   efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
 488   efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
 489   efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
 490   efficiencyDproj->SetStats(0);
 491   efficiencyDproj->SetTitle("");
 492   efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
 493   efficiencyDproj->SetMarkerStyle(26);
 494   efficiencyDproj->SetMarkerColor(kBlue);
 495   efficiencyDproj->SetLineColor(kBlue);
 496   efficiencyDproj->Sumw2();
 497   efficiencyDproj->Draw("P");
 498   cunfolding->cd(4);
 499   TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
 500   projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
 501   projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
 502   projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
 503   projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
 504   projectioncorr->SetStats(0);
 505   projectioncorr->SetTitle("");
 506   projectioncorr->Draw("colz");
 507
 508   if(fWriteToFile){
 509     cunfolding->SaveAs("Unfolding.png");
 510   }
 511
 512 }
 513 //____________________________________________________________
 514 void AliHFEInclusiveSpectrumQA::DrawResult()
 515 {
 516   //
 517   // Draw Results
 518   //
 519   TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
 520   TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt( kFinalResultDirectEfficiency);
 521   if(!correctedspectrumD || !alltogetherspectrumD) return;
 522
 523   SetStyle();
 524
 525   TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
 526   ccorrected->Divide(2,1);
 527   ccorrected->cd(1);
 528   gPad->SetLogy();
 529   gPad->SetTicks();
 530   correctedspectrumD->SetTitle("");
 531   correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
 532   correctedspectrumD->GetYaxis()->SetRangeUser(0.000001,100.0);
 533   correctedspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
 534   correctedspectrumD->SetMarkerStyle(26);
 535   correctedspectrumD->SetMarkerColor(kBlue);
 536   correctedspectrumD->SetLineColor(kBlue);
 537   correctedspectrumD->Draw("AP");
 538   alltogetherspectrumD->SetTitle("");
 539   alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
 540   alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 541   alltogetherspectrumD->SetMarkerStyle(25);
 542   alltogetherspectrumD->SetMarkerColor(kBlack);
 543   alltogetherspectrumD->SetLineColor(kBlack);
 544   alltogetherspectrumD->Draw("P");
 545   TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
 546   legcorrected->AddEntry(correctedspectrumD,"Unfolded","p");
 547   legcorrected->AddEntry(alltogetherspectrumD,"Direct corrected","p");
 548   legcorrected->SetFillStyle(0);
 549   legcorrected->SetLineStyle(0);
 550   legcorrected->SetLineColor(0);
 551   legcorrected->Draw("same");
 552   ccorrected->cd(2);
 553   gPad->SetTicks();
 554   TH1D* ratiocorrected = DivideSpectra(correctedspectrumD,alltogetherspectrumD);
 555   ratiocorrected->SetName("ratiocorrected");
 556   ratiocorrected->SetTitle("");
 557   ratiocorrected->GetYaxis()->SetTitleOffset(1.5);
 558   ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 559   ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 560   ratiocorrected->GetXaxis()->SetRangeUser(0.0,fPtMax);
 561   ratiocorrected->GetYaxis()->SetRangeUser(0.4,1.4);
 562   ratiocorrected->SetStats(0);
 563   ratiocorrected->Draw();
 564   if(fWriteToFile)ccorrected->SaveAs("CorrectedResults.png");
 565
 566 }
 567 //____________________________________________________________
 568 TH1D *AliHFEInclusiveSpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb)
 569 {
 570   //
 571   // Divide Spectra
 572   //
 573
 574   TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
 575   if(!afterE) return 0x0;
 576
 577   TH1D *histoB = (TH1D*) afterE->Clone();
 578   histoB->Sumw2();
 579   histoB->SetName("ratio");
 580   TH1D *histoa = (TH1D*) afterE->Clone();
 581   histoa->Sumw2();
 582   histoa->SetName("a");
 583   TH1D *histob = (TH1D*) afterE->Clone();
 584   histob->Sumw2();
 585   histob->SetName("b");
 586
 587   double xa,ya,xb,yb,eya,eyb;
 588   Int_t npointsa = ga->GetN();
 589   Int_t npointsb = gb->GetN();
 590   if(npointsa != npointsb) {
 591     printf("Problem the two spectra have not the same number of points\n");
 592     return 0x0;
 593   }
 594   for(Int_t k = 0; k < npointsa; k++){
 595     ga->GetPoint(k,xa,ya);
 596     gb->GetPoint(k,xb,yb);
 597     //
 598     Double_t centerhisto = histoa->GetBinCenter(k+1);
 599     //
 600     if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
 601       printf("Problem not the same x axis\n");
 602       return 0x0;
 603     }
 604     histoa->SetBinContent(k+1,ya);
 605     histob->SetBinContent(k+1,yb);
 606     //
 607     eya = ga->GetErrorY(k);
 608     eyb = gb->GetErrorY(k);
 609     //
 610     histoa->SetBinError(k+1,eya);
 611     histob->SetBinError(k+1,eyb);
 612
 613   }
 614
 615   histoB->Sumw2();
 616   histoB->Divide(histoa,histob,1.0,1.0,"B");
 617
 618   return histoB;
 619
 620 }
 621 //__________________________________________
 622 void AliHFEInclusiveSpectrumQA::SetStyle() const
 623 {
 624   //
 625   // Set style
 626   //
 627
 628   gStyle->SetPalette(1);
 629   gStyle->SetOptStat(1111);
 630   gStyle->SetPadBorderMode(0);
 631   gStyle->SetCanvasColor(10);
 632   gStyle->SetPadLeftMargin(0.13);
 633   gStyle->SetPadRightMargin(0.13);
 634
 635 }