]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEspectrum.cxx
Updates for TRD HFE analysis
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEspectrum.cxx
index 61af42b2a4827b3862bc6acc69a0bde93be5def4..bbd1ab5654cb1cb3de5051a6474f9c9badac67c8 100644 (file)
@@ -80,6 +80,7 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   fIPParameterizedEff(kFALSE),
   fNonHFEsyst(0),
   fBeauty2ndMethod(kFALSE),
+  fIPEffCombinedSamples(kTRUE),
   fNbDimensions(1),
   fNMCEvents(0),
   fStepMC(-1),
@@ -91,12 +92,19 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   fNumberOfIterations(5),
   fChargeChoosen(kAllCharge),
   fNCentralityBinAtTheEnd(0),
+  fTestCentralityLow(-1),
+  fTestCentralityHigh(-1),
+  fFillMoreCorrelationMatrix(kFALSE),
   fHadronEffbyIPcut(NULL),
+  fConversionEffbgc(NULL),
+  fNonHFEEffbgc(NULL),      
   fBSpectrum2ndMethod(NULL),
   fkBeauty2ndMethodfilename(""),
   fBeamType(0),
+  fEtaSyst(kTRUE),
   fDebugLevel(0),
-  fWriteToFile(kFALSE)
+  fWriteToFile(kFALSE),
+  fUnfoldBG(kFALSE)
 {
   //
   // Default constructor
@@ -110,8 +118,10 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
       if(k<kCentrality)
       {
           fEfficiencyTOFPIDD[k] = 0;
+          fEfficiencyesdTOFPIDD[k] = 0;
           fEfficiencyIPCharmD[k] = 0;     
           fEfficiencyIPBeautyD[k] = 0;    
+          fEfficiencyIPBeautyesdD[k] = 0;
           fEfficiencyIPConversionD[k] = 0;
           fEfficiencyIPNonhfeD[k] = 0;   
 
@@ -121,6 +131,7 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
          fBeautyEff[k] = 0;
          fEfficiencyCharmSigD[k] = 0;
           fEfficiencyBeautySigD[k] = 0;
+          fEfficiencyBeautySigesdD[k] = 0;
       }
   }
   memset(fEtaRange, 0, sizeof(Double_t) * 2);
@@ -154,7 +165,7 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   // and the appropriate correlation matrix
   //
 
-
+  
   if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
 
   Int_t kNdim = 3;
@@ -217,8 +228,8 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
   if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
 
-  AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen);
-  AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen);
+  AliCFContainer *datacontainerD = GetSlicedContainer(datacontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  AliCFContainer *contaminationcontainerD = GetSlicedContainer(contaminationcontainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
   if((!datacontainerD) || (!contaminationcontainerD)) return kFALSE;
 
   SetContainer(datacontainerD,AliHFEspectrum::kDataContainer);
@@ -226,6 +237,7 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
 
   // MC container: ESD/MC efficiency maps + MC/MC efficiency maps 
   AliCFContainer *mccontaineresd = 0x0;
+  AliCFContainer *mccontaineresdbg = 0x0;
   AliCFContainer *mccontainermc = 0x0;
   AliCFContainer *mccontainermcbg = 0x0;
   AliCFContainer *nonHFEweightedContainer = 0x0;
@@ -239,6 +251,7 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   }
   else {
     mccontaineresd = mchfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
+    mccontaineresdbg = bghfecontainer->MakeMergedCFContainer("sumesd","sumesd","MCTrackCont:recTrackContReco:recTrackContDEReco");
     mccontainermc = mchfecontainer->MakeMergedCFContainer("summc","summc","MCTrackCont:recTrackContMC:recTrackContDEMC");
     mccontainermcbg = bghfecontainer->MakeMergedCFContainer("summcbg","summcbg","MCTrackCont:recTrackContMC:recTrackContDEMC");
 
@@ -258,26 +271,28 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
              }
              if(fBeamType==1)
              {
-                 fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr);
-                 fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr);
+              
+               fConvSourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
+               fNonHFESourceContainer[iSource][iLevel][icentr] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,icentr,icentr);
              }
-             if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
+//           if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
          }
+          if(fBeamType == 1)break;
        }
       }
     }
-    else{      
+    // else{      
       nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
       convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
       if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE;  
-    }
+      //}
   }
   if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
   
   Int_t source = -1;
   if(!fInclusiveSpectrum) source = 1; //beauty
-  AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen);
-  AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen);
+  AliCFContainer *mccontaineresdD = GetSlicedContainer(mccontaineresd, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  AliCFContainer *mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
   if((!mccontaineresdD) || (!mccontainermcD)) return kFALSE;
   SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerMC);
   SetContainer(mccontaineresdD,AliHFEspectrum::kMCContainerESD);
@@ -288,20 +303,21 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
    mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen);
    SetContainer(mccontainermcD,AliHFEspectrum::kMCContainerCharmMC);
 
-   SetParameterizedEff(mccontainermc, mccontainermcbg, dims);
-
-   if(!fNonHFEsyst){
+   //if(!fNonHFEsyst){
      AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
      SetContainer(nonHFEweightedContainerD,AliHFEspectrum::kMCWeightedContainerNonHFEESD);
      AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen);
      SetContainer(convweightedContainerD,AliHFEspectrum::kMCWeightedContainerConversionESD);
-   }
+     //}
+
+     SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
+
   }
   // MC container: correlation matrix
   THnSparseF *mccorrelation = 0x0;
   if(fInclusiveSpectrum) {
-    if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
-    else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
+    if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 2)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
+    else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD + 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
     else if(fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1)) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
     else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
@@ -311,7 +327,13 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
   //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
   if(!mccorrelation) return kFALSE;
-  THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims);
+  Int_t testCentralityLow = fTestCentralityLow;
+  Int_t testCentralityHigh = fTestCentralityHigh;
+  if(fFillMoreCorrelationMatrix) {
+    testCentralityLow = fTestCentralityLow-1;
+    testCentralityHigh = fTestCentralityHigh+1;
+  }
+  THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims,testCentralityLow,testCentralityHigh);
   if(!mccorrelationD) {
     printf("No correlation\n");
     return kFALSE;
@@ -322,7 +344,7 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
   if(v0hfecontainer) {
     AliCFContainer *containerV0 = v0hfecontainer->GetCFContainer("taggedTrackContainerReco");
     if(!containerV0) return kFALSE;
-    AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron);
+    AliCFContainer *containerV0Electron = GetSlicedContainer(containerV0, fNbDimensions, dims, AliPID::kElectron,fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
     if(!containerV0Electron) return kFALSE;
     SetContainer(containerV0Electron,AliHFEspectrum::kDataContainerV0);
   }
@@ -368,9 +390,17 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
         ptcorrelation->Draw("colz");
       }
     }
-    if(fWriteToFile) ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
+    if(fWriteToFile){ 
+      ccontaminationspectrum->SaveAs("contaminationspectrum.eps");
+      ccorrelation->SaveAs("correlationMatrix.eps");
+    }
   }
 
+  TFile *file = TFile::Open("tests.root","recreate");
+  datacontainerD->Write("data");
+  mccontainermcD->Write("mcefficiency");
+  mccorrelationD->Write("correlationmatrix");
+  file->Close();
 
   return kTRUE;
 }
@@ -877,32 +907,59 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
       alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
       alltogetherCorrection->Write();
       //
+      if(unnormalizedRawSpectrum) {
       unnormalizedRawSpectrum->SetName("beautyAfterIP");
       unnormalizedRawSpectrum->Write();
-      
+      }
+
       if(gNormalizedRawSpectrum){
         gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
         gNormalizedRawSpectrum->Write();
       }
 
+      if(fBeamType==0) {
+          Int_t countpp=0;
+         fEfficiencyCharmSigD[countpp]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",countpp));
+         fEfficiencyCharmSigD[countpp]->SetName(Form("IPEfficiencyForCharmSigCent%i",countpp));
+         fEfficiencyCharmSigD[countpp]->Write();
+         fEfficiencyBeautySigD[countpp]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",countpp));
+         fEfficiencyBeautySigD[countpp]->SetName(Form("IPEfficiencyForBeautySigCent%i",countpp));
+         fEfficiencyBeautySigD[countpp]->Write();
+         fCharmEff[countpp]->SetTitle(Form("IPEfficiencyForCharmCent%i",countpp));
+         fCharmEff[countpp]->SetName(Form("IPEfficiencyForCharmCent%i",countpp));
+         fCharmEff[countpp]->Write();
+         fBeautyEff[countpp]->SetTitle(Form("IPEfficiencyForBeautyCent%i",countpp));
+         fBeautyEff[countpp]->SetName(Form("IPEfficiencyForBeautyCent%i",countpp));
+         fBeautyEff[countpp]->Write();
+         fConversionEff[countpp]->SetTitle(Form("IPEfficiencyForConversionCent%i",countpp));
+         fConversionEff[countpp]->SetName(Form("IPEfficiencyForConversionCent%i",countpp));
+         fConversionEff[countpp]->Write();
+         fNonHFEEff[countpp]->SetTitle(Form("IPEfficiencyForNonHFECent%i",countpp));
+         fNonHFEEff[countpp]->SetName(Form("IPEfficiencyForNonHFECent%i",countpp));
+         fNonHFEEff[countpp]->Write();
+      }
+
       if(fBeamType==1) {
 
          TGraphErrors* correctedspectrumDc[kCentrality];
           TGraphErrors* alltogetherspectrumDc[kCentrality];
-         for(Int_t i=0;i<kCentrality-1;i++)
+         for(Int_t i=0;i<kCentrality-2;i++)
          {
              correctedspectrum->GetAxis(0)->SetRange(i+1,i+1);
              correctedspectrumDc[i] = Normalize(correctedspectrum,i);
-             correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
-             correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
+              if(correctedspectrumDc[i]){
+                correctedspectrumDc[i]->SetTitle(Form("UnfoldingCorrectedSpectrum_%i",i));
+                correctedspectrumDc[i]->SetName(Form("UnfoldingCorrectedSpectrum_%i",i));
+                correctedspectrumDc[i]->Write();
+              }
              alltogetherCorrection->GetAxis(0)->SetRange(i+1,i+1);
              alltogetherspectrumDc[i] = Normalize(alltogetherCorrection,i);
-             alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
-             alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
-             correctedspectrumDc[i]->Write();
-             alltogetherspectrumDc[i]->Write();
-             
-
+              if(alltogetherspectrumDc[i]){
+                alltogetherspectrumDc[i]->SetTitle(Form("AlltogetherSpectrum_%i",i));
+                alltogetherspectrumDc[i]->SetName(Form("AlltogetherSpectrum_%i",i));
+                alltogetherspectrumDc[i]->Write();
+              }
+           
              TH1D *centrcrosscheck = correctedspectrum->Projection(0);
              centrcrosscheck->SetTitle(Form("centrality_%i",i));
              centrcrosscheck->SetName(Form("centrality_%i",i));
@@ -922,11 +979,11 @@ Bool_t AliHFEspectrum::CorrectBeauty(Bool_t subtractcontamination){
              ratiocorrectedc->SetName(Form("RatioUnfoldingAlltogetherSpectrum_%i",i));
               ratiocorrectedc->Write();
 
-              fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
-             fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
+              fEfficiencyCharmSigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
+             fEfficiencyCharmSigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
               fEfficiencyCharmSigD[i]->Write();
-             fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForCharmSigCent%i",i));
-             fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForCharmSigCent%i",i));
+             fEfficiencyBeautySigD[i]->SetTitle(Form("IPEfficiencyForBeautySigCent%i",i));
+             fEfficiencyBeautySigD[i]->SetName(Form("IPEfficiencyForBeautySigCent%i",i));
               fEfficiencyBeautySigD[i]->Write();
              fCharmEff[i]->SetTitle(Form("IPEfficiencyForCharmCent%i",i));
              fCharmEff[i]->SetName(Form("IPEfficiencyForCharmCent%i",i));
@@ -983,6 +1040,7 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
   AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
   dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
  
+
   // Background Estimate
   AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
   if(!backgroundContainer){
@@ -995,17 +1053,33 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
 
   if(!fInclusiveSpectrum){
     //Background subtraction for IP analysis
+
+    TH1D *incElecCent[kCentrality-1];
+    TH1D *charmCent[kCentrality-1];
+    TH1D *convCent[kCentrality-1];
+    TH1D *nonHFECent[kCentrality-1]; 
+    TH1D *subtractedCent[kCentrality-1];
     TH1D *measuredTH1Draw = (TH1D *) dataspectrumbeforesubstraction->Project(ptpr);
     CorrectFromTheWidth(measuredTH1Draw);
-    TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",600,500);
+    if(fBeamType==1){
+      THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();
+      for(Int_t icent =  1; icent < kCentrality-1; icent++){
+        sparseIncElec->GetAxis(0)->SetRange(icent,icent);
+        incElecCent[icent-1] = (TH1D *) sparseIncElec->Projection(ptpr);
+        CorrectFromTheWidth(incElecCent[icent-1]);
+      }
+    }
+    TCanvas *rawspectra = new TCanvas("rawspectra","rawspectra",500,400);
     rawspectra->cd();
     rawspectra->SetLogy();
-    TLegend *lRaw = new TLegend(0.65,0.65,0.95,0.95);
+    gStyle->SetOptStat(0);
+    TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
     measuredTH1Draw->SetMarkerStyle(20);
     measuredTH1Draw->Draw();
+    measuredTH1Draw->GetXaxis()->SetRangeUser(0.0,7.9);
     lRaw->AddEntry(measuredTH1Draw,"measured raw spectrum");
     TH1D* htemp;
-    Int_t* bins=new Int_t[1];
+    Int_t* bins=new Int_t[2];
     if(fIPanaHadronBgSubtract){
       // Hadron background
        printf("Hadron background for IP analysis subtracted!\n");
@@ -1017,7 +1091,6 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
        }
        if(fBeamType==1)
        {
-           bins=new Int_t[2];
            htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
            bins[0]=htemp->GetNbinsX();
            htemp  = (TH1D *) fHadronEffbyIPcut->Projection(1);
@@ -1051,6 +1124,14 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       lRaw->AddEntry(charmbg,"charm elecs");
       // subtract charm background
       spectrumSubtracted->Add(charmbgContainer,-1.0);
+      if(fBeamType==1){
+        THnSparseF* sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();
+        for(Int_t icent =  1; icent < kCentrality-1; icent++){ 
+          sparseCharmElec->GetAxis(0)->SetRange(icent,icent);
+          charmCent[icent-1] = (TH1D *) sparseCharmElec->Projection(ptpr);
+          CorrectFromTheWidth(charmCent[icent-1]);
+        }
+      }
     }
     if(fIPanaConversionBgSubtract){
       // Conversion background
@@ -1065,7 +1146,14 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       lRaw->AddEntry(conversionbg,"conversion elecs");
       // subtract conversion background
       spectrumSubtracted->Add(conversionbgContainer,-1.0);
-      printf("Conversion background subtraction is preliminary!\n");
+      if(fBeamType==1){
+        THnSparseF* sparseconvElec = (THnSparseF *) conversionbgContainer->GetGrid();
+        for(Int_t icent =  1; icent < kCentrality-1; icent++){
+          sparseconvElec->GetAxis(0)->SetRange(icent,icent);
+          convCent[icent-1] = (TH1D *) sparseconvElec->Projection(ptpr);
+          CorrectFromTheWidth(convCent[icent-1]);
+        }
+      }
     }
     if(fIPanaNonHFEBgSubtract){
       // NonHFE background
@@ -1080,8 +1168,16 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
       lRaw->AddEntry(nonhfebg,"non-HF elecs");
       // subtract Dalitz/dielectron background
       spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
-      printf("Non HFE background subtraction is preliminary!\n");
+      if(fBeamType==1){
+        THnSparseF* sparseNonHFEElec = (THnSparseF *) nonHFEbgContainer->GetGrid();
+        for(Int_t icent =  1; icent < kCentrality-1; icent++){
+          sparseNonHFEElec->GetAxis(0)->SetRange(icent,icent);
+          nonHFECent[icent-1] = (TH1D *) sparseNonHFEElec->Projection(ptpr);
+          CorrectFromTheWidth(nonHFECent[icent-1]);
+        }
+      }
     }
+
     TH1D *rawbgsubtracted = (TH1D *) spectrumSubtracted->Project(ptpr);
     CorrectFromTheWidth(rawbgsubtracted);
     rawbgsubtracted->SetMarkerStyle(24);
@@ -1089,6 +1185,52 @@ AliCFDataGrid* AliHFEspectrum::SubtractBackground(Bool_t setBackground){
     lRaw->AddEntry(rawbgsubtracted,"subtracted raw spectrum");
     rawbgsubtracted->Draw("samep");
     lRaw->Draw("SAME");
+    gPad->SetGrid();
+    //rawspectra->SaveAs("rawspectra.eps");
+    
+    if(fBeamType==1){
+      THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
+      for(Int_t icent =  1; icent < kCentrality-1; icent++){
+        sparseSubtracted->GetAxis(0)->SetRange(icent,icent);
+        subtractedCent[icent-1] = (TH1D *) sparseSubtracted->Projection(ptpr);
+        CorrectFromTheWidth(subtractedCent[icent-1]);
+      }
+    
+      TLegend *lCentRaw = new TLegend(0.55,0.55,0.85,0.85);
+      TCanvas *centRaw = new TCanvas("centRaw","centRaw",1000,800);
+      centRaw->Divide(3,3);
+      for(Int_t icent = 1; icent < kCentrality-1; icent++){
+        centRaw->cd(icent);
+        gPad->SetLogx();
+        gPad->SetLogy();
+        incElecCent[icent-1]->GetXaxis()->SetRangeUser(0.4,8.);
+        incElecCent[icent-1]->Draw("p");
+        incElecCent[icent-1]->SetMarkerColor(1);
+        incElecCent[icent-1]->SetMarkerStyle(20);
+        charmCent[icent-1]->Draw("samep");
+        charmCent[icent-1]->SetMarkerColor(3);
+        charmCent[icent-1]->SetMarkerStyle(20);
+        convCent[icent-1]->Draw("samep");
+        convCent[icent-1]->SetMarkerColor(4);
+        convCent[icent-1]->SetMarkerStyle(20);
+        nonHFECent[icent-1]->Draw("samep");
+        nonHFECent[icent-1]->SetMarkerColor(6);
+        nonHFECent[icent-1]->SetMarkerStyle(20);
+        subtractedCent[icent-1]->Draw("samep");
+        subtractedCent[icent-1]->SetMarkerStyle(24);
+        if(icent == 1){
+          lCentRaw->AddEntry(incElecCent[0],"inclusive electron spectrum");
+          lCentRaw->AddEntry(charmCent[0],"charm elecs");
+          lCentRaw->AddEntry(convCent[0],"conversion elecs");
+          lCentRaw->AddEntry(nonHFECent[0],"non-HF elecs");
+          lCentRaw->AddEntry(subtractedCent[0],"subtracted electron spectrum");
+          lCentRaw->Draw("SAME");
+        }
+      }
+    }
+
+    delete[] bins; 
+
   }
   else{
     // Subtract 
@@ -1298,14 +1440,12 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   charmBackgroundGrid = new AliCFDataGrid("charmBackgroundGrid","charmBackgroundGrid",*mcContainer, fStepMC-1); // use MC eff. up to right before PID
 
   TH1D *charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
-  Int_t* bins=new Int_t[1];
+  Int_t* bins=new Int_t[2];
   bins[0]=charmbgaftertofpid->GetNbinsX();
 
   if(fBeamType==1)
   {
-      charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(0);
-      bins=new Int_t[2];
-      bins[0]=charmbgaftertofpid->GetNbinsX();
+      bins[0]=12;
       charmbgaftertofpid = (TH1D *) charmBackgroundGrid->Project(1);
       bins[1]=charmbgaftertofpid->GetNbinsX();
 
@@ -1315,24 +1455,41 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
   TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(ptpr);
 
-  if(fBeamType==0)charmBackgroundGrid->Multiply(ipWeightedCharmContainer,evtnorm);
-  Double_t contents[2];
-  AliCFDataGrid *eventTemp = new AliCFDataGrid("eventTemp","eventTemp",nDim,bins);
-  if(fBeamType==1)
-  {
-      for(Int_t kCentr=0;kCentr<bins[0];kCentr++)
+  charmBackgroundGrid->Multiply(ipWeightedCharmContainer,1.);
+
+  Int_t *nBinpp=new Int_t[1];
+  Int_t* binspp=new Int_t[1];
+  binspp[0]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
+
+  Int_t *nBinPbPb=new Int_t[2];
+  Int_t* binsPbPb=new Int_t[2];
+  binsPbPb[1]=charmbgaftertofpid->GetNbinsX();  // number of pt bins
+  binsPbPb[0]=12;
+
+  Int_t looppt=binspp[0];
+  if(fBeamType==1) looppt=binsPbPb[1];
+
+  for(Long_t iBin=1; iBin<= looppt;iBin++){
+      if(fBeamType==0)
       {
-         for(Int_t kpt=0;kpt<bins[1];kpt++)
-         {
-             Double_t evtnormPbPb=0;
-             if(fNMCbgEvents[kCentr]) evtnormPbPb= double(fNEvents[kCentr])/double(fNMCbgEvents[kCentr]);
-             contents[0]=kCentr;
-             contents[1]=kpt;
-              eventTemp->Fill(contents,evtnormPbPb);
-         }
+          nBinpp[0]=iBin;
+          charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
+          charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
+      }
+      if(fBeamType==1)
+      {
+          // loop over centrality
+          for(Long_t iiBin=1; iiBin<= binsPbPb[0];iiBin++){
+              nBinPbPb[0]=iiBin;
+              nBinPbPb[1]=iBin;
+              Double_t evtnormPbPb=0;
+              if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
+              charmBackgroundGrid->SetElementError(nBinPbPb, charmBackgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
+              charmBackgroundGrid->SetElement(nBinPbPb,charmBackgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
+          }
       }
-   charmBackgroundGrid->Multiply(eventTemp,1);
   }
+
   TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(ptpr);
 
   AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
@@ -1424,6 +1581,14 @@ AliCFDataGrid* AliHFEspectrum::GetCharmBackground(){
   CorrectStatErr(charmBackgroundGrid);
   if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
 
+  delete[] bins;
+  delete[] nBinpp;
+  delete[] binspp;
+  delete[] nBinPbPb;
+  delete[] binsPbPb;
+
+  if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
+
   return charmBackgroundGrid;
 }
 
@@ -1455,7 +1620,6 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
   }
   
   Int_t stepbackground = 3;
-  if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
 
   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ConversionBgGrid","ConversionBgGrid",*backgroundContainer,stepbackground);
   Int_t *nBinpp=new Int_t[1];
@@ -1484,7 +1648,7 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
              nBinPbPb[0]=iiBin;
              nBinPbPb[1]=iBin;
               Double_t evtnormPbPb=0;
-              if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]);
+              if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
              backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
              backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
          }
@@ -1497,7 +1661,12 @@ AliCFDataGrid* AliHFEspectrum::GetConversionBackground(){
   else weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",2,binsPbPb);
   weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
   backgroundGrid->Multiply(weightedConversionContainer,1.0);
-  
+
+  delete[] nBinpp;
+  delete[] binspp;
+  delete[] nBinPbPb;
+  delete[] binsPbPb; 
+
   return backgroundGrid;
 }
 
@@ -1530,7 +1699,6 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
   
   
   Int_t stepbackground = 3;
-  if(TMath::Abs(fEtaRange[0]) < 0.5) stepbackground = 4;
 
   AliCFDataGrid *backgroundGrid = new AliCFDataGrid("NonHFEBgGrid","NonHFEBgGrid",*backgroundContainer,stepbackground);
   Int_t *nBinpp=new Int_t[1];
@@ -1559,7 +1727,7 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
              nBinPbPb[0]=iiBin;
              nBinPbPb[1]=iBin;
              Double_t evtnormPbPb=0;
-              if(fNMCbgEvents[iiBin]) evtnormPbPb= double(fNEvents[iiBin])/double(fNMCbgEvents[iiBin]);
+              if(fNMCbgEvents[iiBin-1]) evtnormPbPb= double(fNEvents[iiBin-1])/double(fNMCbgEvents[iiBin-1]);
              backgroundGrid->SetElementError(nBinPbPb, backgroundGrid->GetElementError(nBinPbPb)*evtnormPbPb);
              backgroundGrid->SetElement(nBinPbPb,backgroundGrid->GetElement(nBinPbPb)*evtnormPbPb);
          }
@@ -1572,6 +1740,11 @@ AliCFDataGrid* AliHFEspectrum::GetNonHFEBackground(){
   weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
   backgroundGrid->Multiply(weightedNonHFEContainer,1.0);
 
+  delete[] nBinpp;
+  delete[] binspp;
+  delete[] nBinPbPb;
+  delete[] binsPbPb; 
+
   return backgroundGrid;
 }
 
@@ -1911,7 +2084,7 @@ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
          bins[0]=35;
 
          AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
-         beffContainer->SetGrid(GetBeautyIPEff());
+         beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
          efficiencyD->Multiply(beffContainer,1);
       }
   }
@@ -1920,7 +2093,6 @@ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
   // Unfold 
   
   AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
-  //unfolding.SetUseCorrelatedErrors();
   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
   if(fSetSmoothing) unfolding.UseSmoothing();
@@ -2015,7 +2187,7 @@ AliCFDataGrid *AliHFEspectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpe
          bins[0]=35;
 
          AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
-         beffContainer->SetGrid(GetBeautyIPEff());
+         beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
          efficiencyD->Multiply(beffContainer,1);
       }
   }
@@ -2161,7 +2333,7 @@ TGraphErrors *AliHFEspectrum::Normalize(THnSparse * const spectrum,Int_t i) cons
   // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
   // Give the final pt spectrum to be compared
   //
+
   if(fNEvents[i] > 0) {
 
     Int_t ptpr = 0;
@@ -2315,7 +2487,7 @@ AliCFContainer *AliHFEspectrum::GetContainer(AliHFEspectrum::CFContainer_t type)
   return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
 }
 //____________________________________________________________
-AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centrality) {
+AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh) {
   //
   // Slice bin for a given source of electron
   // nDim is the number of dimension the corrections are done
@@ -2361,12 +2533,22 @@ AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, In
       AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[0]));
     }
     if(ivar == 5){
-       if((centrality>= 0) && (centrality<container->GetNBins(ivar))) {
-           varMin[ivar] = binLimits[centrality];
-           varMax[ivar] = binLimits[centrality];
+       if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar))) {
+           varMin[ivar] = binLimits[centralitylow];
+           varMax[ivar] = binLimits[centralityhigh];
+
+           TAxis *axistest = container->GetAxis(5,0);
+           printf("GetSlicedContainer: Number of bin in centrality direction %d\n",axistest->GetNbins());
+           printf("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]);
+           Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
+           Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
+           printf("GetSlicedContainer: Low centrality %f and high centrality %f\n",lowcentrality,highcentrality);
+       
        }
+
     }
 
+
     delete[] binLimits;
   }
   
@@ -2379,7 +2561,7 @@ AliCFContainer *AliHFEspectrum::GetSlicedContainer(AliCFContainer *container, In
 }
 
 //_________________________________________________________________________
-THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions) const {
+THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions, Int_t centralitylow, Int_t centralityhigh) const {
   //
   // Slice correlation
   //
@@ -2390,6 +2572,43 @@ THnSparseF *AliHFEspectrum::GetSlicedCorrelation(THnSparseF *correlationmatrix,
     AliError("Problem in the dimensions");
     return NULL;
   }
+  
+  // Cut in centrality is centrality > -1
+  if((centralitylow >=0) && (centralityhigh >=0)) {
+
+    TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
+    TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
+
+    Int_t bins0 = axiscentrality0->GetNbins();
+    Int_t bins1 = axiscentrality1->GetNbins();
+    
+    printf("Number of centrality bins: %d and %d\n",bins0,bins1);
+    if(bins0 != bins1) {
+      AliError("Problem in the dimensions");
+      return NULL;
+    }
+    
+    if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0)) {
+      axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
+      axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
+
+      Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
+      Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
+      Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
+      Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
+      printf("GetCorrelation0: Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0);
+      printf("GetCorrelation1: Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1);
+
+      TCanvas * ctestcorrelation = new TCanvas("testcorrelationprojection","testcorrelationprojection",1000,700);
+      ctestcorrelation->cd(1);
+      TH2D* projection = (TH2D *) correlationmatrix->Projection(5,5+((Int_t)(ndimensions/2.)));
+      projection->Draw("colz");
+
+    }
+    
+  }
+
+
   Int_t ndimensionsContainer = (Int_t) ndimensions/2;
   //printf("Number of dimension %d container\n",ndimensionsContainer);
 
@@ -2517,9 +2736,15 @@ THnSparse* AliHFEspectrum::GetCharmWeights(){
       loopcentr=nBinPbPb[0];
   }
 
+  // Weighting factor for pp
+  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};
+  
+  // Weighting factor for PbPb (0-20%)
+  //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};
+
+  // Weighting factor for PbPb (40-80%)
+  //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};
 
-     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};
-//{1.11865,1.11444,1.13395,1.15751,1.13412,1.14446,1.15372,1.09458,1.13446,1.11707,1.1352,1.10979,1.08887,1.11164,1.10772,1.10727,1.07027,1.07016,1.04826,1.03248,1.0093,0.973534,0.932574,0.869838,0.799316,0.734015,0.643001,0.584319,0.527147,0.495661,0.470002,0.441129,0.431156,0.417242,0.39246,0.379611,0.390845,0.368857,0.34959,0.387186,0.376364,0.384437,0.349585,0.383225};
   //points
   Double_t pt[1];
   Double_t contents[2];
@@ -2562,11 +2787,14 @@ THnSparse* AliHFEspectrum::GetCharmWeights(){
     fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
   }
 
+  delete[] binsvar;
+  delete[] binfill;
+
   return fWeightCharm;
 }
 
 //____________________________________________________________________________
-void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, Int_t *dimensions){
+void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
 
    // TOF PID efficiencies
    Int_t ptpr;
@@ -2580,45 +2808,71 @@ void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContain
      loopcentr=nCentralitybinning;
    }
 
-   TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,20.);
-   TF1 *fipfit2 = new TF1("fipfit2","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,10.0);
-   TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",1.0,7.);
+   TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
+   TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.9,8.);
+   TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.3,10.0);
 
-   TCanvas * cefficiencyParam = new TCanvas("efficiencyParam","efficiencyParam",1200,600);
-   cefficiencyParam->Divide(2,1);
-   cefficiencyParam->cd(1);
+   TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
+   cefficiencyParamtof->cd();
 
    AliCFContainer *mccontainermcD = 0x0;
+   AliCFContainer *mccontaineresdD = 0x0;
    TH1D* efficiencysigTOFPIDD[nCentralitybinning];
    TH1D* efficiencyTOFPIDD[nCentralitybinning];
+   TH1D* efficiencysigesdTOFPIDD[nCentralitybinning];
+   TH1D* efficiencyesdTOFPIDD[nCentralitybinning];
    Int_t source = -1; //get parameterized TOF PID efficiencies
 
    for(int icentr=0; icentr<loopcentr; icentr++) {
       // signal sample
       if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
-      if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
       AliCFEffGrid* efficiencymcsigParamTOFPID= new AliCFEffGrid("efficiencymcsigParamTOFPID","",*mccontainermcD);
       efficiencymcsigParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
 
       // mb sample for double check
       if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
-      if(fBeamType==1)  mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
       AliCFEffGrid* efficiencymcParamTOFPID= new AliCFEffGrid("efficiencymcParamTOFPID","",*mccontainermcD);
       efficiencymcParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
 
+      // mb sample with reconstructed variables
+      if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
+      else mccontainermcD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      AliCFEffGrid* efficiencyesdParamTOFPID= new AliCFEffGrid("efficiencyesdParamTOFPID","",*mccontainermcD);
+      efficiencyesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+      // mb sample with reconstructed variables
+      if(fBeamType==0)  mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
+      else mccontainermcD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,icentr+1);
+      AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
+      efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+
+      //fill histo
       efficiencysigTOFPIDD[icentr] = (TH1D *) efficiencymcsigParamTOFPID->Project(ptpr);
       efficiencyTOFPIDD[icentr] = (TH1D *) efficiencymcParamTOFPID->Project(ptpr);
+      efficiencysigesdTOFPIDD[icentr] = (TH1D *) efficiencysigesdParamTOFPID->Project(ptpr);
+      efficiencyesdTOFPIDD[icentr] = (TH1D *) efficiencyesdParamTOFPID->Project(ptpr);
       efficiencysigTOFPIDD[icentr]->SetName(Form("efficiencysigTOFPIDD%d",icentr));
       efficiencyTOFPIDD[icentr]->SetName(Form("efficiencyTOFPIDD%d",icentr));
+      efficiencysigesdTOFPIDD[icentr]->SetName(Form("efficiencysigesdTOFPIDD%d",icentr));
+      efficiencyesdTOFPIDD[icentr]->SetName(Form("efficiencyesdTOFPIDD%d",icentr));
 
+      //fit (mc enhenced sample)
       fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
-      fittofpid->SetLineColor(3);
       efficiencysigTOFPIDD[icentr]->Fit(fittofpid,"R");
       efficiencysigTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
       fEfficiencyTOFPIDD[icentr] = efficiencysigTOFPIDD[icentr]->GetFunction("fittofpid");
+
+      //fit (esd enhenced sample)
+      efficiencysigesdTOFPIDD[icentr]->Fit(fittofpid,"R");
+      efficiencysigesdTOFPIDD[icentr]->GetYaxis()->SetTitle("Efficiency");
+      fEfficiencyesdTOFPIDD[icentr] = efficiencysigesdTOFPIDD[icentr]->GetFunction("fittofpid");
+
    }
 
    // draw (for PbPb, only 1st bin)
+   //sig mc
    efficiencysigTOFPIDD[0]->SetTitle("");
    efficiencysigTOFPIDD[0]->SetStats(0);
    efficiencysigTOFPIDD[0]->SetMarkerStyle(25);
@@ -2626,61 +2880,150 @@ void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContain
    efficiencysigTOFPIDD[0]->SetLineColor(2);
    efficiencysigTOFPIDD[0]->Draw();
 
-   fEfficiencyTOFPIDD[0]->Draw("same");
-
+   //mb mc
    efficiencyTOFPIDD[0]->SetTitle("");
    efficiencyTOFPIDD[0]->SetStats(0);
    efficiencyTOFPIDD[0]->SetMarkerStyle(24);
+   efficiencyTOFPIDD[0]->SetMarkerColor(4);
+   efficiencyTOFPIDD[0]->SetLineColor(4);
    efficiencyTOFPIDD[0]->Draw("same");
 
+   //sig esd
+   efficiencysigesdTOFPIDD[0]->SetTitle("");
+   efficiencysigesdTOFPIDD[0]->SetStats(0);
+   efficiencysigesdTOFPIDD[0]->SetMarkerStyle(25);
+   efficiencysigesdTOFPIDD[0]->SetMarkerColor(3);
+   efficiencysigesdTOFPIDD[0]->SetLineColor(3);
+   efficiencysigesdTOFPIDD[0]->Draw("same");
+
+   //mb esd
+   efficiencyesdTOFPIDD[0]->SetTitle("");
+   efficiencyesdTOFPIDD[0]->SetStats(0);
+   efficiencyesdTOFPIDD[0]->SetMarkerStyle(25);
+   efficiencyesdTOFPIDD[0]->SetMarkerColor(1);
+   efficiencyesdTOFPIDD[0]->SetLineColor(1);
+   efficiencyesdTOFPIDD[0]->Draw("same");
+
+   //signal mc fit
+   if(fEfficiencyTOFPIDD[0]){
+     fEfficiencyTOFPIDD[0]->SetLineColor(2);
+     fEfficiencyTOFPIDD[0]->Draw("same");
+   }
+   //mb esd fit
+   if(fEfficiencyesdTOFPIDD[0]){
+       fEfficiencyesdTOFPIDD[0]->SetLineColor(3);
+       fEfficiencyesdTOFPIDD[0]->Draw("same");
+     }
 
-   cefficiencyParam->cd(2);
+   TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
+   legtofeff->AddEntry(efficiencysigTOFPIDD[0],"TOF PID Step Efficiency","");
+   legtofeff->AddEntry(efficiencysigTOFPIDD[0],"vs MC p_{t} for enhenced samples","p");
+   legtofeff->AddEntry(efficiencyTOFPIDD[0],"vs MC p_{t} for mb samples","p");
+   legtofeff->AddEntry(efficiencysigesdTOFPIDD[0],"vs esd p_{t} for enhenced samples","p");
+   legtofeff->AddEntry(efficiencyesdTOFPIDD[0],"vs esd p_{t} for mb samples","p");
+   legtofeff->Draw("same");
+
+
+   TCanvas * cefficiencyParamIP = new TCanvas("efficiencyParamIP","efficiencyParamIP",500,500);
+   cefficiencyParamIP->cd();
+   gStyle->SetOptStat(0);
 
    // IP cut efficiencies
-  
    for(int icentr=0; icentr<loopcentr; icentr++)  {
-     source = 0; //charm
+
+     AliCFContainer *charmCombined = 0x0; 
+     AliCFContainer *beautyCombined = 0x0;
+     AliCFContainer *beautyCombinedesd = 0x0;
+
+     printf("centrality printing %i \n",icentr);
+
+     source = 0; //charm enhenced
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
      efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
 
-     source = 1; //beauty
+     charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");  
+
+     source = 1; //beauty enhenced
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
      efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
 
-     source = 0; //charm
+     beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined"); 
+
+     if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen);
+     else mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyBeautySigesd = new AliCFEffGrid("efficiencyBeautySigesd","",*mccontaineresdD);
+     efficiencyBeautySigesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+     beautyCombinedesd = (AliCFContainer*)mccontaineresdD->Clone("beautyCombinedesd");
+
+     source = 0; //charm mb
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
      efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
 
-     source = 1; //beauty
+     charmCombined->Add(mccontainermcD); 
+     AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined); 
+     efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
+
+     source = 1; //beauty mb
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
      efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
 
-     source = 2; //conversion
+     beautyCombined->Add(mccontainermcD);
+     AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined); 
+     efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
+
+     if(fBeamType==0) mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen);
+     else mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     AliCFEffGrid* efficiencyBeautyesd = new AliCFEffGrid("efficiencyBeautyesd","",*mccontaineresdD);
+     efficiencyBeautyesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+     beautyCombinedesd->Add(mccontaineresdD);
+     AliCFEffGrid* efficiencyBeautyCombinedesd = new AliCFEffGrid("efficiencyBeautyCombinedesd","",*beautyCombinedesd);
+     efficiencyBeautyCombinedesd->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1);
+
+     source = 2; //conversion mb
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyConv = new AliCFEffGrid("efficiencyConv","",*mccontainermcD);
      efficiencyConv->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
 
-     source = 3; //non HFE except for the conversion
+     source = 3; //non HFE except for the conversion mb
      if(fBeamType==0) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen);
-     if(fBeamType==1) mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
+     else mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen, icentr+1);
      AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
      efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
 
-     fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr);
-     fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr);
-     fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr);
-     fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr);
-     fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr);
-     fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr);
+     if(fIPEffCombinedSamples){
+       fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmCombined->Project(ptpr); //signal enhenced + mb 
+       fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautyCombined->Project(ptpr); //signal enhenced + mb
+       fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautyCombinedesd->Project(ptpr); //signal enhenced + mb
+     }
+     else{
+       fEfficiencyCharmSigD[icentr] = (TH1D*)efficiencyCharmSig->Project(ptpr); //signal enhenced only
+       fEfficiencyBeautySigD[icentr] = (TH1D*)efficiencyBeautySig->Project(ptpr); //signal enhenced only
+       fEfficiencyBeautySigesdD[icentr] = (TH1D*)efficiencyBeautySigesd->Project(ptpr); //signal enhenced only
+     }
+     fCharmEff[icentr] = (TH1D*)efficiencyCharm->Project(ptpr); //mb only
+     fBeautyEff[icentr] = (TH1D*)efficiencyBeauty->Project(ptpr); //mb only
+     fConversionEff[icentr] = (TH1D*)efficiencyConv->Project(ptpr); //mb only
+     fNonHFEEff[icentr] = (TH1D*)efficiencyNonhfe->Project(ptpr); //mb only
+
+   }
+
+   if(fBeamType==0){
+     AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
+     fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
+
+     AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
+     fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
    }
 
    for(int icentr=0; icentr<loopcentr; icentr++)  {
@@ -2688,38 +3031,48 @@ void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContain
      fipfit->SetLineColor(2);
      fEfficiencyBeautySigD[icentr]->Fit(fipfit,"R");
      fEfficiencyBeautySigD[icentr]->GetYaxis()->SetTitle("Efficiency");
-     if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit");
+     if(fBeauty2ndMethod)fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[0]->GetFunction("fipfit"); //why do we need this line?
      else fEfficiencyIPBeautyD[icentr] = fEfficiencyBeautySigD[icentr]->GetFunction("fipfit");
 
+     fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+     fipfit->SetLineColor(6);
+     fEfficiencyBeautySigesdD[icentr]->Fit(fipfit,"R");
+     fEfficiencyBeautySigesdD[icentr]->GetYaxis()->SetTitle("Efficiency");
+     if(fBeauty2ndMethod)fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[0]->GetFunction("fipfit"); //why do we need this line?
+     else fEfficiencyIPBeautyesdD[icentr] = fEfficiencyBeautySigesdD[icentr]->GetFunction("fipfit");
+
+     fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+     fipfit->SetLineColor(1);
+     fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
+     fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
+     fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
+     
      if(fIPParameterizedEff){
-       fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
-       fipfit->SetLineColor(1);
-       fEfficiencyCharmSigD[icentr]->Fit(fipfit,"R");
-       fEfficiencyCharmSigD[icentr]->GetYaxis()->SetTitle("Efficiency");
-       fEfficiencyIPCharmD[icentr] = fEfficiencyCharmSigD[icentr]->GetFunction("fipfit");
-
-       fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
-       fipfit2->SetLineColor(3);
-       fConversionEff[icentr]->Fit(fipfit2,"R");
+       fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fipfitnonhfe->SetLineColor(3);
+       fConversionEff[icentr]->Fit(fipfitnonhfe,"R");
        fConversionEff[icentr]->GetYaxis()->SetTitle("Efficiency");
-       fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfit2");
+       fEfficiencyIPConversionD[icentr] = fConversionEff[icentr]->GetFunction("fipfitnonhfe");
 
-       fipfit2->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
-       fipfit2->SetLineColor(4);
-       fNonHFEEff[icentr]->Fit(fipfit2,"R");
+       fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fipfitnonhfe->SetLineColor(4);
+       fNonHFEEff[icentr]->Fit(fipfitnonhfe,"R");
        fNonHFEEff[icentr]->GetYaxis()->SetTitle("Efficiency");
-       fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfit2");
+       fEfficiencyIPNonhfeD[icentr] = fNonHFEEff[icentr]->GetFunction("fipfitnonhfe");
      }
    }
 
    // draw (for PbPb, only 1st bin)
-
    fEfficiencyCharmSigD[0]->SetMarkerStyle(21);
    fEfficiencyCharmSigD[0]->SetMarkerColor(1);
    fEfficiencyCharmSigD[0]->SetLineColor(1);
    fEfficiencyBeautySigD[0]->SetMarkerStyle(21);
    fEfficiencyBeautySigD[0]->SetMarkerColor(2);
    fEfficiencyBeautySigD[0]->SetLineColor(2);
+   fEfficiencyBeautySigesdD[0]->SetStats(0);
+   fEfficiencyBeautySigesdD[0]->SetMarkerStyle(21);
+   fEfficiencyBeautySigesdD[0]->SetMarkerColor(6);
+   fEfficiencyBeautySigesdD[0]->SetLineColor(6);
    fCharmEff[0]->SetMarkerStyle(24);
    fCharmEff[0]->SetMarkerColor(1);
    fCharmEff[0]->SetLineColor(1);
@@ -2734,22 +3087,54 @@ void AliHFEspectrum::SetParameterizedEff(AliCFContainer *container, AliCFContain
    fNonHFEEff[0]->SetLineColor(4);
 
    fEfficiencyCharmSigD[0]->Draw();
-   fEfficiencyBeautySigD[0]->Draw("same");
-   fCharmEff[0]->Draw("same");
-   fBeautyEff[0]->Draw("same");
-   fConversionEff[0]->Draw("same");
-   fNonHFEEff[0]->Draw("same");
+   fEfficiencyCharmSigD[0]->GetXaxis()->SetRangeUser(0.0,7.9);
+   fEfficiencyCharmSigD[0]->GetYaxis()->SetRangeUser(0.0,0.5);
 
-   fEfficiencyIPBeautyD[0]->Draw("same");
-   if(fIPParameterizedEff){
+   fEfficiencyBeautySigD[0]->Draw("same");
+   fEfficiencyBeautySigesdD[0]->Draw("same");
+   //fCharmEff[0]->Draw("same");
+   //fBeautyEff[0]->Draw("same");
+
+   if(fBeamType==0){
+     fConversionEffbgc->SetMarkerStyle(25);
+     fConversionEffbgc->SetMarkerColor(3);
+     fConversionEffbgc->SetLineColor(3);
+     fNonHFEEffbgc->SetMarkerStyle(25);
+     fNonHFEEffbgc->SetMarkerColor(4);
+     fNonHFEEffbgc->SetLineColor(4);
+     fConversionEffbgc->Draw("same");
+     fNonHFEEffbgc->Draw("same");
+   }
+   else{
+     fConversionEff[0]->Draw("same");
+     fNonHFEEff[0]->Draw("same");
+   }
+   if(fEfficiencyIPBeautyD[0])
+      fEfficiencyIPBeautyD[0]->Draw("same");
+   if(fEfficiencyIPBeautyesdD[0])
+     fEfficiencyIPBeautyesdD[0]->Draw("same");
+   if( fEfficiencyIPCharmD[0])
      fEfficiencyIPCharmD[0]->Draw("same");
+   if(fIPParameterizedEff){
      fEfficiencyIPConversionD[0]->Draw("same");
      fEfficiencyIPNonhfeD[0]->Draw("same");
    }
+   TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
+   legipeff->AddEntry(fEfficiencyBeautySigD[0],"IP Step Efficiency","");
+   legipeff->AddEntry(fEfficiencyBeautySigD[0],"beauty e","p");
+   legipeff->AddEntry(fEfficiencyBeautySigesdD[0],"beauty e(esd pt)","p");
+   legipeff->AddEntry(fEfficiencyCharmSigD[0],"charm e","p");
+   legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
+   legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
+   //legipeff->AddEntry(fConversionEff[0],"conversion e","p");
+   //legipeff->AddEntry(fNonHFEEff[0],"Dalitz e","p");
+   legipeff->Draw("same");
+   gPad->SetGrid();
+   //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
 }
 
 //____________________________________________________________________________
-THnSparse* AliHFEspectrum::GetBeautyIPEff(){
+THnSparse* AliHFEspectrum::GetBeautyIPEff(Bool_t isMCpt){
   //
   // Return beauty electron IP cut efficiency
   //
@@ -2784,16 +3169,19 @@ THnSparse* AliHFEspectrum::GetBeautyIPEff(){
   }
   Double_t pt[1];
   Double_t weight[nCentralitybinning];
+  Double_t weightErr[nCentralitybinning];
   Double_t contents[2];
 
   for(Int_t a=0;a<11;a++)
   {
       weight[a] = 1.0;
+      weightErr[a] = 1.0;
   }
 
 
   Int_t looppt=nBin[0];
   Int_t loopcentr=1;
+  Int_t ibin[2];
   if(fBeamType==1)
   {
       loopcentr=nBinPbPb[0];
@@ -2805,25 +3193,44 @@ THnSparse* AliHFEspectrum::GetBeautyIPEff(){
       for(int i=0; i<looppt; i++)
       {
          pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
-          if(fEfficiencyIPBeautyD[icentr])
-            weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
+          if(isMCpt){
+            if(fEfficiencyIPBeautyD[icentr]){
+              weight[icentr]=fEfficiencyIPBeautyD[icentr]->Eval(pt[0]);
+              weightErr[icentr] = 0;
+            }
+            else{
+              printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
+              weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); 
+              weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
+            }
+          }
           else{
-            printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
-            weight[icentr] = fEfficiencyBeautySigD[icentr]->GetBinContent(i+1); 
+            if(fEfficiencyIPBeautyesdD[icentr]){
+              weight[icentr]=fEfficiencyIPBeautyesdD[icentr]->Eval(pt[0]);
+              weightErr[icentr] = 0;
+            }
+            else{
+              printf("Fit failed on beauty IP cut efficiency for centrality %d. Contents in histo used!\n",icentr);
+              weight[icentr] = fEfficiencyBeautySigesdD[icentr]->GetBinContent(i+1);
+              weightErr[icentr] = fEfficiencyBeautySigD[icentr]->GetBinError(i+1);
+            }
           }
 
           if(fBeamType==1){
               contents[0]=icentr;
               contents[1]=pt[0];
+              ibin[0]=icentr;
+              ibin[1]=i+1;
           }
           if(fBeamType==0){
               contents[0]=pt[0];
+              ibin[0]=i+1;
           }
           ipcut->Fill(contents,weight[icentr]);
+          ipcut->SetBinError(ibin,weightErr[icentr]);
       }
   } 
 
-
   Int_t nDimSparse = ipcut->GetNdimensions();
   Int_t* binsvar = new Int_t[nDimSparse]; // number of bins for each variable
   Long_t nBins = 1; // used to calculate the total number of bins in the THnSparse
@@ -2841,9 +3248,12 @@ THnSparse* AliHFEspectrum::GetBeautyIPEff(){
       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
       bintmp /= binsvar[iVar] ;
     }
-    ipcut->SetBinError(binfill,0.); // put 0 everywhere
+    //ipcut->SetBinError(binfill,0.); // put 0 everywhere
   }
 
+  delete[] binsvar;
+  delete[] binfill;
+
   return ipcut;
 }
 
@@ -2882,15 +3292,18 @@ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
 
   Double_t pt[1];
   Double_t weight[nCentralitybinning];
+  Double_t weightErr[nCentralitybinning];
   Double_t contents[2];
 
-  for(Int_t a=0;a<11;a++)
+  for(Int_t a=0;a<nCentralitybinning;a++)
   {
       weight[a] = 1.0;
+      weightErr[a] = 1.0;
   }
 
   Int_t looppt=nBin[0];
   Int_t loopcentr=1;
+  Int_t ibin[2];
   if(fBeamType==1)
   {
       loopcentr=nBinPbPb[0];
@@ -2904,54 +3317,101 @@ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
          pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
 
           Double_t tofpideff = 0.;
+          Double_t tofpideffesd = 0.;
           if(fEfficiencyTOFPIDD[icentr])
-            tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]);
+            tofpideff = fEfficiencyTOFPIDD[icentr]->Eval(pt[0]); 
           else{
             printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
           }  
+          if(fEfficiencyesdTOFPIDD[icentr])
+            tofpideffesd = fEfficiencyesdTOFPIDD[icentr]->Eval(pt[0]);
+          else{
+            printf("TOF PID fit failed on conversion for centrality %d. The result is wrong!\n",icentr);
+          }
 
           //tof pid eff x tpc pid eff x ip cut eff
           if(fIPParameterizedEff){
             if(source==0) {
-              if(fEfficiencyIPCharmD[icentr])
+              if(fEfficiencyIPCharmD[icentr]){
                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
+                weightErr[icentr] = 0; 
+              }
               else{
                 printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
                 weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
+                weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1); 
               }
             } 
            else if(source==2) {
-              if(fEfficiencyIPConversionD[icentr])
-                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); 
+              if(fEfficiencyIPConversionD[icentr]){
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD[icentr]->Eval(pt[0]); 
+                weightErr[icentr] = 0; 
+              }
               else{
                 printf("Fit failed on conversion IP cut efficiency for centrality %d\n",icentr);
-                weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1);
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
               }
             }
            else if(source==3) {
-              if(fEfficiencyIPNonhfeD[icentr])
-                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); 
+              if(fEfficiencyIPNonhfeD[icentr]){
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD[icentr]->Eval(pt[0]); 
+                weightErr[icentr] = 0; 
+              }
               else{
                 printf("Fit failed on dalitz IP cut efficiency for centrality %d\n",icentr);
-                weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1);
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
               }  
             }
           }
           else{
-            if(source==0) weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1); //charm
-           else if(source==2) weight[icentr] = tofpideff*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
-           else if(source==3) weight[icentr] = tofpideff*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
+            if(source==0){ 
+              if(fEfficiencyIPCharmD[icentr]){
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD[icentr]->Eval(pt[0]);
+                weightErr[icentr] = 0;
+              }
+              else{
+                printf("Fit failed on charm IP cut efficiency for centrality %d\n",icentr);
+                weight[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinContent(i+1);
+                weightErr[icentr] = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD[icentr]->GetBinError(i+1);
+              }
+            }
+           else if(source==2){
+              if(fBeamType==0){
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
+              }
+              else{
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinContent(i+1); // conversion
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fConversionEff[icentr]->GetBinError(i+1);
+              }
+            }
+           else if(source==3){
+              if(fBeamType==0){
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
+              }
+              else{ 
+                weight[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinContent(i+1); // Dalitz
+                weightErr[icentr] = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff[icentr]->GetBinError(i+1);
+              }
+            }
           }
 
          if(fBeamType==1){
              contents[0]=icentr;
              contents[1]=pt[0];
+              ibin[0]=icentr;
+              ibin[1]=i+1;
          }
          if(fBeamType==0){
              contents[0]=pt[0];
+              ibin[0]=i+1;
          }
 
          pideff->Fill(contents,weight[icentr]);
+          pideff->SetBinError(ibin,weightErr[icentr]);
       }
   }
 
@@ -2972,10 +3432,10 @@ THnSparse* AliHFEspectrum::GetPIDxIPEff(Int_t source){
       binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
       bintmp /= binsvar[iVar] ;
     }
-    pideff->SetBinError(binfill,0.); // put 0 everywhere
   }
 
-
+  delete[] binsvar;
+  delete[] binfill;
 
 
   return pideff;
@@ -2998,13 +3458,13 @@ AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
        nDim=2;
     }
 
-    const Int_t nPtbinning1 = 35;//number of pt bins, according to new binning
+    const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
     const Int_t nCentralitybinning=11;//number of centrality bins
-    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
+    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};
+   
     Double_t kCentralityRange[nCentralitybinning+1] = {0.,1.,2., 3., 4., 5., 6., 7.,8.,9., 10., 11.};
     Int_t nBin[1] = {nPtbinning1};
     Int_t nBinPbPb[2] = {nCentralitybinning,nPtbinning1};
-    //fBSpectrum2ndMethod->GetNbinsX()
 
     AliCFDataGrid *rawBeautyContainer;
     if(fBeamType==0)  rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
@@ -3085,7 +3545,10 @@ AliCFDataGrid *AliHFEspectrum::GetRawBspectra2ndMethod(){
     fBSpectrum2ndMethod->Draw("p");
     hRawBeautySpectra->SetMarkerStyle(25);
     hRawBeautySpectra->Draw("samep");
-    
+
+    delete[] binfill;
+    delete[] binsvar; 
+
     return rawBeautyContainer;
 }
 
@@ -3100,17 +3563,18 @@ void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
     return;
 
   Double_t evtnorm[1] = {0.0};
-  if(fNMCbgEvents) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
+  if(fNMCbgEvents[0]>0) evtnorm[0]= double(fNEvents[0])/double(fNMCbgEvents[0]);
   
   AliCFDataGrid *convSourceGrid[kElecBgSources][kBgLevels];
   AliCFDataGrid *nonHFESourceGrid[kElecBgSources][kBgLevels];
 
-  AliCFDataGrid *bgLevelGrid[kBgLevels];
+  AliCFDataGrid *bgLevelGrid[2][kBgLevels];//for pi0 and eta based errors
   AliCFDataGrid *bgNonHFEGrid[kBgLevels];
   AliCFDataGrid *bgConvGrid[kBgLevels];
 
   Int_t stepbackground = 3;
   Int_t* bins=new Int_t[1];
+  const Char_t *bgBase[2] = {"pi0","eta"};
  
   bins[0]=fConversionEff[centrality]->GetNbinsX();
    
@@ -3129,98 +3593,132 @@ void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
     }
     
     bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
-    for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+    for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){
       bgConvGrid[iLevel]->Add(convSourceGrid[iSource][iLevel]);
     }
+    if(!fEtaSyst)
+      bgConvGrid[iLevel]->Add(convSourceGrid[1][iLevel]);
     
     bgNonHFEGrid[iLevel] = (AliCFDataGrid*)nonHFESourceGrid[0][iLevel]->Clone(); 
-    for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){//add other sources to get overall background from all meson decays
+    for(Int_t iSource = 2; iSource < kElecBgSources; iSource++){//add other sources to pi0, to get overall background from all meson decays, exception: eta (independent error calculation)
       bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[iSource][iLevel]);
     }
+    if(!fEtaSyst)
+      bgNonHFEGrid[iLevel]->Add(nonHFESourceGrid[1][iLevel]);
     
-    bgLevelGrid[iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
-    bgLevelGrid[iLevel]->Add(bgNonHFEGrid[iLevel]);
+    bgLevelGrid[0][iLevel] = (AliCFDataGrid*)bgConvGrid[iLevel]->Clone();
+    bgLevelGrid[0][iLevel]->Add(bgNonHFEGrid[iLevel]);
+    if(fEtaSyst){
+      bgLevelGrid[1][iLevel] = (AliCFDataGrid*)nonHFESourceGrid[1][iLevel]->Clone();//background for eta source
+      bgLevelGrid[1][iLevel]->Add(convSourceGrid[1][iLevel]);
+    }
   }
   
-  
-  //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)
-  AliCFDataGrid *bgErrorGrid[2];
-  bgErrorGrid[0] = (AliCFDataGrid*)bgLevelGrid[1]->Clone();
-  bgErrorGrid[1] = (AliCFDataGrid*)bgLevelGrid[2]->Clone();
-  bgErrorGrid[0]->Add(bgLevelGrid[0],-1.);
-  bgErrorGrid[1]->Add(bgLevelGrid[0],-1.);
+  //Now subtract the mean from upper, and lower from mean container to get the error based on the pion yield uncertainty (-> this error sums linearly, since its contribution to all meson yields is correlated; exception: eta errors in pp 7 TeV sum with others the gaussian way, as they are independent from pi0) 
+  AliCFDataGrid *bgErrorGrid[2][2];//for pions/eta error base, for lower/upper
+  TH1D* hBaseErrors[2][2];//pi0/eta and lower/upper
+  for(Int_t iErr = 0; iErr < 2; iErr++){//errors for pi0 and eta base
+    bgErrorGrid[iErr][0] = (AliCFDataGrid*)bgLevelGrid[iErr][1]->Clone();
+    bgErrorGrid[iErr][0]->Add(bgLevelGrid[iErr][0],-1.);
+    bgErrorGrid[iErr][1] = (AliCFDataGrid*)bgLevelGrid[iErr][2]->Clone();    
+    bgErrorGrid[iErr][1]->Add(bgLevelGrid[iErr][0],-1.);
+
+  //plot absolute differences between limit yields (upper/lower limit, based on pi0 and eta errors) and best estimate
  
-  //plot absolute differences between limit yields (upper/lower limit, based on pi0 errors) and best estimate
-  TH1D* hpiErrors[2];
-  hpiErrors[0] = (TH1D*)bgErrorGrid[0]->Project(0);
-  hpiErrors[0]->Scale(-1.);
-  hpiErrors[0]->SetTitle("Absolute systematic errors from non-HF meson decays and conversions");
-  hpiErrors[1] = (TH1D*)bgErrorGrid[1]->Project(0);
-
+    hBaseErrors[iErr][0] = (TH1D*)bgErrorGrid[iErr][0]->Project(0);
+    hBaseErrors[iErr][0]->Scale(-1.);
+    hBaseErrors[iErr][0]->SetTitle(Form("Absolute %s-based systematic errors from non-HF meson decays and conversions",bgBase[iErr]));
+    hBaseErrors[iErr][1] = (TH1D*)bgErrorGrid[iErr][1]->Project(0);
+    if(!fEtaSyst)break;
+  }
   
-
-   //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
+  //Calculate the scaling errors for electrons from all mesons except for pions (and in pp 7 TeV case eta): square sum of (0.3 * best yield estimate), where 0.3 is the error generally assumed for m_t scaling
   TH1D *hSpeciesErrors[kElecBgSources-1];
   for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+    if(fEtaSyst && (iSource == 1))continue;
     hSpeciesErrors[iSource-1] = (TH1D*)convSourceGrid[iSource][0]->Project(0);
     TH1D *hNonHFEtemp = (TH1D*)nonHFESourceGrid[iSource][0]->Project(0);
     hSpeciesErrors[iSource-1]->Add(hNonHFEtemp);
     hSpeciesErrors[iSource-1]->Scale(0.3);   
   }
-
-  TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[0]->Clone();
-  TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[0]->Clone();
-  TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[0]->Clone();
+  
+  //Int_t firstBgSource = 0;//if eta systematics are not from scaling
+  //if(fEtaSyst){firstBgSource = 1;}//source 0 histograms are not filled if eta errors are independently determined!
+  TH1D *hOverallSystErrLow = (TH1D*)hSpeciesErrors[1]->Clone();
+  TH1D *hOverallSystErrUp = (TH1D*)hSpeciesErrors[1]->Clone();
+  TH1D *hScalingErrors = (TH1D*)hSpeciesErrors[1]->Clone();
 
   TH1D *hOverallBinScaledErrsUp = (TH1D*)hOverallSystErrUp->Clone();
   TH1D *hOverallBinScaledErrsLow = (TH1D*)hOverallSystErrLow->Clone();
 
   for(Int_t iBin = 1; iBin <= kBgPtBins; iBin++){
-    Double_t pi0basedErrLow = hpiErrors[0]->GetBinContent(iBin); 
-    Double_t pi0basedErrUp = hpiErrors[1]->GetBinContent(iBin);
-    
+    Double_t pi0basedErrLow,pi0basedErrUp,etaErrLow,etaErrUp;    
+    pi0basedErrLow = hBaseErrors[0][0]->GetBinContent(iBin); 
+    pi0basedErrUp = hBaseErrors[0][1]->GetBinContent(iBin);
+    if(fEtaSyst){
+      etaErrLow = hBaseErrors[1][0]->GetBinContent(iBin); 
+      etaErrUp = hBaseErrors[1][1]->GetBinContent(iBin);
+    }
+    else{ etaErrLow = etaErrUp = 0.;}
+
     Double_t sqrsumErrs= 0;
     for(Int_t iSource = 1; iSource < kElecBgSources; iSource++){
+      if(fEtaSyst && (iSource == 1))continue;
       Double_t scalingErr=hSpeciesErrors[iSource-1]->GetBinContent(iBin);
       sqrsumErrs+=(scalingErr*scalingErr);
     }
     for(Int_t iErr = 0; iErr < 2; iErr++){
-      hpiErrors[iErr]->SetBinContent(iBin,hpiErrors[iErr]->GetBinContent(iBin)/hpiErrors[iErr]->GetBinWidth(iBin));
+      for(Int_t iLevel = 0; iLevel < 2; iLevel++){
+        hBaseErrors[iErr][iLevel]->SetBinContent(iBin,hBaseErrors[iErr][iLevel]->GetBinContent(iBin)/hBaseErrors[iErr][iLevel]->GetBinWidth(iBin));
+      }
+      if(!fEtaSyst)break;
     }
-    hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+sqrsumErrs));
-    hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+sqrsumErrs));
+    hOverallSystErrUp->SetBinContent(iBin, TMath::Sqrt((pi0basedErrUp*pi0basedErrUp)+(etaErrUp*etaErrUp)+sqrsumErrs));
+    hOverallSystErrLow->SetBinContent(iBin, TMath::Sqrt((pi0basedErrLow*pi0basedErrLow)+(etaErrLow*etaErrLow)+sqrsumErrs));
     hScalingErrors->SetBinContent(iBin, TMath::Sqrt(sqrsumErrs)/hScalingErrors->GetBinWidth(iBin));
 
     hOverallBinScaledErrsUp->SetBinContent(iBin,hOverallSystErrUp->GetBinContent(iBin)/hOverallBinScaledErrsUp->GetBinWidth(iBin));
     hOverallBinScaledErrsLow->SetBinContent(iBin,hOverallSystErrLow->GetBinContent(iBin)/hOverallBinScaledErrsLow->GetBinWidth(iBin));           
   }
    
-
-   // /hOverallSystErrUp->GetBinWidth(iBin))
   
   TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
   cPiErrors->cd();
   cPiErrors->SetLogx();
   cPiErrors->SetLogy();
-  hpiErrors[0]->Draw();
-  hpiErrors[1]->SetMarkerColor(kBlack);
-  hpiErrors[1]->SetLineColor(kBlack);
-  hpiErrors[1]->Draw("SAME");
-  hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
-  hOverallBinScaledErrsUp->SetLineColor(kBlue);
-  hOverallBinScaledErrsUp->Draw("SAME");
+  hBaseErrors[0][0]->Draw();
+  //hBaseErrors[0][1]->SetMarkerColor(kBlack);
+  //hBaseErrors[0][1]->SetLineColor(kBlack);
+  //hBaseErrors[0][1]->Draw("SAME");
+  if(fEtaSyst){
+    hBaseErrors[1][0]->Draw("SAME");
+    hBaseErrors[1][0]->SetMarkerColor(kBlack);
+    hBaseErrors[1][0]->SetLineColor(kBlack);
+  //hBaseErrors[1][1]->SetMarkerColor(13);
+  //hBaseErrors[1][1]->SetLineColor(13);
+  //hBaseErrors[1][1]->Draw("SAME");
+  }
+  //hOverallBinScaledErrsUp->SetMarkerColor(kBlue);
+  //hOverallBinScaledErrsUp->SetLineColor(kBlue);
+  //hOverallBinScaledErrsUp->Draw("SAME");
   hOverallBinScaledErrsLow->SetMarkerColor(kGreen);
   hOverallBinScaledErrsLow->SetLineColor(kGreen);
   hOverallBinScaledErrsLow->Draw("SAME");
-  hScalingErrors->SetLineColor(11);
+  hScalingErrors->SetLineColor(kBlue);
   hScalingErrors->Draw("SAME");
 
   TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
-  lPiErr->AddEntry(hpiErrors[0],"Lower error from pion error");
-  lPiErr->AddEntry(hpiErrors[1],"Upper error from pion error");
+  lPiErr->AddEntry(hBaseErrors[0][0],"Lower error from pion error");
+  //lPiErr->AddEntry(hBaseErrors[0][1],"Upper error from pion error");
+  if(fEtaSyst){
+  lPiErr->AddEntry(hBaseErrors[1][0],"Lower error from eta error");
+  //lPiErr->AddEntry(hBaseErrors[1][1],"Upper error from eta error");
+  }
   lPiErr->AddEntry(hScalingErrors, "scaling error");
   lPiErr->AddEntry(hOverallBinScaledErrsLow, "overall lower systematics");
-  lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
+  //lPiErr->AddEntry(hOverallBinScaledErrsUp, "overall upper systematics");
   lPiErr->Draw("SAME");
 
   //Normalize errors
@@ -3255,7 +3753,10 @@ void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
 
 
   AliCFDataGrid *bgYieldGrid;
-  bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0]->Clone();
+  if(fEtaSyst){
+    bgLevelGrid[0][0]->Add(bgLevelGrid[1][0]);//Addition of the eta background best estimate to the rest. Needed to be separated for error treatment - now overall background necessary! If no separate eta systematics exist, the corresponding grid has already been added before.
+  }
+  bgYieldGrid = (AliCFDataGrid*)bgLevelGrid[0][0]->Clone();
 
   TH1D *hBgYield = (TH1D*)bgYieldGrid->Project(0);
   TH1D* hRelErrUp = (TH1D*)hOverallSystErrUp->Clone();
@@ -3303,3 +3804,96 @@ void AliHFEspectrum::CalculateNonHFEsyst(Int_t centrality){
   
 }
 
+//____________________________________________________________
+void AliHFEspectrum::UnfoldBG(AliCFDataGrid* const bgsubpectrum){
+
+  //
+  // Unfold backgrounds to check its sanity
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
+  //AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return;
+  }
+
+  if(!fCorrelation){
+    AliError("No Correlation map available");
+    return;
+  }
+
+  // Data 
+  AliCFDataGrid *dataGrid = 0x0;
+  dataGrid = bgsubpectrum;
+
+  // Guessed
+  AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
+  THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(fStepMC+2,fStepTrue);
+
+  // Unfold background spectra
+  Int_t nDim=1;
+  if(fBeamType==0)nDim = 1;
+  if(fBeamType==1)nDim = 2;
+  AliCFUnfolding unfolding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse);
+  if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
+  unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
+  if(fSetSmoothing) unfolding.UseSmoothing();
+  unfolding.Unfold();
+
+  // Results
+  THnSparse* result = unfolding.GetUnfolded();
+  TCanvas *ctest = new TCanvas("yvonnetest","yvonnetest",1000,600);
+  if(fBeamType==1)
+  {
+      ctest->Divide(2,2);
+      ctest->cd(1);
+      result->GetAxis(0)->SetRange(1,1);
+      TH1D* htest1=(TH1D*)result->Projection(0);
+      htest1->Draw();
+      ctest->cd(2);
+      result->GetAxis(0)->SetRange(1,1);
+      TH1D* htest2=(TH1D*)result->Projection(1);
+      htest2->Draw();
+      ctest->cd(3);
+      result->GetAxis(0)->SetRange(6,6);
+      TH1D* htest3=(TH1D*)result->Projection(0);
+      htest3->Draw();
+      ctest->cd(4);
+      result->GetAxis(0)->SetRange(6,6);
+      TH1D* htest4=(TH1D*)result->Projection(1);
+      htest4->Draw();
+
+  }
+
+
+
+
+
+  TGraphErrors* unfoldedbgspectrumD = Normalize(result);
+  if(!unfoldedbgspectrumD) {
+    AliError("Unfolded background spectrum doesn't exist");
+  }
+  else{
+    TFile *file = TFile::Open("unfoldedbgspectrum.root","recreate");
+    if(fBeamType==0)unfoldedbgspectrumD->Write("unfoldedbgspectrum");
+
+    if(fBeamType==1)
+    {
+        Int_t centr=1;
+       result->GetAxis(0)->SetRange(centr,centr);
+       unfoldedbgspectrumD = Normalize(result,centr-1);
+       unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr0_20");
+        centr=6;
+       result->GetAxis(0)->SetRange(centr,centr);
+       unfoldedbgspectrumD = Normalize(result,centr-1);
+        unfoldedbgspectrumD->Write("unfoldedbgspectrum_centr40_80");
+    }
+
+    file->Close();
+  }
+}