PWGHF/hfe/AliHFEspectrum.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 //            Markus Fasel <M.Fasel@gsi.de>
  29 //
  30
  31 #include <TArrayD.h>
  32 #include <TH1.h>
  33 #include <TList.h>
  34 #include <TObjArray.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 "AliHFEspectrum.h"
  56 #include "AliHFEcuts.h"
  57 #include "AliHFEcontainer.h"
  58 #include "AliHFEtools.h"
  59
  60 ClassImp(AliHFEspectrum)
  61
  62 //____________________________________________________________
  63 AliHFEspectrum::AliHFEspectrum(const char *name):
  64   TNamed(name, ""),
  65   fCFContainers(new TObjArray(kDataContainerV0+1)),
  66   fTemporaryObjects(NULL),
  67   fCorrelation(NULL),
  68   fBackground(NULL),
  69   fEfficiencyFunction(NULL),
  70   fWeightCharm(NULL),
  71   fInclusiveSpectrum(kTRUE),
  72   fDumpToFile(kFALSE),
  73   fEtaSelected(kFALSE),
  74   fUnSetCorrelatedErrors(kTRUE),
  75   fSetSmoothing(kFALSE),
  76   fIPanaHadronBgSubtract(kFALSE),
  77   fIPanaCharmBgSubtract(kFALSE),
  78   fIPanaConversionBgSubtract(kFALSE),
  79   fIPanaNonHFEBgSubtract(kFALSE),
  80   fIPParameterizedEff(kFALSE),
  81   fNonHFEsyst(0),
  82   fBeauty2ndMethod(kFALSE),
  83   fIPEffCombinedSamples(kTRUE),
  84   fNbDimensions(1),
  85   fNMCEvents(0),
  86   fStepMC(-1),
  87   fStepTrue(-1),
  88   fStepData(-1),
  89   fStepBeforeCutsV0(-1),
  90   fStepAfterCutsV0(-1),
  91   fStepGuessedUnfolding(-1),
  92   fNumberOfIterations(5),
  93   fNRandomIter(0),
  94   fChargeChoosen(kAllCharge),
  95   fNCentralityBinAtTheEnd(0),
  96   fTestCentralityLow(-1),
  97   fTestCentralityHigh(-1),
  98   fFillMoreCorrelationMatrix(kFALSE),
  99   fHadronEffbyIPcut(NULL),
 100   fConversionEffbgc(NULL),
 101   fNonHFEEffbgc(NULL),
 102   fBSpectrum2ndMethod(NULL),
 103   fkBeauty2ndMethodfilename(""),
 104   fBeamType(0),
 105   fEtaSyst(kTRUE),
 106   fDebugLevel(0),
 107   fWriteToFile(kFALSE),
 108   fUnfoldBG(kFALSE)
 109 {
 110   //
 111   // Default constructor
 112   //
 113
 114   for(Int_t k = 0; k < 20; k++){
 115       fNEvents[k] = 0;
 116       fNMCbgEvents[k] = 0;
 117       fLowBoundaryCentralityBinAtTheEnd[k] = 0;
 118       fHighBoundaryCentralityBinAtTheEnd[k] = 0;
 119       if(k<kCentrality)
 120       {
 121           fEfficiencyTOFPIDD[k] = 0;
 122           fEfficiencyesdTOFPIDD[k] = 0;
 123           fEfficiencyIPCharmD[k] = 0;
 124           fEfficiencyIPBeautyD[k] = 0;
 125           fEfficiencyIPBeautyesdD[k] = 0;
 126           fEfficiencyIPConversionD[k] = 0;
 127           fEfficiencyIPNonhfeD[k] = 0;
 128
 129           fConversionEff[k] = 0;
 130           fNonHFEEff[k] = 0;
 131           fCharmEff[k] = 0;
 132           fBeautyEff[k] = 0;
 133           fEfficiencyCharmSigD[k] = 0;
 134           fEfficiencyBeautySigD[k] = 0;
 135           fEfficiencyBeautySigesdD[k] = 0;
 136       }
 137   }
 138   memset(fEtaRange, 0, sizeof(Double_t) * 2);
 139   memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
 140   memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
 141   memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
 142 }
 143
 144 //____________________________________________________________
 145 AliHFEspectrum::~AliHFEspectrum(){
 146   //
 147   // Destructor
 148   //
 149   if(fCFContainers) delete fCFContainers;
 150   if(fTemporaryObjects){
 151     fTemporaryObjects->Clear();
 152     delete fTemporaryObjects;
 153   }
 154 }
 155 //____________________________________________________________
 156 Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *v0hfecontainer, const AliHFEcontainer *bghfecontainer){
 157   //
 158   // Init what we need for the correction:
 159   //
 160   // Raw spectrum, hadron contamination
 161   // MC efficiency maps, correlation matrix
 162   // V0 efficiency if wanted
 163   //
 164   // This for a given dimension.
 165   // If no inclusive spectrum, then take only efficiency map for beauty electron
 166   // and the appropriate correlation matrix
 167   //
 168
 169
 170   if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
 171
 172   Int_t kNdim = 3;
 173   Int_t kNcentr =1;
 174   Int_t ptpr =0;
 175   if(fBeamType==0) kNdim=3;
 176   if(fBeamType==1)
 177   {
 178       kNdim=4;
 179       kNcentr=11;
 180       ptpr=1;
 181   }
 182
 183   Int_t dims[kNdim];
 184   // Get the requested format
 185   if(fBeamType==0){
 186     // pp analysis
 187     switch(fNbDimensions){
 188       case 1:   dims[0] = 0;
 189                 break;
 190       case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
 191                 break;
 192       case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
 193                 break;
 194       default:
 195                 AliError("Container with this number of dimensions not foreseen (yet)");
 196                 return kFALSE;
 197     };
 198   }
 199
 200   if(fBeamType==1){
 201     // PbPb analysis; centrality as first dimension
 202     Int_t nbDimensions = fNbDimensions;
 203     fNbDimensions = fNbDimensions + 1;
 204     switch(nbDimensions){
 205       case 1: dims[0] = 5;
 206               dims[1] = 0;
 207               break;
 208       case 2: dims[0] = 5;
 209               for(Int_t i = 0; i < 2; i++) dims[(i+1)] = i;
 210               break;
 211       case 3: dims[0] = 5;
 212               for(Int_t i = 0; i < 3; i++) dims[(i+1)] = i;
 213               break;
 214       default:
 215               AliError("Container with this number of dimensions not foreseen (yet)");
 216               return kFALSE;
 217     };
 218   }
 219
 220   // Data container: raw spectrum + hadron contamination
 221   AliCFContainer *datacontainer = 0x0;
 222   if(fInclusiveSpectrum) {
 223     datacontainer = datahfecontainer->GetCFContainer("recTrackContReco");
 224   }
 225   else{
 226
 227       datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
 228   }
 229   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
 230   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
 231
 232   AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 233   AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 234   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
 235
 236   SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
 237   SetContainer(contaminationcontainerD,AliHFEspectrum::kBackgroundData);
 238
 239   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps
 240   AliCFContainer *mccontaineresd = 0x0;
 241   AliCFContainer *mccontaineresdbg = 0x0;
 242   AliCFContainer *mccontainermc = 0x0;
 243   AliCFContainer *mccontainermcbg = 0x0;
 244   AliCFContainer *nonHFEweightedContainer = 0x0;
 245   AliCFContainer *convweightedContainer = 0x0;
 246   AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
 247   AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
 248
 249   if(fInclusiveSpectrum) {
 250     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco");
 251     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC");
 252   }
 253   else {
 254     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
 255     mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
 256     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
 257     mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
 258
 259     if(fNonHFEsyst){
 260       const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
 261       const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
 262       for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
 263         for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
 264           nonHFEtempContainer =  bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
 265           convtempContainer =  bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
 266           for(Int_t icentr=0;icentr<kNcentr;icentr++)
 267           {
 268               if(fBeamType==0)
 269               {
 270                   fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
 271                   fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen);
 272               }
 273               if(fBeamType==1)
 274               {
 275
 276                 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
 277                 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
 278               }
 279 //            if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
 280           }
 281           if(fBeamType == 1)break;
 282         }
 283       }
 284     }
 285     // else{
 286       nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
 287       convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
 288       if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;
 289       //}
 290   }
 291   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;
 292
 293   Int_t source = -1;
 294   if(!fInclusiveSpectrum) source = 1; //beauty
 295   AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 296   AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 297   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
 298   SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
 299   SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
 300
 301   // set charm, nonHFE container to estimate BG
 302   if(!fInclusiveSpectrum){
 303    source = 0; //charm
 304    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
 305    SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
 306
 307    //if(!fNonHFEsyst){
 308      AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
 309      SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
 310      AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
 311      SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
 312      //}
 313
 314      SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
 315
 316   }
 317   // MC container: correlation matrix
 318   THnSparseF *mccorrelation = 0x0;
 319   if(fInclusiveSpectrum) {
 320     if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 321     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 322     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 323     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
 324     else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
 325
 326     if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
 327   }
 328   else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
 329   //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
 330   if(!mccorrelation) return kFALSE;
 331   Int_t testCentralityLow = fTestCentralityLow;
 332   Int_t testCentralityHigh = fTestCentralityHigh;
 333   if(fFillMoreCorrelationMatrix) {
 334     testCentralityLow = fTestCentralityLow-1;
 335     testCentralityHigh = fTestCentralityHigh+1;
 336   }
 337   THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
 338   if(!mccorrelationD) {
 339     printf("No correlation\n");
 340     return kFALSE;
 341   }
 342   SetCorrelation(mccorrelationD);
 343
 344   // V0 container Electron, pt eta phi
 345   if(v0hfecontainer) {
 346     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
 347     if(!containerV0) return kFALSE;
 348     AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
 349     if(!containerV0Electron) return kFALSE;
 350     SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
 351   }
 352
 353
 354   if(fDebugLevel>0){
 355
 356     AliCFDataGrid *contaminationspectrum = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(contaminationcontainer,1))->Clone();
 357     contaminationspectrum->SetName("contaminationspectrum");
 358     TCanvas * ccontaminationspectrum = new TCanvas("contaminationspectrum","contaminationspectrum",1000,700);
 359     ccontaminationspectrum->Divide(3,1);
 360     ccontaminationspectrum->cd(1);
 361     TH2D * contaminationspectrum2dpteta = (TH2D *) contaminationspectrum->Project(1,0);
 362     TH2D * contaminationspectrum2dptphi = (TH2D *) contaminationspectrum->Project(2,0);
 363     TH2D * contaminationspectrum2detaphi = (TH2D *) contaminationspectrum->Project(1,2);
 364     contaminationspectrum2dpteta->SetStats(0);
 365     contaminationspectrum2dpteta->SetTitle("");
 366     contaminationspectrum2dpteta->GetXaxis()->SetTitle("#eta");
 367     contaminationspectrum2dpteta->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 368     contaminationspectrum2dptphi->SetStats(0);
 369     contaminationspectrum2dptphi->SetTitle("");
 370     contaminationspectrum2dptphi->GetXaxis()->SetTitle("#phi [rad]");
 371     contaminationspectrum2dptphi->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 372     contaminationspectrum2detaphi->SetStats(0);
 373     contaminationspectrum2detaphi->SetTitle("");
 374     contaminationspectrum2detaphi->GetXaxis()->SetTitle("#eta");
 375     contaminationspectrum2detaphi->GetYaxis()->SetTitle("#phi [rad]");
 376     contaminationspectrum2dptphi->Draw("colz");
 377     ccontaminationspectrum->cd(2);
 378     contaminationspectrum2dpteta->Draw("colz");
 379     ccontaminationspectrum->cd(3);
 380     contaminationspectrum2detaphi->Draw("colz");
 381
 382     TCanvas * ccorrelation = new TCanvas("ccorrelation","ccorrelation",1000,700);
 383     ccorrelation->cd(1);
 384     if(mccorrelationD) {
 385       if(fBeamType==0){
 386         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions,0);
 387         ptcorrelation->Draw("colz");
 388       }
 389       if(fBeamType==1){
 390         TH2D * ptcorrelation = (TH2D *) mccorrelationD->Projection(fNbDimensions+1,1);
 391         ptcorrelation->Draw("colz");
 392       }
 393     }
 394     if(fWriteToFile){
 395       ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
 396       ccorrelation->SaveAs("correlationMatrix.eps");
 397     }
 398   }
 399
 400   TFile *file = TFile::Open("tests.root","recreate");
 401   datacontainerD->Write("data");
 402   mccontainermcD->Write("mcefficiency");
 403   mccorrelationD->Write("correlationmatrix");
 404   file->Close();
 405
 406   return kTRUE;
 407 }
 408
 409
 410 //____________________________________________________________
 411 void AliHFEspectrum::CallInputFileForBeauty2ndMethod(){
 412   //
 413   // get spectrum for beauty 2nd method
 414   //
 415   //
 416     TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
 417     fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
 418 }
 419
 420
 421 //____________________________________________________________
 422 Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
 423   //
 424   // Correct the spectrum for efficiency and unfolding
 425   // with both method and compare
 426   //
 427
 428   gStyle->SetPalette(1);
 429   gStyle->SetOptStat(1111);
 430   gStyle->SetPadBorderMode(0);
 431   gStyle->SetCanvasColor(10);
 432   gStyle->SetPadLeftMargin(0.13);
 433   gStyle->SetPadRightMargin(0.13);
 434
 435   printf("Steps are: stepdata %d, stepMC %d, steptrue %d, stepV0after %d, stepV0before %d\n",fStepData,fStepMC,fStepTrue,fStepAfterCutsV0,fStepBeforeCutsV0);
 436
 437   ///////////////////////////
 438   // Check initialization
 439   ///////////////////////////
 440
 441   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
 442     AliInfo("You have to init before");
 443     return kFALSE;
 444   }
 445
 446   if((fStepTrue < 0) && (fStepMC < 0) && (fStepData < 0)) {
 447     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
 448     return kFALSE;
 449   }
 450
 451   SetNumberOfIteration(10);
 452   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 453
 454   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
 455   //////////////////////////////////
 456   // Subtract hadron background
 457   /////////////////////////////////
 458   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
 459   if(subtractcontamination) {
 460     dataspectrumaftersubstraction = SubtractBackground(kTRUE);
 461     dataGridAfterFirstSteps = dataspectrumaftersubstraction;
 462   }
 463
 464   printf("cloning spectrum\n");
 465   AliCFDataGrid *rawsave(NULL);
 466   if(dataspectrumaftersubstraction){
 467      rawsave = (AliCFDataGrid *)dataspectrumaftersubstraction->Clone("rawdata");
 468   } else {
 469     AliCFContainer *dataContainer = GetContainer(kDataContainer);
 470     if(!dataContainer){
 471       AliError("Data Container not available");
 472     }
 473     rawsave = new AliCFDataGrid("rawsave", "raw spectrum after subtraction",*dataContainer, fStepData);
 474   }
 475   printf("cloned: %p\n", rawsave);
 476
 477   ////////////////////////////////////////////////
 478   // Correct for TPC efficiency from V0
 479   ///////////////////////////////////////////////
 480   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
 481   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
 482   if(dataContainerV0){printf("Got the V0 container\n");
 483     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
 484     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
 485   }
 486
 487   //////////////////////////////////////////////////////////////////////////////
 488   // Correct for efficiency parametrized (if TPC efficiency is parametrized)
 489   /////////////////////////////////////////////////////////////////////////////
 490   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
 491   if(fEfficiencyFunction){
 492     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
 493     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
 494   }
 495
 496   ///////////////
 497   // Unfold
 498   //////////////
 499   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
 500   if(!listunfolded){
 501     printf("Unfolded failed\n");
 502     return kFALSE;
 503   }
 504   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
 505   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
 506   if(!correctedspectrum){
 507     AliError("No corrected spectrum\n");
 508     return kFALSE;
 509   }
 510   if(!residualspectrum){
 511     AliError("No residul spectrum\n");
 512     return kFALSE;
 513   }
 514
 515   /////////////////////
 516   // Simply correct
 517   ////////////////////
 518   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
 519
 520
 521   //////////
 522   // Plot
 523   //////////
 524   if(fDebugLevel > 0.0) {
 525
 526     Int_t ptpr = 0;
 527     if(fBeamType==0) ptpr=0;
 528     if(fBeamType==1) ptpr=1;
 529
 530     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
 531     ccorrected->Divide(2,1);
 532     ccorrected->cd(1);
 533     gPad->SetLogy();
 534     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 535     correctedspectrumD->SetTitle("");
 536     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
 537     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 538     correctedspectrumD->SetMarkerStyle(26);
 539     correctedspectrumD->SetMarkerColor(kBlue);
 540     correctedspectrumD->SetLineColor(kBlue);
 541     correctedspectrumD->Draw("AP");
 542     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 543     alltogetherspectrumD->SetTitle("");
 544     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
 545     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 546     alltogetherspectrumD->SetMarkerStyle(25);
 547     alltogetherspectrumD->SetMarkerColor(kBlack);
 548     alltogetherspectrumD->SetLineColor(kBlack);
 549     alltogetherspectrumD->Draw("P");
 550     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
 551     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
 552     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
 553     legcorrected->Draw("same");
 554     ccorrected->cd(2);
 555     TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
 556     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
 557     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
 558     ratiocorrected->SetName("ratiocorrected");
 559     ratiocorrected->SetTitle("");
 560     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 561     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 562     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
 563     ratiocorrected->SetStats(0);
 564     ratiocorrected->Draw();
 565     if(fWriteToFile)ccorrected->SaveAs("CorrectedPbPb.eps");
 566
 567     //TH1D unfoldingspectrac[fNCentralityBinAtTheEnd];
 568     //TGraphErrors unfoldingspectracn[fNCentralityBinAtTheEnd];
 569     //TH1D correctedspectrac[fNCentralityBinAtTheEnd];
 570     //TGraphErrors correctedspectracn[fNCentralityBinAtTheEnd];
 571
 572     TH1D *unfoldingspectrac = new TH1D[fNCentralityBinAtTheEnd];
 573     TGraphErrors *unfoldingspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
 574     TH1D *correctedspectrac = new TH1D[fNCentralityBinAtTheEnd];
 575     TGraphErrors *correctedspectracn = new TGraphErrors[fNCentralityBinAtTheEnd];
 576
 577     if(fBeamType==1) {
 578
 579       TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
 580       ccorrectedallspectra->Divide(2,1);
 581       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
 582       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
 583
 584       THnSparseF* sparsesu = (THnSparseF *) correctedspectrum;
 585       TAxis *cenaxisa = sparsesu->GetAxis(0);
 586       THnSparseF* sparsed = (THnSparseF *) alltogetherCorrection->GetGrid();
 587       TAxis *cenaxisb = sparsed->GetAxis(0);
 588       Int_t nbbin = cenaxisb->GetNbins();
 589       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
 590       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
 591       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
 592         TString titlee("corrected_centrality_bin_");
 593         titlee += "[";
 594         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
 595         titlee += "_";
 596         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
 597         titlee += "[";
 598         TString titleec("corrected_check_projection_bin_");
 599         titleec += "[";
 600         titleec += fLowBoundaryCentralityBinAtTheEnd[binc];
 601         titleec += "_";
 602         titleec += fHighBoundaryCentralityBinAtTheEnd[binc];
 603         titleec += "[";
 604         TString titleenameunotnormalized("Unfolded_Notnormalized_centrality_bin_");
 605         titleenameunotnormalized += "[";
 606         titleenameunotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
 607         titleenameunotnormalized += "_";
 608         titleenameunotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
 609         titleenameunotnormalized += "[";
 610         TString titleenameunormalized("Unfolded_normalized_centrality_bin_");
 611         titleenameunormalized += "[";
 612         titleenameunormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
 613         titleenameunormalized += "_";
 614         titleenameunormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
 615         titleenameunormalized += "[";
 616         TString titleenamednotnormalized("Dirrectcorrected_Notnormalized_centrality_bin_");
 617         titleenamednotnormalized += "[";
 618         titleenamednotnormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
 619         titleenamednotnormalized += "_";
 620         titleenamednotnormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
 621         titleenamednotnormalized += "[";
 622         TString titleenamednormalized("Dirrectedcorrected_normalized_centrality_bin_");
 623         titleenamednormalized += "[";
 624         titleenamednormalized += fLowBoundaryCentralityBinAtTheEnd[binc];
 625         titleenamednormalized += "_";
 626         titleenamednormalized += fHighBoundaryCentralityBinAtTheEnd[binc];
 627         titleenamednormalized += "[";
 628         Int_t nbEvents = 0;
 629         for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
 630           printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
 631           nbEvents += fNEvents[k];
 632         }
 633         Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
 634         Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
 635         printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
 636         Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
 637         Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
 638         printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
 639         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
 640         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
 641         TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
 642         ccorrectedcheck->cd(1);
 643         TH1D *aftersuc = (TH1D *) sparsesu->Projection(0);
 644         TH1D *aftersdc = (TH1D *) sparsed->Projection(0);
 645         aftersuc->Draw();
 646         aftersdc->Draw("same");
 647         TCanvas * ccorrectede = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
 648         ccorrectede->Divide(2,1);
 649         ccorrectede->cd(1);
 650         gPad->SetLogy();
 651         TH1D *aftersu = (TH1D *) sparsesu->Projection(1);
 652         CorrectFromTheWidth(aftersu);
 653         aftersu->SetName((const char*)titleenameunotnormalized);
 654         unfoldingspectrac[binc] = *aftersu;
 655         ccorrectede->cd(1);
 656         TGraphErrors* aftersun = NormalizeTH1N(aftersu,nbEvents);
 657         aftersun->SetTitle("");
 658         aftersun->GetYaxis()->SetTitleOffset(1.5);
 659         aftersun->GetYaxis()->SetRangeUser(0.000000001,1.0);
 660         aftersun->SetMarkerStyle(26);
 661         aftersun->SetMarkerColor(kBlue);
 662         aftersun->SetLineColor(kBlue);
 663         aftersun->Draw("AP");
 664         aftersun->SetName((const char*)titleenameunormalized);
 665         unfoldingspectracn[binc] = *aftersun;
 666         ccorrectede->cd(1);
 667         TH1D *aftersd = (TH1D *) sparsed->Projection(1);
 668         CorrectFromTheWidth(aftersd);
 669         aftersd->SetName((const char*)titleenamednotnormalized);
 670         correctedspectrac[binc] = *aftersd;
 671         ccorrectede->cd(1);
 672         TGraphErrors* aftersdn = NormalizeTH1N(aftersd,nbEvents);
 673         aftersdn->SetTitle("");
 674         aftersdn->GetYaxis()->SetTitleOffset(1.5);
 675         aftersdn->GetYaxis()->SetRangeUser(0.000000001,1.0);
 676         aftersdn->SetMarkerStyle(25);
 677         aftersdn->SetMarkerColor(kBlack);
 678         aftersdn->SetLineColor(kBlack);
 679         aftersdn->Draw("P");
 680         aftersdn->SetName((const char*)titleenamednormalized);
 681         correctedspectracn[binc] = *aftersdn;
 682         TLegend *legcorrectedud = new TLegend(0.4,0.6,0.89,0.89);
 683         legcorrectedud->AddEntry(aftersun,"Corrected","p");
 684         legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
 685         legcorrectedud->Draw("same");
 686         ccorrectedallspectra->cd(1);
 687         gPad->SetLogy();
 688         TH1D *aftersunn = (TH1D *) aftersun->Clone();
 689         aftersunn->SetMarkerStyle(stylee[binc]);
 690         aftersunn->SetMarkerColor(colorr[binc]);
 691         if(binc==0) aftersunn->Draw("AP");
 692         else aftersunn->Draw("P");
 693         legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
 694         ccorrectedallspectra->cd(2);
 695         gPad->SetLogy();
 696         TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
 697         aftersdnn->SetMarkerStyle(stylee[binc]);
 698         aftersdnn->SetMarkerColor(colorr[binc]);
 699         if(binc==0) aftersdnn->Draw("AP");
 700         else aftersdnn->Draw("P");
 701         legtotalg->AddEntry(aftersdnn,(const char*) titlee,"p");
 702         ccorrectede->cd(2);
 703         TH1D* ratiocorrectedbinc = (TH1D*)aftersu->Clone();
 704         TString titleee("ratiocorrected_bin_");
 705         titleee += binc;
 706         ratiocorrectedbinc->SetName((const char*) titleee);
 707         ratiocorrectedbinc->SetTitle("");
 708         ratiocorrectedbinc->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 709         ratiocorrectedbinc->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 710         ratiocorrectedbinc->Divide(aftersu,aftersd,1,1);
 711         ratiocorrectedbinc->SetStats(0);
 712         ratiocorrectedbinc->Draw();
 713       }
 714
 715       ccorrectedallspectra->cd(1);
 716       legtotal->Draw("same");
 717       ccorrectedallspectra->cd(2);
 718       legtotalg->Draw("same");
 719
 720       cenaxisa->SetRange(0,nbbin);
 721       cenaxisb->SetRange(0,nbbin);
 722       if(fWriteToFile) ccorrectedallspectra->SaveAs("CorrectedPbPb.eps");
 723     }
 724
 725     // Dump to file if needed
 726     if(fDumpToFile) {
 727       TFile *out = new TFile("finalSpectrum.root","recreate");
 728       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
 729       correctedspectrumD->Write();
 730       alltogetherspectrumD->SetName("AlltogetherSpectrum");
 731       alltogetherspectrumD->Write();
 732       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
 733       ratiocorrected->Write();
 734       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
 735       correctedspectrum->Write();
 736       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
 737       alltogetherCorrection->Write();
 738       rawsave->Write();
 739       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
 740               unfoldingspectrac[binc].Write();
 741               unfoldingspectracn[binc].Write();
 742               correctedspectrac[binc].Write();
 743               correctedspectracn[binc].Write();
 744       }
 745       out->Close(); delete out;
 746     }
 747
 748     if (unfoldingspectrac) delete[] unfoldingspectrac;
 749     if (unfoldingspectracn)  delete[] unfoldingspectracn;
 750     if (correctedspectrac) delete[] correctedspectrac;
 751     if (correctedspectracn) delete[] correctedspectracn;
 752
 753   }
 754
 755   return kTRUE;
 756 }
 757
 758 //____________________________________________________________
 759 Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
 760   //
 761   // Correct the spectrum for efficiency and unfolding for beauty analysis
 762   // with both method and compare
 763   //
 764
 765   gStyle->SetPalette(1);
 766   gStyle->SetOptStat(1111);
 767   gStyle->SetPadBorderMode(0);
 768   gStyle->SetCanvasColor(10);
 769   gStyle->SetPadLeftMargin(0.13);
 770   gStyle->SetPadRightMargin(0.13);
 771
 772   ///////////////////////////
 773   // Check initialization
 774   ///////////////////////////
 775
 776   if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
 777     AliInfo("You have to init before");
 778     return kFALSE;
 779   }
 780
 781   if((fStepTrue == 0) && (fStepMC == 0) && (fStepData == 0)) {
 782     AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
 783     return kFALSE;
 784   }
 785
 786   SetNumberOfIteration(10);
 787   SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 788
 789   AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
 790   //////////////////////////////////
 791   // Subtract hadron background
 792   /////////////////////////////////
 793   AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
 794   AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
 795   TGraphErrors *gNormalizedRawSpectrum = 0x0;
 796   if(subtractcontamination) {
 797       if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
 798       else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
 799       unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
 800       dataGridAfterFirstSteps = dataspectrumaftersubstraction;
 801       gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
 802   }
 803
 804   printf("after normalize getting IP \n");
 805
 806   /////////////////////////////////////////////////////////////////////////////////////////
 807   // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
 808   /////////////////////////////////////////////////////////////////////////////////////////
 809
 810   AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
 811   AliCFDataGrid *dataspectrumafterV0efficiencycorrection = 0x0;
 812   AliCFContainer *dataContainerV0 = GetContainer(kDataContainerV0);
 813
 814   if(fEfficiencyFunction){
 815     dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
 816     dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
 817   }
 818   else if(dataContainerV0){
 819     dataspectrumafterV0efficiencycorrection = CorrectV0Efficiency(dataspectrumaftersubstraction);
 820     dataGridAfterFirstSteps = dataspectrumafterV0efficiencycorrection;
 821   }
 822
 823
 824
 825   ///////////////
 826   // Unfold
 827   //////////////
 828   TList *listunfolded = Unfold(dataGridAfterFirstSteps);
 829   if(!listunfolded){
 830     printf("Unfolded failed\n");
 831     return kFALSE;
 832   }
 833   THnSparse *correctedspectrum = (THnSparse *) listunfolded->At(0);
 834   THnSparse *residualspectrum = (THnSparse *) listunfolded->At(1);
 835   if(!correctedspectrum){
 836     AliError("No corrected spectrum\n");
 837     return kFALSE;
 838   }
 839   if(!residualspectrum){
 840     AliError("No residual spectrum\n");
 841     return kFALSE;
 842   }
 843
 844   /////////////////////
 845   // Simply correct
 846   ////////////////////
 847
 848   AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
 849
 850
 851   //////////
 852   // Plot
 853   //////////
 854
 855   if(fDebugLevel > 0.0) {
 856
 857     Int_t ptpr = 0;
 858     if(fBeamType==0) ptpr=0;
 859     if(fBeamType==1) ptpr=1;
 860
 861     TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
 862     ccorrected->Divide(2,1);
 863     ccorrected->cd(1);
 864     gPad->SetLogy();
 865     TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
 866     correctedspectrumD->SetTitle("");
 867     correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
 868     correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 869     correctedspectrumD->SetMarkerStyle(26);
 870     correctedspectrumD->SetMarkerColor(kBlue);
 871     correctedspectrumD->SetLineColor(kBlue);
 872     correctedspectrumD->Draw("AP");
 873     TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
 874     alltogetherspectrumD->SetTitle("");
 875     alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
 876     alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
 877     alltogetherspectrumD->SetMarkerStyle(25);
 878     alltogetherspectrumD->SetMarkerColor(kBlack);
 879     alltogetherspectrumD->SetLineColor(kBlack);
 880     alltogetherspectrumD->Draw("P");
 881     TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
 882     legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
 883     legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
 884     legcorrected->Draw("same");
 885     ccorrected->cd(2);
 886     TH1D *correctedTH1D = correctedspectrum->Projection(ptpr);
 887     TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(ptpr);
 888     TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
 889     ratiocorrected->SetName("ratiocorrected");
 890     ratiocorrected->SetTitle("");
 891     ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
 892     ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 893     ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
 894     ratiocorrected->SetStats(0);
 895     ratiocorrected->Draw();
 896     if(fWriteToFile) ccorrected->SaveAs("CorrectedBeauty.eps");
 897
 898     if(fBeamType == 0){
 899         if(fNonHFEsyst){
 900             CalculateNonHFEsyst(0);
 901         }
 902     }
 903
 904     // Dump to file if needed
 905
 906     if(fDumpToFile) {
 907         // to do centrality dependent
 908
 909       TFile *out;
 910       out = new TFile("finalSpectrum.root","recreate");
 911       out->cd();
 912       //
 913       correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
 914       correctedspectrumD->Write();
 915       alltogetherspectrumD->SetName("AlltogetherSpectrum");
 916       alltogetherspectrumD->Write();
 917       ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
 918       ratiocorrected->Write();
 919       //
 920       correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
 921       correctedspectrum->Write();
 922       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
 923       alltogetherCorrection->Write();
 924       //
 925       if(unnormalizedRawSpectrum) {
 926       unnormalizedRawSpectrum->SetName("beautyAfterIP");
 927       unnormalizedRawSpectrum->Write();
 928       }
 929
 930       if(gNormalizedRawSpectrum){
 931         gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
 932         gNormalizedRawSpectrum->Write();
 933       }
 934
 935       if(fBeamType==0) {
 936           Int_t countpp=0;
 937           fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
 938           fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
 939           fEfficiencyCharmSigD[countpp]->Write();
 940           fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
 941           fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
 942           fEfficiencyBeautySigD[countpp]->Write();
 943           fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
 944           fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
 945           fCharmEff[countpp]->Write();
 946           fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
 947           fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
 948           fBeautyEff[countpp]->Write();
 949           fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
 950           fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
 951           fConversionEff[countpp]->Write();
 952           fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
 953           fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
 954           fNonHFEEff[countpp]->Write();
 955       }
 956
 957       if(fBeamType==1) {
 958
 959           TGraphErrors* correctedspectrumDc[kCentrality];
 960           TGraphErrors* alltogetherspectrumDc[kCentrality];
 961           for(Int_t i=0;i<kCentrality-2;i++)
 962           {
 963               correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
 964               correctedspectrumDc[i] = Normalize(correctedspectrum,i);
 965               if(correctedspectrumDc[i]){
 966                 correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
 967                 correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
 968                 correctedspectrumDc[i]->Write();
 969               }
 970               alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
 971               alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
 972               if(alltogetherspectrumDc[i]){
 973                 alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
 974                 alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
 975                 alltogetherspectrumDc[i]->Write();
 976               }
 977
 978               TH1D *centrcrosscheck = correctedspectrum->Projection(0);
 979               centrcrosscheck->SetTitle(Form("centrality_%i",i));
 980               centrcrosscheck->SetName(Form("centrality_%i",i));
 981               centrcrosscheck->Write();
 982
 983               TH1D *correctedTH1Dc = correctedspectrum->Projection(ptpr);
 984               TH1D *alltogetherTH1Dc = (TH1D *) alltogetherCorrection->Project(ptpr);
 985
 986               TH1D *centrcrosscheck2 =  (TH1D *) alltogetherCorrection->Project(0);
 987               centrcrosscheck2->SetTitle(Form("centrality2_%i",i));
 988               centrcrosscheck2->SetName(Form("centrality2_%i",i));
 989               centrcrosscheck2->Write();
 990
 991               TH1D* ratiocorrectedc = (TH1D*)correctedTH1D->Clone();
 992               ratiocorrectedc->Divide(correctedTH1Dc,alltogetherTH1Dc,1,1);
 993               ratiocorrectedc->SetTitle(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
 994               ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
 995               ratiocorrectedc->Write();
 996
 997               fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
 998               fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
 999               fEfficiencyCharmSigD[i]->Write();
1000               fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
1001               fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
1002               fEfficiencyBeautySigD[i]->Write();
1003               fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
1004               fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
1005               fCharmEff[i]->Write();
1006               fBeautyEff[i]->SetTitle(Form("IPEfficiencyForBeautyCent%i",i));
1007               fBeautyEff[i]->SetName(Form("IPEfficiencyForBeautyCent%i",i));
1008               fBeautyEff[i]->Write();
1009               fConversionEff[i]->SetTitle(Form("IPEfficiencyForConversionCent%i",i));
1010               fConversionEff[i]->SetName(Form("IPEfficiencyForConversionCent%i",i));
1011               fConversionEff[i]->Write();
1012               fNonHFEEff[i]->SetTitle(Form("IPEfficiencyForNonHFECent%i",i));
1013               fNonHFEEff[i]->SetName(Form("IPEfficiencyForNonHFECent%i",i));
1014               fNonHFEEff[i]->Write();
1015           }
1016
1017       }
1018
1019       out->Close();
1020       delete out;
1021     }
1022   }
1023
1024   return kTRUE;
1025 }
1026
1027 //____________________________________________________________
1028 AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
1029   //
1030   // Apply background subtraction
1031   //
1032
1033     Int_t ptpr = 0;
1034     Int_t nbins=1;
1035     if(fBeamType==0)
1036     {
1037         ptpr=0;
1038         nbins=1;
1039     }
1040     if(fBeamType==1)
1041     {
1042         ptpr=1;
1043         nbins=2;
1044     }
1045
1046   // Raw spectrum
1047   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1048   if(!dataContainer){
1049     AliError("Data Container not available");
1050     return NULL;
1051   }
1052   printf("Step data: %d\n",fStepData);
1053   AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
1054
1055   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
1056   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction");
1057
1058
1059   // Background Estimate
1060   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
1061   if(!backgroundContainer){
1062     AliError("MC background container not found");
1063     return NULL;
1064   }
1065
1066   Int_t stepbackground = 1; // 2 for !fInclusiveSpectrum analysis(old method)
1067   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
1068
1069   if(!fInclusiveSpectrum){
1070     //Background subtraction for IP analysis
1071
1072     TH1D *incElecCent[kCentrality-1];
1073     TH1D *charmCent[kCentrality-1];
1074     TH1D *convCent[kCentrality-1];
1075     TH1D *nonHFECent[kCentrality-1];
1076     TH1D *subtractedCent[kCentrality-1];
1077     TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
1078     CorrectFromTheWidth(measuredTH1Draw);
1079     if(fBeamType==1){
1080       THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1081       for(Int_t icent =  1; icent < kCentrality-1; icent++){
1082         sparseIncElec->GetAxis(0)->SetRange(icent,icent);
1083         incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
1084         CorrectFromTheWidth(incElecCent[icent-1]);
1085       }
1086     }
1087     TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
1088     rawspectra->cd();
1089     rawspectra->SetLogy();
1090     gStyle->SetOptStat(0);
1091     TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
1092     measuredTH1Draw->SetMarkerStyle(20);
1093     measuredTH1Draw->Draw();
1094     measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
1095     lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
1096     TH1D* htemp;
1097     Int_t* bins=new Int_t[2];
1098     if(fIPanaHadronBgSubtract){
1099       // Hadron background
1100         printf("Hadron background for IP analysis subtracted!\n");
1101         if(fBeamType==0)
1102         {
1103
1104             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
1105             bins[0]=htemp->GetNbinsX();
1106         }
1107         if(fBeamType==1)
1108         {
1109             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
1110             bins[0]=htemp->GetNbinsX();
1111             htemp  = (TH1D *) fHadronEffbyIPcut->Projection(1);
1112             bins[1]=htemp->GetNbinsX();
1113         }
1114       AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
1115       hbgContainer->SetGrid(fHadronEffbyIPcut);
1116       backgroundGrid->Multiply(hbgContainer,1);
1117       // draw raw hadron bg spectra
1118       TH1D *hadronbg= (TH1D *) backgroundGrid->Project(ptpr);
1119       CorrectFromTheWidth(hadronbg);
1120       hadronbg->SetMarkerColor(7);
1121       hadronbg->SetMarkerStyle(20);
1122       rawspectra->cd();
1123       hadronbg->Draw("samep");
1124       lRaw->AddEntry(hadronbg,"hadrons");
1125       // subtract hadron contamination
1126       spectrumSubtracted->Add(backgroundGrid,-1.0);
1127     }
1128     if(fIPanaCharmBgSubtract){
1129       // Charm background
1130       printf("Charm background for IP analysis subtracted!\n");
1131       AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
1132       // draw charm bg spectra
1133       TH1D *charmbg= (TH1D *) charmbgContainer->Project(ptpr);
1134       CorrectFromTheWidth(charmbg);
1135       charmbg->SetMarkerColor(3);
1136       charmbg->SetMarkerStyle(20);
1137       rawspectra->cd();
1138       charmbg->Draw("samep");
1139       lRaw->AddEntry(charmbg,"charm elecs");
1140       // subtract charm background
1141       spectrumSubtracted->Add(charmbgContainer,-1.0);
1142       if(fBeamType==1){
1143         THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
1144         for(Int_t icent =  1; icent < kCentrality-1; icent++){
1145           sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
1146           charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
1147           CorrectFromTheWidth(charmCent[icent-1]);
1148         }
1149       }
1150     }
1151     if(fIPanaConversionBgSubtract){
1152       // Conversion background
1153       AliCFDataGrid *conversionbgContainer = (AliCFDataGrid *) GetConversionBackground();
1154       // draw conversion bg spectra
1155       TH1D *conversionbg= (TH1D *) conversionbgContainer->Project(ptpr);
1156       CorrectFromTheWidth(conversionbg);
1157       conversionbg->SetMarkerColor(4);
1158       conversionbg->SetMarkerStyle(20);
1159       rawspectra->cd();
1160       conversionbg->Draw("samep");
1161       lRaw->AddEntry(conversionbg,"conversion elecs");
1162       // subtract conversion background
1163       spectrumSubtracted->Add(conversionbgContainer,-1.0);
1164       if(fBeamType==1){
1165         THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
1166         for(Int_t icent =  1; icent < kCentrality-1; icent++){
1167           sparseconvElec->GetAxis(0)->SetRange(icent,icent);
1168           convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
1169           CorrectFromTheWidth(convCent[icent-1]);
1170         }
1171       }
1172     }
1173     if(fIPanaNonHFEBgSubtract){
1174       // NonHFE background
1175       AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground();
1176       // draw Dalitz/dielectron bg spectra
1177       TH1D *nonhfebg= (TH1D *) nonHFEbgContainer->Project(ptpr);
1178       CorrectFromTheWidth(nonhfebg);
1179       nonhfebg->SetMarkerColor(6);
1180       nonhfebg->SetMarkerStyle(20);
1181       rawspectra->cd();
1182       nonhfebg->Draw("samep");
1183       lRaw->AddEntry(nonhfebg,"non-HF elecs");
1184       // subtract Dalitz/dielectron background
1185       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
1186       if(fBeamType==1){
1187         THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
1188         for(Int_t icent =  1; icent < kCentrality-1; icent++){
1189           sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
1190           nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
1191           CorrectFromTheWidth(nonHFECent[icent-1]);
1192         }
1193       }
1194     }
1195
1196     TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
1197     CorrectFromTheWidth(rawbgsubtracted);
1198     rawbgsubtracted->SetMarkerStyle(24);
1199     rawspectra->cd();
1200     lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
1201     rawbgsubtracted->Draw("samep");
1202     lRaw->Draw("SAME");
1203     gPad->SetGrid();
1204     //rawspectra->SaveAs("rawspectra.eps");
1205
1206     if(fBeamType==1){
1207       THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1208       for(Int_t icent =  1; icent < kCentrality-1; icent++){
1209         sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
1210         subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
1211         CorrectFromTheWidth(subtractedCent[icent-1]);
1212       }
1213
1214       TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
1215       TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
1216       centRaw->Divide(3,3);
1217       for(Int_t icent = 1; icent < kCentrality-1; icent++){
1218         centRaw->cd(icent);
1219         gPad->SetLogx();
1220         gPad->SetLogy();
1221         incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
1222         incElecCent[icent-1]->Draw("p");
1223         incElecCent[icent-1]->SetMarkerColor(1);
1224         incElecCent[icent-1]->SetMarkerStyle(20);
1225         charmCent[icent-1]->Draw("samep");
1226         charmCent[icent-1]->SetMarkerColor(3);
1227         charmCent[icent-1]->SetMarkerStyle(20);
1228         convCent[icent-1]->Draw("samep");
1229         convCent[icent-1]->SetMarkerColor(4);
1230         convCent[icent-1]->SetMarkerStyle(20);
1231         nonHFECent[icent-1]->Draw("samep");
1232         nonHFECent[icent-1]->SetMarkerColor(6);
1233         nonHFECent[icent-1]->SetMarkerStyle(20);
1234         subtractedCent[icent-1]->Draw("samep");
1235         subtractedCent[icent-1]->SetMarkerStyle(24);
1236         if(icent == 1){
1237           lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
1238           lCentRaw->AddEntry(charmCent[0],"charm elecs");
1239           lCentRaw->AddEntry(convCent[0],"conversion elecs");
1240           lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
1241           lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
1242           lCentRaw->Draw("SAME");
1243         }
1244       }
1245     }
1246
1247     delete[] bins;
1248
1249   }
1250   else{
1251     // Subtract
1252     spectrumSubtracted->Add(backgroundGrid,-1.0);
1253   }
1254
1255   if(setBackground){
1256     if(fBackground) delete fBackground;
1257     fBackground = backgroundGrid;
1258   } else delete backgroundGrid;
1259
1260
1261   if(fDebugLevel > 0) {
1262
1263     Int_t ptprd;
1264     if(fBeamType==0) ptprd=0;
1265     if(fBeamType==1) ptprd=1;
1266
1267     TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
1268     cbackgroundsubtraction->Divide(3,1);
1269     cbackgroundsubtraction->cd(1);
1270     gPad->SetLogy();
1271     TH1D *measuredTH1Daftersubstraction = (TH1D *) spectrumSubtracted->Project(ptprd);
1272     TH1D *measuredTH1Dbeforesubstraction = (TH1D *) dataspectrumbeforesubstraction->Project(ptprd);
1273     CorrectFromTheWidth(measuredTH1Daftersubstraction);
1274     CorrectFromTheWidth(measuredTH1Dbeforesubstraction);
1275     measuredTH1Daftersubstraction->SetStats(0);
1276     measuredTH1Daftersubstraction->SetTitle("");
1277     measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1278     measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1279     measuredTH1Daftersubstraction->SetMarkerStyle(25);
1280     measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
1281     measuredTH1Daftersubstraction->SetLineColor(kBlack);
1282     measuredTH1Dbeforesubstraction->SetStats(0);
1283     measuredTH1Dbeforesubstraction->SetTitle("");
1284     measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1285     measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1286     measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
1287     measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
1288     measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
1289     measuredTH1Daftersubstraction->Draw();
1290     measuredTH1Dbeforesubstraction->Draw("same");
1291     TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
1292     legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
1293     legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
1294     legsubstraction->Draw("same");
1295     cbackgroundsubtraction->cd(2);
1296     gPad->SetLogy();
1297     TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
1298     ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
1299     ratiomeasuredcontamination->SetTitle("");
1300     ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1301     ratiomeasuredcontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1302     ratiomeasuredcontamination->Sumw2();
1303     ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
1304     ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
1305     ratiomeasuredcontamination->SetStats(0);
1306     ratiomeasuredcontamination->SetMarkerStyle(26);
1307     ratiomeasuredcontamination->SetMarkerColor(kBlack);
1308     ratiomeasuredcontamination->SetLineColor(kBlack);
1309     for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
1310       ratiomeasuredcontamination->SetBinError(k+1,0.0);
1311     }
1312     ratiomeasuredcontamination->Draw("P");
1313     cbackgroundsubtraction->cd(3);
1314     TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(ptprd);
1315     CorrectFromTheWidth(measuredTH1background);
1316     measuredTH1background->SetStats(0);
1317     measuredTH1background->SetTitle("");
1318     measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1319     measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1320     measuredTH1background->SetMarkerStyle(26);
1321     measuredTH1background->SetMarkerColor(kBlack);
1322     measuredTH1background->SetLineColor(kBlack);
1323     measuredTH1background->Draw();
1324     if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.eps");
1325
1326     if(fBeamType==1) {
1327
1328       TCanvas * cbackgrounde = new TCanvas("BackgroundSubtraction_allspectra","BackgroundSubtraction_allspectra",1000,700);
1329       cbackgrounde->Divide(2,1);
1330       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1331       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1332
1333       THnSparseF* sparsesubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
1334       TAxis *cenaxisa = sparsesubtracted->GetAxis(0);
1335       THnSparseF* sparsebefore = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
1336       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1337       Int_t nbbin = cenaxisb->GetNbins();
1338       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1339       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1340       for(Int_t binc = 0; binc < nbbin; binc++){
1341         TString titlee("BackgroundSubtraction_centrality_bin_");
1342         titlee += binc;
1343         TCanvas * cbackground = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1344         cbackground->Divide(2,1);
1345         cbackground->cd(1);
1346         gPad->SetLogy();
1347         cenaxisa->SetRange(binc+1,binc+1);
1348         cenaxisb->SetRange(binc+1,binc+1);
1349         TH1D *aftersubstraction = (TH1D *) sparsesubtracted->Projection(1);
1350         TH1D *beforesubstraction = (TH1D *) sparsebefore->Projection(1);
1351         CorrectFromTheWidth(aftersubstraction);
1352         CorrectFromTheWidth(beforesubstraction);
1353         aftersubstraction->SetStats(0);
1354         aftersubstraction->SetTitle((const char*)titlee);
1355         aftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1356         aftersubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1357         aftersubstraction->SetMarkerStyle(25);
1358         aftersubstraction->SetMarkerColor(kBlack);
1359         aftersubstraction->SetLineColor(kBlack);
1360         beforesubstraction->SetStats(0);
1361         beforesubstraction->SetTitle((const char*)titlee);
1362         beforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1363         beforesubstraction->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1364         beforesubstraction->SetMarkerStyle(24);
1365         beforesubstraction->SetMarkerColor(kBlue);
1366         beforesubstraction->SetLineColor(kBlue);
1367         aftersubstraction->Draw();
1368         beforesubstraction->Draw("same");
1369         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1370         lega->AddEntry(beforesubstraction,"With hadron contamination","p");
1371         lega->AddEntry(aftersubstraction,"Without hadron contamination ","p");
1372         lega->Draw("same");
1373         cbackgrounde->cd(1);
1374         gPad->SetLogy();
1375         TH1D *aftersubtractionn = (TH1D *) aftersubstraction->Clone();
1376         aftersubtractionn->SetMarkerStyle(stylee[binc]);
1377         aftersubtractionn->SetMarkerColor(colorr[binc]);
1378         if(binc==0) aftersubtractionn->Draw();
1379         else aftersubtractionn->Draw("same");
1380         legtotal->AddEntry(aftersubtractionn,(const char*) titlee,"p");
1381         cbackgrounde->cd(2);
1382         gPad->SetLogy();
1383         TH1D *aftersubtractionng = (TH1D *) aftersubstraction->Clone();
1384         aftersubtractionng->SetMarkerStyle(stylee[binc]);
1385         aftersubtractionng->SetMarkerColor(colorr[binc]);
1386         if(fNEvents[binc] > 0.0) aftersubtractionng->Scale(1/(Double_t)fNEvents[binc]);
1387         if(binc==0) aftersubtractionng->Draw();
1388         else aftersubtractionng->Draw("same");
1389         legtotalg->AddEntry(aftersubtractionng,(const char*) titlee,"p");
1390         cbackground->cd(2);
1391         TH1D* ratiocontamination = (TH1D*)beforesubstraction->Clone();
1392         ratiocontamination->SetName("ratiocontamination");
1393         ratiocontamination->SetTitle((const char*)titlee);
1394         ratiocontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
1395         ratiocontamination->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1396         ratiocontamination->Add(aftersubstraction,-1.0);
1397         ratiocontamination->Divide(beforesubstraction);
1398         Int_t totalbin = ratiocontamination->GetXaxis()->GetNbins();
1399         for(Int_t nbinpt = 0; nbinpt < totalbin; nbinpt++) {
1400           ratiocontamination->SetBinError(nbinpt+1,0.0);
1401         }
1402         ratiocontamination->SetStats(0);
1403         ratiocontamination->SetMarkerStyle(26);
1404         ratiocontamination->SetMarkerColor(kBlack);
1405         ratiocontamination->SetLineColor(kBlack);
1406         ratiocontamination->Draw("P");
1407       }
1408
1409       cbackgrounde->cd(1);
1410       legtotal->Draw("same");
1411       cbackgrounde->cd(2);
1412       legtotalg->Draw("same");
1413
1414       cenaxisa->SetRange(0,nbbin);
1415       cenaxisb->SetRange(0,nbbin);
1416       if(fWriteToFile) cbackgrounde->SaveAs("BackgroundSubtractedPbPb.eps");
1417     }
1418   }
1419
1420   return spectrumSubtracted;
1421 }
1422
1423 //____________________________________________________________
1424 AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
1425   //
1426   // calculate charm background
1427   //
1428     Int_t ptpr = 0;
1429     Int_t nDim = 1;
1430     if(fBeamType==0)
1431     {
1432         ptpr=0;
1433     }
1434     if(fBeamType==1)
1435     {
1436         ptpr=1;
1437         nDim=2;
1438     }
1439
1440   Double_t evtnorm=0;
1441   if(fNMCbgEvents[0]) evtnorm= double(fNEvents[0])/double(fNMCbgEvents[0]);
1442
1443   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
1444   if(!mcContainer){
1445     AliError("MC Container not available");
1446     return NULL;
1447   }
1448
1449   if(!fCorrelation){
1450     AliError("No Correlation map available");
1451     return NULL;
1452   }
1453
1454   AliCFDataGrid *charmBackgroundGrid= 0x0;
1455   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
1456
1457   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
1458   Int_t* bins=new Int_t[2];
1459   bins[0]=charmbgaftertofpid->GetNbinsX();
1460
1461   if(fBeamType==1)
1462   {
1463       bins[0]=12;
1464       charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
1465       bins[1]=charmbgaftertofpid->GetNbinsX();
1466
1467   }
1468
1469   AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
1470   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
1471   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
1472
1473   charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
1474
1475   Int_t *nBinpp=new Int_t[1];
1476   Int_t* binspp=new Int_t[1];
1477   binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
1478
1479   Int_t *nBinPbPb=new Int_t[2];
1480   Int_t* binsPbPb=new Int_t[2];
1481   binsPbPb[1]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
1482   binsPbPb[0]=12;
1483
1484   Int_t looppt=binspp[0];
1485   if(fBeamType==1) looppt=binsPbPb[1];
1486
1487   for(Long_t iBin=1; iBin<= looppt;iBin++){
1488       if(fBeamType==0)
1489       {
1490           nBinpp[0]=iBin;
1491           charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
1492           charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
1493       }
1494       if(fBeamType==1)
1495       {
1496           // loop over centrality
1497           for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1498               nBinPbPb[0]=iiBin;
1499               nBinPbPb[1]=iBin;
1500               Double_t evtnormPbPb=0;
1501               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1502               charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1503               charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1504           }
1505       }
1506   }
1507
1508   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
1509
1510   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
1511   weightedCharmContainer->SetGrid(GetCharmWeights()); // get charm weighting factors
1512   TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(ptpr);
1513   charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
1514   TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(ptpr);
1515
1516   // Efficiency (set efficiency to 1 for only folding)
1517   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
1518   efficiencyD->CalculateEfficiency(0,0);
1519
1520   // Folding
1521   if(fBeamType==0)nDim = 1;
1522   if(fBeamType==1)nDim = 2;
1523   AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
1524   folding.SetMaxNumberOfIterations(1);
1525   folding.Unfold();
1526
1527   // Results
1528   THnSparse* result1= folding.GetEstMeasured(); // folded spectra
1529   THnSparse* result=(THnSparse*)result1->Clone();
1530   charmBackgroundGrid->SetGrid(result);
1531   TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(ptpr);
1532
1533   //Charm background evaluation plots
1534
1535   TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
1536   cCharmBgEval->Divide(3,1);
1537
1538   cCharmBgEval->cd(1);
1539
1540   if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
1541   if(fBeamType==1)
1542   {
1543       Double_t evtnormPbPb=0;
1544       for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
1545       {
1546           if(fNMCbgEvents[kCentr]) evtnormPbPb= evtnormPbPb+double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
1547       }
1548       charmbgaftertofpid->Scale(evtnormPbPb);
1549   }
1550
1551   CorrectFromTheWidth(charmbgaftertofpid);
1552   charmbgaftertofpid->SetMarkerStyle(25);
1553   charmbgaftertofpid->Draw("p");
1554   charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
1555   charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1556   gPad->SetLogy();
1557
1558   CorrectFromTheWidth(charmbgafteripcut);
1559   charmbgafteripcut->SetMarkerStyle(24);
1560   charmbgafteripcut->Draw("samep");
1561
1562   CorrectFromTheWidth(charmbgafterweight);
1563   charmbgafterweight->SetMarkerStyle(24);
1564   charmbgafterweight->SetMarkerColor(4);
1565   charmbgafterweight->Draw("samep");
1566
1567   CorrectFromTheWidth(charmbgafterfolding);
1568   charmbgafterfolding->SetMarkerStyle(24);
1569   charmbgafterfolding->SetMarkerColor(2);
1570   charmbgafterfolding->Draw("samep");
1571
1572   cCharmBgEval->cd(2);
1573   parametrizedcharmpidipeff->SetMarkerStyle(24);
1574   parametrizedcharmpidipeff->Draw("p");
1575   parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1576
1577   cCharmBgEval->cd(3);
1578   charmweightingfc->SetMarkerStyle(24);
1579   charmweightingfc->Draw("p");
1580   charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
1581   charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1582
1583   cCharmBgEval->cd(1);
1584   TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
1585   legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
1586   legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
1587   legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
1588   legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
1589   legcharmbg->Draw("same");
1590
1591   cCharmBgEval->cd(2);
1592   TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
1593   legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
1594   legcharmbg2->Draw("same");
1595
1596   CorrectStatErr(charmBackgroundGrid);
1597   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
1598
1599   delete[] bins;
1600   delete[] nBinpp;
1601   delete[] binspp;
1602   delete[] nBinPbPb;
1603   delete[] binsPbPb;
1604
1605   if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
1606
1607   return charmBackgroundGrid;
1608 }
1609
1610 //____________________________________________________________
1611 AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
1612   //
1613   // calculate conversion background
1614   //
1615
1616   Double_t evtnorm[1] = {0.0};
1617   if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1618   printf("check event!!! %lf \n",evtnorm[0]);
1619
1620   AliCFContainer *backgroundContainer = 0x0;
1621
1622   if(fNonHFEsyst){
1623     backgroundContainer = (AliCFContainer*)fConvSourceContainer[0][0][0]->Clone();
1624     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1625       backgroundContainer->Add(fConvSourceContainer[iSource][0][0]); // make centrality dependent
1626     }
1627   }
1628   else{
1629     // Background Estimate
1630     backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
1631   }
1632   if(!backgroundContainer){
1633     AliError("MC background container not found");
1634     return NULL;
1635   }
1636
1637   Int_t stepbackground = 3;
1638
1639   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
1640   Int_t *nBinpp=new Int_t[1];
1641   Int_t* binspp=new Int_t[1];
1642   binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1643
1644   Int_t *nBinPbPb=new Int_t[2];
1645   Int_t* binsPbPb=new Int_t[2];
1646   binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1647   binsPbPb[0]=12;
1648
1649   Int_t looppt=binspp[0];
1650   if(fBeamType==1) looppt=binsPbPb[1];
1651
1652   for(Long_t iBin=1; iBin<= looppt;iBin++){
1653       if(fBeamType==0)
1654       {
1655           nBinpp[0]=iBin;
1656           backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1657           backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1658       }
1659       if(fBeamType==1)
1660       {
1661           // loop over centrality
1662           for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
1663               nBinPbPb[0]=iiBin;
1664               nBinPbPb[1]=iBin;
1665               Double_t evtnormPbPb=0;
1666               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1667               backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1668               backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1669           }
1670       }
1671   }
1672   //end of workaround for statistical errors
1673
1674   AliCFDataGrid *weightedConversionContainer;
1675   if(fBeamType==0) weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,binspp);
1676   else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
1677   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
1678   backgroundGrid->Multiply(weightedConversionContainer,1.0);
1679
1680   delete[] nBinpp;
1681   delete[] binspp;
1682   delete[] nBinPbPb;
1683   delete[] binsPbPb;
1684
1685   return backgroundGrid;
1686 }
1687
1688
1689 //____________________________________________________________
1690 AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
1691   //
1692   // calculate non-HFE background
1693   //
1694
1695   Double_t evtnorm[1] = {0.0};
1696   if(fNMCbgEvents[0]) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
1697   printf("check event!!! %lf \n",evtnorm[0]);
1698
1699   AliCFContainer *backgroundContainer = 0x0;
1700   if(fNonHFEsyst){
1701     backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[0][0][0]->Clone();
1702     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
1703       if(iSource == 1)
1704         backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0],1.41);//correction for the eta Dalitz decay branching ratio in PYTHIA
1705       else
1706         backgroundContainer->Add(fNonHFESourceContainer[iSource][0][0]);
1707     }
1708   }
1709   else{
1710     // Background Estimate
1711     backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
1712   }
1713   if(!backgroundContainer){
1714     AliError("MC background container not found");
1715     return NULL;
1716   }
1717
1718
1719   Int_t stepbackground = 3;
1720
1721   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
1722   Int_t *nBinpp=new Int_t[1];
1723   Int_t* binspp=new Int_t[1];
1724   binspp[0]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1725
1726   Int_t *nBinPbPb=new Int_t[2];
1727   Int_t* binsPbPb=new Int_t[2];
1728   binsPbPb[1]=fConversionEff[0]->GetNbinsX();  // number of pt bins
1729   binsPbPb[0]=12;
1730
1731   Int_t looppt=binspp[0];
1732   if(fBeamType==1) looppt=binsPbPb[1];
1733
1734
1735   for(Long_t iBin=1; iBin<= looppt;iBin++){
1736       if(fBeamType==0)
1737       {
1738           nBinpp[0]=iBin;
1739           backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
1740           backgroundGrid->SetElement(nBinpp,backgroundGrid->GetElement(nBinpp)*evtnorm[0]);
1741       }
1742       if(fBeamType==1)
1743       {
1744           for(Long_t iiBin=1; iiBin<=binsPbPb[0];iiBin++){
1745               nBinPbPb[0]=iiBin;
1746               nBinPbPb[1]=iBin;
1747               Double_t evtnormPbPb=0;
1748               if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
1749               backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
1750               backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
1751           }
1752       }
1753   }
1754   //end of workaround for statistical errors
1755   AliCFDataGrid *weightedNonHFEContainer;
1756   if(fBeamType==0) weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,binspp);
1757   else weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",2,binsPbPb);
1758   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
1759   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
1760
1761   delete[] nBinpp;
1762   delete[] binspp;
1763   delete[] nBinPbPb;
1764   delete[] binsPbPb;
1765
1766   return backgroundGrid;
1767 }
1768
1769 //____________________________________________________________
1770 AliCFDataGrid *AliHFEspectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){
1771
1772   //
1773   // Apply TPC pid efficiency correction from parametrisation
1774   //
1775
1776   // Data in the right format
1777   AliCFDataGrid *dataGrid = 0x0;
1778   if(bgsubpectrum) {
1779     dataGrid = bgsubpectrum;
1780   }
1781   else {
1782
1783     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1784     if(!dataContainer){
1785       AliError("Data Container not available");
1786       return NULL;
1787     }
1788     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1789   }
1790   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1791   result->SetName("ParametrizedEfficiencyBefore");
1792   THnSparse *h = result->GetGrid();
1793   Int_t nbdimensions = h->GetNdimensions();
1794   //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
1795   AliCFContainer *dataContainer = GetContainer(kDataContainer);
1796   if(!dataContainer){
1797     AliError("Data Container not available");
1798     return NULL;
1799   }
1800   AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
1801   dataContainerbis->Add(dataContainerbis,-1.0);
1802
1803
1804   Int_t* coord = new Int_t[nbdimensions];
1805   memset(coord, 0, sizeof(Int_t) * nbdimensions);
1806   Double_t* points = new Double_t[nbdimensions];
1807
1808   ULong64_t nEntries = h->GetNbins();
1809   for (ULong64_t i = 0; i < nEntries; ++i) {
1810
1811     Double_t value = h->GetBinContent(i, coord);
1812     //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
1813     //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
1814
1815     // Get the bin co-ordinates given an coord
1816     for (Int_t j = 0; j < nbdimensions; ++j)
1817       points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
1818
1819     if (!fEfficiencyFunction->IsInside(points))
1820          continue;
1821     TF1::RejectPoint(kFALSE);
1822
1823     // Evaulate function at points
1824     Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
1825     //printf("Value efficiency is %f\n",valueEfficiency);
1826
1827     if(valueEfficiency > 0.0) {
1828       h->SetBinContent(coord,value/valueEfficiency);
1829       dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
1830     }
1831     Double_t error = h->GetBinError(i);
1832     h->SetBinError(coord,error/valueEfficiency);
1833     dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);
1834
1835
1836   }
1837
1838   delete[] coord;
1839   delete[] points;
1840
1841   AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
1842
1843   if(fDebugLevel > 0) {
1844
1845     TCanvas * cEfficiencyParametrized = new TCanvas("EfficiencyParametrized","EfficiencyParametrized",1000,700);
1846     cEfficiencyParametrized->Divide(2,1);
1847     cEfficiencyParametrized->cd(1);
1848     TH1D *afterE = (TH1D *) resultt->Project(0);
1849     TH1D *beforeE = (TH1D *) dataGrid->Project(0);
1850     CorrectFromTheWidth(afterE);
1851     CorrectFromTheWidth(beforeE);
1852     afterE->SetStats(0);
1853     afterE->SetTitle("");
1854     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1855     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1856     afterE->SetMarkerStyle(25);
1857     afterE->SetMarkerColor(kBlack);
1858     afterE->SetLineColor(kBlack);
1859     beforeE->SetStats(0);
1860     beforeE->SetTitle("");
1861     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1862     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1863     beforeE->SetMarkerStyle(24);
1864     beforeE->SetMarkerColor(kBlue);
1865     beforeE->SetLineColor(kBlue);
1866     gPad->SetLogy();
1867     afterE->Draw();
1868     beforeE->Draw("same");
1869     TLegend *legefficiencyparametrized = new TLegend(0.4,0.6,0.89,0.89);
1870     legefficiencyparametrized->AddEntry(beforeE,"Before Efficiency correction","p");
1871     legefficiencyparametrized->AddEntry(afterE,"After Efficiency correction","p");
1872     legefficiencyparametrized->Draw("same");
1873     cEfficiencyParametrized->cd(2);
1874     fEfficiencyFunction->Draw();
1875     //cEfficiencyParametrized->cd(3);
1876     //TH1D *ratioefficiency = (TH1D *) beforeE->Clone();
1877     //ratioefficiency->Divide(afterE);
1878     //ratioefficiency->Draw();
1879
1880     if(fWriteToFile) cEfficiencyParametrized->SaveAs("efficiency.eps");
1881
1882   }
1883
1884   return resultt;
1885
1886 }
1887 //____________________________________________________________
1888 AliCFDataGrid *AliHFEspectrum::CorrectV0Efficiency(AliCFDataGrid* const bgsubpectrum){
1889
1890   //
1891   // Apply TPC pid efficiency correction from V0
1892   //
1893
1894   AliCFContainer *v0Container = GetContainer(kDataContainerV0);
1895   if(!v0Container){
1896     AliError("V0 Container not available");
1897     return NULL;
1898   }
1899
1900   // Efficiency
1901   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*v0Container);
1902   efficiencyD->CalculateEfficiency(fStepAfterCutsV0,fStepBeforeCutsV0);
1903
1904   // Data in the right format
1905   AliCFDataGrid *dataGrid = 0x0;
1906   if(bgsubpectrum) {
1907     dataGrid = bgsubpectrum;
1908   }
1909   else {
1910
1911     AliCFContainer *dataContainer = GetContainer(kDataContainer);
1912     if(!dataContainer){
1913       AliError("Data Container not available");
1914       return NULL;
1915     }
1916
1917     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
1918   }
1919
1920   // Correct
1921   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
1922   result->ApplyEffCorrection(*efficiencyD);
1923
1924   if(fDebugLevel > 0) {
1925
1926     Int_t ptpr;
1927     if(fBeamType==0) ptpr=0;
1928     if(fBeamType==1) ptpr=1;
1929
1930     TCanvas * cV0Efficiency = new TCanvas("V0Efficiency","V0Efficiency",1000,700);
1931     cV0Efficiency->Divide(2,1);
1932     cV0Efficiency->cd(1);
1933     TH1D *afterE = (TH1D *) result->Project(ptpr);
1934     TH1D *beforeE = (TH1D *) dataGrid->Project(ptpr);
1935     afterE->SetStats(0);
1936     afterE->SetTitle("");
1937     afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1938     afterE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1939     afterE->SetMarkerStyle(25);
1940     afterE->SetMarkerColor(kBlack);
1941     afterE->SetLineColor(kBlack);
1942     beforeE->SetStats(0);
1943     beforeE->SetTitle("");
1944     beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1945     beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1946     beforeE->SetMarkerStyle(24);
1947     beforeE->SetMarkerColor(kBlue);
1948     beforeE->SetLineColor(kBlue);
1949     afterE->Draw();
1950     beforeE->Draw("same");
1951     TLegend *legV0efficiency = new TLegend(0.4,0.6,0.89,0.89);
1952     legV0efficiency->AddEntry(beforeE,"Before Efficiency correction","p");
1953     legV0efficiency->AddEntry(afterE,"After Efficiency correction","p");
1954     legV0efficiency->Draw("same");
1955     cV0Efficiency->cd(2);
1956     TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(ptpr);
1957     efficiencyDproj->SetTitle("");
1958     efficiencyDproj->SetStats(0);
1959     efficiencyDproj->SetMarkerStyle(25);
1960     efficiencyDproj->Draw();
1961
1962     if(fBeamType==1) {
1963
1964       TCanvas * cV0Efficiencye = new TCanvas("V0Efficiency_allspectra","V0Efficiency_allspectra",1000,700);
1965       cV0Efficiencye->Divide(2,1);
1966       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
1967       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
1968
1969       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
1970       TAxis *cenaxisa = sparseafter->GetAxis(0);
1971       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
1972       TAxis *cenaxisb = sparsebefore->GetAxis(0);
1973       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
1974       TAxis *cenaxisc = efficiencya->GetAxis(0);
1975       Int_t nbbin = cenaxisb->GetNbins();
1976       Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
1977       Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
1978       for(Int_t binc = 0; binc < nbbin; binc++){
1979         TString titlee("V0Efficiency_centrality_bin_");
1980         titlee += binc;
1981         TCanvas * ccV0Efficiency = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
1982         ccV0Efficiency->Divide(2,1);
1983         ccV0Efficiency->cd(1);
1984         gPad->SetLogy();
1985         cenaxisa->SetRange(binc+1,binc+1);
1986         cenaxisb->SetRange(binc+1,binc+1);
1987         cenaxisc->SetRange(binc+1,binc+1);
1988         TH1D *aftere = (TH1D *) sparseafter->Projection(1);
1989         TH1D *beforee = (TH1D *) sparsebefore->Projection(1);
1990         CorrectFromTheWidth(aftere);
1991         CorrectFromTheWidth(beforee);
1992         aftere->SetStats(0);
1993         aftere->SetTitle((const char*)titlee);
1994         aftere->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
1995         aftere->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1996         aftere->SetMarkerStyle(25);
1997         aftere->SetMarkerColor(kBlack);
1998         aftere->SetLineColor(kBlack);
1999         beforee->SetStats(0);
2000         beforee->SetTitle((const char*)titlee);
2001         beforee->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2002         beforee->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2003         beforee->SetMarkerStyle(24);
2004         beforee->SetMarkerColor(kBlue);
2005         beforee->SetLineColor(kBlue);
2006         aftere->Draw();
2007         beforee->Draw("same");
2008         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2009         lega->AddEntry(beforee,"Before correction","p");
2010         lega->AddEntry(aftere,"After correction","p");
2011         lega->Draw("same");
2012         cV0Efficiencye->cd(1);
2013         gPad->SetLogy();
2014         TH1D *afteree = (TH1D *) aftere->Clone();
2015         afteree->SetMarkerStyle(stylee[binc]);
2016         afteree->SetMarkerColor(colorr[binc]);
2017         if(binc==0) afteree->Draw();
2018         else afteree->Draw("same");
2019         legtotal->AddEntry(afteree,(const char*) titlee,"p");
2020         cV0Efficiencye->cd(2);
2021         gPad->SetLogy();
2022         TH1D *aftereeu = (TH1D *) aftere->Clone();
2023         aftereeu->SetMarkerStyle(stylee[binc]);
2024         aftereeu->SetMarkerColor(colorr[binc]);
2025         if(fNEvents[binc] > 0.0) aftereeu->Scale(1/(Double_t)fNEvents[binc]);
2026         if(binc==0) aftereeu->Draw();
2027         else aftereeu->Draw("same");
2028         legtotalg->AddEntry(aftereeu,(const char*) titlee,"p");
2029         ccV0Efficiency->cd(2);
2030         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(1);
2031         efficiencyDDproj->SetTitle("");
2032         efficiencyDDproj->SetStats(0);
2033         efficiencyDDproj->SetMarkerStyle(25);
2034         efficiencyDDproj->Draw();
2035       }
2036
2037       cV0Efficiencye->cd(1);
2038       legtotal->Draw("same");
2039       cV0Efficiencye->cd(2);
2040       legtotalg->Draw("same");
2041
2042       cenaxisa->SetRange(0,nbbin);
2043       cenaxisb->SetRange(0,nbbin);
2044       cenaxisc->SetRange(0,nbbin);
2045
2046       if(fWriteToFile) cV0Efficiencye->SaveAs("V0efficiency.eps");
2047     }
2048
2049   }
2050
2051
2052   return result;
2053
2054 }
2055 //____________________________________________________________
2056 TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
2057
2058   //
2059   // Unfold and eventually correct for efficiency the bgsubspectrum
2060   //
2061
2062   AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
2063   if(!mcContainer){
2064     AliError("MC Container not available");
2065     return NULL;
2066   }
2067
2068   if(!fCorrelation){
2069     AliError("No Correlation map available");
2070     return NULL;
2071   }
2072
2073   // Data
2074   AliCFDataGrid *dataGrid = 0x0;
2075   if(bgsubpectrum) {
2076     dataGrid = bgsubpectrum;
2077   }
2078   else {
2079
2080     AliCFContainer *dataContainer = GetContainer(kDataContainer);
2081     if(!dataContainer){
2082       AliError("Data Container not available");
2083       return NULL;
2084     }
2085
2086     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2087   }
2088
2089   // Guessed
2090   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
2091   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
2092
2093   // Efficiency
2094   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2095   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2096
2097   if(!fBeauty2ndMethod)
2098   {
2099       // Consider parameterized IP cut efficiency
2100       if(!fInclusiveSpectrum){
2101           Int_t* bins=new Int_t[1];
2102           bins[0]=35;
2103
2104           AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2105           beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
2106           efficiencyD->Multiply(beffContainer,1);
2107       }
2108   }
2109
2110
2111   // Unfold
2112
2113   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
2114   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
2115   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
2116   if(fNRandomIter > 0) unfolding.SetNRandomIterations(fNRandomIter);
2117   if(fSetSmoothing) unfolding.UseSmoothing();
2118   unfolding.Unfold();
2119
2120   // Results
2121   THnSparse* result = unfolding.GetUnfolded();
2122   THnSparse* residual = unfolding.GetEstMeasured();
2123   TList *listofresults = new TList;
2124   listofresults->SetOwner();
2125   listofresults->AddAt((THnSparse*)result->Clone(),0);
2126   listofresults->AddAt((THnSparse*)residual->Clone(),1);
2127
2128   if(fDebugLevel > 0) {
2129
2130     Int_t ptpr;
2131     if(fBeamType==0) ptpr=0;
2132     if(fBeamType==1) ptpr=1;
2133
2134     TCanvas * cresidual = new TCanvas("residual","residual",1000,700);
2135     cresidual->Divide(2,1);
2136     cresidual->cd(1);
2137     gPad->SetLogy();
2138     TGraphErrors* residualspectrumD = Normalize(residual);
2139     if(!residualspectrumD) {
2140       AliError("Number of Events not set for the normalization");
2141       return NULL;
2142     }
2143     residualspectrumD->SetTitle("");
2144     residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
2145     residualspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2146     residualspectrumD->SetMarkerStyle(26);
2147     residualspectrumD->SetMarkerColor(kBlue);
2148     residualspectrumD->SetLineColor(kBlue);
2149     residualspectrumD->Draw("AP");
2150     AliCFDataGrid *dataGridBis = (AliCFDataGrid *) (dataGrid->Clone());
2151     dataGridBis->SetName("dataGridBis");
2152     TGraphErrors* measuredspectrumD = Normalize(dataGridBis);
2153     measuredspectrumD->SetTitle("");
2154     measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
2155     measuredspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
2156     measuredspectrumD->SetMarkerStyle(25);
2157     measuredspectrumD->SetMarkerColor(kBlack);
2158     measuredspectrumD->SetLineColor(kBlack);
2159     measuredspectrumD->Draw("P");
2160     TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
2161     legres->AddEntry(residualspectrumD,"Residual","p");
2162     legres->AddEntry(measuredspectrumD,"Measured","p");
2163     legres->Draw("same");
2164     cresidual->cd(2);
2165     TH1D *residualTH1D = residual->Projection(ptpr);
2166     TH1D *measuredTH1D = (TH1D *) dataGridBis->Project(ptpr);
2167     TH1D* ratioresidual = (TH1D*)residualTH1D->Clone();
2168     ratioresidual->SetName("ratioresidual");
2169     ratioresidual->SetTitle("");
2170     ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
2171     ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2172     ratioresidual->Divide(residualTH1D,measuredTH1D,1,1);
2173     ratioresidual->SetStats(0);
2174     ratioresidual->Draw();
2175
2176     if(fWriteToFile) cresidual->SaveAs("Unfolding.eps");
2177   }
2178
2179   return listofresults;
2180
2181 }
2182
2183 //____________________________________________________________
2184 AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
2185
2186   //
2187   // Apply unfolding and efficiency correction together to bgsubspectrum
2188   //
2189
2190   AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
2191   if(!mcContainer){
2192     AliError("MC Container not available");
2193     return NULL;
2194   }
2195
2196   // Efficiency
2197   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
2198   efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
2199
2200
2201   if(!fBeauty2ndMethod)
2202   {
2203       // Consider parameterized IP cut efficiency
2204       if(!fInclusiveSpectrum){
2205           Int_t* bins=new Int_t[1];
2206           bins[0]=35;
2207
2208           AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
2209           beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
2210           efficiencyD->Multiply(beffContainer,1);
2211       }
2212   }
2213
2214   // Data in the right format
2215   AliCFDataGrid *dataGrid = 0x0;
2216   if(bgsubpectrum) {
2217     dataGrid = bgsubpectrum;
2218   }
2219   else {
2220
2221     AliCFContainer *dataContainer = GetContainer(kDataContainer);
2222     if(!dataContainer){
2223       AliError("Data Container not available");
2224       return NULL;
2225     }
2226
2227     dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
2228   }
2229
2230   // Correct
2231   AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
2232   result->ApplyEffCorrection(*efficiencyD);
2233
2234   if(fDebugLevel > 0) {
2235
2236     Int_t ptpr;
2237     if(fBeamType==0) ptpr=0;
2238     if(fBeamType==1) ptpr=1;
2239
2240     printf("Step MC: %d\n",fStepMC);
2241     printf("Step tracking: %d\n",AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2242     printf("Step MC true: %d\n",fStepTrue);
2243     AliCFEffGrid  *efficiencymcPID = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
2244     AliCFEffGrid  *efficiencymctrackinggeo = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack,fStepTrue);
2245     AliCFEffGrid  *efficiencymcall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerMC),fStepMC,fStepTrue);
2246
2247     AliCFEffGrid  *efficiencyesdall = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCContainerESD),fStepMC,fStepTrue);
2248
2249     TCanvas * cefficiency = new TCanvas("efficiency","efficiency",1000,700);
2250     cefficiency->cd(1);
2251     TH1D* efficiencymcPIDD = (TH1D *) efficiencymcPID->Project(ptpr);
2252     efficiencymcPIDD->SetTitle("");
2253     efficiencymcPIDD->SetStats(0);
2254     efficiencymcPIDD->SetMarkerStyle(25);
2255     efficiencymcPIDD->Draw();
2256     TH1D* efficiencymctrackinggeoD = (TH1D *) efficiencymctrackinggeo->Project(ptpr);
2257     efficiencymctrackinggeoD->SetTitle("");
2258     efficiencymctrackinggeoD->SetStats(0);
2259     efficiencymctrackinggeoD->SetMarkerStyle(26);
2260     efficiencymctrackinggeoD->Draw("same");
2261     TH1D* efficiencymcallD = (TH1D *) efficiencymcall->Project(ptpr);
2262     efficiencymcallD->SetTitle("");
2263     efficiencymcallD->SetStats(0);
2264     efficiencymcallD->SetMarkerStyle(27);
2265     efficiencymcallD->Draw("same");
2266     TH1D* efficiencyesdallD = (TH1D *) efficiencyesdall->Project(ptpr);
2267     efficiencyesdallD->SetTitle("");
2268     efficiencyesdallD->SetStats(0);
2269     efficiencyesdallD->SetMarkerStyle(24);
2270     efficiencyesdallD->Draw("same");
2271     TLegend *legeff = new TLegend(0.4,0.6,0.89,0.89);
2272     legeff->AddEntry(efficiencymcPIDD,"PID efficiency","p");
2273     legeff->AddEntry(efficiencymctrackinggeoD,"Tracking geometry efficiency","p");
2274     legeff->AddEntry(efficiencymcallD,"Overall efficiency","p");
2275     legeff->AddEntry(efficiencyesdallD,"Overall efficiency ESD","p");
2276     legeff->Draw("same");
2277
2278     if(fBeamType==1) {
2279
2280       THnSparseF* sparseafter = (THnSparseF *) result->GetGrid();
2281       TAxis *cenaxisa = sparseafter->GetAxis(0);
2282       THnSparseF* sparsebefore = (THnSparseF *) dataGrid->GetGrid();
2283       TAxis *cenaxisb = sparsebefore->GetAxis(0);
2284       THnSparseF* efficiencya = (THnSparseF *) efficiencyD->GetGrid();
2285       TAxis *cenaxisc = efficiencya->GetAxis(0);
2286       //Int_t nbbin = cenaxisb->GetNbins();
2287       //Int_t stylee[20] = {20,21,22,23,24,25,26,27,28,30,4,5,7,29,29,29,29,29,29,29};
2288       //Int_t colorr[20] = {2,3,4,5,6,7,8,9,46,38,29,30,31,32,33,34,35,37,38,20};
2289       for(Int_t binc = 0; binc < fNCentralityBinAtTheEnd; binc++){
2290         TString titlee("Efficiency_centrality_bin_");
2291         titlee += fLowBoundaryCentralityBinAtTheEnd[binc];
2292         titlee += "_";
2293         titlee += fHighBoundaryCentralityBinAtTheEnd[binc];
2294         TCanvas * cefficiencye = new TCanvas((const char*) titlee,(const char*) titlee,1000,700);
2295         cefficiencye->Divide(2,1);
2296         cefficiencye->cd(1);
2297         gPad->SetLogy();
2298         cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2299         cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
2300         TH1D *afterefficiency = (TH1D *) sparseafter->Projection(ptpr);
2301         TH1D *beforeefficiency = (TH1D *) sparsebefore->Projection(ptpr);
2302         CorrectFromTheWidth(afterefficiency);
2303         CorrectFromTheWidth(beforeefficiency);
2304         afterefficiency->SetStats(0);
2305         afterefficiency->SetTitle((const char*)titlee);
2306         afterefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2307         afterefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2308         afterefficiency->SetMarkerStyle(25);
2309         afterefficiency->SetMarkerColor(kBlack);
2310         afterefficiency->SetLineColor(kBlack);
2311         beforeefficiency->SetStats(0);
2312         beforeefficiency->SetTitle((const char*)titlee);
2313         beforeefficiency->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
2314         beforeefficiency->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2315         beforeefficiency->SetMarkerStyle(24);
2316         beforeefficiency->SetMarkerColor(kBlue);
2317         beforeefficiency->SetLineColor(kBlue);
2318         afterefficiency->Draw();
2319         beforeefficiency->Draw("same");
2320         TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
2321         lega->AddEntry(beforeefficiency,"Before efficiency correction","p");
2322         lega->AddEntry(afterefficiency,"After efficiency correction","p");
2323         lega->Draw("same");
2324         cefficiencye->cd(2);
2325         cenaxisc->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fLowBoundaryCentralityBinAtTheEnd[binc]+1);
2326         TH1D* efficiencyDDproj = (TH1D *) efficiencya->Projection(ptpr);
2327         efficiencyDDproj->SetTitle("");
2328         efficiencyDDproj->SetStats(0);
2329         efficiencyDDproj->SetMarkerStyle(25);
2330         efficiencyDDproj->SetMarkerColor(2);
2331         efficiencyDDproj->Draw();
2332         cenaxisc->SetRange(fHighBoundaryCentralityBinAtTheEnd[binc],fHighBoundaryCentralityBinAtTheEnd[binc]);
2333         TH1D* efficiencyDDproja = (TH1D *) efficiencya->Projection(ptpr);
2334         efficiencyDDproja->SetTitle("");
2335         efficiencyDDproja->SetStats(0);
2336         efficiencyDDproja->SetMarkerStyle(26);
2337         efficiencyDDproja->SetMarkerColor(4);
2338         efficiencyDDproja->Draw("same");
2339       }
2340     }
2341
2342     if(fWriteToFile) cefficiency->SaveAs("efficiencyCorrected.eps");
2343   }
2344
2345   return result;
2346
2347 }
2348
2349 //__________________________________________________________________________________
2350 TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) const {
2351   //
2352   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2353   // Give the final pt spectrum to be compared
2354   //
2355
2356   if(fNEvents[i] > 0) {
2357
2358     Int_t ptpr = 0;
2359     if(fBeamType==0) ptpr=0;
2360     if(fBeamType==1) ptpr=1;
2361
2362     TH1D* projection = spectrum->Projection(ptpr);
2363     CorrectFromTheWidth(projection);
2364     TGraphErrors *graphError = NormalizeTH1(projection,i);
2365     return graphError;
2366
2367   }
2368
2369   return 0x0;
2370
2371
2372 }
2373 //__________________________________________________________________________________
2374 TGraphErrors *AliHFEspectrum::Normalize(AliCFDataGrid * const spectrum,Int_t i) const {
2375   //
2376   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2377   // Give the final pt spectrum to be compared
2378   //
2379   if(fNEvents[i] > 0) {
2380
2381     Int_t ptpr=0;
2382     if(fBeamType==0) ptpr=0;
2383     if(fBeamType==1) ptpr=1;
2384
2385     TH1D* projection = (TH1D *) spectrum->Project(ptpr);
2386     CorrectFromTheWidth(projection);
2387     TGraphErrors *graphError = NormalizeTH1(projection,i);
2388
2389     return graphError;
2390
2391   }
2392
2393   return 0x0;
2394
2395
2396 }
2397 //__________________________________________________________________________________
2398 TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
2399   //
2400   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2401   // Give the final pt spectrum to be compared
2402   //
2403   Double_t chargecoefficient = 0.5;
2404   if(fChargeChoosen != 0) chargecoefficient = 1.0;
2405
2406   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2407   printf("Normalizing Eta Range %f\n", etarange);
2408   if(fNEvents[i] > 0) {
2409
2410     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2411     Double_t p = 0, dp = 0; Int_t point = 1;
2412     Double_t n = 0, dN = 0;
2413     Double_t nCorr = 0, dNcorr = 0;
2414     //Double_t errdN = 0, errdp = 0;
2415     Double_t errdN = 0;
2416     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2417       point = ibin - input->GetXaxis()->GetFirst();
2418       p = input->GetXaxis()->GetBinCenter(ibin);
2419       dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2420       n = input->GetBinContent(ibin);
2421       dN = input->GetBinError(ibin);
2422       // New point
2423       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
2424       errdN = 1./(2. * TMath::Pi() * p);
2425       //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2426       //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2427       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN);
2428
2429       spectrumNormalized->SetPoint(point, p, nCorr);
2430       spectrumNormalized->SetPointError(point, dp, dNcorr);
2431     }
2432     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2433     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2434     spectrumNormalized->SetMarkerStyle(22);
2435     spectrumNormalized->SetMarkerColor(kBlue);
2436     spectrumNormalized->SetLineColor(kBlue);
2437
2438     return spectrumNormalized;
2439
2440   }
2441
2442   return 0x0;
2443
2444
2445 }
2446 //__________________________________________________________________________________
2447 TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) const {
2448   //
2449   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
2450   // Give the final pt spectrum to be compared
2451   //
2452   Double_t chargecoefficient = 0.5;
2453   if(fChargeChoosen != kAllCharge) chargecoefficient = 1.0;
2454
2455   Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
2456   printf("Normalizing Eta Range %f\n", etarange);
2457   if(normalization > 0) {
2458
2459     TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
2460     Double_t p = 0, dp = 0; Int_t point = 1;
2461     Double_t n = 0, dN = 0;
2462     Double_t nCorr = 0, dNcorr = 0;
2463     //Double_t errdN = 0, errdp = 0;
2464     Double_t errdN = 0;
2465     for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
2466       point = ibin - input->GetXaxis()->GetFirst();
2467       p = input->GetXaxis()->GetBinCenter(ibin);
2468       //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
2469       n = input->GetBinContent(ibin);
2470       dN = input->GetBinError(ibin);
2471       // New point
2472       nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
2473       errdN = 1./(2. * TMath::Pi() * p);
2474       //errdp = 1./(2. * TMath::Pi() * p*p) * n;
2475       //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
2476       dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN);
2477
2478       spectrumNormalized->SetPoint(point, p, nCorr);
2479       spectrumNormalized->SetPointError(point, dp, dNcorr);
2480     }
2481     spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
2482     spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
2483     spectrumNormalized->SetMarkerStyle(22);
2484     spectrumNormalized->SetMarkerColor(kBlue);
2485     spectrumNormalized->SetLineColor(kBlue);
2486
2487     return spectrumNormalized;
2488
2489   }
2490
2491   return 0x0;
2492
2493
2494 }
2495 //____________________________________________________________
2496 void AliHFEspectrum::SetContainer(AliCFContainer *cont, AliHFEspectrum::CFContainer_t type){
2497   //
2498   // Set the container for a given step to the
2499   //
2500   if(!fCFContainers) fCFContainers = new TObjArray(kDataContainerV0+1);
2501   fCFContainers->AddAt(cont, type);
2502 }
2503
2504 //____________________________________________________________
2505 AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type){
2506   //
2507   // Get Correction Framework Container for given type
2508   //
2509   if(!fCFContainers) return NULL;
2510   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
2511 }
2512 //____________________________________________________________
2513 AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
2514   //
2515   // Slice bin for a given source of electron
2516   // nDim is the number of dimension the corrections are done
2517   // dimensions are the definition of the dimensions
2518   // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
2519   // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
2520   // centrality (-1 means we do not cut on centrality)
2521   //
2522
2523   Double_t *varMin = new Double_t[container->GetNVar()],
2524            *varMax = new Double_t[container->GetNVar()];
2525
2526   Double_t *binLimits;
2527   for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
2528
2529     binLimits = new Double_t[container->GetNBins(ivar)+1];
2530     container->GetBinLimits(ivar,binLimits);
2531     varMin[ivar] = binLimits[0];
2532     varMax[ivar] = binLimits[container->GetNBins(ivar)];
2533     // source
2534     if(ivar == 4){
2535       if((source>= 0) && (source<container->GetNBins(ivar))) {
2536               varMin[ivar] = container->GetAxis(4,0)->GetBinLowEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2537               varMax[ivar] = container->GetAxis(4,0)->GetBinUpEdge(container->GetAxis(4,0)->FindBin(binLimits[source]));
2538       }
2539     }
2540     // charge
2541     if(ivar == 3) {
2542       if(charge != kAllCharge){
2543         varMin[ivar] = container->GetAxis(3,0)->GetBinLowEdge(container->GetAxis(3,0)->FindBin(charge));
2544         varMax[ivar] = container->GetAxis(3,0)->GetBinUpEdge(container->GetAxis(3,0)->FindBin(charge));
2545       }
2546     }
2547     // eta
2548     if(ivar == 1){
2549       for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
2550         AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
2551       if(fEtaSelected){
2552         varMin[ivar] = fEtaRange[0];
2553         varMax[ivar] = fEtaRange[1];
2554       }
2555     }
2556     if(fEtaSelected){
2557       fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
2558       fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
2559       AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
2560     }
2561     if(ivar == 5){
2562         if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
2563             varMin[ivar] = binLimits[centralitylow];
2564             varMax[ivar] = binLimits[centralityhigh];
2565
2566             TAxis *axistest = container->GetAxis(5,0);
2567             printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
2568             printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
2569             Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
2570             Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
2571             printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
2572
2573         }
2574
2575     }
2576
2577
2578     delete[] binLimits;
2579   }
2580
2581   AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
2582   AddTemporaryObject(k);
2583   delete[] varMin; delete[] varMax;
2584
2585   return k;
2586
2587 }
2588
2589 //_________________________________________________________________________
2590 THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
2591   //
2592   // Slice correlation
2593   //
2594
2595   Int_t ndimensions = correlationmatrix->GetNdimensions();
2596   //printf("Number of dimension %d correlation map\n",ndimensions);
2597   if(ndimensions < (2*nDim)) {
2598     AliError("Problem in the dimensions");
2599     return NULL;
2600   }
2601
2602   // Cut in centrality is centrality > -1
2603   if((centralitylow >=0) && (centralityhigh >=0)) {
2604
2605     TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
2606     TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
2607
2608     Int_t bins0 = axiscentrality0->GetNbins();
2609     Int_t bins1 = axiscentrality1->GetNbins();
2610
2611     printf("Number of centrality bins: %d and %d\n",bins0,bins1);
2612     if(bins0 != bins1) {
2613       AliError("Problem in the dimensions");
2614       return NULL;
2615     }
2616
2617     if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
2618       axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
2619       axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
2620
2621       Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
2622       Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
2623       Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
2624       Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
2625       printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
2626       printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
2627
2628       TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
2629       ctestcorrelation->cd(1);
2630       TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
2631       projection->Draw("colz");
2632
2633     }
2634
2635   }
2636
2637
2638   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
2639   //printf("Number of dimension %d container\n",ndimensionsContainer);
2640
2641   Int_t *dim = new Int_t[nDim*2];
2642   for(Int_t iter=0; iter < nDim; iter++){
2643     dim[iter] = dimensions[iter];
2644     dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
2645     //printf("For iter %d: %d and iter+nDim %d: %d\n",iter,dimensions[iter],iter+nDim,ndimensionsContainer + dimensions[iter]);
2646   }
2647
2648   THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);
2649
2650   delete[] dim;
2651   return k;
2652
2653 }
2654 //___________________________________________________________________________
2655 void AliHFEspectrum::CorrectFromTheWidth(TH1D *h1) const {
2656   //
2657   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
2658   //
2659
2660   TAxis *axis = h1->GetXaxis();
2661   Int_t nbinX = h1->GetNbinsX();
2662
2663   for(Int_t i = 1; i <= nbinX; i++) {
2664
2665     Double_t width = axis->GetBinWidth(i);
2666     Double_t content = h1->GetBinContent(i);
2667     Double_t error = h1->GetBinError(i);
2668     h1->SetBinContent(i,content/width);
2669     h1->SetBinError(i,error/width);
2670   }
2671
2672 }
2673
2674 //___________________________________________________________________________
2675 void AliHFEspectrum::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
2676   //
2677   // Correct statistical error
2678   //
2679
2680   TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
2681   Int_t nbinX = h1->GetNbinsX();
2682   Int_t bins[1];
2683   for(Long_t i = 1; i <= nbinX; i++) {
2684     bins[0] = i;
2685     Float_t content = h1->GetBinContent(i);
2686     if(content>0){
2687       Float_t error = TMath::Sqrt(content);
2688       backgroundGrid->SetElementError(bins, error);
2689     }
2690   }
2691 }
2692
2693 //____________________________________________________________
2694 void AliHFEspectrum::AddTemporaryObject(TObject *o){
2695   //
2696   // Emulate garbage collection on class level
2697   //
2698   if(!fTemporaryObjects) {
2699     fTemporaryObjects= new TList;
2700     fTemporaryObjects->SetOwner();
2701   }
2702   fTemporaryObjects->Add(o);
2703 }
2704
2705 //____________________________________________________________
2706 void AliHFEspectrum::ClearObject(TObject *o){
2707   //
2708   // Do a safe deletion
2709   //
2710   if(fTemporaryObjects){
2711     if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
2712     delete o;
2713   } else{
2714     // just delete
2715     delete o;
2716   }
2717 }
2718 //_________________________________________________________________________
2719 TObject* AliHFEspectrum::GetSpectrum(const AliCFContainer * const c, Int_t step) {
2720   AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
2721   return data;
2722 }
2723 //_________________________________________________________________________
2724 TObject* AliHFEspectrum::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
2725   //
2726   // Create efficiency grid and calculate efficiency
2727   // of step to step0
2728   //
2729   TString name("eff");
2730   name += step;
2731   name+= step0;
2732   AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
2733   eff->CalculateEfficiency(step,step0);
2734   return eff;
2735 }
2736
2737 //____________________________________________________________________________
2738 THnSparse* AliHFEspectrum::GetCharmWeights(){
2739
2740   //
2741   // Measured D->e based weighting factors
2742   //
2743
2744   const Int_t nDimpp=1;
2745   Int_t nBinpp[nDimpp] = {35};
2746   Double_t ptbinning1[36] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};
2747   const Int_t nDimPbPb=2;
2748   Int_t nBinPbPb[nDimPbPb] = {11,35};
2749   Double_t kCentralityRange[12] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
2750   Int_t loopcentr=1;
2751   Int_t looppt=nBinpp[0];
2752   if(fBeamType==0)
2753   {
2754       fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
2755       fWeightCharm->SetBinEdges(0, ptbinning1);
2756   }
2757   if(fBeamType==1)
2758   {
2759       fWeightCharm = new THnSparseF("weightHisto", "weighting factor; centrality bin; pt[GeV/c]", nDimPbPb, nBinPbPb);
2760       fWeightCharm->SetBinEdges(1, ptbinning1);
2761       fWeightCharm->SetBinEdges(0, kCentralityRange);
2762       loopcentr=nBinPbPb[0];
2763   }
2764
2765   // Weighting factor for pp
2766   Double_t weight[35]={0.859260, 0.872552, 0.847475, 0.823631, 0.839386, 0.874024, 0.916755, 0.942801, 0.965856, 0.933905, 0.933414, 0.931936, 0.847826, 0.810902, 0.796608, 0.727002, 0.659227, 0.583610, 0.549956, 0.512633, 0.472254, 0.412364, 0.353191, 0.319145, 0.305412, 0.290334, 0.269863, 0.254646, 0.230245, 0.200859, 0.275953, 0.276271, 0.227332, 0.197004, 0.474385};
2767
2768   // Weighting factor for PbPb (0-20%)
2769   //Double_t weight[35]={0.641897,  0.640472,  0.615228,  0.650469,  0.737762,  0.847867,  1.009317,  1.158594,  1.307482,  1.476973,  1.551131,  1.677131,  1.785478,  1.888933,  2.017957,  2.074757,  1.926700,  1.869495,  1.546558,  1.222873,  1.160313,  0.903375,  0.799642,  0.706244,  0.705449,  0.599947,  0.719570,  0.499422,  0.703978,  0.477452,  0.325057,  0.093391,  0.096675,  0.000000,  0.000000};
2770
2771   // Weighting factor for PbPb (40-80%)
2772   //Double_t weight[35]={0.181953,  0.173248,  0.166799,  0.182558,  0.206581,  0.236955,  0.279390,  0.329129,  0.365260,  0.423059,  0.452057,  0.482726,  0.462627,  0.537770,  0.584663,  0.579452,  0.587194,  0.499498,  0.443299,  0.398596,  0.376695,  0.322331,  0.260890,  0.374834,  0.249114,  0.310330,  0.287326,  0.243174,  0.758945,  0.138867,  0.170576,  0.107797,  0.011390,  0.000000,  0.000000};
2773
2774   //points
2775   Double_t pt[1];
2776   Double_t contents[2];
2777
2778   for(int icentr=0; icentr<loopcentr; icentr++)
2779   {
2780       for(int i=0; i<looppt; i++){
2781           pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
2782           if(fBeamType==1)
2783           {
2784               contents[0]=icentr;
2785               contents[1]=pt[0];
2786           }
2787           if(fBeamType==0)
2788           {
2789               contents[0]=pt[0];
2790           }
2791
2792           fWeightCharm->Fill(contents,weight[i]);
2793       }
2794   }
2795
2796   Int_t nDimSparse = fWeightCharm->GetNdimensions();
2797   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
2798   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
2799
2800   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2801       binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
2802       nBins *= binsvar[iVar];
2803   }
2804
2805   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
2806   // loop that sets 0 error in each bin
2807   for (Long_t iBin=0; iBin<nBins; iBin++) {
2808     Long_t bintmp = iBin ;
2809     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
2810       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
2811       bintmp /= binsvar[iVar] ;
2812     }
2813     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
2814   }
2815
2816   delete[] binsvar;
2817   delete[] binfill;
2818
2819   return fWeightCharm;
2820 }
2821
2822 //____________________________________________________________________________
2823 void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
2824
2825    // TOF PID efficiencies
2826    Int_t ptpr;
2827    if(fBeamType==0) ptpr=0;
2828    if(fBeamType==1) ptpr=1;
2829
2830    Int_t loopcentr=1;
2831    const Int_t nCentralitybinning=11; //number of centrality bins
2832    if(fBeamType==1)
2833    {
2834      loopcentr=nCentralitybinning;
2835    }
2836
2837    TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2838    TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
2839    TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
2840
2841    TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
2842    cefficiencyParamtof->cd();
2843
2844    AliCFContainer *mccontainermcD = 0x0;
2845    AliCFContainer *mccontaineresdD = 0x0;
2846    TH1D* efficiencysigTOFPIDD[nCentralitybinning];
2847    TH1D* efficiencyTOFPIDD[nCentralitybinning];
2848    TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
2849    TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
2850    Int_t source = -1; //get parameterized TOF PID efficiencies
2851
2852    for(int icentr=0; icentr<loopcentr; icentr++) {
2853       // signal sample
2854       if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2855       else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2856       AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
2857       efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2858
2859       // mb sample for double check
2860       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2861       else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2862       AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
2863       efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2864
2865       // mb sample with reconstructed variables
2866       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
2867       else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2868       AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
2869       efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2870
2871       // mb sample with reconstructed variables
2872       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2873       else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
2874       AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
2875       efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
2876
2877       //fill histo
2878       efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
2879       efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
2880       efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
2881       efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
2882       efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
2883       efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
2884       efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
2885       efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
2886
2887       //fit (mc enhenced sample)
2888       fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
2889       efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
2890       efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2891       fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
2892
2893       //fit (esd enhenced sample)
2894       efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
2895       efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
2896       fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
2897
2898    }
2899
2900    // draw (for PbPb, only 1st bin)
2901    //sig mc
2902    efficiencysigTOFPIDD[0]->SetTitle("");
2903    efficiencysigTOFPIDD[0]->SetStats(0);
2904    efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
2905    efficiencysigTOFPIDD[0]->SetMarkerColor(2);
2906    efficiencysigTOFPIDD[0]->SetLineColor(2);
2907    efficiencysigTOFPIDD[0]->Draw();
2908
2909    //mb mc
2910    efficiencyTOFPIDD[0]->SetTitle("");
2911    efficiencyTOFPIDD[0]->SetStats(0);
2912    efficiencyTOFPIDD[0]->SetMarkerStyle(24);
2913    efficiencyTOFPIDD[0]->SetMarkerColor(4);
2914    efficiencyTOFPIDD[0]->SetLineColor(4);
2915    efficiencyTOFPIDD[0]->Draw("same");
2916
2917    //sig esd
2918    efficiencysigesdTOFPIDD[0]->SetTitle("");
2919    efficiencysigesdTOFPIDD[0]->SetStats(0);
2920    efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
2921    efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
2922    efficiencysigesdTOFPIDD[0]->SetLineColor(3);
2923    efficiencysigesdTOFPIDD[0]->Draw("same");
2924
2925    //mb esd
2926    efficiencyesdTOFPIDD[0]->SetTitle("");
2927    efficiencyesdTOFPIDD[0]->SetStats(0);
2928    efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
2929    efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
2930    efficiencyesdTOFPIDD[0]->SetLineColor(1);
2931    efficiencyesdTOFPIDD[0]->Draw("same");
2932
2933    //signal mc fit
2934    if(fEfficiencyTOFPIDD[0]){
2935      fEfficiencyTOFPIDD[0]->SetLineColor(2);
2936      fEfficiencyTOFPIDD[0]->Draw("same");
2937    }
2938    //mb esd fit
2939    if(fEfficiencyesdTOFPIDD[0]){
2940        fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
2941        fEfficiencyesdTOFPIDD[0]->Draw("same");
2942      }
2943
2944    TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
2945    legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
2946    legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
2947    legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
2948    legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
2949    legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
2950    legtofeff->Draw("same");
2951
2952
2953    TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
2954    cefficiencyParamIP->cd();
2955    gStyle->SetOptStat(0);
2956
2957    // IP cut efficiencies
2958    for(int icentr=0; icentr<loopcentr; icentr++)  {
2959
2960      AliCFContainer *charmCombined = 0x0;
2961      AliCFContainer *beautyCombined = 0x0;
2962      AliCFContainer *beautyCombinedesd = 0x0;
2963
2964      printf("centrality printing %i \n",icentr);
2965
2966      source = 0; //charm enhenced
2967      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2968      else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2969      AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
2970      efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2971
2972      charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");
2973
2974      source = 1; //beauty enhenced
2975      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
2976      else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2977      AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
2978      efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2979
2980      beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined");
2981
2982      if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
2983      else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2984      AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
2985      efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2986
2987      beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
2988
2989      source = 0; //charm mb
2990      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
2991      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
2992      AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
2993      efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
2994
2995      charmCombined->Add(mccontainermcD);
2996      AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined);
2997      efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
2998
2999      source = 1; //beauty mb
3000      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3001      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3002      AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
3003      efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3004
3005      beautyCombined->Add(mccontainermcD);
3006      AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined);
3007      efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3008
3009      if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
3010      else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3011      AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
3012      efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3013
3014      beautyCombinedesd->Add(mccontaineresdD);
3015      AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
3016      efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
3017
3018      source = 2; //conversion mb
3019      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3020      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3021      AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
3022      efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3023
3024      source = 3; //non HFE except for the conversion mb
3025      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
3026      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
3027      AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
3028      efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
3029
3030      if(fIPEffCombinedSamples){
3031        fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb
3032        fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
3033        fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
3034      }
3035      else{
3036        fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
3037        fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
3038        fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
3039      }
3040      fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
3041      fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
3042      fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
3043      fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
3044
3045    }
3046
3047    if(fBeamType==0){
3048      AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
3049      fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
3050
3051      AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
3052      fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
3053    }
3054
3055    for(int icentr=0; icentr<loopcentr; icentr++)  {
3056      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3057      fipfit->SetLineColor(2);
3058      fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
3059      fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3060      if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
3061      else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
3062
3063      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3064      fipfit->SetLineColor(6);
3065      fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
3066      fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
3067      if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
3068      else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
3069
3070      fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3071      fipfit->SetLineColor(1);
3072      fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
3073      fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
3074      fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
3075
3076      if(fIPParameterizedEff){
3077        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3078        fipfitnonhfe->SetLineColor(3);
3079        fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
3080        fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3081        fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
3082
3083        fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
3084        fipfitnonhfe->SetLineColor(4);
3085        fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
3086        fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
3087        fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
3088      }
3089    }
3090
3091    // draw (for PbPb, only 1st bin)
3092    fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
3093    fEfficiencyCharmSigD[0]->SetMarkerColor(1);
3094    fEfficiencyCharmSigD[0]->SetLineColor(1);
3095    fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
3096    fEfficiencyBeautySigD[0]->SetMarkerColor(2);
3097    fEfficiencyBeautySigD[0]->SetLineColor(2);
3098    fEfficiencyBeautySigesdD[0]->SetStats(0);
3099    fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
3100    fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
3101    fEfficiencyBeautySigesdD[0]->SetLineColor(6);
3102    fCharmEff[0]->SetMarkerStyle(24);
3103    fCharmEff[0]->SetMarkerColor(1);
3104    fCharmEff[0]->SetLineColor(1);
3105    fBeautyEff[0]->SetMarkerStyle(24);
3106    fBeautyEff[0]->SetMarkerColor(2);
3107    fBeautyEff[0]->SetLineColor(2);
3108    fConversionEff[0]->SetMarkerStyle(24);
3109    fConversionEff[0]->SetMarkerColor(3);
3110    fConversionEff[0]->SetLineColor(3);
3111    fNonHFEEff[0]->SetMarkerStyle(24);
3112    fNonHFEEff[0]->SetMarkerColor(4);
3113    fNonHFEEff[0]->SetLineColor(4);
3114
3115    fEfficiencyCharmSigD[0]->Draw();
3116    fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
3117    fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
3118
3119    fEfficiencyBeautySigD[0]->Draw("same");
3120    fEfficiencyBeautySigesdD[0]->Draw("same");
3121    //fCharmEff[0]->Draw("same");
3122    //fBeautyEff[0]->Draw("same");
3123
3124    if(fBeamType==0){
3125      fConversionEffbgc->SetMarkerStyle(25);
3126      fConversionEffbgc->SetMarkerColor(3);
3127      fConversionEffbgc->SetLineColor(3);
3128      fNonHFEEffbgc->SetMarkerStyle(25);
3129      fNonHFEEffbgc->SetMarkerColor(4);
3130      fNonHFEEffbgc->SetLineColor(4);
3131      fConversionEffbgc->Draw("same");
3132      fNonHFEEffbgc->Draw("same");
3133    }
3134    else{
3135      fConversionEff[0]->Draw("same");
3136      fNonHFEEff[0]->Draw("same");
3137    }
3138    if(fEfficiencyIPBeautyD[0])
3139       fEfficiencyIPBeautyD[0]->Draw("same");
3140    if(fEfficiencyIPBeautyesdD[0])
3141      fEfficiencyIPBeautyesdD[0]->Draw("same");
3142    if( fEfficiencyIPCharmD[0])
3143      fEfficiencyIPCharmD[0]->Draw("same");
3144    if(fIPParameterizedEff){
3145      fEfficiencyIPConversionD[0]->Draw("same");
3146      fEfficiencyIPNonhfeD[0]->Draw("same");
3147    }
3148    TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
3149    legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
3150    legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
3151    legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
3152    legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
3153    legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
3154    legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
3155    //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
3156    //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
3157    legipeff->Draw("same");
3158    gPad->SetGrid();
3159    //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
3160 }
3161
3162 //____________________________________________________________________________
3163 THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
3164   //
3165   // Return beauty electron IP cut efficiency
3166   //
3167
3168   const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3169   const Int_t nCentralitybinning=11;//number of centrality bins
3170   Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3171   Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3172   Int_t ptpr = 0;
3173   Int_t nDim=1;  //dimensions of the efficiency weighting grid
3174   if(fBeamType==1)
3175   {
3176     ptpr=1;
3177     nDim=2; //dimensions of the efficiency weighting grid
3178   }
3179   Int_t nBin[1] = {nPtbinning1};
3180   Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3181
3182
3183   THnSparseF *ipcut;
3184   if(fBeamType==0) ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
3185   else ipcut = new THnSparseF("beff", "b IP efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3186
3187   for(Int_t idim = 0; idim < nDim; idim++)
3188   {
3189     if(nDim==1) ipcut->SetBinEdges(idim, kPtRange);
3190     if(nDim==2)
3191       {
3192         ipcut->SetBinEdges(0, kCentralityRange);
3193         ipcut->SetBinEdges(1, kPtRange);
3194       }
3195   }
3196   Double_t pt[1];
3197   Double_t weight[nCentralitybinning];
3198   Double_t weightErr[nCentralitybinning];
3199   Double_t contents[2];
3200
3201   for(Int_t a=0;a<11;a++)
3202   {
3203       weight[a] = 1.0;
3204       weightErr[a] = 1.0;
3205   }
3206
3207
3208   Int_t looppt=nBin[0];
3209   Int_t loopcentr=1;
3210   Int_t ibin[2];
3211   if(fBeamType==1)
3212   {
3213       loopcentr=nBinPbPb[0];
3214   }
3215
3216
3217   for(int icentr=0; icentr<loopcentr; icentr++)
3218   {
3219       for(int i=0; i<looppt; i++)
3220       {
3221           pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3222           if(isMCpt){
3223             if(fEfficiencyIPBeautyD[icentr]){
3224               weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
3225               weightErr[icentr] = 0;
3226             }
3227             else{
3228               printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3229               weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1);
3230               weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3231             }
3232           }
3233           else{
3234             if(fEfficiencyIPBeautyesdD[icentr]){
3235               weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
3236               weightErr[icentr] = 0;
3237             }
3238             else{
3239               printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
3240               weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
3241               weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
3242             }
3243           }
3244
3245           if(fBeamType==1){
3246               contents[0]=icentr;
3247               contents[1]=pt[0];
3248               ibin[0]=icentr;
3249               ibin[1]=i+1;
3250           }
3251           if(fBeamType==0){
3252               contents[0]=pt[0];
3253               ibin[0]=i+1;
3254           }
3255           ipcut->Fill(contents,weight[icentr]);
3256           ipcut->SetBinError(ibin,weightErr[icentr]);
3257       }
3258   }
3259
3260   Int_t nDimSparse = ipcut->GetNdimensions();
3261   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3262   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3263
3264   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3265       binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
3266       nBins *= binsvar[iVar];
3267   }
3268
3269   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3270   // loop that sets 0 error in each bin
3271   for (Long_t iBin=0; iBin<nBins; iBin++) {
3272     Long_t bintmp = iBin ;
3273     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3274       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3275       bintmp /= binsvar[iVar] ;
3276     }
3277     //ipcut->SetBinError(binfill,0.); // put 0 everywhere
3278   }
3279
3280   delete[] binsvar;
3281   delete[] binfill;
3282
3283   return ipcut;
3284 }
3285
3286 //____________________________________________________________________________
3287 THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
3288   //
3289   // Return PID x IP cut efficiency
3290   //
3291     const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
3292     const Int_t nCentralitybinning=11;//number of centrality bins
3293     Double_t kPtRange[nPtbinning1+1] = { 0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 4.5, 5., 5.5, 6., 7., 8., 10., 12., 14., 16., 18., 20.};//pt bin limits
3294     Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3295     Int_t ptpr = 0;
3296     Int_t nDim=1;  //dimensions of the efficiency weighting grid
3297     if(fBeamType==1)
3298     {
3299         ptpr=1;
3300         nDim=2; //dimensions of the efficiency weighting grid
3301     }
3302     Int_t nBin[1] = {nPtbinning1};
3303     Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3304
3305     THnSparseF *pideff;
3306     if(fBeamType==0) pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
3307     else pideff = new THnSparseF("pideff", "PID efficiency; centrality bin; p_{t}(GeV/c)", nDim, nBinPbPb);
3308     for(Int_t idim = 0; idim < nDim; idim++)
3309     {
3310
3311         if(nDim==1) pideff->SetBinEdges(idim, kPtRange);
3312         if(nDim==2)
3313         {
3314             pideff->SetBinEdges(0, kCentralityRange);
3315             pideff->SetBinEdges(1, kPtRange);
3316         }
3317     }
3318
3319   Double_t pt[1];
3320   Double_t weight[nCentralitybinning];
3321   Double_t weightErr[nCentralitybinning];
3322   Double_t contents[2];
3323
3324   for(Int_t a=0;a<nCentralitybinning;a++)
3325   {
3326       weight[a] = 1.0;
3327       weightErr[a] = 1.0;
3328   }
3329
3330   Int_t looppt=nBin[0];
3331   Int_t loopcentr=1;
3332   Int_t ibin[2];
3333   if(fBeamType==1)
3334   {
3335       loopcentr=nBinPbPb[0];
3336   }
3337
3338   for(int icentr=0; icentr<loopcentr; icentr++)
3339   {
3340       Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
3341       for(int i=0; i<looppt; i++)
3342       {
3343           pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3344
3345           Double_t tofpideff = 0.;
3346           Double_t tofpideffesd = 0.;
3347           if(fEfficiencyTOFPIDD[icentr])
3348             tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
3349           else{
3350             printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3351           }
3352           if(fEfficiencyesdTOFPIDD[icentr])
3353             tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
3354           else{
3355             printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
3356           }
3357
3358           //tof pid eff x tpc pid eff x ip cut eff
3359           if(fIPParameterizedEff){
3360             if(source==0) {
3361               if(fEfficiencyIPCharmD[icentr]){
3362                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3363                 weightErr[icentr] = 0;
3364               }
3365               else{
3366                 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3367                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3368                 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3369               }
3370             }
3371             else if(source==2) {
3372               if(fEfficiencyIPConversionD[icentr]){
3373                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]);
3374                 weightErr[icentr] = 0;
3375               }
3376               else{
3377                 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
3378                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
3379                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3380               }
3381             }
3382             else if(source==3) {
3383               if(fEfficiencyIPNonhfeD[icentr]){
3384                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]);
3385                 weightErr[icentr] = 0;
3386               }
3387               else{
3388                 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
3389                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
3390                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3391               }
3392             }
3393           }
3394           else{
3395             if(source==0){
3396               if(fEfficiencyIPCharmD[icentr]){
3397                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
3398                 weightErr[icentr] = 0;
3399               }
3400               else{
3401                 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
3402                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
3403                 weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
3404               }
3405             }
3406             else if(source==2){
3407               if(fBeamType==0){
3408                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
3409                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
3410               }
3411               else{
3412                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
3413                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
3414               }
3415             }
3416             else if(source==3){
3417               if(fBeamType==0){
3418                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
3419                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
3420               }
3421               else{
3422                 weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
3423                 weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
3424               }
3425             }
3426           }
3427
3428           if(fBeamType==1){
3429               contents[0]=icentr;
3430               contents[1]=pt[0];
3431               ibin[0]=icentr;
3432               ibin[1]=i+1;
3433           }
3434           if(fBeamType==0){
3435               contents[0]=pt[0];
3436               ibin[0]=i+1;
3437           }
3438
3439           pideff->Fill(contents,weight[icentr]);
3440           pideff->SetBinError(ibin,weightErr[icentr]);
3441       }
3442   }
3443
3444   Int_t nDimSparse = pideff->GetNdimensions();
3445   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3446   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3447
3448   for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3449       binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
3450       nBins *= binsvar[iVar];
3451   }
3452
3453   Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3454   // loop that sets 0 error in each bin
3455   for (Long_t iBin=0; iBin<nBins; iBin++) {
3456     Long_t bintmp = iBin ;
3457     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3458       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3459       bintmp /= binsvar[iVar] ;
3460     }
3461   }
3462
3463   delete[] binsvar;
3464   delete[] binfill;
3465
3466
3467   return pideff;
3468 }
3469
3470 //__________________________________________________________________________
3471 AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
3472  //
3473  // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
3474     //
3475     Int_t ptpr = 0;
3476     Int_t nDim = 1;
3477     if(fBeamType==0)
3478     {
3479         ptpr=0;
3480     }
3481     if(fBeamType==1)
3482     {
3483         ptpr=1;
3484         nDim=2;
3485     }
3486
3487     const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
3488     const Int_t nCentralitybinning=11;//number of centrality bins
3489     Double_t kPtRange[19] = {0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 2, 2.5, 3, 4, 5, 6, 8, 12, 16, 20};
3490
3491     Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
3492     Int_t nBin[1] = {nPtbinning1};
3493     Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
3494
3495     AliCFDataGrid *rawBeautyContainer;
3496     if(fBeamType==0)  rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
3497     else rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBinPbPb);
3498     //  printf("number of bins= %d\n",bins[0]);
3499
3500
3501
3502
3503     THnSparseF *brawspectra;
3504     if(fBeamType==0) brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
3505     else brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBinPbPb);
3506     if(fBeamType==0) brawspectra->SetBinEdges(0, kPtRange);
3507     if(fBeamType==1)
3508       {
3509         //      brawspectra->SetBinEdges(0, centralityBins);
3510         brawspectra->SetBinEdges(0, kCentralityRange);
3511         brawspectra->SetBinEdges(1, kPtRange);
3512       }
3513
3514     Double_t pt[1];
3515     Double_t yields= 0.;
3516     Double_t valuesb[2];
3517
3518     //Int_t looppt=nBin[0];
3519     Int_t loopcentr=1;
3520     if(fBeamType==1)
3521       {
3522         loopcentr=nBinPbPb[0];
3523       }
3524
3525     for(int icentr=0; icentr<loopcentr; icentr++)
3526       {
3527
3528         for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
3529           pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
3530
3531           yields = fBSpectrum2ndMethod->GetBinContent(i+1);
3532
3533           if(fBeamType==1)
3534             {
3535               valuesb[0]=icentr;
3536               valuesb[1]=pt[0];
3537             }
3538           if(fBeamType==0) valuesb[0]=pt[0];
3539           brawspectra->Fill(valuesb,yields);
3540         }
3541       }
3542
3543
3544
3545     Int_t nDimSparse = brawspectra->GetNdimensions();
3546     Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
3547     Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
3548
3549     for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3550       binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
3551       nBins *= binsvar[iVar];
3552     }
3553
3554     Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
3555     // loop that sets 0 error in each bin
3556     for (Long_t iBin=0; iBin<nBins; iBin++) {
3557       Long_t bintmp = iBin ;
3558       for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
3559         binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
3560         bintmp /= binsvar[iVar] ;
3561       }
3562       brawspectra->SetBinError(binfill,0.); // put 0 everywhere
3563     }
3564
3565
3566     rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
3567     TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(ptpr);
3568
3569     new TCanvas;
3570     fBSpectrum2ndMethod->SetMarkerStyle(24);
3571     fBSpectrum2ndMethod->Draw("p");
3572     hRawBeautySpectra->SetMarkerStyle(25);
3573     hRawBeautySpectra->Draw("samep");
3574
3575     delete[] binfill;
3576     delete[] binsvar;
3577
3578     return rawBeautyContainer;
3579 }
3580
3581 //__________________________________________________________________________
3582 void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
3583   //
3584   // Calculate non HFE sys
3585   //
3586   //
3587
3588   if(!fNonHFEsyst)
3589     return;
3590
3591   Double_t evtnorm[1] = {0.0};
3592   if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
3593
3594   AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
3595   AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
3596
3597   AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
3598   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
3599   AliCFDataGrid *bgConvGrid[kBgLevels];
3600
3601   Int_t stepbackground = 3;
3602   Int_t* bins=new Int_t[1];
3603   const Char_t *bgBase[2] = {"pi0","eta"};
3604
3605   bins[0]=fConversionEff[centrality]->GetNbinsX();
3606
3607   AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
3608   AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
3609
3610   for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
3611     for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
3612       convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),Form("convGrid_%d_%d_%d",iSource,iLevel,centrality),*fConvSourceContainer[iSource][iLevel][centrality],stepbackground);
3613       weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
3614       convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
3615
3616       nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),Form("nonHFEGrid_%d_%d_%d",iSource,iLevel,centrality),*fNonHFESourceContainer[iSource][iLevel][centrality],stepbackground);
3617       weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
3618       nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
3619     }
3620
3621     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
3622     for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
3623       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
3624     }
3625     if(!fEtaSyst)
3626       bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
3627
3628     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone();
3629     for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
3630       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
3631     }
3632     if(!fEtaSyst)
3633       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
3634
3635     bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
3636     bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
3637     if(fEtaSyst){
3638       bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
3639       bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
3640     }
3641   }
3642
3643
3644   //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0)
3645   AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
3646   TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
3647   for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
3648     bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
3649     bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
3650     bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();
3651     bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
3652
3653   //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
3654
3655     hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
3656     hBaseErrors[iErr][0]->Scale(-1.);
3657     hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
3658     hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
3659     if(!fEtaSyst)break;
3660   }
3661
3662
3663   //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
3664   TH1D *hSpeciesErrors[kElecBgSources-1];
3665   for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3666     if(fEtaSyst && (iSource == 1))continue;
3667     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
3668     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
3669     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
3670     hSpeciesErrors[iSource-1]->Scale(0.3);
3671   }
3672
3673   //Int_t firstBgSource = 0;//if eta systematics are not from scaling
3674   //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
3675   TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
3676   TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
3677   TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
3678
3679   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
3680   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
3681
3682   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
3683     Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;
3684     pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin);
3685     pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
3686     if(fEtaSyst){
3687       etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin);
3688       etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
3689     }
3690     else{ etaErrLow = etaErrUp = 0.;}
3691
3692     Double_t sqrsumErrs= 0;
3693     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
3694       if(fEtaSyst && (iSource == 1))continue;
3695       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
3696       sqrsumErrs+=(scalingErr*scalingErr);
3697     }
3698     for(Int_t iErr = 0; iErr < 2; iErr++){
3699       for(Int_t iLevel = 0; iLevel < 2; iLevel++){
3700         hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
3701       }
3702       if(!fEtaSyst)break;
3703     }
3704     hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
3705     hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
3706     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
3707
3708     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
3709     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));
3710   }
3711
3712
3713   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
3714   cPiErrors->cd();
3715   cPiErrors->SetLogx();
3716   cPiErrors->SetLogy();
3717   hBaseErrors[0][0]->Draw();
3718   //hBaseErrors[0][1]->SetMarkerColor(kBlack);
3719   //hBaseErrors[0][1]->SetLineColor(kBlack);
3720   //hBaseErrors[0][1]->Draw("SAME");
3721   if(fEtaSyst){
3722     hBaseErrors[1][0]->Draw("SAME");
3723     hBaseErrors[1][0]->SetMarkerColor(kBlack);
3724     hBaseErrors[1][0]->SetLineColor(kBlack);
3725   //hBaseErrors[1][1]->SetMarkerColor(13);
3726   //hBaseErrors[1][1]->SetLineColor(13);
3727   //hBaseErrors[1][1]->Draw("SAME");
3728   }
3729   //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
3730   //hOverallBinScaledErrsUp->SetLineColor(kBlue);
3731   //hOverallBinScaledErrsUp->Draw("SAME");
3732   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
3733   hOverallBinScaledErrsLow->SetLineColor(kGreen);
3734   hOverallBinScaledErrsLow->Draw("SAME");
3735   hScalingErrors->SetLineColor(kBlue);
3736   hScalingErrors->Draw("SAME");
3737
3738   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
3739   lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
3740   //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
3741   if(fEtaSyst){
3742   lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
3743   //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
3744   }
3745   lPiErr->AddEntry(hScalingErrors, "scaling error");
3746   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
3747   //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
3748   lPiErr->Draw("SAME");
3749
3750   //Normalize errors
3751   TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
3752   TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
3753   hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
3754   hLowSystScaled->Scale(evtnorm[0]);
3755   TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
3756   TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
3757   //histograms to be normalized to TGraphErrors
3758   CorrectFromTheWidth(hNormAllSystErrUp);
3759   CorrectFromTheWidth(hNormAllSystErrLow);
3760
3761   TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
3762   cNormOvErrs->cd();
3763   cNormOvErrs->SetLogx();
3764   cNormOvErrs->SetLogy();
3765
3766   TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
3767   TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
3768   gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
3769   gOverallSystErrUp->SetMarkerColor(kBlack);
3770   gOverallSystErrUp->SetLineColor(kBlack);
3771   gOverallSystErrLow->SetMarkerColor(kRed);
3772   gOverallSystErrLow->SetLineColor(kRed);
3773   gOverallSystErrUp->Draw("AP");
3774   gOverallSystErrLow->Draw("PSAME");
3775   TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
3776   lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
3777   lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
3778   lAllSys->Draw("same");
3779
3780
3781   AliCFDataGrid *bgYieldGrid;
3782   if(fEtaSyst){
3783     bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
3784   }
3785   bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
3786
3787   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
3788   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
3789   hRelErrUp->Divide(hBgYield);
3790   TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
3791   hRelErrLow->Divide(hBgYield);
3792
3793   TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
3794   cRelErrs->cd();
3795   cRelErrs->SetLogx();
3796   hRelErrUp->SetTitle("Relative error of non-HFE background yield");
3797   hRelErrUp->Draw();
3798   hRelErrLow->SetLineColor(kBlack);
3799   hRelErrLow->Draw("SAME");
3800
3801   TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
3802   lRel->AddEntry(hRelErrUp, "upper");
3803   lRel->AddEntry(hRelErrLow, "lower");
3804   lRel->Draw("SAME");
3805
3806   //CorrectFromTheWidth(hBgYield);
3807   //hBgYield->Scale(evtnorm[0]);
3808
3809
3810   //write histograms/TGraphs to file
3811   TFile *output = new TFile("systHists.root","recreate");
3812
3813   hBgYield->SetName("hBgYield");
3814   hBgYield->Write();
3815   hRelErrUp->SetName("hRelErrUp");
3816   hRelErrUp->Write();
3817   hRelErrLow->SetName("hRelErrLow");
3818   hRelErrLow->Write();
3819   hUpSystScaled->SetName("hOverallSystErrUp");
3820   hUpSystScaled->Write();
3821   hLowSystScaled->SetName("hOverallSystErrLow");
3822   hLowSystScaled->Write();
3823   gOverallSystErrUp->SetName("gOverallSystErrUp");
3824   gOverallSystErrUp->Write();
3825   gOverallSystErrLow->SetName("gOverallSystErrLow");
3826   gOverallSystErrLow->Write();
3827
3828   output->Close();
3829   delete output;
3830
3831 }
3832
3833 //____________________________________________________________
3834 void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
3835
3836   //
3837   // Unfold backgrounds to check its sanity
3838   //
3839
3840   AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
3841   //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
3842   if(!mcContainer){
3843     AliError("MC Container not available");
3844     return;
3845   }
3846
3847   if(!fCorrelation){
3848     AliError("No Correlation map available");
3849     return;
3850   }
3851
3852   // Data
3853   AliCFDataGrid *dataGrid = 0x0;
3854   dataGrid = bgsubpectrum;
3855
3856   // Guessed
3857   AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
3858   THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
3859
3860   // Efficiency
3861   AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
3862   efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
3863
3864   // Unfold background spectra
3865   Int_t nDim=1;
3866   if(fBeamType==0)nDim = 1;
3867   if(fBeamType==1)nDim = 2;
3868   AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
3869   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
3870   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
3871   if(fSetSmoothing) unfolding.UseSmoothing();
3872   unfolding.Unfold();
3873
3874   // Results
3875   THnSparse* result = unfolding.GetUnfolded();
3876   TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
3877   if(fBeamType==1)
3878   {
3879       ctest->Divide(2,2);
3880       ctest->cd(1);
3881       result->GetAxis(0)->SetRange(1,1);
3882       TH1D* htest1=(TH1D*)result->Projection(0);
3883       htest1->Draw();
3884       ctest->cd(2);
3885       result->GetAxis(0)->SetRange(1,1);
3886       TH1D* htest2=(TH1D*)result->Projection(1);
3887       htest2->Draw();
3888       ctest->cd(3);
3889       result->GetAxis(0)->SetRange(6,6);
3890       TH1D* htest3=(TH1D*)result->Projection(0);
3891       htest3->Draw();
3892       ctest->cd(4);
3893       result->GetAxis(0)->SetRange(6,6);
3894       TH1D* htest4=(TH1D*)result->Projection(1);
3895       htest4->Draw();
3896
3897   }
3898
3899
3900
3901
3902
3903   TGraphErrors* unfoldedbgspectrumD = Normalize(result);
3904   if(!unfoldedbgspectrumD) {
3905     AliError("Unfolded background spectrum doesn't exist");
3906   }
3907   else{
3908     TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
3909     if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
3910
3911     if(fBeamType==1)
3912     {
3913         Int_t centr=1;
3914         result->GetAxis(0)->SetRange(centr,centr);
3915         unfoldedbgspectrumD = Normalize(result,centr-1);
3916         unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
3917         centr=6;
3918         result->GetAxis(0)->SetRange(centr,centr);
3919         unfoldedbgspectrumD = Normalize(result,centr-1);
3920         unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");
3921     }
3922
3923     file->Close();
3924   }
3925 }