PWG3/hfe/AliHFEspectrum.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 // Class for spectrum correction
  17 // Subtraction of hadronic background, Unfolding of the data and
  18 // Renormalization done here
  19 // The following containers have to be set:
  20 //  - Correction framework container for real data
  21 //  - Correction framework container for MC (Efficiency Map)
  22 //  - Correction framework container for background coming from data
  23 //  - Correction framework container for background coming from MC
  24 //
  25 //  Author:
  26 //            Raphaelle Bailhache <R.Bailhache@gsi.de>
  27 //            Markus Fasel <M.Fasel@gsi.de>
  28 //
  29
  30 #include <TArrayD.h>
  31 #include <TH1.h>
  32 #include <TList.h>
  33 #include <TObjArray.h>
  34 #include <TROOT.h>
  35 #include <TCanvas.h>
  36 #include <TLegend.h>
  37 #include <TStyle.h>
  38 #include <TMath.h>
  39 #include <TAxis.h>
  40 #include <TGraphErrors.h>
  41 #include <TFile.h>
  42 #include <TPad.h>
  43
  44
  45 #include "AliCFContainer.h"
  46 #include "AliCFDataGrid.h"
  47 #include "AliCFEffGrid.h"
  48 #include "AliCFGridSparse.h"
  49 #include "AliCFUnfolding.h"
  50 #include "AliLog.h"
  51
  52 #include "AliHFEspectrum.h"
  53
  54 ClassImp(AliHFEspectrum)
  55
  56 //____________________________________________________________
  57 AliHFEspectrum::AliHFEspectrum(const char *name):
  58   TNamed(name, ""),
  59   fCFContainers(NULL),
  60   fTemporaryObjects(NULL),
  61   fCorrelation(NULL),
  62   fBackground(NULL),
  63   fBackgroundSource(kMCbackground),
  64   fDumpToFile(kFALSE),
  65   fNEvents(0),
  66   fStepMC(-1),
  67   fStepTrue(-1),
  68   fStepData(-1),
  69   fStepGuessedUnfolding(-1),
  70   fNumberOfIterations(5)
  71 {
  72   //
  73   // Default constructor
  74   //
  75 }
  76
  77 //____________________________________________________________
  78 AliHFEspectrum::~AliHFEspectrum(){
  79   //
  80   // Destructor
  81   //
  82   if(fCFContainers) delete fCFContainers;
  83   if(fTemporaryObjects){
  84     fTemporaryObjects->Clear();
  85     delete fTemporaryObjects;
  86   }
  87 }
  88
  89 //____________________________________________________________
  90 void AliHFEspectrum::Correct(AliCFContainer *datacontainer,AliCFContainer *mccontainer,THnSparseF *mccorrelation){
  91   //
  92   // Correct the spectrum for efficiency and unfolding
  93   // with both method and compare
  94   //
  95
  96   gStyle->SetPalette(1);
  97   gStyle->SetOptStat(1111);
  98   gStyle->SetPadBorderMode(0);
  99   gStyle->SetCanvasColor(10);
 100   gStyle->SetPadLeftMargin(0.13);
 101   gStyle->SetPadRightMargin(0.13);
 102
 103   //////////////////
 104   // Take only pt
 105   /////////////////
 106
 107   AliCFDataGrid *dataspectrum = (AliCFDataGrid *)GetSpectrum(datacontainer,20);
 108
 109   //
 110
 111   SetCorrelation(mccorrelation);
 112   SetContainer(datacontainer,AliHFEspectrum::kDataContainer);
 113   SetContainer(mccontainer,AliHFEspectrum::kMCContainer);
 114
 115   SetMCEffStep(8);
 116   SetMCTruthStep(0);
 117   SetNumberOfIteration(5);
 118   SetStepGuessedUnfolding(16);
 119   SetNumberOfIteration(5);
 120   SetStepToCorrect(20);
 121
 122   // Unfold
 123
 124   TList *listunfolded = Unfold(1);
 125   if(!listunfolded){
 126     printf("Unfolded failed\n");
 127     return;
 128   }
 129   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
 130   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
 131   if(!correctedspectrum){
 132     AliError("No corrected spectrum\n");
 133     return;
 134   }
 135  if(!residualspectrum){
 136     AliError("No residul spectrum\n");
 137     return;
 138   }
 139
 140   // Simply correct
 141  SetMCEffStep(20);
 142  AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(1);
 143
 144  //////////
 145   // Plot
 146   //////////
 147
 148   TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
 149   cresidual->Divide(2,1);
 150   cresidual->cd(1);
 151   gPad->SetLogy();
 152   TGraphErrors* residualspectrumD = Normalize(residualspectrum);
 153   if(!residualspectrumD) {
 154     AliError("Number of Events not set for the normalization");
 155     return;
 156   }
 157   residualspectrumD->SetTitle("");
 158   residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
 159   residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 160   residualspectrumD->SetMarkerStyle(26);
 161   residualspectrumD->SetMarkerColor(kBlue);
 162   residualspectrumD->SetLineColor(kBlue);
 163   residualspectrumD->Draw("AP");
 164   TGraphErrors* measuredspectrumD = Normalize(dataspectrum);
 165   measuredspectrumD->SetTitle("");
 166   measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
 167   measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 168   measuredspectrumD->SetMarkerStyle(25);
 169   measuredspectrumD->SetMarkerColor(kBlack);
 170   measuredspectrumD->SetLineColor(kBlack);
 171   measuredspectrumD->Draw("P");
 172   TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
 173   legres->AddEntry(residualspectrumD,"Residual","p");
 174   legres->AddEntry(measuredspectrumD,"Measured","p");
 175   legres->Draw("same");
 176   cresidual->cd(2);
 177   TH1D *residualTH1D = residualspectrum->Projection(0);
 178   TH1D *measuredTH1D = dataspectrum->Project(0);
 179   TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
 180   ratioresidual->SetName("ratioresidual");
 181   ratioresidual->SetTitle("");
 182   ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
 183   ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 184   ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
 185   ratioresidual->SetStats(0);
 186   ratioresidual->Draw();
 187
 188   //
 189
 190
 191   TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
 192   ccorrected->Divide(2,1);
 193   ccorrected->cd(1);
 194   gPad->SetLogy();
 195   TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 196   correctedspectrumD->SetTitle("");
 197   correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
 198   correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 199   correctedspectrumD->SetMarkerStyle(26);
 200   correctedspectrumD->SetMarkerColor(kBlue);
 201   correctedspectrumD->SetLineColor(kBlue);
 202   correctedspectrumD->Draw("AP");
 203   TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 204   alltogetherspectrumD->SetTitle("");
 205   alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
 206   alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 207   alltogetherspectrumD->SetMarkerStyle(25);
 208   alltogetherspectrumD->SetMarkerColor(kBlack);
 209   alltogetherspectrumD->SetLineColor(kBlack);
 210   alltogetherspectrumD->Draw("P");
 211   TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
 212   legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
 213   legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
 214   legcorrected->Draw("same");
 215   ccorrected->cd(2);
 216   TH1D *correctedTH1D = correctedspectrum->Projection(0);
 217   TH1D *alltogetherTH1D = alltogetherCorrection->Project(0);
 218   TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
 219   ratiocorrected->SetName("ratiocorrected");
 220   ratiocorrected->SetTitle("");
 221   ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 222   ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 223   ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
 224   ratiocorrected->SetStats(0);
 225   ratiocorrected->Draw();
 226
 227   // Efficiency correction
 228
 229   AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(mccontainer,8,7);
 230   AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(mccontainer,7,0);
 231   AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(mccontainer,8,0);
 232
 233   AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(mccontainer,20,0);
 234
 235   TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
 236   cefficiency->cd(1);
 237   TH1D* efficiencymcPIDD = efficiencymcPID->Project(0);
 238   efficiencymcPIDD->SetTitle("");
 239   efficiencymcPIDD->SetStats(0);
 240   efficiencymcPIDD->SetMarkerStyle(25);
 241   efficiencymcPIDD->Draw();
 242   TH1D* efficiencymctrackinggeoD = efficiencymctrackinggeo->Project(0);
 243   efficiencymctrackinggeoD->SetTitle("");
 244   efficiencymctrackinggeoD->SetStats(0);
 245   efficiencymctrackinggeoD->SetMarkerStyle(26);
 246   efficiencymctrackinggeoD->Draw("same");
 247   TH1D* efficiencymcallD = efficiencymcall->Project(0);
 248   efficiencymcallD->SetTitle("");
 249   efficiencymcallD->SetStats(0);
 250   efficiencymcallD->SetMarkerStyle(27);
 251   efficiencymcallD->Draw("same");
 252   TH1D* efficiencyesdallD = efficiencyesdall->Project(0);
 253   efficiencyesdallD->SetTitle("");
 254   efficiencyesdallD->SetStats(0);
 255   efficiencyesdallD->SetMarkerStyle(24);
 256   efficiencyesdallD->Draw("same");
 257   TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
 258   legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
 259   legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
 260   legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
 261   legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
 262   legeff->Draw("same");
 263
 264   // Dump to file if needed
 265
 266   if(fDumpToFile) {
 267
 268     TFile *out = new TFile("finalSpectrum.root","recreate");
 269     out->cd();
 270     //
 271     residualspectrumD->SetName("UnfoldingResidualSpectrum");
 272     residualspectrumD->Write();
 273     measuredspectrumD->SetName("MeasuredSpectrum");
 274     measuredspectrumD->Write();
 275     ratioresidual->SetName("RatioResidualSpectrum");
 276     ratioresidual->Write();
 277     //
 278     correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
 279     correctedspectrumD->Write();
 280     alltogetherspectrumD->SetName("AlltogetherSpectrum");
 281     alltogetherspectrumD->Write();
 282     ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
 283     ratiocorrected->Write();
 284     //
 285     correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
 286     correctedspectrum->Write();
 287     alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
 288     alltogetherCorrection->Write();
 289     //
 290     out->Close(); delete out;
 291
 292   }
 293
 294
 295 }
 296
 297 //____________________________________________________________
 298 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Int_t dimensions, Bool_t setBackground){
 299   //
 300   // Apply background subtraction
 301   //
 302
 303   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 304   if(!dataContainer){
 305     AliError("Data Container not available");
 306     return NULL;
 307   }
 308
 309   // Get the data grid in the requested format
 310   Bool_t initContainer = kFALSE;
 311   Int_t dims[3];
 312   switch(dimensions){
 313     case 1:   dims[0] = 0;
 314               initContainer = kTRUE;
 315               break;
 316     case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 317               initContainer = kTRUE;
 318               break;
 319     default:
 320               AliError("Container with this number of dimensions not foreseen (yet)");
 321   };
 322   if(!initContainer){
 323     AliError("Creation of the sliced containers failed. Background estimation not possible");
 324     return NULL;
 325   }
 326   AliCFContainer *spectrumSliced = GetSlicedContainer(dataContainer, dimensions, dims);
 327
 328   // Make Background Estimate
 329   AliCFDataGrid *backgroundGrid = NULL;
 330   if(fBackgroundSource == kDataBackground)
 331     backgroundGrid = MakeBackgroundEstimateFromData(dimensions);
 332   else
 333     backgroundGrid = MakeBackgroundEstimateFromMC(dimensions);
 334   if(!backgroundGrid){
 335     AliError("Background Estimate not created");
 336     ClearObject(spectrumSliced);
 337     return NULL;
 338   }
 339
 340
 341   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *spectrumSliced);
 342
 343   spectrumSubtracted->ApplyBGCorrection(*backgroundGrid);
 344   if(setBackground){
 345     if(fBackground) delete fBackground;
 346     fBackground = backgroundGrid;
 347   } else delete backgroundGrid;
 348
 349   return spectrumSubtracted;
 350 }
 351
 352 //____________________________________________________________
 353 TList *AliHFEspectrum::Unfold(Int_t dimensions,AliCFDataGrid* const bgsubpectrum){
 354
 355   //
 356   // Unfold and eventually correct for efficiency the bgsubspectrum
 357   //
 358
 359   AliCFContainer *mcContainer = GetContainer(kMCContainer);
 360   if(!mcContainer){
 361     AliError("MC Container not available");
 362     return NULL;
 363   }
 364
 365   if(!fCorrelation){
 366     AliError("No Correlation map available");
 367     return NULL;
 368   }
 369
 370   // Get the mc grid in the requested format
 371   Bool_t initContainer = kFALSE;
 372   Int_t dims[3];
 373   switch(dimensions){
 374     case 1:   dims[0] = 0;
 375               initContainer = kTRUE;
 376               break;
 377     case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 378               initContainer = kTRUE;
 379               break;
 380     default:
 381               AliError("Container with this number of dimensions not foreseen (yet)");
 382   };
 383   if(!initContainer){
 384     AliError("Creation of the sliced containers failed. Background estimation not possible");
 385     return NULL;
 386   }
 387   AliCFContainer *mcContainerD = GetSlicedContainer(mcContainer, dimensions, dims);
 388   THnSparse *correlationD = GetSlicedCorrelation(dimensions, dims);
 389
 390   // Data in the right format
 391   AliCFDataGrid *dataGrid = 0x0;
 392   if(bgsubpectrum) {
 393     if(bgsubpectrum->GetNVar()!= dimensions) {
 394       AliError("Not the expected number of dimensions for the AliCFDataGrid\n");
 395       return NULL;
 396     }
 397     dataGrid = bgsubpectrum;
 398
 399   }
 400   else {
 401
 402     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 403     if(!dataContainer){
 404       AliError("Data Container not available");
 405       return NULL;
 406     }
 407
 408     AliCFContainer *dataContainerD = GetSlicedContainer(dataContainer, dimensions, dims);
 409     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainerD, fStepData);
 410
 411   }
 412
 413
 414   // Guessed
 415   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainerD, fStepGuessedUnfolding);
 416   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
 417
 418   // Efficiency
 419   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainerD);
 420   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 421
 422   // Unfold
 423
 424   AliCFUnfolding unfolding("unfolding","",dimensions,correlationD,efficiencyD->GetGrid(),((AliCFGridSparse*)dataGrid->GetData())->GetGrid(),guessedTHnSparse);
 425   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
 426   unfolding.UseSmoothing();
 427   unfolding.Unfold();
 428
 429   // Results
 430   THnSparse* result = unfolding.GetUnfolded();
 431   THnSparse* residual = unfolding.GetEstMeasured();
 432   TList *listofresults = new TList;
 433   listofresults->SetOwner();
 434   listofresults->AddAt((THnSparse*)result->Clone(),0);
 435   listofresults->AddAt((THnSparse*)residual->Clone(),1);
 436
 437   return listofresults;
 438
 439 }
 440
 441 //____________________________________________________________
 442 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(Int_t dimensions,AliCFDataGrid* const bgsubpectrum){
 443
 444   //
 445   // Apply unfolding and efficiency correction together to bgsubspectrum
 446   //
 447
 448   AliCFContainer *mcContainer = GetContainer(kMCContainer);
 449   if(!mcContainer){
 450     AliError("MC Container not available");
 451     return NULL;
 452   }
 453
 454   // Get the data grid in the requested format
 455   Bool_t initContainer = kFALSE;
 456   Int_t dims[3];
 457   switch(dimensions){
 458     case 1:   dims[0] = 0;
 459               initContainer = kTRUE;
 460               break;
 461     case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 462               initContainer = kTRUE;
 463               break;
 464     default:
 465               AliError("Container with this number of dimensions not foreseen (yet)");
 466   };
 467   if(!initContainer){
 468     AliError("Creation of the sliced containers failed. Background estimation not possible");
 469     return NULL;
 470   }
 471   AliCFContainer *mcContainerD = GetSlicedContainer(mcContainer, dimensions, dims);
 472
 473   // Efficiency
 474   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainerD);
 475   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
 476
 477   // Data in the right format
 478   AliCFDataGrid *dataGrid = 0x0;
 479   if(bgsubpectrum) {
 480     if(bgsubpectrum->GetNVar()!= dimensions) {
 481       AliError("Not the expected number of dimensions for the AliCFDataGrid\n");
 482       return NULL;
 483     }
 484     dataGrid = bgsubpectrum;
 485
 486   }
 487   else {
 488
 489     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 490     if(!dataContainer){
 491       AliError("Data Container not available");
 492       return NULL;
 493     }
 494
 495     AliCFContainer *dataContainerD = GetSlicedContainer(dataContainer, dimensions, dims);
 496     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainerD, fStepData);
 497
 498   }
 499
 500   // Correct
 501   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
 502   result->ApplyEffCorrection(*efficiencyD);
 503
 504
 505   return result;
 506
 507 }
 508 //__________________________________________________________________________________
 509 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum) const {
 510   //
 511   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
 512   // Give the final pt spectrum to be compared
 513   //
 514
 515   if(fNEvents > 0) {
 516
 517     TH1D* projection = spectrum->Projection(0);
 518     CorrectFromTheWidth(projection);
 519     TGraphErrors *graphError = NormalizeTH1(projection);
 520     return graphError;
 521
 522   }
 523
 524   return 0x0;
 525
 526
 527 }
 528 //__________________________________________________________________________________
 529 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum) const {
 530   //
 531   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
 532   // Give the final pt spectrum to be compared
 533   //
 534   if(fNEvents > 0) {
 535
 536     TH1D* projection = spectrum->Project(0);
 537     CorrectFromTheWidth(projection);
 538     TGraphErrors *graphError = NormalizeTH1(projection);
 539     return graphError;
 540
 541   }
 542
 543   return 0x0;
 544
 545
 546 }
 547 //__________________________________________________________________________________
 548 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input) const {
 549   //
 550   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
 551   // Give the final pt spectrum to be compared
 552   //
 553   if(fNEvents > 0) {
 554
 555     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
 556     Double_t p = 0, dp = 0; Int_t point = 1;
 557     Double_t n = 0, dN = 0;
 558     Double_t nCorr = 0, dNcorr = 0;
 559     Double_t errdN = 0, errdp = 0;
 560     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
 561       point = ibin - input->GetXaxis()->GetFirst();
 562       p = input->GetXaxis()->GetBinCenter(ibin);
 563       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
 564       n = input->GetBinContent(ibin);
 565       dN = input->GetBinError(ibin);
 566
 567       // New point
 568       nCorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * 1/(2. * TMath::Pi() * p) * n;
 569       errdN = 1./(2. * TMath::Pi() * p);
 570       errdp = 1./(2. * TMath::Pi() * p*p) * n;
 571       dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
 572
 573       spectrumNormalized->SetPoint(point, p, nCorr);
 574       spectrumNormalized->SetPointError(point, dp, dNcorr);
 575     }
 576     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 577     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
 578     spectrumNormalized->SetMarkerStyle(22);
 579     spectrumNormalized->SetMarkerColor(kBlue);
 580     spectrumNormalized->SetLineColor(kBlue);
 581
 582     return spectrumNormalized;
 583
 584   }
 585
 586   return 0x0;
 587
 588
 589 }
 590 //____________________________________________________________
 591 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
 592   //
 593   // Set the container for a given step to the
 594   //
 595   if(!fCFContainers) fCFContainers = new TList;
 596   fCFContainers->AddAt(cont, type);
 597 }
 598
 599 //____________________________________________________________
 600 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
 601   //
 602   // Get Correction Framework Container for given type
 603   //
 604   if(!fCFContainers) return NULL;
 605   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
 606 }
 607
 608 //____________________________________________________________
 609 AliCFDataGrid *AliHFEspectrum::MakeBackgroundEstimateFromData(Int_t nDim){
 610   //
 611   // Make Background Estimate from Data
 612   // Container having background source from
 613   // Correction has to be done for background selection Efficiency
 614   //
 615
 616   AliCFContainer *backgroundContainer = GetContainer(kBackgroundMC);
 617   AliCFContainer *dataContainer = GetContainer(kBackgroundData);
 618   if(!backgroundContainer){
 619     AliError("MC background container not found");
 620     return NULL;
 621   }
 622   if(!dataContainer){
 623     AliError("Data container not found");
 624     return NULL;
 625   }
 626   AliInfo(Form("Making background Estimate from Data in %d Dimensions", nDim));
 627
 628   Bool_t initContainer = kFALSE;
 629   Int_t dimensions[3];
 630   switch(nDim){
 631     case 1:   dimensions[0] = 1;
 632               initContainer = kTRUE;
 633               break;
 634     case 3:   for(Int_t i = 0; i < 3; i++) dimensions[i] = i;
 635               initContainer = kTRUE;
 636               break;
 637     default:
 638               AliError("Container with this number of dimensions not foreseen (yet)");
 639   };
 640   if(!initContainer){
 641     AliError("Creation of the sliced containers failed. Background estimation not possible");
 642     return NULL;
 643   }
 644   AliCFContainer *slicedBackground = GetSlicedContainer(backgroundContainer, nDim, dimensions);
 645   AliCFContainer *slicedData = GetSlicedContainer(dataContainer, nDim, dimensions);
 646
 647   // Create Efficiency Grid and data grid
 648   AliCFEffGrid backgroundRatio("backgroundRatio", "BackgroundRatio", *slicedBackground);
 649   backgroundRatio.CalculateEfficiency(2, 1);
 650   AliCFDataGrid *backgroundEstimate = new AliCFDataGrid("backgroundEstimate", "Grid for Background Estimate", *slicedData, fStepData);
 651   backgroundEstimate->ApplyEffCorrection(backgroundRatio);
 652
 653   return backgroundEstimate;
 654 }
 655
 656 //____________________________________________________________
 657 AliCFDataGrid *AliHFEspectrum::MakeBackgroundEstimateFromMC(Int_t nDim){
 658   //
 659   // Make Background Estimate using MC
 660   // Calculate ratio of hadronic background from MC and
 661   // apply this on data
 662   //
 663   AliCFContainer *backgroundContainer = GetContainer(kBackgroundMC);
 664   AliCFContainer *dataContainer = GetContainer(kDataContainer);
 665   if(!backgroundContainer){
 666     AliError("MC background container not found");
 667     return NULL;
 668   }
 669   if(!dataContainer){
 670     AliError("Data container not found");
 671     return NULL;
 672   }
 673   AliInfo(Form("Making background Estimate from MC in %d Dimensions", nDim));
 674
 675   Bool_t initContainer = kFALSE;
 676   Int_t dimensions[3];
 677   switch(nDim){
 678     case 1:   dimensions[0] = 1;
 679               initContainer = kTRUE;
 680               break;
 681     case 3:   for(Int_t i = 0; i < 3; i++) dimensions[i] = i;
 682               initContainer = kTRUE;
 683               break;
 684     default:
 685               AliError("Container with this number of dimensions not foreseen (yet)");
 686   };
 687   if(!initContainer){
 688     AliError("Creation of the sliced containers failed. Background estimation not possible");
 689     return NULL;
 690   }
 691   AliCFContainer *slicedBackground = GetSlicedContainer(backgroundContainer, nDim, dimensions);
 692   AliCFContainer *slicedData = GetSlicedContainer(dataContainer, nDim, dimensions);
 693
 694   // Create Efficiency Grid and data grid
 695   AliCFEffGrid backgroundRatio("backgroundRatio", "BackgroundRatio", *slicedBackground);
 696   backgroundRatio.CalculateEfficiency(1, 0);
 697   AliCFDataGrid *backgroundEstimate = new AliCFDataGrid("backgroundEstimate", "Grid for Background Estimate", *slicedData, fStepData);
 698   backgroundEstimate->ApplyEffCorrection(backgroundRatio);
 699
 700   return backgroundEstimate;
 701 }
 702
 703 //____________________________________________________________
 704 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions) {
 705   //
 706   // Slice Pt bin
 707   //
 708
 709   Double_t *varMin = new Double_t[container->GetNVar()],
 710            *varMax = new Double_t[container->GetNVar()];
 711
 712   Double_t *binLimits;
 713   for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
 714
 715     binLimits = new Double_t[container->GetNBins(ivar)+1];
 716     container->GetBinLimits(ivar,binLimits);
 717     varMin[ivar] = binLimits[0];
 718     varMax[ivar] = binLimits[container->GetNBins(ivar)];
 719     delete[] binLimits;
 720   }
 721
 722   AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
 723   AddTemporaryObject(k);
 724   delete[] varMin; delete[] varMax;
 725
 726   return k;
 727
 728 }
 729
 730 //_________________________________________________________________________
 731 THnSparse *AliHFEspectrum::GetSlicedCorrelation(Int_t nDim, Int_t *dimensions) const {
 732   //
 733   // Slice Pt correlation
 734   //
 735
 736   Int_t ndimensions = fCorrelation->GetNdimensions();
 737   printf("Number of dimension %d correlation map\n",ndimensions);
 738   if(ndimensions < (2*nDim)) {
 739     AliError("Problem in the dimensions");
 740     return NULL;
 741   }
 742   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
 743   printf("Number of dimension %d container\n",ndimensionsContainer);
 744
 745   Int_t *dim = new Int_t[nDim*2];
 746   for(Int_t iter=0; iter < nDim; iter++){
 747     dim[iter] = dimensions[iter];
 748     dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
 749     printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
 750   }
 751
 752   THnSparse *k = fCorrelation->Projection(nDim*2,dim);
 753
 754   delete[] dim;
 755   return k;
 756
 757 }
 758 //___________________________________________________________________________
 759 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
 760   //
 761   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
 762   //
 763
 764   TAxis *axis = h1->GetXaxis();
 765   Int_t nbinX = h1->GetNbinsX();
 766
 767   for(Int_t i = 1; i <= nbinX; i++) {
 768
 769     Double_t width = axis->GetBinWidth(i);
 770     Double_t content = h1->GetBinContent(i);
 771     Double_t error = h1->GetBinError(i);
 772     h1->SetBinContent(i,content/width);
 773     h1->SetBinError(i,error/width);
 774   }
 775
 776 }
 777 //____________________________________________________________
 778 void AliHFEspectrum::AddTemporaryObject(TObject *o){
 779   //
 780   // Emulate garbage collection on class level
 781   //
 782   if(!fTemporaryObjects) {
 783     fTemporaryObjects= new TList;
 784     fTemporaryObjects->SetOwner();
 785   }
 786   fTemporaryObjects->Add(o);
 787 }
 788
 789 //____________________________________________________________
 790 void AliHFEspectrum::ClearObject(TObject *o){
 791   //
 792   // Do a safe deletion
 793   //
 794   if(fTemporaryObjects){
 795     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
 796     delete o;
 797   } else{
 798     // just delete
 799     delete o;
 800   }
 801 }
 802 //_________________________________________________________________________
 803 TObject* AliHFEspectrum::GetSpectrum(AliCFContainer * const c, Int_t step) {
 804   AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
 805   return data;
 806 }
 807 //_________________________________________________________________________
 808 TObject* AliHFEspectrum::GetEfficiency(AliCFContainer * const c, Int_t step, Int_t step0) {
 809   //
 810   // Create efficiency grid and calculate efficiency
 811   // of step to step0
 812   //
 813   TString name("eff");
 814   name += step;
 815   name+= step0;
 816   AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
 817   eff->CalculateEfficiency(step,step0);
 818   return eff;
 819 }