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