Add Classes (Thanks Constantin)
authorrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Feb 2013 08:20:30 +0000 (08:20 +0000)
committerrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Feb 2013 08:20:30 +0000 (08:20 +0000)
PWGHF/hfe/AliHFEBeautySpectrum.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEBeautySpectrum.h [new file with mode: 0644]
PWGHF/hfe/AliHFEBeautySpectrumQA.cxx [new file with mode: 0644]
PWGHF/hfe/AliHFEBeautySpectrumQA.h [new file with mode: 0644]

diff --git a/PWGHF/hfe/AliHFEBeautySpectrum.cxx b/PWGHF/hfe/AliHFEBeautySpectrum.cxx
new file mode 100644 (file)
index 0000000..f105eb6
--- /dev/null
@@ -0,0 +1,2239 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: 
+//            Raphaelle Bailhache <R.Bailhache@gsi.de>
+//            Markus Fasel <M.Fasel@gsi.de>
+//
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+#include <TH2D.h>
+#include <TF1.h>
+
+#include "AliPID.h"
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFEBeautySpectrum.h"
+#include "AliHFEBeautySpectrumQA.h"
+#include "AliHFEcuts.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEtools.h"
+
+ClassImp(AliHFEBeautySpectrum)
+
+//____________________________________________________________
+AliHFEBeautySpectrum::AliHFEBeautySpectrum(const char *name):
+AliHFECorrectSpectrumBase(name),
+  fQA(NULL),
+  fTemporaryObjects(NULL),
+  fBackground(NULL),
+  fEfficiencyTOFPIDD(NULL),
+  fEfficiencyesdTOFPIDD(NULL),
+  fEfficiencyIPCharmD(NULL),
+  fEfficiencyIPBeautyD(NULL),
+  fEfficiencyIPBeautyesdD(NULL),
+  fEfficiencyIPConversionD(NULL),
+  fEfficiencyIPNonhfeD(NULL), 
+  fNonHFEbg(NULL),
+  fWeightCharm(NULL),
+  fInclusiveSpectrum(kFALSE),
+  fDumpToFile(kFALSE),
+  fUnSetCorrelatedErrors(kTRUE),
+  fIPanaHadronBgSubtract(kFALSE),
+  fIPanaCharmBgSubtract(kFALSE),
+  fIPanaNonHFEBgSubtract(kFALSE),
+  fIPParameterizedEff(kFALSE),
+  fNonHFEsyst(0),
+  fBeauty2ndMethod(kFALSE),
+  fIPEffCombinedSamples(kTRUE),
+  fNMCEvents(0),
+  fNMCbgEvents(0),
+  fNCentralityBinAtTheEnd(0),
+  fFillMoreCorrelationMatrix(kFALSE),
+  fHadronEffbyIPcut(NULL),
+  fEfficiencyCharmSigD(NULL),
+  fEfficiencyBeautySigD(NULL),
+  fEfficiencyBeautySigesdD(NULL),
+  fConversionEff(NULL),
+  fNonHFEEff(NULL),
+  fCharmEff(NULL),
+  fBeautyEff(NULL),
+  fConversionEffbgc(NULL),
+  fNonHFEEffbgc(NULL),      
+  fBSpectrum2ndMethod(NULL),
+  fkBeauty2ndMethodfilename(""),
+  fBeamType(0),
+  fEtaSyst(kTRUE),
+  fDebugLevel(0),
+  fWriteToFile(kFALSE),
+  fUnfoldBG(kFALSE)
+{
+  //
+  // Default constructor
+  //
+
+  fQA = new AliHFEBeautySpectrumQA("AliHFEBeautySpectrumQA");
+
+  for(Int_t k = 0; k < 20; k++){
+    fLowBoundaryCentralityBinAtTheEnd[k] = 0;
+    fHighBoundaryCentralityBinAtTheEnd[k] = 0;
+  }
+  fNMCbgEvents = 0;  
+  fEfficiencyTOFPIDD = 0;
+  fEfficiencyesdTOFPIDD = 0;
+  fEfficiencyIPCharmD = 0;     
+  fEfficiencyIPBeautyD = 0;    
+  fEfficiencyIPBeautyesdD = 0;
+  fEfficiencyIPConversionD = 0;
+  fEfficiencyIPNonhfeD = 0;   
+  
+  fConversionEff = 0;
+  fNonHFEEff = 0;
+  fCharmEff = 0;
+  fBeautyEff = 0;
+  fEfficiencyCharmSigD = 0;
+  fEfficiencyBeautySigD = 0;
+  fEfficiencyBeautySigesdD = 0;
+  
+  
+  memset(fConvSourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
+  memset(fNonHFESourceContainer, 0, sizeof(AliCFContainer*) * kElecBgSources * kBgLevels);
+}
+//____________________________________________________________
+AliHFEBeautySpectrum::AliHFEBeautySpectrum(const AliHFEBeautySpectrum &ref):
+  AliHFECorrectSpectrumBase(ref),
+  fQA(ref.fQA),
+  fTemporaryObjects(ref.fTemporaryObjects),
+  fBackground(ref.fBackground),
+  fEfficiencyTOFPIDD(ref.fEfficiencyTOFPIDD),
+  fEfficiencyesdTOFPIDD(ref.fEfficiencyesdTOFPIDD),
+  fEfficiencyIPCharmD(ref.fEfficiencyIPCharmD),
+  fEfficiencyIPBeautyD(ref.fEfficiencyIPBeautyD),
+  fEfficiencyIPBeautyesdD(ref.fEfficiencyIPBeautyesdD),
+  fEfficiencyIPConversionD(ref.fEfficiencyIPConversionD),
+  fEfficiencyIPNonhfeD(ref.fEfficiencyIPNonhfeD), 
+  fNonHFEbg(ref.fNonHFEbg),
+  fWeightCharm(ref.fWeightCharm),
+  fInclusiveSpectrum(ref.fInclusiveSpectrum),
+  fDumpToFile(ref.fDumpToFile),
+  fUnSetCorrelatedErrors(ref.fUnSetCorrelatedErrors),
+  fIPanaHadronBgSubtract(ref.fIPanaHadronBgSubtract),
+  fIPanaCharmBgSubtract(ref.fIPanaCharmBgSubtract),
+  fIPanaNonHFEBgSubtract(ref.fIPanaNonHFEBgSubtract),
+  fIPParameterizedEff(ref.fIPParameterizedEff),
+  fNonHFEsyst(ref.fNonHFEsyst),
+  fBeauty2ndMethod(ref.fBeauty2ndMethod),
+  fIPEffCombinedSamples(ref.fIPEffCombinedSamples),
+  fNMCEvents(ref.fNMCEvents),
+  fNMCbgEvents(ref.fNMCbgEvents),
+  fNCentralityBinAtTheEnd(ref.fNCentralityBinAtTheEnd),
+  fFillMoreCorrelationMatrix(ref.fFillMoreCorrelationMatrix),
+  fHadronEffbyIPcut(ref.fHadronEffbyIPcut),
+  fEfficiencyCharmSigD(ref.fEfficiencyCharmSigD),
+  fEfficiencyBeautySigD(ref.fEfficiencyBeautySigD),
+  fEfficiencyBeautySigesdD(ref.fEfficiencyBeautySigesdD),
+  fConversionEff(ref.fConversionEff),
+  fNonHFEEff(ref.fNonHFEEff),
+  fCharmEff(ref.fCharmEff),
+  fBeautyEff(ref.fBeautyEff),
+  fConversionEffbgc(ref.fConversionEffbgc),
+  fNonHFEEffbgc(ref.fNonHFEEffbgc),      
+  fBSpectrum2ndMethod(ref.fBSpectrum2ndMethod),
+  fkBeauty2ndMethodfilename(ref.fkBeauty2ndMethodfilename),
+  fBeamType(ref.fBeamType),
+  fEtaSyst(ref.fEtaSyst),
+  fDebugLevel(ref.fDebugLevel),
+  fWriteToFile(ref.fWriteToFile),
+  fUnfoldBG(ref.fUnfoldBG)
+{
+  //
+  // Copy constructor
+  //
+
+  ref.Copy(*this);  
+}
+//____________________________________________________________
+AliHFEBeautySpectrum &AliHFEBeautySpectrum::operator=(const AliHFEBeautySpectrum &ref){
+ //
+  // Assignment operator
+  //
+  if(this == &ref) 
+    ref.Copy(*this);
+  return *this;
+}
+//____________________________________________________________
+void AliHFEBeautySpectrum::Copy(TObject &o) const {
+  // 
+  // Copy into object o
+  //
+  AliHFEBeautySpectrum &target = dynamic_cast<AliHFEBeautySpectrum &>(o);
+  target.fQA = fQA;
+
+  target.fQA = fQA;
+  target.fTemporaryObjects = fTemporaryObjects;
+  target.fBackground = fBackground;
+  target.fWeightCharm = fWeightCharm;
+  target.fDumpToFile = fDumpToFile;
+  target.fUnSetCorrelatedErrors = fUnSetCorrelatedErrors;
+  target.fIPanaHadronBgSubtract = fIPanaHadronBgSubtract;
+  target.fIPanaCharmBgSubtract = fIPanaCharmBgSubtract;
+  target.fIPanaNonHFEBgSubtract = fIPanaNonHFEBgSubtract;
+  target.fIPParameterizedEff = fIPParameterizedEff;
+  target.fNonHFEsyst = fNonHFEsyst;
+  target.fBeauty2ndMethod = fBeauty2ndMethod;
+  target.fIPEffCombinedSamples = fIPEffCombinedSamples;
+  target.fNMCEvents = fNMCEvents;
+  target.fNCentralityBinAtTheEnd = fNCentralityBinAtTheEnd;
+  target.fFillMoreCorrelationMatrix = fFillMoreCorrelationMatrix;
+  target.fHadronEffbyIPcut = fHadronEffbyIPcut;
+  target.fConversionEffbgc = fConversionEffbgc;
+  target.fNonHFEEffbgc = fNonHFEEffbgc;      
+  target.fBSpectrum2ndMethod = fBSpectrum2ndMethod;
+  target.fkBeauty2ndMethodfilename = fkBeauty2ndMethodfilename;
+  target.fBeamType = fBeamType;
+  target.fEtaSyst = fEtaSyst;
+  target.fDebugLevel = fDebugLevel;
+  target.fWriteToFile = fWriteToFile;
+  target.fUnfoldBG = fUnfoldBG;
+}
+//____________________________________________________________
+AliHFEBeautySpectrum::~AliHFEBeautySpectrum(){
+  //
+  // Destructor
+  //
+  if(fTemporaryObjects){
+    fTemporaryObjects->Clear();
+    delete fTemporaryObjects;
+  }
+  if(fQA) delete fQA;
+}
+//____________________________________________________________
+Bool_t AliHFEBeautySpectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *bghfecontainer, const AliHFEcontainer */*v0hfecontainer*/,AliCFContainer */*photoniccontainerD*/){
+  //
+  // Init what we need for the correction:
+  //
+  // Raw spectrum, hadron contamination
+  // MC efficiency maps, correlation matrix
+  //
+  // This for a given dimension.
+  // If no inclusive spectrum, then take only efficiency map for beauty electron
+  // and the appropriate correlation matrix
+  //
+
+  
+  if(fBeauty2ndMethod) CallInputFileForBeauty2ndMethod();
+
+  Int_t kNdim = 3;
+  Int_t dims[kNdim];
+  
+  switch(fNbDimensions){
+  case 1:   dims[0] = 0;
+    break;
+  case 2:   for(Int_t i = 0; i < 2; i++) dims[i] = i;
+    break;
+  case 3:   for(Int_t i = 0; i < 3; i++) dims[i] = i;
+    break;
+  default:
+    AliError("Container with this number of dimensions not foreseen (yet)");
+    return kFALSE;
+  };
+  
+  //
+  // Data container: raw spectrum + hadron contamination  
+  //
+  AliCFContainer *datacontainer = datahfecontainer->MakeMergedCFContainer("sumreco","sumreco","recTrackContReco:recTrackContDEReco");
+  AliCFContainer *contaminationcontainer = datahfecontainer->GetCFContainer("hadronicBackground");
+  if((!datacontainer) || (!contaminationcontainer)) return kFALSE;
+  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,AliHFECorrectSpectrumBase::kDataContainer);
+  SetContainer(contaminationcontainerD,AliHFECorrectSpectrumBase::kBackgroundData);
+
+  // QA : see with MJ if it is necessary for Beauty
+  Int_t dimqa = datacontainer->GetNVar();
+  Int_t dimsqa[dimqa];
+  for(Int_t i = 0; i < dimqa; i++) dimsqa[i] = i;
+  AliCFContainer *datacontainerDQA = GetSlicedContainer(datacontainer, dimqa, dimsqa, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  fQA->AddResultAt(datacontainerDQA,AliHFEBeautySpectrumQA::kDataProjection);
+
+  //
+  // 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;
+  AliCFContainer *convweightedContainer = 0x0;
+  AliCFContainer *nonHFEweightedContainerSig = 0x0;//mjmj
+  AliCFContainer *convweightedContainerSig = 0x0;//mjmj
+  AliCFContainer *nonHFEtempContainer = 0x0;//temporary container to be sliced for the fnonHFESourceContainers
+  AliCFContainer *convtempContainer = 0x0;//temporary container to be sliced for the fConvSourceContainers
+   
+  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");
+  
+  if(fNonHFEsyst){   
+    const Char_t *sourceName[kElecBgSources]={"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
+    const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
+    for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
+      for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
+        if(!(nonHFEtempContainer = bghfecontainer->GetCFContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]))))continue;
+        convtempContainer = bghfecontainer->GetCFContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));     
+        if(!(fConvSourceContainer[iSource][iLevel] = GetSlicedContainer(convtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
+        if(!(fNonHFESourceContainer[iSource][iLevel] = GetSlicedContainer(nonHFEtempContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh)))continue;
+        //           if((!fConvSourceContainer[iSource][iLevel][icentr])||(!fNonHFESourceContainer[iSource][iLevel])) return kFALSE;
+        if(fBeamType == 1)break;//no systematics yet for PbPb non-HFE backgrounds
+      }
+    }
+  }
+  // else{      
+  nonHFEweightedContainer = bghfecontainer->GetCFContainer("mesonElecs");
+  convweightedContainer = bghfecontainer->GetCFContainer("conversionElecs");
+  if((!convweightedContainer)||(!nonHFEweightedContainer)) return kFALSE; 
+  nonHFEweightedContainerSig = mchfecontainer->GetCFContainer("mesonElecs");//mjmj
+  convweightedContainerSig = mchfecontainer->GetCFContainer("conversionElecs");//mjmj
+  if((!convweightedContainerSig)||(!nonHFEweightedContainerSig)) return kFALSE;
+  //}
+  
+  if((!mccontaineresd) || (!mccontainermc)) return kFALSE;  
+  
+  Int_t source = 1; //beauty
+  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,AliHFECorrectSpectrumBase::kMCContainerMC);
+  SetContainer(mccontaineresdD,AliHFECorrectSpectrumBase::kMCContainerESD);
+
+  // set charm, nonHFE container to estimate BG
+  source = 0; //charm
+  mccontainermcD = GetSlicedContainer(mccontainermc, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh); //mjmjmj
+  //mccontainermcD = GetSlicedContainer(mccontainermcbg, fNbDimensions, dims, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  SetContainer(mccontainermcD,AliHFECorrectSpectrumBase::kMCContainerCharmMC);
+
+  //if(!fNonHFEsyst){
+  AliCFContainer *nonHFEweightedContainerD = GetSlicedContainer(nonHFEweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  SetContainer(nonHFEweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESD);
+  AliCFContainer *convweightedContainerD = GetSlicedContainer(convweightedContainer, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  SetContainer(convweightedContainerD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESD);
+  
+  //mjmj
+  AliCFContainer *nonHFEweightedContainerSigD = GetSlicedContainer(nonHFEweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  SetContainer(nonHFEweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig);
+  if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig)){printf("merged non-HFE container not found!\n");return kFALSE;};
+  AliCFContainer *convweightedContainerSigD = GetSlicedContainer(convweightedContainerSig, fNbDimensions, dims, -1, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+  SetContainer(convweightedContainerSigD,AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig); 
+  if(!GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig)){printf(" merged conversion container not found!\n");return kFALSE;}; 
+  //}
+  
+  SetParameterizedEff(mccontainermc, mccontainermcbg, mccontaineresd, mccontaineresdbg, dims);
+  
+  //
+  // MC container: correlation matrix
+  //
+  THnSparseF *mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
+  //mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
+  if(!mccorrelation) return kFALSE;
+  Int_t testCentralityLow = fTestCentralityLow;
+  Int_t testCentralityHigh = fTestCentralityHigh;
+  if(fFillMoreCorrelationMatrix) {
+    testCentralityLow = fTestCentralityLow-1;
+    testCentralityHigh = fTestCentralityHigh+1;
+  }
+  THnSparseF *mccorrelationD = GetSlicedCorrelation(mccorrelation, fNbDimensions, dims, fChargeChoosen, testCentralityLow,testCentralityHigh);
+  if(!mccorrelationD) {
+    printf("No correlation\n");
+    return kFALSE;
+  }
+  SetCorrelation(mccorrelationD);
+
+  // QA
+  fQA->AddResultAt(mccorrelation,AliHFEBeautySpectrumQA::kCMProjection); 
+  //
+  fQA->DrawProjections();
+
+
+  return kTRUE;
+}
+
+
+//____________________________________________________________
+void AliHFEBeautySpectrum::CallInputFileForBeauty2ndMethod(){
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+    TFile *inputfile2ndmethod=TFile::Open(fkBeauty2ndMethodfilename);
+    fBSpectrum2ndMethod = new TH1D(*static_cast<TH1D*>(inputfile2ndmethod->Get("BSpectrum")));
+}
+
+//____________________________________________________________
+Bool_t AliHFEBeautySpectrum::Correct(Bool_t subtractcontamination, Bool_t /*subtractphotonic*/){
+  //
+  // Correct the spectrum for efficiency and unfolding for beauty analysis
+  // with both method and compare
+  //
+  
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(1111);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadLeftMargin(0.13);
+  gStyle->SetPadRightMargin(0.13);
+
+  printf("Steps are: stepdata %d, stepMC %d, steptrue %d\n",fStepData,fStepMC,fStepTrue);
+
+  ///////////////////////////
+  // Check initialization
+  ///////////////////////////
+
+  if((!GetContainer(kDataContainer)) || (!GetContainer(kMCContainerMC)) || (!GetContainer(kMCContainerESD))){
+    AliInfo("You have to init before");
+    return kFALSE;
+  }
+  
+  if((fStepTrue <= 0) && (fStepMC <= 0) && (fStepData <= 0)) {
+    AliInfo("You have to set the steps before: SetMCTruthStep, SetMCEffStep, SetStepToCorrect");
+    return kFALSE;
+  }
+  SetStepGuessedUnfolding(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
+    
+  AliCFDataGrid *dataGridAfterFirstSteps = 0x0;
+  //////////////////////////////////
+  // Subtract hadron background
+  /////////////////////////////////
+  AliCFDataGrid *dataspectrumaftersubstraction = 0x0;
+  AliCFDataGrid *unnormalizedRawSpectrum = 0x0;
+  TGraphErrors *gNormalizedRawSpectrum = 0x0;
+  if(subtractcontamination) {
+      if(!fBeauty2ndMethod) dataspectrumaftersubstraction = SubtractBackground(kTRUE);
+      else dataspectrumaftersubstraction = GetRawBspectra2ndMethod();
+      unnormalizedRawSpectrum = (AliCFDataGrid*)dataspectrumaftersubstraction->Clone();
+      dataGridAfterFirstSteps = dataspectrumaftersubstraction;
+      gNormalizedRawSpectrum = Normalize(unnormalizedRawSpectrum);
+  }
+
+  /////////////////////////////////////////////////////////////////////////////////////////
+  // Correct for IP efficiency for beauty electrons after subtracting all the backgrounds
+  /////////////////////////////////////////////////////////////////////////////////////////
+
+  AliCFDataGrid *dataspectrumafterefficiencyparametrizedcorrection = 0x0;
+  
+  if(fEfficiencyFunction){
+    dataspectrumafterefficiencyparametrizedcorrection = CorrectParametrizedEfficiency(dataGridAfterFirstSteps);
+    dataGridAfterFirstSteps = dataspectrumafterefficiencyparametrizedcorrection;
+  }
+  ///////////////
+  // Unfold
+  //////////////
+  THnSparse *correctedspectrum = Unfold(dataGridAfterFirstSteps);
+  if(!correctedspectrum){
+    AliError("No corrected spectrum\n");
+    return kFALSE;
+  }
+  
+  /////////////////////
+  // Simply correct
+  ////////////////////
+
+  AliCFDataGrid *alltogetherCorrection = CorrectForEfficiency(dataGridAfterFirstSteps);
+  
+
+  //////////
+  // Plot
+  //////////
+
+  // QA final results
+  TGraphErrors* correctedspectrumD = Normalize(correctedspectrum);
+  TGraphErrors* alltogetherspectrumD = Normalize(alltogetherCorrection);
+  fQA->AddResultAt(correctedspectrumD,AliHFEBeautySpectrumQA::kFinalResultUnfolded);
+  fQA->AddResultAt(alltogetherspectrumD,AliHFEBeautySpectrumQA::kFinalResultDirectEfficiency);
+  fQA->AddResultAt(correctedspectrum,AliHFEBeautySpectrumQA::kFinalResultUnfSparse);
+  fQA->AddResultAt(alltogetherCorrection,AliHFEBeautySpectrumQA::kFinalResultDirectEffSparse);
+  fQA->DrawResult();
+
+  if(fBeamType == 0){
+    if(fNonHFEsyst){
+      CalculateNonHFEsyst();
+    }
+  }
+
+  // Dump to file if needed
+  
+  if(fDumpToFile) {
+    TFile *out;
+    out = new TFile("finalSpectrum.root","update");
+    out->cd();
+    // to do centrality dependent
+    correctedspectrum->SetName("UnfoldingCorrectedNotNormalizedSpectrum");
+    correctedspectrum->Write();
+    alltogetherCorrection->SetName("AlltogetherCorrectedNotNormalizedSpectrum");
+    alltogetherCorrection->Write();
+    //
+    if(unnormalizedRawSpectrum) {
+      unnormalizedRawSpectrum->SetName("beautyAfterIP");
+      unnormalizedRawSpectrum->Write();
+    }
+    
+    if(gNormalizedRawSpectrum){
+      gNormalizedRawSpectrum->SetName("normalizedBeautyAfterIP");
+      gNormalizedRawSpectrum->Write();
+    }
+        
+    fEfficiencyCharmSigD->SetTitle("IPEfficiencyForCharmSig");
+    fEfficiencyCharmSigD->SetName("IPEfficiencyForCharmSig");
+    fEfficiencyCharmSigD->Write();
+    fEfficiencyBeautySigD->SetTitle("IPEfficiencyForBeautySig");
+    fEfficiencyBeautySigD->SetName("IPEfficiencyForBeautySig");
+    fEfficiencyBeautySigD->Write();
+    fCharmEff->SetTitle("IPEfficiencyForCharm");
+    fCharmEff->SetName("IPEfficiencyForCharm");
+    fCharmEff->Write();
+    fBeautyEff->SetTitle("IPEfficiencyForBeauty");
+    fBeautyEff->SetName("IPEfficiencyForBeauty");
+    fBeautyEff->Write();
+    fConversionEff->SetTitle("IPEfficiencyForConversion");
+    fConversionEff->SetName("IPEfficiencyForConversion");
+    fConversionEff->Write();
+    fNonHFEEff->SetTitle("IPEfficiencyForNonHFE");
+    fNonHFEEff->SetName("IPEfficiencyForNonHFE");
+    fNonHFEEff->Write();
+    
+   
+    out->Close(); 
+    delete out;
+   }
+  return kTRUE;
+}
+//____________________________________________________________
+AliCFDataGrid* AliHFEBeautySpectrum::SubtractBackground(Bool_t setBackground){
+  //
+  // Apply background subtraction for IP analysis
+  //
+  
+  Int_t nbins=1;
+  printf("subtracting background\n");
+  // Raw spectrum
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  if(!dataContainer){
+    AliError("Data Container not available");
+    return NULL;
+  }
+  printf("Step data: %d\n",fStepData);
+  AliCFDataGrid *spectrumSubtracted = new AliCFDataGrid("spectrumSubtracted", "Data Grid for spectrum after Background subtraction", *dataContainer,fStepData);
+  
+  AliCFDataGrid *dataspectrumbeforesubstraction = (AliCFDataGrid *) ((AliCFDataGrid *)GetSpectrum(GetContainer(kDataContainer),fStepData))->Clone();
+  dataspectrumbeforesubstraction->SetName("dataspectrumbeforesubstraction"); 
+  
+  // Background Estimate
+  AliCFContainer *backgroundContainer = GetContainer(kBackgroundData);
+  if(!backgroundContainer){
+    AliError("MC background container not found");
+    return NULL;
+  }
+    
+  Int_t stepbackground = 1;
+  AliCFDataGrid *backgroundGrid = new AliCFDataGrid("ContaminationGrid","ContaminationGrid",*backgroundContainer,stepbackground);
+  
+  const Int_t bgPlots = 8;
+  TH1D *incElec = 0x0;
+  TH1D *hadron = 0x0;
+  TH1D *charm = 0x0;
+  TH1D *nonHFE[bgPlots];
+  TH1D *subtracted = 0x0;
+  
+  THnSparseF* sparseIncElec = (THnSparseF *) dataspectrumbeforesubstraction->GetGrid();  
+  incElec = (TH1D *) sparseIncElec->Projection(0);
+  CorrectFromTheWidth(incElec);
+
+  TH1D* htemp;
+  Int_t* bins=new Int_t[2];
+
+  if(fIPanaHadronBgSubtract){
+    // Hadron background
+    printf("Hadron background for IP analysis subtracted!\n");
+    htemp  = (TH1D *) fHadronEffbyIPcut->Projection(0);
+    bins[0]=htemp->GetNbinsX();
+    
+    AliCFDataGrid *hbgContainer = new AliCFDataGrid("hbgContainer","hadron bg after IP cut",nbins,bins);
+    hbgContainer->SetGrid(fHadronEffbyIPcut);
+    backgroundGrid->Multiply(hbgContainer,1);
+    // subtract hadron contamination
+    spectrumSubtracted->Add(backgroundGrid,-1.0);
+    // draw raw hadron bg spectra     
+    THnSparseF* sparseHbg = (THnSparseF *) hbgContainer->GetGrid();
+    hadron = (TH1D *) sparseHbg->Projection(0);
+    CorrectFromTheWidth(hadron);
+  }
+
+  if(fIPanaCharmBgSubtract){
+    // Charm background
+    printf("Charm background for IP analysis subtracted!\n");
+    AliCFDataGrid *charmbgContainer = (AliCFDataGrid *) GetCharmBackground();
+    spectrumSubtracted->Add(charmbgContainer,-1.0);
+    THnSparseF *sparseCharmElec = (THnSparseF *) charmbgContainer->GetGrid();  
+    charm = (TH1D *) sparseCharmElec->Projection(0);
+    CorrectFromTheWidth(charm); 
+  }
+
+  const Char_t *sourceName[bgPlots]={"Conversion","Dalitz","K0s secondary pion Dalitz","K0s secondary pion conversions","Other secondary pion Dalitz","Other secondary pion Dalitz conversions","Other Dalitz", "Other conversions"};
+  const Int_t style[2] = {21,22};
+  const Int_t color[3] = {7,12,8};
+  if(fIPanaNonHFEBgSubtract){ 
+    Int_t iPlot = 0;
+    for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
+      if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
+      for(Int_t decay = 0; decay < 2; decay++){
+        AliCFDataGrid *nonHFEbgContainer = (AliCFDataGrid *) GetNonHFEBackground(decay,iSource);//conversions/Dalitz decays from different sources
+        if(iSource == 0)
+          spectrumSubtracted->Add(nonHFEbgContainer,-1.0);
+        THnSparseF* sparseNonHFEelecs = (THnSparseF *) nonHFEbgContainer->GetGrid();
+        nonHFE[iPlot] = (TH1D *) sparseNonHFEelecs->Projection(0);
+        CorrectFromTheWidth(nonHFE[iPlot]);
+        iPlot++;
+      }
+      if(!(fNonHFEsyst && (fBeamType == 1)))break;
+    } 
+  }
+
+  TLegend *lRaw = new TLegend(0.55,0.55,0.85,0.85);
+  TCanvas *cRaw = new TCanvas("cRaw","cRaw",1000,800);     
+  cRaw->cd();
+  gPad->SetLogx();
+  gPad->SetLogy();
+  incElec->GetXaxis()->SetRangeUser(0.4,8.);
+  incElec->Draw("p");
+  incElec->SetMarkerColor(1);
+  incElec->SetMarkerStyle(20);
+  lRaw->AddEntry(incElec,"inclusive electron spectrum");
+  
+  if(fIPanaCharmBgSubtract){
+    charm->Draw("samep");
+    charm->SetMarkerColor(3);
+    charm->SetMarkerStyle(20);
+    charm->GetXaxis()->SetRangeUser(0.4,8.);
+    lRaw->AddEntry(charm,"charm elecs");
+  }
+  
+  if(fIPanaNonHFEBgSubtract){
+    Int_t iPlot = 0;
+    for(Int_t iSource = 0; iSource < kElecBgSources; iSource++){
+      if((iSource>0)&&(iSource<6))continue;//standard sources are summed up in case of iSource == 0; not treated separately in plots
+      for(Int_t decay = 0; decay < 2; decay++){
+        nonHFE[iPlot]->Draw("samep");
+        if(iSource == 0){
+          nonHFE[iPlot]->SetMarkerStyle(20);
+          if(decay == 0)
+            nonHFE[iPlot]->SetMarkerColor(4);
+          else
+            nonHFE[iPlot]->SetMarkerColor(6);
+        }
+        else{     
+          nonHFE[iPlot]->SetMarkerColor(color[iSource-6]);      
+          nonHFE[iPlot]->SetMarkerStyle(style[decay]);
+        }
+        nonHFE[iPlot]->GetXaxis()->SetRangeUser(0.4,8.);
+        lRaw->AddEntry(nonHFE[iPlot],sourceName[iPlot]);
+        iPlot++;
+      }
+      if(!(fNonHFEsyst && (fBeamType == 1)))break;//only draw standard sources
+    }
+  }
+
+  THnSparseF* sparseSubtracted = (THnSparseF *) spectrumSubtracted->GetGrid();
+  subtracted = (TH1D *) sparseSubtracted->Projection(0);
+  CorrectFromTheWidth(subtracted);
+  subtracted->Draw("samep");
+  subtracted->SetMarkerStyle(24);
+  lRaw->AddEntry(subtracted,"subtracted electron spectrum");
+  lRaw->Draw("SAME");
+  delete[] bins; 
+
+  TH1D *measuredTH1background = (TH1D *) backgroundGrid->Project(0);
+  CorrectFromTheWidth(measuredTH1background);
+
+  if(setBackground){
+    if(fBackground) delete fBackground;
+    fBackground = backgroundGrid;
+  } else delete backgroundGrid;
+  
+  fQA->AddResultAt(subtracted,AliHFEBeautySpectrumQA::kAfterSC);
+  fQA->AddResultAt(incElec,AliHFEBeautySpectrumQA::kBeforeSC);
+  fQA->AddResultAt(measuredTH1background, AliHFEBeautySpectrumQA::kMeasBG);
+  fQA->DrawSubtractContamination();
+  
+  return spectrumSubtracted;
+}
+//____________________________________________________________________________________________
+AliCFDataGrid* AliHFEBeautySpectrum::GetNonHFEBackground(Int_t decay, Int_t source){
+  //
+  // calculate non-HF electron background
+  // Arguments: decay = 0 for conversions, decay = 1 for meson decays to electrons
+  // source: gives mother either of the electron (for meson->electron) or of the gamma (conversions);
+  // source numbers as given by array in AliAnalysisTaskHFE: 
+  // {"Pion","Eta","Omega","Phi","EtaPrime","Rho","K0sSec","OtherSecPi0","Else"};
+  //
+  
+  Double_t evtnorm[1] = {1.};
+  if(fNMCbgEvents) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
+  
+  Int_t stepbackground = 1;//for non-standard background, step after IP cut is taken directly
+  if(source == 0) stepbackground = 3;//for standard sources, take step before PID -> correction of PID x IP efficiencies from enhanced samples
+  AliCFContainer *backgroundContainer = 0x0;
+  if(fNonHFEsyst){
+    if(decay == 0){//conversions
+      backgroundContainer = (AliCFContainer*)fConvSourceContainer[source][0]->Clone();
+      if(source == 0)//standard conversion containers
+        for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){//add all standard sources
+          backgroundContainer->Add(fConvSourceContainer[iSource][0]);     
+        }
+    }
+    else if(decay == 1){//Dalitz decays, same procedure as above
+      backgroundContainer = (AliCFContainer*)fNonHFESourceContainer[source][0]->Clone();
+      if(source == 0)
+        for(Int_t iSource = 1; iSource < kElecBgSources-3; iSource++){
+          backgroundContainer->Add(fNonHFESourceContainer[iSource][0]);   
+        }
+    }   
+  }
+  else{//careful: these containers include the non-standard electron sources (K0s, etc.). If necessary, use SetRange in species axis to fix it?
+    if(source == 0){
+      // Background Estimate
+      if(decay == 0)
+        backgroundContainer = GetContainer(kMCWeightedContainerConversionESD);
+      else if(decay == 1)
+        backgroundContainer = GetContainer(kMCWeightedContainerNonHFEESD);
+    } 
+  }
+  if(!backgroundContainer){
+    AliError("MC background container not found");
+    return NULL;
+  }
+  AliCFDataGrid *backgroundGrid = new AliCFDataGrid("backgroundGrid","backgroundGrid",*backgroundContainer,stepbackground);
+  Double_t rangelow = 1.;
+  Double_t rangeup = 6.;
+  if(decay == 1) rangelow = 0.9;
+
+
+  const Char_t *dmode[2]={"Conversions","Dalitz decays"};
+  TF1 *fitHagedorn = new TF1("fitHagedorn", "[0]/TMath::Power(([1]+x/[2]), [3])", rangelow, rangeup);
+  fNonHFEbg = 0x0;
+  TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
+  TAxis *axis = h1->GetXaxis();
+  CorrectFromTheWidth(h1);
+  if(source == 0){
+    fitHagedorn->SetParameter(0, 0.15);
+    fitHagedorn->SetParameter(1, 0.09);
+    fitHagedorn->SetParameter(2, 8.4);
+    fitHagedorn->SetParameter(3, 6.3);
+    // TCanvas *ch1conversion = new TCanvas("ch1conversion","ch1conversion",500,400);
+    // ch1conversion->cd();
+    //fHagedorn->SetLineColor(2);
+    h1->Fit(fitHagedorn,"R");
+    // h1->Draw();
+    if(!(fNonHFEbg = h1->GetFunction("fitHagedorn")))printf("electron background fit failed for %s\n",dmode[decay]);
+  }
+  
+  Int_t *nBinpp=new Int_t[1];
+  Int_t *binspp=new Int_t[1];
+  binspp[0]=kSignalPtBins;// number of pt bins
+  
+  Int_t looppt=binspp[0];
+  
+  for(Long_t iBin=1; iBin<= looppt;iBin++){       
+    Double_t iipt= h1->GetBinCenter(iBin);
+    Double_t iiwidth= axis->GetBinWidth(iBin);
+    nBinpp[0]=iBin;
+    Double_t fitFactor = backgroundGrid->GetElement(nBinpp);//if no fit available, just take bin-by-bin information
+    if(fNonHFEbg)fitFactor = fNonHFEbg->Eval(iipt)*iiwidth;
+    backgroundGrid->SetElementError(nBinpp, backgroundGrid->GetElementError(nBinpp)*evtnorm[0]);
+    backgroundGrid->SetElement(nBinpp,fitFactor*evtnorm[0]);
+  }
+  delete[] nBinpp;
+  delete[] binspp;
+  //end of workaround for statistical errors
+  if(source == 0){  
+    AliCFDataGrid *weightedBGContainer = 0x0;   
+    weightedBGContainer = new AliCFDataGrid("weightedBGContainer","weightedBGContainer",1,binspp);
+    if(decay == 0)
+      weightedBGContainer->SetGrid(GetPIDxIPEff(2));
+    else if(decay == 1)
+      weightedBGContainer->SetGrid(GetPIDxIPEff(3));
+    if(stepbackground == 3){    
+      backgroundGrid->Multiply(weightedBGContainer,1.0);
+    }  
+  }  
+  return backgroundGrid;
+}
+//____________________________________________________________
+AliCFDataGrid* AliHFEBeautySpectrum::GetCharmBackground(){
+  //
+  // calculate charm background; when using centrality binning 2: centBin = 0 for 0-20%, centBin = 5 for 40-80%
+  //
+
+  Int_t nDim = 1;
+  Double_t evtnorm=0;
+  if(fNMCEvents) evtnorm= double(fNEvents)/double(fNMCEvents); //mjmjmj
+  //if(fNMCbgEvents) evtnorm= double(fNEvents)/double(fNMCbgEvents);
+  printf("events for data: %d",fNEvents);
+  printf("events for MC: %d",fNMCEvents);
+  printf("normalization factor for charm: %f",evtnorm);
+  
+  AliCFContainer *mcContainer = GetContainer(kMCContainerCharmMC);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  if(!fCorrelation){
+    AliError("No Correlation map available");
+    return NULL;
+  }
+  AliCFDataGrid *charmBackgroundGrid= 0x0;
+  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[2];
+  bins[0]=charmbgaftertofpid->GetNbinsX();
+   
+  AliCFDataGrid *ipWeightedCharmContainer = new AliCFDataGrid("ipWeightedCharmContainer","ipWeightedCharmContainer",nDim,bins);
+  ipWeightedCharmContainer->SetGrid(GetPIDxIPEff(0)); // get charm efficiency
+  TH1D* parametrizedcharmpidipeff = (TH1D*)ipWeightedCharmContainer->Project(0);
+  
+  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 looppt=binspp[0];
+  
+  for(Long_t iBin=1; iBin<= looppt;iBin++){
+    nBinpp[0]=iBin;
+    charmBackgroundGrid->SetElementError(nBinpp, charmBackgroundGrid->GetElementError(nBinpp)*evtnorm);
+    charmBackgroundGrid->SetElement(nBinpp,charmBackgroundGrid->GetElement(nBinpp)*evtnorm);
+  }
+  
+  TH1D* charmbgafteripcut = (TH1D*)charmBackgroundGrid->Project(0);
+  
+  AliCFDataGrid *weightedCharmContainer = new AliCFDataGrid("weightedCharmContainer","weightedCharmContainer",nDim,bins);
+  weightedCharmContainer->SetGrid(GetCharmWeights(fTestCentralityLow)); // get charm weighting factors
+  TH1D* charmweightingfc = (TH1D*)weightedCharmContainer->Project(0);
+  charmBackgroundGrid->Multiply(weightedCharmContainer,1.);
+  TH1D* charmbgafterweight = (TH1D*)charmBackgroundGrid->Project(0);
+  
+  // Efficiency (set efficiency to 1 for only folding) 
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(0,0);
+  
+  // Folding
+  if(fBeamType==0)nDim = 1;
+  AliCFUnfolding folding("unfolding","",nDim,fCorrelation,efficiencyD->GetGrid(),charmBackgroundGrid->GetGrid(),charmBackgroundGrid->GetGrid());
+  folding.SetMaxNumberOfIterations(1);
+  folding.Unfold();
+  
+  // Results
+  THnSparse* result1= folding.GetEstMeasured(); // folded spectra
+  THnSparse* result=(THnSparse*)result1->Clone();
+  charmBackgroundGrid->SetGrid(result);
+  TH1D* charmbgafterfolding = (TH1D*)charmBackgroundGrid->Project(0);
+
+  //Charm background evaluation plots
+  
+  TCanvas *cCharmBgEval = new TCanvas("cCharmBgEval","cCharmBgEval",1000,600);
+  cCharmBgEval->Divide(3,1);
+  
+  cCharmBgEval->cd(1);
+  
+  if(fBeamType==0)charmbgaftertofpid->Scale(evtnorm);
+    
+  CorrectFromTheWidth(charmbgaftertofpid);
+  charmbgaftertofpid->SetMarkerStyle(25);
+  charmbgaftertofpid->Draw("p");
+  charmbgaftertofpid->GetYaxis()->SetTitle("yield normalized by # of data events");
+  charmbgaftertofpid->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+  gPad->SetLogy();
+  
+  CorrectFromTheWidth(charmbgafteripcut);
+  charmbgafteripcut->SetMarkerStyle(24);
+  charmbgafteripcut->Draw("samep");
+  
+  CorrectFromTheWidth(charmbgafterweight);
+  charmbgafterweight->SetMarkerStyle(24);
+  charmbgafterweight->SetMarkerColor(4);
+  charmbgafterweight->Draw("samep");
+  
+  CorrectFromTheWidth(charmbgafterfolding);
+  charmbgafterfolding->SetMarkerStyle(24);
+  charmbgafterfolding->SetMarkerColor(2);
+  charmbgafterfolding->Draw("samep");
+  
+  cCharmBgEval->cd(2);
+  parametrizedcharmpidipeff->SetMarkerStyle(24);
+  parametrizedcharmpidipeff->Draw("p");
+  parametrizedcharmpidipeff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+  cCharmBgEval->cd(3);
+  charmweightingfc->SetMarkerStyle(24);
+  charmweightingfc->Draw("p");
+  charmweightingfc->GetYaxis()->SetTitle("weighting factor for charm electron");
+  charmweightingfc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+  cCharmBgEval->cd(1);
+  TLegend *legcharmbg = new TLegend(0.3,0.7,0.89,0.89);
+  legcharmbg->AddEntry(charmbgaftertofpid,"After TOF PID","p");
+  legcharmbg->AddEntry(charmbgafteripcut,"After IP cut","p");
+  legcharmbg->AddEntry(charmbgafterweight,"After Weighting","p");
+  legcharmbg->AddEntry(charmbgafterfolding,"After Folding","p");
+  legcharmbg->Draw("same");
+
+  cCharmBgEval->cd(2);
+  TLegend *legcharmbg2 = new TLegend(0.3,0.7,0.89,0.89);
+  legcharmbg2->AddEntry(parametrizedcharmpidipeff,"PID + IP cut eff.","p");
+  legcharmbg2->Draw("same");
+
+  CorrectStatErr(charmBackgroundGrid);
+  if(fWriteToFile) cCharmBgEval->SaveAs("CharmBackground.eps");
+
+  delete[] bins;
+  delete[] nBinpp;
+  delete[] binspp;
+  
+  if(fUnfoldBG) UnfoldBG(charmBackgroundGrid);
+
+  return charmBackgroundGrid;
+}
+//____________________________________________________________
+AliCFDataGrid *AliHFEBeautySpectrum::CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum){  
+  //
+  // Apply TPC pid efficiency correction from parametrisation
+  //
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+    
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+  AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
+  result->SetName("ParametrizedEfficiencyBefore");
+  THnSparse *h = result->GetGrid();
+  Int_t nbdimensions = h->GetNdimensions();
+  //printf("CorrectParametrizedEfficiency::We have dimensions %d\n",nbdimensions);
+  AliCFContainer *dataContainer = GetContainer(kDataContainer);
+  if(!dataContainer){
+    AliError("Data Container not available");
+    return NULL;
+  }
+  AliCFContainer *dataContainerbis = (AliCFContainer *) dataContainer->Clone();
+  dataContainerbis->Add(dataContainerbis,-1.0);
+
+  Int_t* coord = new Int_t[nbdimensions];
+  memset(coord, 0, sizeof(Int_t) * nbdimensions);
+  Double_t* points = new Double_t[nbdimensions];
+
+  ULong64_t nEntries = h->GetNbins();
+  for (ULong64_t i = 0; i < nEntries; ++i) {   
+    Double_t value = h->GetBinContent(i, coord);
+    //Double_t valuecontainer = dataContainerbis->GetBinContent(coord,fStepData);
+    //printf("Value %f, and valuecontainer %f\n",value,valuecontainer);
+    
+    // Get the bin co-ordinates given an coord
+    for (Int_t j = 0; j < nbdimensions; ++j)
+      points[j] = h->GetAxis(j)->GetBinCenter(coord[j]);
+
+    if (!fEfficiencyFunction->IsInside(points))
+         continue;
+    TF1::RejectPoint(kFALSE);
+
+    // Evaulate function at points
+    Double_t valueEfficiency = fEfficiencyFunction->EvalPar(points, NULL);
+    //printf("Value efficiency is %f\n",valueEfficiency);
+
+    if(valueEfficiency > 0.0) {
+      h->SetBinContent(coord,value/valueEfficiency);
+      dataContainerbis->SetBinContent(coord,fStepData,value/valueEfficiency);
+    }
+    Double_t error = h->GetBinError(i);
+    h->SetBinError(coord,error/valueEfficiency);
+    dataContainerbis->SetBinError(coord,fStepData,error/valueEfficiency);  
+  } 
+
+  delete[] coord;
+  delete[] points;
+
+  AliCFDataGrid *resultt = new AliCFDataGrid("spectrumEfficiencyParametrized", "Data Grid for spectrum after Efficiency parametrized", *dataContainerbis,fStepData);
+
+// QA
+  TH1D *afterE = (TH1D *) resultt->Project(0);
+  CorrectFromTheWidth(afterE);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  CorrectFromTheWidth(beforeE);
+  fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterPE);
+  fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforePE);
+  fQA->AddResultAt(fEfficiencyFunction,AliHFEBeautySpectrumQA::kPEfficiency);
+  fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kParametrized);
+  
+  return resultt;
+}
+//____________________________________________________________
+THnSparse *AliHFEBeautySpectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Unfold and eventually correct for efficiency the bgsubspectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainerMC);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  if(!fCorrelation){
+    AliError("No Correlation map available");
+    return NULL;
+  }
+
+  // Data 
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+  
+  // Guessed
+  AliCFDataGrid* guessedGrid = new AliCFDataGrid("guessed","",*mcContainer, fStepGuessedUnfolding);
+  THnSparse* guessedTHnSparse = ((AliCFGridSparse*)guessedGrid->GetData())->GetGrid();
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+  if(!fBeauty2ndMethod)
+  {
+    // Consider parameterized IP cut efficiency
+    Int_t* bins=new Int_t[1];
+    bins[0]=kSignalPtBins;
+    
+    AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+    beffContainer->SetGrid(GetBeautyIPEff(kTRUE));
+    efficiencyD->Multiply(beffContainer,1);
+  }
+  
+
+  // Unfold 
+  
+  AliCFUnfolding unfolding("unfolding","",fNbDimensions,fCorrelation,efficiencyD->GetGrid(),dataGrid->GetGrid(),guessedTHnSparse,1.e-06,0,fNumberOfIterations);//check with MinJung if last arguments are correct (taken from inclusive analysis...
+  if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
+  unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
+  if(fSetSmoothing) unfolding.UseSmoothing();
+  unfolding.Unfold();
+
+  // Results
+  THnSparse* result = unfolding.GetUnfolded();
+  THnSparse* residual = unfolding.GetEstMeasured();
+
+ // QA
+  TH1D *residualh = (TH1D *) residual->Projection(0);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
+  TH1D *afterE = (TH1D *) result->Projection(0);
+
+  CorrectFromTheWidth(residualh);
+  CorrectFromTheWidth(beforeE);
+  CorrectFromTheWidth(afterE);
+  fQA->AddResultAt(residualh,AliHFEBeautySpectrumQA::kResidualU);
+  fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterU);
+  fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeU);
+  fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kUEfficiency);
+  fQA->DrawUnfolding();
+
+  return (THnSparse *) result->Clone();
+}
+
+//____________________________________________________________
+AliCFDataGrid *AliHFEBeautySpectrum::CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum){
+  
+  //
+  // Apply unfolding and efficiency correction together to bgsubspectrum
+  //
+
+  AliCFContainer *mcContainer = GetContainer(kMCContainerESD);
+  if(!mcContainer){
+    AliError("MC Container not available");
+    return NULL;
+  }
+
+  // Efficiency
+  AliCFEffGrid* efficiencyD = new AliCFEffGrid("efficiency","",*mcContainer);
+  efficiencyD->CalculateEfficiency(fStepMC,fStepTrue);
+
+
+  if(!fBeauty2ndMethod)
+  {
+    // Consider parameterized IP cut efficiency
+    Int_t* bins=new Int_t[1];
+    bins[0]=kSignalPtBins;
+    
+    AliCFEffGrid *beffContainer = new AliCFEffGrid("beffContainer","beffContainer",1,bins);
+    beffContainer->SetGrid(GetBeautyIPEff(kFALSE));
+    efficiencyD->Multiply(beffContainer,1);
+  }
+
+  // Data in the right format
+  AliCFDataGrid *dataGrid = 0x0;  
+  if(bgsubpectrum) {
+    dataGrid = bgsubpectrum;
+  }
+  else {
+    
+    AliCFContainer *dataContainer = GetContainer(kDataContainer);
+    if(!dataContainer){
+      AliError("Data Container not available");
+      return NULL;
+    }
+
+    dataGrid = new AliCFDataGrid("dataGrid","dataGrid",*dataContainer, fStepData);
+  } 
+
+  // Correct
+  AliCFDataGrid *result = (AliCFDataGrid *) dataGrid->Clone();
+  result->ApplyEffCorrection(*efficiencyD);
+
+  // QA
+  TH1D *afterE = (TH1D *) result->Project(0);
+  CorrectFromTheWidth(afterE);
+  TH1D *beforeE = (TH1D *) dataGrid->Project(0);
+  CorrectFromTheWidth(beforeE);
+  TH1D* efficiencyDproj = (TH1D *) efficiencyD->Project(0);
+  fQA->AddResultAt(afterE,AliHFEBeautySpectrumQA::kAfterMCE);
+  fQA->AddResultAt(beforeE,AliHFEBeautySpectrumQA::kBeforeMCE);
+  fQA->AddResultAt(efficiencyDproj,AliHFEBeautySpectrumQA::kMCEfficiency);
+  fQA->DrawCorrectWithEfficiency(AliHFEBeautySpectrumQA::kMC);
+  
+  return result;
+
+}
+//____________________________________________________________
+void AliHFEBeautySpectrum::AddTemporaryObject(TObject *o){
+  // 
+  // Emulate garbage collection on class level
+  // 
+  if(!fTemporaryObjects) {
+    fTemporaryObjects= new TList;
+    fTemporaryObjects->SetOwner();
+  }
+  fTemporaryObjects->Add(o);
+}
+
+//____________________________________________________________
+void AliHFEBeautySpectrum::ClearObject(TObject *o){
+  //
+  // Do a safe deletion
+  //
+  if(fTemporaryObjects){
+    if(fTemporaryObjects->FindObject(o)) fTemporaryObjects->Remove(o);
+    delete o;
+  } else{
+    // just delete
+    delete o;
+  }
+}
+//_________________________________________________________________________
+THnSparse* AliHFEBeautySpectrum::GetCharmWeights(Int_t centBin){
+  //
+  // Measured D->e based weighting factors
+  //
+
+  const Int_t nDimpp=1;
+  Int_t nBinpp[nDimpp] = {kSignalPtBins};
+  Double_t ptbinning1[kSignalPtBins+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.};
+  Int_t looppt=nBinpp[0];
+  
+  fWeightCharm = new THnSparseF("weightHisto", "weighting factor; pt[GeV/c]", nDimpp, nBinpp);
+  fWeightCharm->SetBinEdges(0, ptbinning1);
+  // //if(fBeamType == 0){// Weighting factor for pp
+  //Double_t weightpp[kSignalPtBins]={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};
+//TPC+TOF standard cut 4800
+  Double_t weightpp[kSignalPtBins]={0.050326, 0.045826, 0.042043, 0.039641, 0.039872, 0.041796, 0.044320, 0.046273, 0.047376, 0.047657, 0.047973, 0.048307, 0.045325, 0.044067, 0.043367, 0.040417, 0.037048, 0.033695, 0.032192, 0.029270, 0.027270, 0.024451, 0.020846, 0.019032, 0.018210, 0.017554, 0.015604, 0.015194, 0.013542, 0.013447, 0.015160, 0.015507, 0.014989, 0.012533, 0.025603}; 
+  //}
+  //else{
+  //if(centBin == 0){
+      // Weighting factor for PbPb (0-20%)
+  Double_t weightPbPb1[kSignalPtBins]={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};
+  //}
+  //else{
+  // Weighting factor for PbPb (40-80%)
+  Double_t weightPbPb2[kSignalPtBins]={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[kSignalPtBins];
+  for(Int_t ipt = 0; ipt < kSignalPtBins; ipt++){
+    if(fBeamType == 0)weight[ipt] = weightpp[ipt];
+    else if(centBin == 0)weight[ipt] = weightPbPb1[ipt];
+    else weight[ipt] = weightPbPb2[ipt];
+  }  
+
+  //points
+  Double_t pt[1];
+  Double_t contents[2];
+  
+  for(int i=0; i<looppt; i++){
+    pt[0]=(ptbinning1[i]+ptbinning1[i+1])/2.;
+    contents[0]=pt[0];
+    fWeightCharm->Fill(contents,weight[i]);
+  }
+  
+  
+  Int_t nDimSparse = fWeightCharm->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
+
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = fWeightCharm->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    fWeightCharm->SetBinError(binfill,0.); // put 0 everywhere
+  }
+  delete[] binsvar;
+  delete[] binfill;
+
+  return fWeightCharm;
+}
+//____________________________________________________________________________
+void AliHFEBeautySpectrum::SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions){
+
+   // TOF PID efficiencies
+  
+   TF1 *fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
+   TF1 *fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,6.);
+   TF1 *fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.0);
+
+   if(fBeamType == 1){//not in use yet - adapt as soon as possible!
+     fittofpid = new TF1("fittofpid","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
+     fipfit = new TF1("fipfit","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.95,8.);
+     fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
+   }
+
+   TCanvas * cefficiencyParamtof = new TCanvas("efficiencyParamtof","efficiencyParamtof",600,600);
+   cefficiencyParamtof->cd();
+
+   AliCFContainer *mccontainermcD = 0x0;
+   AliCFContainer *mccontaineresdD = 0x0;
+   TH1D* efficiencysigTOFPIDD;
+   TH1D* efficiencyTOFPIDD;
+   TH1D* efficiencysigesdTOFPIDD;
+   TH1D* efficiencyesdTOFPIDD;
+   Int_t source = -1; //get parameterized TOF PID efficiencies
+
+   // signal sample
+   mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   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,fTestCentralityLow,fTestCentralityHigh);
+   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,fTestCentralityLow,fTestCentralityHigh);
+   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,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencysigesdParamTOFPID= new AliCFEffGrid("efficiencysigesdParamTOFPID","",*mccontainermcD);
+   efficiencysigesdParamTOFPID->CalculateEfficiency(fStepMC,fStepMC-1); // TOF PID efficiencies
+   
+   //fill histo
+   efficiencysigTOFPIDD = (TH1D *) efficiencymcsigParamTOFPID->Project(0);
+   efficiencyTOFPIDD = (TH1D *) efficiencymcParamTOFPID->Project(0);
+   efficiencysigesdTOFPIDD = (TH1D *) efficiencysigesdParamTOFPID->Project(0);
+   efficiencyesdTOFPIDD = (TH1D *) efficiencyesdParamTOFPID->Project(0);
+   efficiencysigTOFPIDD->SetName("efficiencysigTOFPIDD");
+   efficiencyTOFPIDD->SetName("efficiencyTOFPIDD");
+   efficiencysigesdTOFPIDD->SetName("efficiencysigesdTOFPIDD");
+   efficiencyesdTOFPIDD->SetName("efficiencyesdTOFPIDD");
+   
+   //fit (mc enhenced sample)
+   fittofpid->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+   efficiencysigTOFPIDD->Fit(fittofpid,"R");
+   efficiencysigTOFPIDD->GetYaxis()->SetTitle("Efficiency");
+   fEfficiencyTOFPIDD = efficiencysigTOFPIDD->GetFunction("fittofpid");
+   
+   //fit (esd enhenced sample)
+   efficiencysigesdTOFPIDD->Fit(fittofpid,"R");
+   efficiencysigesdTOFPIDD->GetYaxis()->SetTitle("Efficiency");
+   fEfficiencyesdTOFPIDD = efficiencysigesdTOFPIDD->GetFunction("fittofpid");
+   
+   // draw (for PbPb, only 1st bin)
+   //sig mc
+   efficiencysigTOFPIDD->SetTitle("");
+   efficiencysigTOFPIDD->SetStats(0);
+   efficiencysigTOFPIDD->SetMarkerStyle(25);
+   efficiencysigTOFPIDD->SetMarkerColor(2);
+   efficiencysigTOFPIDD->SetLineColor(2);
+   efficiencysigTOFPIDD->Draw();
+
+   //mb mc
+   efficiencyTOFPIDD->SetTitle("");
+   efficiencyTOFPIDD->SetStats(0);
+   efficiencyTOFPIDD->SetMarkerStyle(24);
+   efficiencyTOFPIDD->SetMarkerColor(4);
+   efficiencyTOFPIDD->SetLineColor(4);
+   efficiencyTOFPIDD->Draw("same");
+
+   //sig esd
+   efficiencysigesdTOFPIDD->SetTitle("");
+   efficiencysigesdTOFPIDD->SetStats(0);
+   efficiencysigesdTOFPIDD->SetMarkerStyle(25);
+   efficiencysigesdTOFPIDD->SetMarkerColor(3);
+   efficiencysigesdTOFPIDD->SetLineColor(3);
+   efficiencysigesdTOFPIDD->Draw("same");
+
+   //mb esd
+   efficiencyesdTOFPIDD->SetTitle("");
+   efficiencyesdTOFPIDD->SetStats(0);
+   efficiencyesdTOFPIDD->SetMarkerStyle(25);
+   efficiencyesdTOFPIDD->SetMarkerColor(1);
+   efficiencyesdTOFPIDD->SetLineColor(1);
+   efficiencyesdTOFPIDD->Draw("same");
+
+   //signal mc fit
+   if(fEfficiencyTOFPIDD){
+     fEfficiencyTOFPIDD->SetLineColor(2);
+     fEfficiencyTOFPIDD->Draw("same");
+   }
+   //mb esd fit
+   if(fEfficiencyesdTOFPIDD){
+       fEfficiencyesdTOFPIDD->SetLineColor(3);
+       fEfficiencyesdTOFPIDD->Draw("same");
+     }
+
+   TLegend *legtofeff = new TLegend(0.3,0.15,0.79,0.44);
+   legtofeff->AddEntry(efficiencysigTOFPIDD,"TOF PID Step Efficiency","");
+   legtofeff->AddEntry(efficiencysigTOFPIDD,"vs MC p_{t} for enhenced samples","p");
+   legtofeff->AddEntry(efficiencyTOFPIDD,"vs MC p_{t} for mb samples","p");
+   legtofeff->AddEntry(efficiencysigesdTOFPIDD,"vs esd p_{t} for enhenced samples","p");
+   legtofeff->AddEntry(efficiencyesdTOFPIDD,"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
+   AliCFContainer *charmCombined = 0x0; 
+   AliCFContainer *beautyCombined = 0x0;
+   AliCFContainer *beautyCombinedesd = 0x0;
+
+   source = 0; //charm enhenced
+   mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencyCharmSig = new AliCFEffGrid("efficiencyCharmSig","",*mccontainermcD);
+   efficiencyCharmSig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+   charmCombined= (AliCFContainer*)mccontainermcD->Clone("charmCombined");  
+
+   source = 1; //beauty enhenced
+   mccontainermcD = GetSlicedContainer(container, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencyBeautySig = new AliCFEffGrid("efficiencyBeautySig","",*mccontainermcD);
+   efficiencyBeautySig->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+
+   beautyCombined = (AliCFContainer*)mccontainermcD->Clone("beautyCombined"); 
+
+   mccontaineresdD = GetSlicedContainer(containeresd, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   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
+   mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencyCharm = new AliCFEffGrid("efficiencyCharm","",*mccontainermcD);
+   efficiencyCharm->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+   
+   charmCombined->Add(mccontainermcD); 
+   AliCFEffGrid* efficiencyCharmCombined = new AliCFEffGrid("efficiencyCharmCombined","",*charmCombined); 
+   efficiencyCharmCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
+
+   source = 1; //beauty mb
+   mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencyBeauty = new AliCFEffGrid("efficiencyBeauty","",*mccontainermcD);
+   efficiencyBeauty->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency. 
+   
+   beautyCombined->Add(mccontainermcD);
+   AliCFEffGrid* efficiencyBeautyCombined = new AliCFEffGrid("efficiencyBeautyCombined","",*beautyCombined); 
+   efficiencyBeautyCombined->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); 
+   
+   mccontaineresdD = GetSlicedContainer(containeresdmb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   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
+   mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   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 mb
+   mccontainermcD = GetSlicedContainer(containermb, fNbDimensions, dimensions, source, fChargeChoosen,fTestCentralityLow,fTestCentralityHigh);
+   AliCFEffGrid* efficiencyNonhfe= new AliCFEffGrid("efficiencyNonhfe","",*mccontainermcD);
+   efficiencyNonhfe->CalculateEfficiency(AliHFEcuts::kNcutStepsMCTrack + fStepData,AliHFEcuts::kNcutStepsMCTrack + fStepData-1); // ip cut efficiency.
+
+   if(fIPEffCombinedSamples){printf("combined samples taken for beauty and charm\n");
+     fEfficiencyCharmSigD = (TH1D*)efficiencyCharmCombined->Project(0); //signal enhenced + mb 
+     fEfficiencyBeautySigD = (TH1D*)efficiencyBeautyCombined->Project(0); //signal enhenced + mb
+     fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautyCombinedesd->Project(0); //signal enhenced + mb
+     }
+     else{printf("signal enhanced samples taken for beauty and charm\n");
+       fEfficiencyCharmSigD = (TH1D*)efficiencyCharmSig->Project(0); //signal enhenced only
+       fEfficiencyBeautySigD = (TH1D*)efficiencyBeautySig->Project(0); //signal enhenced only
+       fEfficiencyBeautySigesdD = (TH1D*)efficiencyBeautySigesd->Project(0); //signal enhenced only
+     }
+     fCharmEff = (TH1D*)efficiencyCharm->Project(0); //mb only
+     fBeautyEff = (TH1D*)efficiencyBeauty->Project(0); //mb only
+     fConversionEff = (TH1D*)efficiencyConv->Project(0); //mb only
+     fNonHFEEff = (TH1D*)efficiencyNonhfe->Project(0); //mb only
+
+   if(fBeamType==0){
+     //AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerNonHFEESD),1,0);
+     AliCFEffGrid  *nonHFEEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerNonHFEESDSig),1,0); //mjmj
+     fNonHFEEffbgc = (TH1D *) nonHFEEffGrid->Project(0);
+     
+     //AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(kMCWeightedContainerConversionESD),1,0);
+     AliCFEffGrid  *conversionEffGrid = (AliCFEffGrid*)  GetEfficiency(GetContainer(AliHFECorrectSpectrumBase::kMCWeightedContainerConversionESDSig),1,0); //mjmj
+     fConversionEffbgc = (TH1D *) conversionEffGrid->Project(0);
+     
+     //MHMH
+     if(fBeamType == 0){
+       fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+       fNonHFEEffbgc->Fit(fipfitnonhfe,"R"); 
+       fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
+       
+       fipfitnonhfe = new TF1("fipfitnonhfe","[0]*([1]+[2]*log(x)+[3]*log(x)*log(x))*tanh([4]*x-[5])",0.5,8.);
+       fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);      
+       fConversionEffbgc->Fit(fipfitnonhfe,"R");
+       fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
+     }
+     //MHMH
+   }
+   
+   fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+   fipfit->SetLineColor(2);
+   fEfficiencyBeautySigD->Fit(fipfit,"R");
+   fEfficiencyBeautySigD->GetYaxis()->SetTitle("Efficiency");
+   fEfficiencyIPBeautyD = fEfficiencyBeautySigD->GetFunction("fipfit");
+   
+   fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+   fipfit->SetLineColor(6);
+   fEfficiencyBeautySigesdD->Fit(fipfit,"R");
+   fEfficiencyBeautySigesdD->GetYaxis()->SetTitle("Efficiency");
+   fEfficiencyIPBeautyesdD = fEfficiencyBeautySigesdD->GetFunction("fipfit");
+   
+   fipfit->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+   fipfit->SetLineColor(1);
+   fEfficiencyCharmSigD->Fit(fipfit,"R");
+   fEfficiencyCharmSigD->GetYaxis()->SetTitle("Efficiency");
+   fEfficiencyIPCharmD = fEfficiencyCharmSigD->GetFunction("fipfit");
+   
+   //if(0){
+   if(fIPParameterizedEff){
+     // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+     fipfitnonhfe->SetLineColor(3);
+     if(fBeamType==1){
+       fConversionEff->Fit(fipfitnonhfe,"R");
+       fConversionEff->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPConversionD = fConversionEff->GetFunction("fipfitnonhfe");
+     }
+     else{
+       fConversionEffbgc->Fit(fipfitnonhfe,"R");
+       fConversionEffbgc->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPConversionD = fConversionEffbgc->GetFunction("fipfitnonhfe");
+     }       
+     // fipfitnonhfe->SetParameters(0.5,0.319,0.0157,0.00664,6.77,2.08);
+     fipfitnonhfe->SetLineColor(4);
+     if(fBeamType==1){
+       fNonHFEEff->Fit(fipfitnonhfe,"R");
+       fNonHFEEff->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPNonhfeD = fNonHFEEff->GetFunction("fipfitnonhfe");
+     }
+     else{
+       fNonHFEEffbgc->Fit(fipfitnonhfe,"R");
+       fNonHFEEffbgc->GetYaxis()->SetTitle("Efficiency");
+       fEfficiencyIPNonhfeD = fNonHFEEffbgc->GetFunction("fipfitnonhfe");
+     }
+   }
+   
+   // draw
+   fEfficiencyCharmSigD->SetMarkerStyle(21);
+   fEfficiencyCharmSigD->SetMarkerColor(1);
+   fEfficiencyCharmSigD->SetLineColor(1);
+   fEfficiencyBeautySigD->SetMarkerStyle(21);
+   fEfficiencyBeautySigD->SetMarkerColor(2);
+   fEfficiencyBeautySigD->SetLineColor(2);
+   fEfficiencyBeautySigesdD->SetStats(0);
+   fEfficiencyBeautySigesdD->SetMarkerStyle(21);
+   fEfficiencyBeautySigesdD->SetMarkerColor(6);
+   fEfficiencyBeautySigesdD->SetLineColor(6);
+   fCharmEff->SetMarkerStyle(24);
+   fCharmEff->SetMarkerColor(1);
+   fCharmEff->SetLineColor(1);
+   fBeautyEff->SetMarkerStyle(24);
+   fBeautyEff->SetMarkerColor(2);
+   fBeautyEff->SetLineColor(2);
+   fConversionEff->SetMarkerStyle(24);
+   fConversionEff->SetMarkerColor(3);
+   fConversionEff->SetLineColor(3);
+   fNonHFEEff->SetMarkerStyle(24);
+   fNonHFEEff->SetMarkerColor(4);
+   fNonHFEEff->SetLineColor(4);
+
+   fEfficiencyCharmSigD->Draw();
+   fEfficiencyCharmSigD->GetXaxis()->SetRangeUser(0.0,7.9);
+   fEfficiencyCharmSigD->GetYaxis()->SetRangeUser(0.0,0.5);
+
+   fEfficiencyBeautySigD->Draw("same");
+   fEfficiencyBeautySigesdD->Draw("same");
+   if(fBeamType == 1){
+     fNonHFEEff->Draw("same");
+     fConversionEff->Draw("same");
+   }
+   //fCharmEff->Draw("same");
+   //fBeautyEff->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");
+   }
+   if(fEfficiencyIPBeautyD)
+      fEfficiencyIPBeautyD->Draw("same");
+   if(fEfficiencyIPBeautyesdD)
+     fEfficiencyIPBeautyesdD->Draw("same");
+   if( fEfficiencyIPCharmD)
+     fEfficiencyIPCharmD->Draw("same");
+   if(fIPParameterizedEff){
+     if(fEfficiencyIPConversionD)
+       fEfficiencyIPConversionD->Draw("same");
+     if(fEfficiencyIPNonhfeD)
+       fEfficiencyIPNonhfeD->Draw("same");
+   }
+   TLegend *legipeff = new TLegend(0.58,0.2,0.88,0.39);
+   legipeff->AddEntry(fEfficiencyBeautySigD,"IP Step Efficiency","");
+   legipeff->AddEntry(fEfficiencyBeautySigD,"beauty e","p");
+   legipeff->AddEntry(fEfficiencyBeautySigesdD,"beauty e(esd pt)","p");
+   legipeff->AddEntry(fEfficiencyCharmSigD,"charm e","p");
+   if(fBeamType == 0){
+     legipeff->AddEntry(fConversionEffbgc,"conversion e(esd pt)","p");
+     legipeff->AddEntry(fNonHFEEffbgc,"Dalitz e(esd pt)","p");
+   }
+   else{
+     legipeff->AddEntry(fConversionEff,"conversion e","p");
+     legipeff->AddEntry(fNonHFEEff,"Dalitz e","p");
+   }
+   legipeff->Draw("same");
+   gPad->SetGrid();
+   //cefficiencyParamIP->SaveAs("efficiencyParamIP.eps");
+}
+
+//____________________________________________________________________________
+THnSparse* AliHFEBeautySpectrum::GetBeautyIPEff(Bool_t isMCpt){
+  //
+  // Return beauty electron IP cut efficiency
+  //
+
+  const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
+  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
+  Int_t nDim=1;  //dimensions of the efficiency weighting grid
+  
+  Int_t nBin[1] = {nPtbinning1};
+  THnSparseF *ipcut = new THnSparseF("beff", "b IP efficiency; p_{t}(GeV/c)", nDim, nBin);
+  
+  ipcut->SetBinEdges(0, kPtRange);
+  
+  Double_t pt[1];
+  Double_t weight;
+  Double_t weightErr;
+  Double_t contents[2];
+
+  weight = 1.0;
+  weightErr = 1.0;
+  
+  Int_t looppt=nBin[0];
+  Int_t ibin[2];
+  
+  for(int i=0; i<looppt; i++)
+    {
+      pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+      if(isMCpt){
+        if(fEfficiencyIPBeautyD){
+          weight=fEfficiencyIPBeautyD->Eval(pt[0]);
+          weightErr = 0;
+        }
+        else{
+          printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
+          weight = fEfficiencyBeautySigD->GetBinContent(i+1); 
+          weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
+        }
+      }
+      else{
+        if(fEfficiencyIPBeautyesdD){
+          weight=fEfficiencyIPBeautyesdD->Eval(pt[0]);
+          weightErr = 0;
+        }
+        else{
+          printf("Fit failed on beauty IP cut efficiency. Contents in histo used!\n");
+          weight = fEfficiencyBeautySigesdD->GetBinContent(i+1);
+          weightErr = fEfficiencyBeautySigD->GetBinError(i+1);
+        }
+      }
+      
+      contents[0]=pt[0];
+      ibin[0]=i+1;
+      
+      ipcut->Fill(contents,weight);
+      ipcut->SetBinError(ibin,weightErr);
+    }
+  
+  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
+
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binsvar[iVar] = ipcut->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    //ipcut->SetBinError(binfill,0.); // put 0 everywhere
+  }
+
+  delete[] binsvar;
+  delete[] binfill;
+
+  return ipcut;
+}
+
+//____________________________________________________________________________
+THnSparse* AliHFEBeautySpectrum::GetPIDxIPEff(Int_t source){
+  //
+  // Return PID x IP cut efficiency
+  //
+  const Int_t nPtbinning1 = kSignalPtBins;//number of pt bins, according to new binning
+  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
+  Int_t nDim=1;  //dimensions of the efficiency weighting grid
+  
+  Int_t nBin[1] = {nPtbinning1};
+  
+  THnSparseF *pideff;
+  pideff = new THnSparseF("pideff", "PID efficiency; p_{t}(GeV/c)", nDim, nBin);
+  pideff->SetBinEdges(0, kPtRange);
+  
+  Double_t pt[1];
+  Double_t weight;
+  Double_t weightErr;
+  Double_t contents[2];
+  
+  weight = 1.0;
+  weightErr = 1.0;
+  
+  Int_t looppt=nBin[0];
+  Int_t ibin[2];
+      
+  Double_t trdtpcPidEfficiency = fEfficiencyFunction->Eval(0); // assume we have constant TRD+TPC PID efficiency
+  for(int i=0; i<looppt; i++)
+    {
+      pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+      
+      Double_t tofpideff = 1.;
+      Double_t tofpideffesd = 1.;
+      if(fEfficiencyTOFPIDD)
+        tofpideff = fEfficiencyTOFPIDD->Eval(pt[0]); 
+      else{
+        printf("TOF PID fit failed on conversion. The result is wrong!\n");
+      }  
+      if(fEfficiencyesdTOFPIDD)
+        tofpideffesd = fEfficiencyesdTOFPIDD->Eval(pt[0]);
+      else{
+        printf("TOF esd PID fit failed on conversion. The result is wrong!\n");
+      }
+      
+      //tof pid eff x tpc pid eff x ip cut eff
+      if(fIPParameterizedEff){
+        if(source==0) {
+          if(fEfficiencyIPCharmD){
+            weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
+            weightErr = 0; 
+          }
+          else{
+            printf("Fit failed on charm IP cut efficiency\n");
+            weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
+            weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1); 
+          }
+        } 
+        else if(source==2) {
+          if(fEfficiencyIPConversionD){
+            weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPConversionD->Eval(pt[0]); 
+            weightErr = 0; 
+          }
+          else{
+            printf("Fit failed on conversion IP cut efficiency\n");
+            weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1);
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
+          }
+        }
+        else if(source==3) {
+          if(fEfficiencyIPNonhfeD){
+            weight = tofpideffesd*trdtpcPidEfficiency*fEfficiencyIPNonhfeD->Eval(pt[0]); 
+            weightErr = 0; 
+          }
+          else{
+            printf("Fit failed on dalitz IP cut efficiency\n");
+            weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1);
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
+          }  
+        }
+      }
+      else{
+        if(source==0){ 
+          if(fEfficiencyIPCharmD){
+            weight = tofpideff*trdtpcPidEfficiency*fEfficiencyIPCharmD->Eval(pt[0]);
+            weightErr = 0;
+          }
+          else{
+            printf("Fit failed on charm IP cut efficiency\n");
+            weight = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinContent(i+1);
+            weightErr = tofpideff*trdtpcPidEfficiency*fEfficiencyCharmSigD->GetBinError(i+1);
+          }
+        }
+        else if(source==2){
+          if(fBeamType==0){
+            weight = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinContent(i+1); // conversion
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEffbgc->GetBinError(i+1);
+          }
+          else{
+            weight = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinContent(i+1); // conversion
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fConversionEff->GetBinError(i+1);
+          }
+        }
+        else if(source==3){
+          if(fBeamType==0){
+            weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinContent(i+1); // conversion
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEffbgc->GetBinError(i+1);
+          }
+          else{ 
+            weight = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinContent(i+1); // Dalitz
+            weightErr = tofpideffesd*trdtpcPidEfficiency*fNonHFEEff->GetBinError(i+1);
+          }
+        }
+      }
+      
+      contents[0]=pt[0];
+      ibin[0]=i+1;
+      
+      pideff->Fill(contents,weight);
+      pideff->SetBinError(ibin,weightErr);
+    }
+  
+  Int_t nDimSparse = pideff->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
+  
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+    binsvar[iVar] = pideff->GetAxis(iVar)->GetNbins();
+      nBins *= binsvar[iVar];
+  }
+
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+  }
+
+  delete[] binsvar;
+  delete[] binfill;
+
+
+  return pideff;
+}
+
+//__________________________________________________________________________
+AliCFDataGrid *AliHFEBeautySpectrum::GetRawBspectra2ndMethod(){
+ //
+ // retrieve AliCFDataGrid for raw beauty spectra obtained from fit method
+    //
+  Int_t nDim = 1;
+  
+  const Int_t nPtbinning1 = 18;//number of pt bins, according to new binning
+  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};
+  
+  Int_t nBin[1] = {nPtbinning1};
+  
+  AliCFDataGrid *rawBeautyContainer = new AliCFDataGrid("rawBeautyContainer","rawBeautyContainer",nDim,nBin);
+      
+  THnSparseF *brawspectra;
+  brawspectra= new THnSparseF("brawspectra", "beauty yields ; p_{t}(GeV/c)", nDim, nBin);
+  
+  brawspectra->SetBinEdges(0, kPtRange);
+      
+  Double_t pt[1];
+  Double_t yields= 0.;
+  Double_t valuesb[2];
+            
+  for(int i=0; i<fBSpectrum2ndMethod->GetNbinsX(); i++){
+    pt[0]=(kPtRange[i]+kPtRange[i+1])/2.;
+    
+    yields = fBSpectrum2ndMethod->GetBinContent(i+1);
+        
+    valuesb[0]=pt[0];
+    brawspectra->Fill(valuesb,yields);
+  }
+  
+  
+  
+  Int_t nDimSparse = brawspectra->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
+    
+  for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+    binsvar[iVar] = brawspectra->GetAxis(iVar)->GetNbins();
+    nBins *= binsvar[iVar];
+  }
+    
+  Int_t *binfill = new Int_t[nDimSparse]; // bin to fill the THnSparse (holding the bin coordinates)
+  // loop that sets 0 error in each bin
+  for (Long_t iBin=0; iBin<nBins; iBin++) {
+    Long_t bintmp = iBin ;
+    for (Int_t iVar=0; iVar<nDimSparse; iVar++) {
+      binfill[iVar] = 1 + bintmp % binsvar[iVar] ;
+      bintmp /= binsvar[iVar] ;
+    }
+    brawspectra->SetBinError(binfill,0.); // put 0 everywhere
+  }
+  
+  
+  rawBeautyContainer->SetGrid(brawspectra); // get charm efficiency
+  TH1D* hRawBeautySpectra = (TH1D*)rawBeautyContainer->Project(0);
+  
+  new TCanvas;
+  fBSpectrum2ndMethod->SetMarkerStyle(24);
+  fBSpectrum2ndMethod->Draw("p");
+  hRawBeautySpectra->SetMarkerStyle(25);
+  hRawBeautySpectra->Draw("samep");
+  
+  delete[] binfill;
+  delete[] binsvar; 
+  
+  return rawBeautyContainer;
+}
+//__________________________________________________________________________
+void AliHFEBeautySpectrum::CalculateNonHFEsyst(){
+  //
+  // Calculate non HFE sys
+  //
+  //
+
+  if(!fNonHFEsyst)
+    return;
+
+  Double_t evtnorm[1] = {0.0};
+  if(fNMCbgEvents>0) evtnorm[0]= double(fNEvents)/double(fNMCbgEvents);
+  
+  AliCFDataGrid *convSourceGrid[kElecBgSources-3][kBgLevels];
+  AliCFDataGrid *nonHFESourceGrid[kElecBgSources-3][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]=kSignalPtBins;//fConversionEff[centrality]->GetNbinsX();
+   
+  AliCFDataGrid *weightedConversionContainer = new AliCFDataGrid("weightedConversionContainer","weightedConversionContainer",1,bins);
+  AliCFDataGrid *weightedNonHFEContainer = new AliCFDataGrid("weightedNonHFEContainer","weightedNonHFEContainer",1,bins);
+
+  for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){   
+    for(Int_t iSource = 0; iSource < kElecBgSources-3; iSource++){
+      convSourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("convGrid_%d_%d",iSource,iLevel),Form("convGrid_%d_%d",iSource,iLevel),*fConvSourceContainer[iSource][iLevel],stepbackground);
+      weightedConversionContainer->SetGrid(GetPIDxIPEff(2));
+      convSourceGrid[iSource][iLevel]->Multiply(weightedConversionContainer,1.0);
+      
+      nonHFESourceGrid[iSource][iLevel] = new AliCFDataGrid(Form("nonHFEGrid_%d_%d",iSource,iLevel),Form("nonHFEGrid_%d_%d",iSource,iLevel),*fNonHFESourceContainer[iSource][iLevel],stepbackground);
+      weightedNonHFEContainer->SetGrid(GetPIDxIPEff(3));
+      nonHFESourceGrid[iSource][iLevel]->Multiply(weightedNonHFEContainer,1.0);
+    }
+    
+    bgConvGrid[iLevel] = (AliCFDataGrid*)convSourceGrid[0][iLevel]->Clone();
+    for(Int_t iSource = 2; iSource < kElecBgSources-3; 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 = 2; iSource < kElecBgSources-3; 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[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; 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
+    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 (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-4];
+  for(Int_t iSource = 1; iSource < kElecBgSources-3; 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);   
+  }
+  
+  //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,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-3; 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++){
+      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)+(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));           
+  }
+     
+  TCanvas *cPiErrors = new TCanvas("cPiErrors","cPiErrors",1000,600);
+  cPiErrors->cd();
+  cPiErrors->SetLogx();
+  cPiErrors->SetLogy();
+  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(kBlue);
+  hScalingErrors->Draw("SAME");
+
+  TLegend *lPiErr = new TLegend(0.6,0.6, 0.95,0.95);
+  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->Draw("SAME");
+
+  //Normalize errors
+  TH1D *hUpSystScaled = (TH1D*)hOverallSystErrUp->Clone();
+  TH1D *hLowSystScaled = (TH1D*)hOverallSystErrLow->Clone();
+  hUpSystScaled->Scale(evtnorm[0]);//scale by N(data)/N(MC), to make data sets comparable to saved subtracted spectrum (calculations in separate macro!)
+  hLowSystScaled->Scale(evtnorm[0]);
+  TH1D *hNormAllSystErrUp = (TH1D*)hUpSystScaled->Clone();
+  TH1D *hNormAllSystErrLow = (TH1D*)hLowSystScaled->Clone();
+  //histograms to be normalized to TGraphErrors
+  CorrectFromTheWidth(hNormAllSystErrUp);
+  CorrectFromTheWidth(hNormAllSystErrLow);
+
+  TCanvas *cNormOvErrs = new TCanvas("cNormOvErrs","cNormOvErrs");
+  cNormOvErrs->cd();
+  cNormOvErrs->SetLogx();
+  cNormOvErrs->SetLogy();
+
+  TGraphErrors* gOverallSystErrUp = NormalizeTH1(hNormAllSystErrUp);
+  TGraphErrors* gOverallSystErrLow = NormalizeTH1(hNormAllSystErrLow);
+  gOverallSystErrUp->SetTitle("Overall Systematic non-HFE Errors");
+  gOverallSystErrUp->SetMarkerColor(kBlack);
+  gOverallSystErrUp->SetLineColor(kBlack);
+  gOverallSystErrLow->SetMarkerColor(kRed);
+  gOverallSystErrLow->SetLineColor(kRed);
+  gOverallSystErrUp->Draw("AP");
+  gOverallSystErrLow->Draw("PSAME");
+  TLegend *lAllSys = new TLegend(0.4,0.6,0.89,0.89);
+  lAllSys->AddEntry(gOverallSystErrLow,"lower","p");
+  lAllSys->AddEntry(gOverallSystErrUp,"upper","p");
+  lAllSys->Draw("same");
+
+
+  AliCFDataGrid *bgYieldGrid;
+  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();
+  hRelErrUp->Divide(hBgYield);
+  TH1D* hRelErrLow = (TH1D*)hOverallSystErrLow->Clone();
+  hRelErrLow->Divide(hBgYield);
+
+  TCanvas *cRelErrs = new TCanvas("cRelErrs","cRelErrs");
+  cRelErrs->cd();
+  cRelErrs->SetLogx();
+  hRelErrUp->SetTitle("Relative error of non-HFE background yield");
+  hRelErrUp->Draw();
+  hRelErrLow->SetLineColor(kBlack);
+  hRelErrLow->Draw("SAME");
+
+  TLegend *lRel = new TLegend(0.6,0.6,0.95,0.95);
+  lRel->AddEntry(hRelErrUp, "upper");
+  lRel->AddEntry(hRelErrLow, "lower");
+  lRel->Draw("SAME");
+
+  //CorrectFromTheWidth(hBgYield);
+  //hBgYield->Scale(evtnorm[0]);
+  //write histograms/TGraphs to file
+  TFile *output = new TFile("systHists.root","recreate");
+
+  hBgYield->SetName("hBgYield");  
+  hBgYield->Write();
+  hRelErrUp->SetName("hRelErrUp");
+  hRelErrUp->Write();
+  hRelErrLow->SetName("hRelErrLow");
+  hRelErrLow->Write();
+  hUpSystScaled->SetName("hOverallSystErrUp");
+  hUpSystScaled->Write();
+  hLowSystScaled->SetName("hOverallSystErrLow");
+  hLowSystScaled->Write();
+  gOverallSystErrUp->SetName("gOverallSystErrUp");
+  gOverallSystErrUp->Write();
+  gOverallSystErrLow->SetName("gOverallSystErrLow");
+  gOverallSystErrLow->Write(); 
+
+  output->Close(); 
+  delete output;  
+}
+//____________________________________________________________
+void AliHFEBeautySpectrum::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;
+  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);
+      ctest->cd(1);
+      TH1D* htest1=(TH1D*)result->Projection(0);
+      htest1->Draw();
+      ctest->cd(2);
+      TH1D* htest2=(TH1D*)result->Projection(1);
+      htest2->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");
+
+    file->Close();
+  }
+}
diff --git a/PWGHF/hfe/AliHFEBeautySpectrum.h b/PWGHF/hfe/AliHFEBeautySpectrum.h
new file mode 100644 (file)
index 0000000..7a208bb
--- /dev/null
@@ -0,0 +1,177 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFEBEAUTYSPECTRUM_H
+#define ALIHFEBEAUTYSPECTRUM_H
+
+#include "AliHFECorrectSpectrumBase.h"
+
+
+
+class TGraphErrors;
+class TObject;
+class TH1;
+class TF1;
+class TList;
+class TObjArray;
+class AliCFContainer;
+class AliHFEcontainer;
+class AliCFDataGrid;
+class AliCFEffGrid;
+class AliHFEInclusiveSpectrumQA;
+class AliHFEBeautySpectrumQA;
+
+class AliHFEBeautySpectrum : public AliHFECorrectSpectrumBase{
+ public:
+  
+  enum{
+    kElecBgSources = 9,
+    kBgLevels = 3,
+    kBgPtBins = 44,
+    kSignalPtBins = 35,
+    kCentrality = 12
+  };
+  
+  AliHFEBeautySpectrum(const char* name);
+  ~AliHFEBeautySpectrum();
+  
+  
+  virtual Bool_t Init(const AliHFEcontainer *datahfecontainer, const AliHFEcontainer *mchfecontainer, const AliHFEcontainer *bghfecontainer=0x0, const AliHFEcontainer */*v0hfecontainer*/ = 0x0,AliCFContainer */*photoniccontainerD*/ = 0x0);
+  virtual Bool_t Correct(Bool_t subtractcontamination=kTRUE, Bool_t /*subtractphotonic*/=kFALSE);
+  
+  AliCFDataGrid *SubtractBackground(Bool_t setBackground = kFALSE);
+  
+  AliCFDataGrid *CorrectParametrizedEfficiency(AliCFDataGrid* const bgsubpectrum = 0x0);
+  
+  THnSparse *Unfold(AliCFDataGrid* const bgsubpectrum = 0x0);
+  void UnfoldBG(AliCFDataGrid* const bgsubpectrum);
+  AliCFDataGrid *CorrectForEfficiency(AliCFDataGrid* const bgsubpectrum = 0x0);
+  
+  void SetPbPbAnalysis(Bool_t isPbPb = kFALSE) { fBeamType=(Char_t) isPbPb; };
+  void SetEtaSyst(Bool_t etaSyst = kTRUE) { fEtaSyst = etaSyst; };
+  
+  void SetParameterizedEff(AliCFContainer *container, AliCFContainer *containermb, AliCFContainer *containeresd, AliCFContainer *containeresdmb, Int_t *dimensions);
+  
+  void SetNumberOfMCEvents(Int_t nEvents) { fNMCEvents = nEvents; };
+  void SetNumberOfMC2Events(Int_t nEvents) { fNMCbgEvents = nEvents; };
+  void SetUnSetCorrelatedErrors(Bool_t unsetcorrelatederrors) {fUnSetCorrelatedErrors = unsetcorrelatederrors;};
+  void SetFillMoreCorrelationMatrix(Bool_t fillMoreCorrelationMatrix) { fFillMoreCorrelationMatrix = fillMoreCorrelationMatrix;}
+  
+  void SetNCentralityBinAtTheEnd(Int_t nCentralityBinAtTheEnd) {fNCentralityBinAtTheEnd = nCentralityBinAtTheEnd; };
+  void SetLowHighBoundaryCentralityBinAtTheEnd(Int_t low, Int_t high, Int_t i) { fLowBoundaryCentralityBinAtTheEnd[i] = low; fHighBoundaryCentralityBinAtTheEnd[i] = high;};
+  
+  void CallInputFileForBeauty2ndMethod();
+  void SetInputFileForBeauty2ndMethod(const char *filenameb = "BSpectrum2ndmethod.root"){fkBeauty2ndMethodfilename = filenameb; };
+  void SetBeautyAnalysis2ndMethod(Bool_t beauty2ndmethod) { fBeauty2ndMethod = beauty2ndmethod; }
+  void SetIPEffCombinedSamples(Bool_t ipEffCombinedSamples) { fIPEffCombinedSamples = ipEffCombinedSamples; }
+  void SetHadronEffbyIPcut(THnSparseF* hsHadronEffbyIPcut) { fHadronEffbyIPcut = hsHadronEffbyIPcut;};
+  void SetNonHFEsyst(Bool_t syst){ fNonHFEsyst = syst; };
+  
+  void SetDumpToFile(Bool_t dumpToFile) { fDumpToFile=dumpToFile; }; 
+  
+  void SetDebugLevel(Int_t debugLevel, Bool_t writeToFile = kFALSE) { fDebugLevel = debugLevel; fWriteToFile = writeToFile; };
+  void SetUnfoldBG() { fUnfoldBG = kTRUE; };
+  
+  
+  AliCFDataGrid* GetRawBspectra2ndMethod();
+  AliCFDataGrid* GetCharmBackground();
+  AliCFDataGrid* GetNonHFEBackground(Int_t decay = 0, Int_t source = 0);
+  THnSparse* GetCharmWeights(Int_t centBin);
+  THnSparse* GetBeautyIPEff(Bool_t isMCpt);
+  THnSparse* GetPIDxIPEff(Int_t source);
+  void CalculateNonHFEsyst();
+  
+  void EnableIPanaHadronBgSubtract() { fIPanaHadronBgSubtract = kTRUE; };
+  void EnableIPanaCharmBgSubtract() { fIPanaCharmBgSubtract = kTRUE; };
+  void EnableIPanaNonHFEBgSubtract() { fIPanaNonHFEBgSubtract = kTRUE; };
+  void EnableIPParameterizedEff() { fIPParameterizedEff = kTRUE; };
+  
+ protected:
+  void AddTemporaryObject(TObject *cont);
+  void ClearObject(TObject *o);
+  
+ private:
+  AliHFEBeautySpectrum(const AliHFEBeautySpectrum &ref);
+  AliHFEBeautySpectrum &operator=(const AliHFEBeautySpectrum &ref);
+  virtual void Copy(TObject &o) const;
+  
+  AliHFEBeautySpectrumQA *fQA; // Beauty QA
+  
+  //shift some of the following to QA or Base class:
+  TList *fTemporaryObjects;     // Emulate garbage collection
+  AliCFDataGrid *fBackground;   // Background Grid
+  TF1 *fEfficiencyTOFPIDD;       // TOF PID efficiency parameterized
+  TF1 *fEfficiencyesdTOFPIDD;    // TOF PID efficiency parameterized
+  TF1 *fEfficiencyIPCharmD;      // IP efficiency parameterized for charm
+  TF1 *fEfficiencyIPBeautyD;     // IP efficiency parameterized for beauty 
+  TF1 *fEfficiencyIPBeautyesdD;  // IP efficiency parameterized for beauty for esd
+  TF1 *fEfficiencyIPConversionD; // IP efficiency parameterized for conversion
+  TF1 *fEfficiencyIPNonhfeD;     // IP efficiency parameterized for nonhfe 
+  TF1 *fNonHFEbg;     // Efficiency Function
+  THnSparseF *fWeightCharm;     // Weight for charm bg
+  
+  AliCFContainer *fConvSourceContainer[kElecBgSources][kBgLevels]; //container for conversion electrons, divided into different photon sources
+  AliCFContainer *fNonHFESourceContainer[kElecBgSources][kBgLevels]; //container for non-HF electrons, divided into different sources
+  
+  Bool_t fInclusiveSpectrum;     // Inclusive Spectrum
+  Bool_t fDumpToFile;           // Write Result in a file
+  
+  Bool_t fUnSetCorrelatedErrors;    // Unset correlated errors
+  
+  Bool_t fIPanaHadronBgSubtract;     // Hadron background subtraction
+  Bool_t fIPanaCharmBgSubtract;      // Charm background subtraction 
+  Bool_t fIPanaNonHFEBgSubtract;     // nonHFE background subtraction
+  Bool_t fIPParameterizedEff;        // switch to use parameterized efficiency for ip analysis
+  Bool_t fNonHFEsyst;            // choose NonHFE background level (upper, lower, central)
+  Bool_t fBeauty2ndMethod;      // 2nd method to get beauty spectrum
+  Bool_t fIPEffCombinedSamples; // flag to combine two different samples
+  
+  Int_t fNMCEvents;             // Number of MC Events
+  Int_t fNMCbgEvents;       // Number of BG MC Events
+  
+  Int_t fNCentralityBinAtTheEnd;// Number of centrality class at the end
+  Int_t fLowBoundaryCentralityBinAtTheEnd[20];  // Boundary of the bins
+  Int_t fHighBoundaryCentralityBinAtTheEnd[20];  // Boundary of the bins
+  Bool_t fFillMoreCorrelationMatrix;             // For low stats to have reasonable errors
+  
+  THnSparseF *fHadronEffbyIPcut;// container for hadron efficiency by IP cut
+  TH1D *fEfficiencyCharmSigD; // charm IP cut eff from signal enhanced MC
+  TH1D *fEfficiencyBeautySigD; // beauty IP cut eff from signal enhanced MC
+  TH1D *fEfficiencyBeautySigesdD; // beauty IP cut eff from signal enhanced MC for esd
+  TH1D *fConversionEff;     // conversion IP cut eff
+  TH1D *fNonHFEEff;         // nonhfe IP cut eff
+  TH1D *fCharmEff;          // charm IP cut eff
+  TH1D *fBeautyEff;         // beauty IP cut eff
+  TH1D *fConversionEffbgc;      // conversion IP cut eff
+  TH1D *fNonHFEEffbgc;          // nonhfe IP cut eff
+  TH1D *fBSpectrum2ndMethod;             // beauty spectrum for 2nd method
+  const char *fkBeauty2ndMethodfilename;      // name of file, which contains beauty spectrum for 2ndmethod
+  Char_t fBeamType;             // beamtype; default -1; pp =0; PbPb=1
+  Bool_t fEtaSyst;              // pp 2.76 TeV (= kTRUE) or 7 TeV (= kFALSE)
+  
+  
+  Int_t fDebugLevel;            // Debug Level
+  Bool_t fWriteToFile;           // Write plots to eps files
+  Bool_t fUnfoldBG;             // flag to unfold backgroud
+  
+  ClassDef(AliHFEBeautySpectrum, 1) 
+};
+#endif
+
diff --git a/PWGHF/hfe/AliHFEBeautySpectrumQA.cxx b/PWGHF/hfe/AliHFEBeautySpectrumQA.cxx
new file mode 100644 (file)
index 0000000..edbae57
--- /dev/null
@@ -0,0 +1,578 @@
+
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// The following containers have to be set:
+//  - Correction framework container for real data
+//  - Correction framework container for MC (Efficiency Map)
+//  - Correction framework container for background coming from data
+//  - Correction framework container for background coming from MC
+//
+//  Author: 
+//            Raphaelle Bailhache <R.Bailhache@gsi.de>
+//
+
+#include <TArrayD.h>
+#include <TH1.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TObject.h>
+#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TAxis.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TPad.h>
+#include <TH2D.h>
+#include <TF1.h>
+
+#include "AliPID.h"
+#include "AliCFContainer.h"
+#include "AliCFDataGrid.h"
+#include "AliCFEffGrid.h"
+#include "AliCFGridSparse.h"
+#include "AliCFUnfolding.h"
+#include "AliLog.h"
+
+#include "AliHFEBeautySpectrumQA.h"
+#include "AliHFECorrectSpectrumBase.h"
+#include "AliHFEcuts.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEtools.h"
+
+ClassImp(AliHFEBeautySpectrumQA)
+
+const Char_t *AliHFEBeautySpectrumQA::fgkNameCanvas[AliHFEBeautySpectrumQA::kNTypeEfficiency] = {
+  "MCEfficiency",
+  "ParametrizedEfficiency"
+};
+
+//____________________________________________________________
+AliHFEBeautySpectrumQA::AliHFEBeautySpectrumQA():
+  TNamed(),
+  fPtMax(8.0),
+  fListOfResult(),
+  fWriteToFile(kTRUE)
+  {
+  //
+  // Default constructor
+  //
+
+  fListOfResult = new TObjArray(kNResults);
+  fListOfResult->SetName("ListOfResults");
+
+}
+//____________________________________________________________
+AliHFEBeautySpectrumQA::AliHFEBeautySpectrumQA(const char *name):
+  TNamed(name, ""),
+  fPtMax(8.0),
+  fListOfResult(),
+  fWriteToFile(kTRUE)
+  {
+  //
+  // Default constructor
+  //
+
+  fListOfResult = new TObjArray(kNResults);
+  fListOfResult->SetName("ListOfResults");
+
+}
+
+//____________________________________________________________
+AliHFEBeautySpectrumQA::~AliHFEBeautySpectrumQA(){
+  //
+  // Destructor
+  //
+  if(fListOfResult) delete fListOfResult;
+}
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::AddResultAt(TObject *obj,Int_t index)
+{
+  //
+  // Init what we need for the correction:
+  //
+
+  if(fListOfResult) fListOfResult->AddAt(obj,index);
+
+}
+//____________________________________________________________
+TObject *AliHFEBeautySpectrumQA::GetResult(Int_t index)
+{
+  //
+  // Get result
+  //
+
+  if(fListOfResult) return fListOfResult->UncheckedAt(index);
+  else return 0x0;
+
+}
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::DrawProjections() const
+{
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+  AliCFContainer *data = (AliCFContainer *) fListOfResult->UncheckedAt(kDataProjection);
+  THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
+  if(!data || !correlation) return;
+
+  Int_t ndimcont = data->GetNVar();
+  Int_t ndimcor = correlation->GetNdimensions();
+  Int_t charge = 3;
+  Int_t centrality = 5;
+  Int_t eta = 1;
+
+  TCanvas * canvas = new TCanvas("Projections","Projections",1000,700);
+  Int_t n = 0;
+  if(charge < ndimcont) n++;
+  if(centrality < ndimcont) n++;
+  if(eta < ndimcont) n++;
+  canvas->Divide(2,n);
+  Int_t counter = 1;
+
+  if(charge < ndimcont) {
+   
+    canvas->cd(counter);
+    TH1 *checkcharge = (TH1 *) data->Project(data->GetNStep()-1,charge);
+    checkcharge->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2F* projectioncharge = (TH2F *) correlation->Projection(charge,charge+((Int_t)(ndimcor/2.)));
+    projectioncharge->Draw("colz");
+    counter++;  
+
+  }
+
+  if(centrality < ndimcont) {
+    canvas->cd(counter);
+    TH1 *checkcentrality = (TH1 *) data->Project(data->GetNStep()-1,centrality);
+    checkcentrality->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2F *projectioncentrality = (TH2F *) correlation->Projection(centrality,centrality+((Int_t)(ndimcor/2.)));
+    projectioncentrality->Draw("colz");
+    counter++;  
+  }
+
+  if(eta < ndimcont) {
+    canvas->cd(counter);
+    TH1 *checketa = (TH1 *) data->Project(data->GetNStep()-1,eta);
+    checketa->Draw();
+    counter++;
+    canvas->cd(counter);
+    TH2D* projectioneta = (TH2D *) correlation->Projection(eta,eta+((Int_t)(ndimcor/2.)));
+    projectioneta->Draw("colz");
+  }
+
+
+}
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::DrawSubtractContamination() const
+{
+  //
+  // get spectrum for beauty 2nd method
+  //
+  //
+  TH1D *measuredTH1Daftersubstraction = (TH1D *) fListOfResult->UncheckedAt(kAfterSC);
+  TH1D *measuredTH1Dbeforesubstraction = (TH1D *) fListOfResult->UncheckedAt(kBeforeSC);
+  TH1D *measuredTH1background = (TH1D *) fListOfResult->UncheckedAt(kMeasBG);
+  if(!measuredTH1Daftersubstraction || !measuredTH1Dbeforesubstraction) return;
+
+  SetStyle();
+
+  TCanvas * cbackgroundsubtraction = new TCanvas("backgroundsubtraction","backgroundsubtraction",1000,700);
+  cbackgroundsubtraction->Divide(3,1);
+  cbackgroundsubtraction->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  measuredTH1Daftersubstraction->SetStats(0);
+  measuredTH1Daftersubstraction->SetTitle("");
+  measuredTH1Daftersubstraction->GetYaxis()->SetTitleOffset(1.5);
+  measuredTH1Daftersubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  measuredTH1Daftersubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  //measuredTH1Daftersubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  //measuredTH1Daftersubstraction->SetMarkerStyle(25);
+  //measuredTH1Daftersubstraction->SetMarkerColor(kBlack);
+  //measuredTH1Daftersubstraction->SetLineColor(kBlack);
+  measuredTH1Dbeforesubstraction->SetStats(0);
+  measuredTH1Dbeforesubstraction->SetTitle("");
+  measuredTH1Dbeforesubstraction->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  measuredTH1Dbeforesubstraction->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  //measuredTH1Dbeforesubstraction->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  //measuredTH1Dbeforesubstraction->SetMarkerStyle(24);
+  //measuredTH1Dbeforesubstraction->SetMarkerColor(kBlue);
+  //measuredTH1Dbeforesubstraction->SetLineColor(kBlue);
+  measuredTH1Daftersubstraction->Draw();
+  measuredTH1Dbeforesubstraction->Draw("same");
+  TLegend *legsubstraction = new TLegend(0.4,0.6,0.89,0.89);
+  legsubstraction->AddEntry(measuredTH1Dbeforesubstraction,"With hadron contamination","p");
+  legsubstraction->AddEntry(measuredTH1Daftersubstraction,"Without hadron contamination ","p");
+  legsubstraction->SetFillStyle(0);
+  legsubstraction->SetLineStyle(0);
+  legsubstraction->SetLineColor(0);
+  legsubstraction->Draw("same");
+  cbackgroundsubtraction->cd(2);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  TH1D* ratiomeasuredcontamination = (TH1D*)measuredTH1Dbeforesubstraction->Clone();
+  ratiomeasuredcontamination->SetName("ratiomeasuredcontamination");
+  ratiomeasuredcontamination->SetTitle("");
+  ratiomeasuredcontamination->GetYaxis()->SetTitleOffset(1.5);
+  ratiomeasuredcontamination->GetYaxis()->SetTitle("(with contamination - without contamination) / with contamination");
+  ratiomeasuredcontamination->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  ratiomeasuredcontamination->GetYaxis()->SetRangeUser(0.8,1.2);
+  ratiomeasuredcontamination->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  ratiomeasuredcontamination->Sumw2();
+  ratiomeasuredcontamination->Add(measuredTH1Daftersubstraction,-1.0);
+  ratiomeasuredcontamination->Divide(measuredTH1Dbeforesubstraction);
+  ratiomeasuredcontamination->SetStats(0);
+  ratiomeasuredcontamination->SetMarkerStyle(26);
+  ratiomeasuredcontamination->SetMarkerColor(kBlack);
+  ratiomeasuredcontamination->SetLineColor(kBlack);
+  for(Int_t k=0; k < ratiomeasuredcontamination->GetNbinsX(); k++){
+    ratiomeasuredcontamination->SetBinError(k+1,0.0);
+  }
+  ratiomeasuredcontamination->Draw("P");
+  cbackgroundsubtraction->cd(3);
+  measuredTH1background->SetStats(0);
+  measuredTH1background->SetTitle("");
+  measuredTH1background->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  measuredTH1background->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  measuredTH1background->SetMarkerStyle(26);
+  measuredTH1background->SetMarkerColor(kBlack);
+  measuredTH1background->SetLineColor(kBlack);
+  measuredTH1background->Draw();
+  if(fWriteToFile) cbackgroundsubtraction->SaveAs("BackgroundSubtracted.png");
+
+}
+
+
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::DrawCorrectWithEfficiency(Int_t typeeff) const
+{
+  //
+  // Correct the spectrum for efficiency and unfolding
+  // with both method and compare
+  //
+  
+  TH1D *afterE = 0x0;
+  TH1D *beforeE = 0x0;
+  TH1D *efficiencyDproj = 0x0;
+  TF1 *efficiencyparametrized = 0x0;
+  
+  if(typeeff== kMC) {
+    afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
+    beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforeMCE);
+    efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kMCEfficiency);
+  }
+  if(typeeff== kParametrized) {
+    afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterPE);
+    beforeE = (TH1D *) fListOfResult->UncheckedAt(kBeforePE);
+    efficiencyparametrized = (TF1 *) fListOfResult->UncheckedAt(kPEfficiency);
+  }
+  
+  if(typeeff==kMC && (!afterE || !beforeE || !efficiencyDproj)) return;
+  if(typeeff==kParametrized && (!afterE || !beforeE || !efficiencyparametrized)) return;
+  
+  SetStyle();
+  
+  TCanvas * cEfficiency = new TCanvas(AliHFEBeautySpectrumQA::fgkNameCanvas[typeeff],AliHFEBeautySpectrumQA::fgkNameCanvas[typeeff],1000,700);
+  cEfficiency->Divide(2,1);
+  cEfficiency->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  afterE->SetStats(0);
+  afterE->SetTitle("");
+  afterE->GetYaxis()->SetTitleOffset(1.5);
+  afterE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  afterE->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  afterE->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  afterE->SetMarkerStyle(25);
+  afterE->SetMarkerColor(kBlack);
+  afterE->SetLineColor(kBlack);
+  beforeE->SetStats(0);
+  beforeE->SetTitle("");
+  beforeE->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  beforeE->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  beforeE->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  beforeE->SetMarkerStyle(24);
+  beforeE->SetMarkerColor(kBlue);
+  beforeE->SetLineColor(kBlue);
+  afterE->Draw();
+  beforeE->Draw("same");
+  TLegend *legefficiency = new TLegend(0.4,0.6,0.89,0.89);
+  legefficiency->AddEntry(beforeE,"Before Efficiency correction","p");
+  legefficiency->AddEntry(afterE,"After Efficiency correction","p");
+  legefficiency->SetFillStyle(0);
+  legefficiency->SetLineStyle(0);
+  legefficiency->SetLineColor(0);
+  legefficiency->Draw("same");
+  cEfficiency->cd(2);
+  gPad->SetTicks();
+  if(typeeff==kMC) {
+    if(efficiencyDproj) {
+      efficiencyDproj->SetTitle("");
+      efficiencyDproj->SetStats(0);
+      efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
+      efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
+      efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
+      efficiencyDproj->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+      efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
+      efficiencyDproj->SetMarkerStyle(25);
+      efficiencyDproj->Draw();
+    }
+  }
+  if(typeeff==kParametrized) {
+    if(efficiencyparametrized) efficiencyparametrized->Draw();
+  }
+  
+  if(fWriteToFile) {
+    if(typeeff==kMC) cEfficiency->SaveAs("EfficiencyMC.png");
+    if(typeeff==kParametrized) cEfficiency->SaveAs("EfficiencyParametrized.png");
+  }
+
+}
+
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::DrawUnfolding() const
+{
+  //
+  // Draw unfolding
+  //
+  TH1D *measuredspectrumD = (TH1D *) fListOfResult->UncheckedAt(kBeforeU);
+  TH1D *residualspectrumD = (TH1D *) fListOfResult->UncheckedAt(kResidualU);
+  TH1D *efficiencyDproj = (TH1D *) fListOfResult->UncheckedAt(kUEfficiency);
+  THnSparseF *correlation = (THnSparseF *) fListOfResult->UncheckedAt(kCMProjection);
+  
+  if(!measuredspectrumD || !residualspectrumD || !efficiencyDproj || !correlation) return;
+  
+  Int_t ndimcor = (Int_t) correlation->GetNdimensions()/2.;
+
+  SetStyle();
+
+  TCanvas * cunfolding = new TCanvas("unfolding","unfolding",1000,700);
+  cunfolding->Divide(2,2);
+  cunfolding->cd(1);
+  gPad->SetLogy();
+  gPad->SetTicks();
+  residualspectrumD->GetYaxis()->SetTitle("dN/dp_{T} [(GeV/c)^{-1}]");
+  residualspectrumD->GetXaxis()->SetTitle("p^{rec}_{T} [GeV/c]");
+  residualspectrumD->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  residualspectrumD->SetStats(0);
+  residualspectrumD->SetTitle("");
+  residualspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  residualspectrumD->SetMarkerStyle(26);
+  residualspectrumD->SetMarkerColor(kBlue);
+  residualspectrumD->SetLineColor(kBlue);
+  residualspectrumD->Sumw2();
+  residualspectrumD->Draw("P");
+  measuredspectrumD->SetStats(0);
+  measuredspectrumD->SetTitle("");  
+  measuredspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  measuredspectrumD->SetMarkerStyle(25);
+  measuredspectrumD->SetMarkerColor(kBlack);
+  measuredspectrumD->SetLineColor(kBlack);
+  measuredspectrumD->Draw("same");
+  TLegend *legres = new TLegend(0.4,0.6,0.89,0.89);
+  legres->AddEntry(residualspectrumD,"Residual","p");
+  legres->AddEntry(measuredspectrumD,"Measured","p");
+  legres->SetFillStyle(0);
+  legres->SetLineStyle(0);
+  legres->SetLineColor(0);
+  legres->Draw("same");
+  cunfolding->cd(2);
+  gPad->SetTicks();
+  TH1D* ratioresidual = (TH1D*)residualspectrumD->Clone();
+  ratioresidual->SetName("ratioresidual");
+  ratioresidual->SetTitle("");
+  ratioresidual->GetYaxis()->SetRangeUser(0.6,1.4);
+  ratioresidual->GetYaxis()->SetTitle("Folded/Measured");
+  ratioresidual->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratioresidual->Divide(measuredspectrumD);
+  ratioresidual->SetStats(0);
+  ratioresidual->Draw();
+  cunfolding->cd(3);
+  gPad->SetTicks();
+  efficiencyDproj->GetYaxis()->SetTitle("Efficiency");
+  efficiencyDproj->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
+  efficiencyDproj->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  efficiencyDproj->GetYaxis()->SetRangeUser(0.0,1.0);
+  efficiencyDproj->SetStats(0);
+  efficiencyDproj->SetTitle("");
+  efficiencyDproj->GetYaxis()->SetTitleOffset(1.5);
+  efficiencyDproj->SetMarkerStyle(26);
+  efficiencyDproj->SetMarkerColor(kBlue);
+  efficiencyDproj->SetLineColor(kBlue);
+  efficiencyDproj->Sumw2();
+  efficiencyDproj->Draw("P");
+  cunfolding->cd(4);
+  TH2F *projectioncorr = (TH2F *) correlation->Projection(0,ndimcor);
+  projectioncorr->GetYaxis()->SetTitle("p^{ESD}_{T} [GeV/c]");
+  projectioncorr->GetXaxis()->SetTitle("p^{MC}_{T} [GeV/c]");
+  projectioncorr->GetXaxis()->SetRangeUser(0.0,fPtMax);
+  projectioncorr->GetYaxis()->SetRangeUser(0.0,fPtMax);
+  projectioncorr->SetStats(0);
+  projectioncorr->SetTitle("");
+  projectioncorr->Draw("colz");
+
+  if(fWriteToFile){
+    cunfolding->SaveAs("Unfolding.png");
+  }
+  
+}
+//____________________________________________________________
+void AliHFEBeautySpectrumQA::DrawResult()
+{
+  //
+  // Draw Results
+  //
+  TGraphErrors* correctedspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultUnfolded);
+  TGraphErrors* alltogetherspectrumD = (TGraphErrors *) fListOfResult->UncheckedAt(kFinalResultDirectEfficiency);
+  THnSparse* correctedspectrum = (THnSparse *) fListOfResult->UncheckedAt(kFinalResultUnfSparse);
+  AliCFDataGrid* alltogetherCorrection = (AliCFDataGrid *) fListOfResult->UncheckedAt( kFinalResultDirectEffSparse);
+  if(!correctedspectrumD || !alltogetherspectrumD) return;
+  
+  SetStyle();
+
+  TCanvas * ccorrected = new TCanvas("corrected","corrected",1000,700);
+  ccorrected->Divide(2,1);
+  ccorrected->cd(1);
+  gPad->SetLogy();
+  correctedspectrumD->SetTitle("");
+  correctedspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  correctedspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  correctedspectrumD->SetMarkerStyle(26);
+  correctedspectrumD->SetMarkerColor(kBlue);
+  correctedspectrumD->SetLineColor(kBlue);
+  correctedspectrumD->Draw("AP");
+  alltogetherspectrumD->SetTitle("");
+  alltogetherspectrumD->GetYaxis()->SetTitleOffset(1.5);
+  alltogetherspectrumD->GetYaxis()->SetRangeUser(0.000000001,1.0);
+  alltogetherspectrumD->SetMarkerStyle(25);
+  alltogetherspectrumD->SetMarkerColor(kBlack);
+  alltogetherspectrumD->SetLineColor(kBlack);
+  alltogetherspectrumD->Draw("P");
+  TLegend *legcorrected = new TLegend(0.4,0.6,0.89,0.89);
+  legcorrected->AddEntry(correctedspectrumD,"Corrected","p");
+  legcorrected->AddEntry(alltogetherspectrumD,"Alltogether","p");
+  legcorrected->Draw("same");
+  ccorrected->cd(2);
+  TH1D *correctedTH1D = correctedspectrum->Projection(0);
+  TH1D *alltogetherTH1D = (TH1D *) alltogetherCorrection->Project(0);
+  TH1D* ratiocorrected = (TH1D*)correctedTH1D->Clone();
+  ratiocorrected->SetName("ratiocorrected");
+  ratiocorrected->SetTitle("");
+  ratiocorrected->GetYaxis()->SetTitle("Unfolded/DirectCorrected");
+  ratiocorrected->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  ratiocorrected->Divide(correctedTH1D,alltogetherTH1D,1,1);
+  ratiocorrected->SetStats(0);
+  ratiocorrected->Draw();
+  if(fWriteToFile){ 
+    ccorrected->SaveAs("CorrectedBeauty.eps");    
+    TFile *out;
+    out = new TFile("finalSpectrum.root","recreate");
+    out->cd();
+    //
+    correctedspectrumD->SetName("UnfoldingCorrectedSpectrum");
+    correctedspectrumD->Write();
+    alltogetherspectrumD->SetName("AlltogetherSpectrum");
+    alltogetherspectrumD->Write();
+    ratiocorrected->SetName("RatioUnfoldingAlltogetherSpectrum");
+    ratiocorrected->Write();
+    //
+    out->Close(); 
+    delete out;
+  }
+}
+//____________________________________________________________
+TH1D *AliHFEBeautySpectrumQA::DivideSpectra(TGraphErrors *ga, TGraphErrors *gb) 
+{
+  //
+  // Divide Spectra
+  //
+
+  TH1D *afterE = (TH1D *) fListOfResult->UncheckedAt(kAfterMCE);
+  if(!afterE) return 0x0;
+
+  TH1D *histoB = (TH1D*) afterE->Clone();
+  histoB->Sumw2();
+  histoB->SetName("ratio");
+  TH1D *histoa = (TH1D*) afterE->Clone();
+  histoa->Sumw2();
+  histoa->SetName("a");
+  TH1D *histob = (TH1D*) afterE->Clone();
+  histob->Sumw2();
+  histob->SetName("b");
+  
+  double xa,ya,xb,yb,eya,eyb;
+  Int_t npointsa = ga->GetN();
+  Int_t npointsb = gb->GetN();
+  if(npointsa != npointsb) {
+    printf("Problem the two spectra have not the same number of points\n");
+    return 0x0;
+  }
+  for(Int_t k = 0; k < npointsa; k++){
+    ga->GetPoint(k,xa,ya);
+    gb->GetPoint(k,xb,yb);
+    //
+    Double_t centerhisto = histoa->GetBinCenter(k+1);
+    //
+    if((TMath::Abs(xa-xb) > 0.0001) || (TMath::Abs(xa-centerhisto) > 0.0001)) {
+      printf("Problem not the same x axis\n");
+      return 0x0;
+    }
+    histoa->SetBinContent(k+1,ya);
+    histob->SetBinContent(k+1,yb);
+    //
+    eya = ga->GetErrorY(k);
+    eyb = gb->GetErrorY(k);
+    //
+    histoa->SetBinError(k+1,eya);
+    histob->SetBinError(k+1,eyb);
+   
+  }
+  
+  histoB->Sumw2();
+  histoB->Divide(histoa,histob,1.0,1.0,"B");
+
+  return histoB;  
+
+}
+//__________________________________________
+void AliHFEBeautySpectrumQA::SetStyle() const
+{
+  //
+  // Set style
+  //
+
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(1111);
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetCanvasColor(10);
+  gStyle->SetPadLeftMargin(0.13);
+  gStyle->SetPadRightMargin(0.13);
+
+}
diff --git a/PWGHF/hfe/AliHFEBeautySpectrumQA.h b/PWGHF/hfe/AliHFEBeautySpectrumQA.h
new file mode 100644 (file)
index 0000000..7bcb734
--- /dev/null
@@ -0,0 +1,104 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+//
+// Class for spectrum correction
+// Subtraction of hadronic background, Unfolding of the data and
+// Renormalization done here
+// For more information see the implementation file
+//
+#ifndef ALIHFEBEAUTYSPECTRUMQA_H
+#define ALIHFEBEAUTYSPECTRUMQA_H
+
+#ifndef ROOT_TNamed
+#include <TNamed.h>
+#endif
+
+
+class TObjArray;
+class TGraphErrors;
+class TObject;
+
+
+class AliHFEBeautySpectrumQA : public TNamed{
+ public:
+
+  enum Results_t{
+    kDataProjection = 0,
+    kCMProjection = 1,
+    kBeforeSC = 2,
+    kAfterSC = 3,
+    kBeforePE = 4,
+    kAfterPE = 5,
+    kPEfficiency = 6,
+    kBeforeMCE = 7,
+    kAfterMCE = 8,
+    kMCEfficiency = 9,
+    kBeforeU = 10,
+    kAfterU = 11,
+    kResidualU = 12,
+    kUEfficiency = 13,
+    kFinalResultUnfolded = 14,
+    kFinalResultDirectEfficiency = 15,
+    kFinalResultUnfSparse = 16,
+    kFinalResultDirectEffSparse = 17,
+    kMeasBG = 18,
+    kNResults = 19
+  };
+
+
+  enum EfficiencyCorrection_t{
+    kMC = 0,
+    kParametrized = 1,
+    kNTypeEfficiency = 2
+  };
+   
+  void AddResultAt(TObject *obj,Int_t index);
+  TObject *GetResult(Int_t index);
+  
+  void DrawProjections() const;
+  void DrawSubtractContamination() const;
+  void DrawCorrectWithEfficiency(Int_t typeeff) const;
+  void DrawUnfolding() const;
+  void DrawResult();
+  
+  void SetStyle() const;
+  void SetWriteToFile(Bool_t writetofile) {fWriteToFile=writetofile; };
+  void SetPtMax(Double_t ptmax) {fPtMax = ptmax; };
+
+
+  TH1D *DivideSpectra(TGraphErrors *ga, TGraphErrors *gb);
+  
+
+  AliHFEBeautySpectrumQA();
+  AliHFEBeautySpectrumQA(const char* name);
+  ~AliHFEBeautySpectrumQA();
+  
+  protected:
+  
+  private:
+
+  static const Char_t* fgkNameCanvas[kNTypeEfficiency];     // Name of canvas
+
+    AliHFEBeautySpectrumQA(const AliHFEBeautySpectrumQA &);
+    AliHFEBeautySpectrumQA &operator=(const AliHFEBeautySpectrumQA &);
+    Double_t   fPtMax;             // Pt max to plot
+    TObjArray *fListOfResult;     // ListOfResults
+    Bool_t fWriteToFile;           // Write plots to eps files
+   
+    ClassDef(AliHFEBeautySpectrumQA, 1) 
+};
+#endif
+