/**************************************************************************
 * 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.                  *
 **************************************************************************/
//
//
// First implementation of a class
// to reject tagged electron coming from conversion, pi0 and eta
// by calculating the e+e- invariant mass 
// of the tagged electron with other tracks
// after looser cuts for the partner.
// PostProcess should extract the background yield
// If running with MC, it can be compared to the expected background yield 
//
// Authors:
//   Raphaelle Bailhache <rbailhache@ikf.uni-frankfurt.de > <R.Bailhache@gsi.de >
//
//

#include <THnSparse.h>
#include <TParticle.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TAxis.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TF1.h>
#include <TString.h>
#include <TLegend.h>
#include <TStyle.h>

#include <AliESDEvent.h>
#include <AliAODEvent.h>
#include <AliESDtrack.h>
#include <AliAODTrack.h>
#include "AliHFEelecbackground.h"
#include <AliMCEvent.h>
#include <AliStack.h>
#include <AliKFParticle.h>
#include "AliCFContainer.h"
#include "AliHFEpid.h"
#include "AliESDpid.h"
#include "AliITSPIDResponse.h"
#include "AliTPCPIDResponse.h"

ClassImp(AliHFEelecbackground)

Bool_t AliHFEelecbackground::fgUseMCPID = kFALSE;
const Double_t    AliHFEelecbackground::fgkMe = 0.00051099892;

//___________________________________________________________________________________________
AliHFEelecbackground::AliHFEelecbackground():
  fESD1(0x0)
  ,fAOD1(0x0)
  ,fMCEvent(0x0)
  ,fBz(0)
  ,fkVertex(0x0)
  ,fPtESD(0.0)
  ,fIndexTrack(0)
  ,fPdg(0)
  ,fLabMother(-1)
  ,fIsFrom(-1)
  ,fMotherGamma(-1)
  ,fMotherPi0(-1)
  ,fMotherC(-1)
  ,fMotherB(-1)
  ,fMotherEta(-1)
  ,fIsPartner(kFALSE)
  ,fIsSplittedTrack(kFALSE)
  ,fOpeningAngleCut(0.35)
  ,fInvMassCut(140.0)
  ,fChi2NdfCut(999999999.0)
  ,fUseAliKFCode(kTRUE)
  ,fSharedClusterCut(kFALSE)
  ,fRequireITSStandalone(0)
  ,fMinNbCls(2)
  ,fMinITSChi2(10.0)
  ,fMinNbClsSDDSPD(2)
  ,fPIDPartner(kFALSE)
  ,fPIDMethodPartner(0x0)
  ,fPIDMethodPartnerITS(0x0)
  ,fDebugLevel(0)
  ,fList(0x0)
  ,fListPostProcess(0x0)
{ 
  //
  // Default constructor
  //
  
}

//_______________________________________________________________________________________________
AliHFEelecbackground::AliHFEelecbackground(const AliHFEelecbackground &p):
  TObject(p)
  ,fESD1(0x0)
  ,fAOD1(0x0)
  ,fMCEvent(0x0)
  ,fBz(p.fBz)
  ,fkVertex(p.fkVertex)  
  ,fPtESD(p.fPtESD)
  ,fIndexTrack(0)
  ,fPdg(0)
  ,fLabMother(-1)
  ,fIsFrom(-1)
  ,fMotherGamma(-1)
  ,fMotherPi0(-1)
  ,fMotherC(-1)
  ,fMotherB(-1)
  ,fMotherEta(-1)
  ,fIsPartner(kFALSE)
  ,fIsSplittedTrack(kFALSE)
  ,fOpeningAngleCut(0.35)
  ,fInvMassCut(140.0)
  ,fChi2NdfCut(999999999.0)
  ,fUseAliKFCode(kTRUE)
  ,fSharedClusterCut(kFALSE)
  ,fRequireITSStandalone(0)
  ,fMinNbCls(2)
  ,fMinITSChi2(10.0)
  ,fMinNbClsSDDSPD(2)
  ,fPIDPartner(kFALSE)
  ,fPIDMethodPartner(0x0)
  ,fPIDMethodPartnerITS(0x0)
  ,fDebugLevel(0)
  ,fList(0x0)  
  ,fListPostProcess(0x0)
{ 
  //
  // Copy constructor
  //
}

//_______________________________________________________________________________________________
AliHFEelecbackground&
AliHFEelecbackground::operator=(const AliHFEelecbackground &)
{
  //
  // Assignment operator
  //

  AliInfo("Not yet implemented.");
  return *this;
}

//_______________________________________________________________________________________________
AliHFEelecbackground::~AliHFEelecbackground()
{
  //
  // Destructor
  //
  if(fPIDMethodPartner) delete fPIDMethodPartner;
  if(fPIDMethodPartnerITS) delete fPIDMethodPartnerITS;

  if(fList){
    fList->Clear();
    delete fList;
  }

  if(fListPostProcess){
    fListPostProcess->Clear();
    delete fListPostProcess;
  }
}
//___________________________________________________________________________________________
Bool_t AliHFEelecbackground::Load(const Char_t * filename)
{
  //
  // Generic container loader
  //

  if(!TFile::Open(filename)){
    return kFALSE;
  }
  TList *o = 0x0;
  if(!(o = (TList*)gFile->Get("Results"))){
    return kFALSE;
  }
  TList *oe = 0x0;
  if(!(oe = (TList*)dynamic_cast<TList *>(o->FindObject("HFEelecbackground")))){
    return kFALSE;
  }
  fList = (TList*)oe->Clone("HFEelecbackground");
  gFile->Close();
  return kTRUE;
}
//___________________________________________________________________________________________
Bool_t AliHFEelecbackground::Load(TList * const outputlist)
{
  //
  // Generic container loader
  //
  if(!outputlist) return kFALSE;
  else   fList = (TList*)outputlist->Clone("HFEelecbackground");
  return kTRUE;
}
//_______________________________________________________________________________________________
void AliHFEelecbackground::Reset()
{
  //
  // Reset variables
  //
  fPtESD = 0.0;
  fIndexTrack = -1;
  fPdg = -1;
  fLabMother = -1;
  fIsFrom = -1;
  fMotherGamma = -1;
  fMotherPi0 = -1;
  fMotherC = -1;
  fMotherB = -1;
  fMotherEta = -1;
  fIsPartner = kFALSE;
  fIsSplittedTrack = kFALSE;
  for(Int_t id = 0; id < 10; id++) {
    fCuts[id] = kFALSE;
  }
 
}
//_______________________________________________________________________________________________
void AliHFEelecbackground::CreateHistograms(TList * const qaList)
{ 
  //
  // create histograms
  //
  if(!qaList) return;

  fList = qaList;
  fList->SetName("HFEelecbackground");  

  //////////
  // bins
  /////////

  Int_t nBinsPt = 25;
  Double_t minPt = 0.01;
  Double_t maxPt = 10.0;
  
  Int_t nBinsPtMore = 100;
  Double_t minPtMore = 0.01;
  Double_t maxPtMore = 10.0;
  
  Int_t nBinsInv = 50;
  Double_t minInv = 0.0;
  Double_t maxInv = 0.2;
  
  Int_t nBinsOp = 50;
  Double_t minOp = 0.0;
  Double_t maxOp = 2;

  Int_t nBinsCh = 4;
  Double_t minCh = 0.0;
  Double_t maxCh = 4.0;
  
  Double_t *binLimLogPt = new Double_t[nBinsPt+1];
  Double_t *binLimPt    = new Double_t[nBinsPt+1];
  for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
  for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);

  Double_t *binLimLogPtMore = new Double_t[nBinsPtMore+1];
  Double_t *binLimPtMore    = new Double_t[nBinsPtMore+1];
  for(Int_t i=0; i<=nBinsPtMore; i++) binLimLogPtMore[i]=(Double_t)TMath::Log10(minPtMore) + (TMath::Log10(maxPtMore)-TMath::Log10(minPtMore))/nBinsPtMore*(Double_t)i ;
  for(Int_t i=0; i<=nBinsPtMore; i++) binLimPtMore[i]=(Double_t)TMath::Power(10,binLimLogPtMore[i]);

  Double_t *binLimInv    = new Double_t[nBinsInv+1];
  for(Int_t i=0; i<=nBinsInv; i++) binLimInv[i]=(Double_t)minInv  + (maxInv-minInv)  /nBinsInv*(Double_t)i ;
  
  Double_t *binLimOp    = new Double_t[nBinsOp+1];
  for(Int_t i=0; i<=nBinsOp; i++) binLimOp[i]=(Double_t)minOp  + (maxOp-minOp) /nBinsOp*(Double_t)i ;
  
  Double_t *binLimCh    = new Double_t[nBinsCh+1];
  for(Int_t i=0; i<=nBinsCh; i++) binLimCh[i]=(Double_t)minCh  + (maxCh-minCh) /nBinsCh*(Double_t)i ;
  
  const Int_t nvarData = 5;
  // pt reconstructed tagged e
  // pt reconstructed mother
  // opening angle
  // invariant mass 
  // Data: charge (0-opposite sign, 1-like sign ++, 2-like sign --, 3-rotated tracks)

  Int_t iBinData[nvarData];
  iBinData[0]=nBinsPt;
  iBinData[1]=nBinsPt;
  iBinData[2]=nBinsOp;
  iBinData[3]=nBinsInv;
  iBinData[4]=nBinsCh;
  
  //
  // Opening angle and invariant mass
  //
  
  THnSparseF *hsSparseData = new THnSparseF("OpeningangleinvmassData","",nvarData,iBinData);
  hsSparseData->SetBinEdges(0,binLimPt);
  hsSparseData->SetBinEdges(1,binLimPt);
  hsSparseData->SetBinEdges(2,binLimOp);
  hsSparseData->SetBinEdges(3,binLimInv);
  hsSparseData->SetBinEdges(4,binLimCh);
  hsSparseData->Sumw2();

  fList->AddAt(hsSparseData,kDatai);

  //
  // Radius, DCA and Chi2Ndf
  //

  TH1F *dataRadiusHisto = new TH1F("DataRadius","", 200, 0.0, 200.0); // recontructed radius from the AliKF of the e+e- pair
  fList->AddAt(dataRadiusHisto,kDatar);

  TH1F *dataDcaHisto = new TH1F("DataDCA","", 100, 0.0, 6.0); // dca distribution
  fList->AddAt(dataDcaHisto,kDatadca); 
  
  TH1F *dataChi2NdfHisto = new TH1F("DataChi2Ndf","", 100, 0.0, 5.0); // chi2Ndf distribution    
  fList->AddAt(dataChi2NdfHisto,kDatachi2Ndf); 
  

  if(HasMCData()) {

    //
    // Opening angle and invariant mass with MC infos
    //

    const Int_t nvarMCo = 6;
    // pt reconstructed tagged e
    // pt reconstructed mother
    // opening angle
    // invariant mass 
    // MC: 0-FromBackground, 1-FromGamma, 2-FromPi0, 3-FromEta, 4-FromC, 5-FromB
    // MCSplitted: 0-not, 1-splittedOs, 2-ksplittedSs
    

    Int_t nBinsMCOrigin = 6;
    Double_t minMCOrigin = 0.0;
    Double_t maxMCOrigin = 6.0;
    
    Double_t *binLimMCOrigin = new Double_t[nBinsMCOrigin+1];
    for(Int_t i=0; i<=nBinsMCOrigin; i++) binLimMCOrigin[i]=(Double_t)minMCOrigin  + (maxMCOrigin-minMCOrigin) /nBinsMCOrigin*(Double_t)i ;

    Int_t nBinsMCSplitted = 3;
    Double_t minMCSplitted = 0.0;
    Double_t maxMCSplitted = 3.0;
    
    Double_t *binLimMCSplitted = new Double_t[nBinsMCSplitted+1];
    for(Int_t i=0; i<=nBinsMCSplitted; i++) binLimMCSplitted[i]=(Double_t)minMCSplitted  + (maxMCSplitted-minMCSplitted) /nBinsMCSplitted*(Double_t)i ;
    
    Int_t iBinMCo[nvarMCo];
    iBinMCo[0]=nBinsPt;
    iBinMCo[1]=nBinsPt;
    iBinMCo[2]=nBinsOp;
    iBinMCo[3]=nBinsInv;
    iBinMCo[4]=nBinsMCOrigin;
    iBinMCo[5]=nBinsMCSplitted;
        
    THnSparseF *hsSparseMCo = new THnSparseF("OpeningangleinvmassMC","",nvarMCo,iBinMCo);
    hsSparseMCo->SetBinEdges(0,binLimPt);
    hsSparseMCo->SetBinEdges(1,binLimPt);
    hsSparseMCo->SetBinEdges(2,binLimOp);
    hsSparseMCo->SetBinEdges(3,binLimInv);
    hsSparseMCo->SetBinEdges(4,binLimMCOrigin);
    hsSparseMCo->SetBinEdges(5,binLimMCSplitted);
    hsSparseMCo->Sumw2();

    fList->AddAt(hsSparseMCo,kMCo);

    //
    // Radius, DCA and Chi2Ndf with MC info
    //

    TH2F *mcRadiusHisto = new TH2F("MCRadius","", 200, 0.0, 200.0,6,-0.5,5.5); // recontructed radius from the AliKF of the e+e- pair
    fList->AddAt(mcRadiusHisto,kMCr);
    
    TH2F *mcDcaHisto = new TH2F("MCDCA","", 100, 0.0, 6.0,6,-0.5,5.5); // dca distribution
    fList->AddAt(mcDcaHisto,kMCdca); 
    
    TH2F *mcChi2NdfHisto = new TH2F("MCChi2Ndf","", 100, 0.0, 5.0,6,-0.5,5.5); // chi2Ndf distribution    
    fList->AddAt(mcChi2NdfHisto,kMCchi2Ndf); 

    //////////////////////////////////////////////////////////
    // if fDebugLevel 1: Rejection efficiencies of the cuts
    //////////////////////////////////////////////////////////

    if(fDebugLevel > 0) {

      if(HasMCData()) {
	
	const Int_t nvarMCe = 3;
	// pt reconstructed tagged e
	// cut passed: 0-all, 1-Partner tracked, 2-Opposite-sign, 3-SingleTrackCutPart, 4-ShareCluster, 5-PID, 6-DCA, 7-chi2Ndf AliKF, 8-Openingangle, 9-Invmass
	// MC: 0-FromBackground, 1-FromGamma, 2-FromPi0, 3-FromEta, 4-FromC, 5-FromB
      
      Int_t nBinsMCCutPassed = 10;
      Double_t minMCCutPassed = -0.5;
      Double_t maxMCCutPassed = 9.5;
      
      Double_t *binLimMCCutPassed = new Double_t[nBinsMCCutPassed+1];
      for(Int_t i=0; i<=nBinsMCCutPassed; i++) binLimMCCutPassed[i]=(Double_t)minMCCutPassed  + (maxMCCutPassed-minMCCutPassed) /nBinsMCCutPassed*(Double_t)i ;
      
      Int_t iBinMCe[nvarMCe];
      iBinMCe[0]=nBinsPt;
      iBinMCe[1]=nBinsMCCutPassed;
      iBinMCe[2]=nBinsMCOrigin;
      
      THnSparseF *hsSparseMCe = new THnSparseF("CutPassedMC","",nvarMCe,iBinMCe);
      hsSparseMCe->SetBinEdges(0,binLimPt);
      hsSparseMCe->SetBinEdges(1,binLimMCCutPassed);
      hsSparseMCe->SetBinEdges(2,binLimMCOrigin);
      hsSparseMCe->Sumw2();
      
      fList->AddAt(hsSparseMCe,kMCe); 
      
      }
    }

    /////////////////////////////////////////////////////////////////
    // if fDebugLevel 1: PIDPartCut and ShareClusters
    /////////////////////////////////////////////////////////////////

    if(fDebugLevel > 1) {

      if(!fRequireITSStandalone){
	
	//
	// TPC 
	//

	TH2F *tpcPartner0 = new TH2F("TPCPartner0","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0); 
	fList->AddAt(tpcPartner0,kMCcutPart0); 
	TH2F *tpcPartner1 = new TH2F("TPCPartner1","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0); 
	fList->AddAt(tpcPartner1,kMCcutPart1); 
      
      }
      else {

	//
	// ITS
	//

	TH2F *itsPartner0 = new TH2F("ITSPartner0","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0); 
	fList->AddAt(itsPartner0,kMCcutPart0); 
	TH2F *itsPartner1 = new TH2F("ITSPartner1","", nBinsPtMore, binLimPtMore, 200, 0.0, 700.0); 
	fList->AddAt(itsPartner1,kMCcutPart1); 

	/////////////////////////////////////////////////////
       	// dEdx of the four layers for the track partner
	/////////////////////////////////////////////////////
	const Int_t nvarITSsignal = 5;
	
	Int_t nBinsITSsignal  = 100;
	Double_t minITSsignal = 0.0;
	Double_t maxITSsignal = 350.0;
	
	Double_t *binLimITSsignal = new Double_t[nBinsITSsignal+1];
	for(Int_t i=0; i<=nBinsITSsignal; i++) binLimITSsignal[i]=(Double_t)minITSsignal  + (maxITSsignal-minITSsignal) /nBinsITSsignal*(Double_t)i ;
	
	Int_t iBinITSsignal[nvarITSsignal];
	iBinITSsignal[0]=nBinsPt;
	iBinITSsignal[1]=nBinsITSsignal;
	iBinITSsignal[2]=nBinsITSsignal;
	iBinITSsignal[3]=nBinsITSsignal;
	iBinITSsignal[4]=nBinsITSsignal;
	
	THnSparseF *hsSparseITSpid = new THnSparseF("SparseITSsignal","",nvarITSsignal,iBinITSsignal);
	hsSparseITSpid->SetBinEdges(0,binLimPt);
	hsSparseITSpid->SetBinEdges(1,binLimITSsignal);
	hsSparseITSpid->SetBinEdges(2,binLimITSsignal);
	hsSparseITSpid->SetBinEdges(3,binLimITSsignal);
	hsSparseITSpid->SetBinEdges(4,binLimITSsignal);
	hsSparseITSpid->Sumw2();
	
	fList->AddAt(hsSparseITSpid,kMCcutPart2); 

	///////////////////////////////////////////////////////////////////////////////////////
 	// dEdx of the four layers for the track partner and track to reject splitted track
	///////////////////////////////////////////////////////////////////////////////////////
	const Int_t nvarITSsignalSplit = 5;

	Int_t nBinsITSSplit  = 2;
	Double_t minITSSplit = 0.0;
	Double_t maxITSSplit = 2.0;
	
	Double_t *binLimITSSplit = new Double_t[nBinsITSSplit+1];
	for(Int_t i=0; i<=nBinsITSSplit; i++) binLimITSSplit[i]=(Double_t)minITSSplit  + (maxITSSplit-minITSSplit) /nBinsITSSplit*(Double_t)i ;


	Int_t nBinsITSsignalSplit  = 50;
	Double_t minITSsignalSplit = -25.0;
	Double_t maxITSsignalSplit = 25.0;
	
	Double_t *binLimITSsignalSplit = new Double_t[nBinsITSsignalSplit+1];
	for(Int_t i=0; i<=nBinsITSsignalSplit; i++) binLimITSsignalSplit[i]=(Double_t)minITSsignalSplit  + (maxITSsignalSplit-minITSsignalSplit) /nBinsITSsignalSplit*(Double_t)i ;
	
	Int_t iBinITSsignalSplit[nvarITSsignalSplit];
	iBinITSsignalSplit[0]=nBinsITSSplit;
	for(Int_t k = 1; k < 5; k++){
	  iBinITSsignalSplit[k]=nBinsITSsignalSplit;
	}
	
	THnSparseF *hsSparseITSpidSplit = new THnSparseF("SparseITSsignalSplit","",nvarITSsignalSplit,iBinITSsignalSplit);
	hsSparseITSpidSplit->SetBinEdges(0,binLimITSSplit);
	for(Int_t k = 1; k < 5; k++) {
	  hsSparseITSpidSplit->SetBinEdges(k,binLimITSsignalSplit);
	}
	hsSparseITSpidSplit->Sumw2();
	
	fList->AddAt(hsSparseITSpidSplit,kMCcutPart3); 
	
      }

    }
    
  }

  //qaList->Add(fList);

}
//_______________________________________________________________________________________________
void AliHFEelecbackground::PairAnalysis(AliESDtrack* const track, AliESDtrack* const trackPart)
{
  //
  // calculate (tagged e-partner) dca, opening angle, invariant mass 
  //

  /////////////////////
  // pt tagged
  //////////////////////
  TVector3 v3Dtagged;
  Double_t *pxyz = new Double_t[3];
  track->PxPyPz(pxyz);
  v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
  fPtESD = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]); 
  
  ////////////////////////
  // Take label
  ////////////////////////
  Int_t indexTrack = TMath::Abs(track->GetLabel());
  Int_t indexTrackPart = TMath::Abs(trackPart->GetLabel());
  
  /////////////////////////
  // If MC data
  ////////////////////////
  
  if(HasMCData()) {
    
    // Take info track if not already done 
    if(indexTrack!= fIndexTrack) {

      for(Int_t id = 0; id < 10; id++) {
	fCuts[id] = kFALSE;
      }
      
      fIsFrom = kElectronFromBackground;
          
      fPdg = GetPdg(indexTrack); 
      fLabMother = GetLabMother(indexTrack);
      
      fMotherGamma = IsMotherGamma(indexTrack);
      if((fMotherGamma != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromGamma;
      fMotherPi0 = IsMotherPi0(indexTrack);
      if((fMotherPi0 != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromPi0;
      fMotherC = IsMotherC(indexTrack);
      if((fMotherC != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromC;
      fMotherB = IsMotherB(indexTrack);
      if((fMotherB != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromB;
      fMotherEta = IsMotherEta(indexTrack);
      if((fMotherEta != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromEta;
      
      fIndexTrack = indexTrack;
      
    }

    // MC PID for tagged
    if(fgUseMCPID) {
      if(TMath::Abs(fPdg) != 11) return;
    }
    
    // Look at trackPart
    fIsPartner = kFALSE;
    Int_t pdgPart = GetPdg(indexTrackPart);
    if(TMath::Abs(pdgPart) == 11) {
      Int_t labMotherPart = GetLabMother(indexTrackPart);
      if((labMotherPart == fLabMother) && (indexTrack != indexTrackPart) && (TMath::Abs(fPdg) == 11) && (fPdg*pdgPart < 0) && (fLabMother >=0) && (fLabMother < (((AliStack *)fMCEvent->Stack())->GetNtrack()))) fIsPartner = kTRUE;
      // special case of c and b
      Int_t motherCPart = IsMotherC(indexTrackPart);
      if((motherCPart != -1) && (fIsFrom == kElectronFromC) && (fPdg*pdgPart < 0)){
	fIsPartner = kTRUE;	
      }
      Int_t motherBPart = IsMotherB(indexTrackPart);
      if((motherBPart != -1) && (fIsFrom == kElectronFromB) && (fPdg*pdgPart < 0)){
	fIsPartner = kTRUE;	
      }
    }

    // Look at splitted tracks
    fIsSplittedTrack = kFALSE;
    if(indexTrackPart == fIndexTrack) fIsSplittedTrack = kTRUE;
    
  }

  //////////////////////
  // Sign
  /////////////////////
  Int_t sign = -1;
  if((track->Charge() > 0.0) && (trackPart->Charge() > 0.0)) sign = kPp; 
  if((track->Charge() < 0.0) && (trackPart->Charge() < 0.0)) sign = kNn; 
  if(((track->Charge() > 0.0) && (trackPart->Charge() < 0.0)) || ((track->Charge() < 0.0) && (trackPart->Charge() > 0.0))) sign = kOs; 
  
  /////////////////////////
  // Cut effects
  ////////////////////////   
  Double_t cuteffect[3];

  if(fDebugLevel > 0) {  
    if(HasMCData()) {
      cuteffect[0] = fPtESD;
      cuteffect[1] = 0.0;
      cuteffect[2] = fIsFrom;
      if(!fCuts[0]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[0] = kTRUE;
      }
    }
  }


  ///////////////////////////////
  // Cut effect: Partner track 
  ///////////////////////////////   
  
  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner) {
      cuteffect[1] = 1.0;
      if(!fCuts[1]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[1] = kTRUE;
      }
    }
  }

  ///////////////////////////////
  // Cut effect: Opposite sign 
  ///////////////////////////////   
  
  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign == kOs)) {
      cuteffect[1] = 2.0;
      if(!fCuts[2]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[2] = kTRUE;
      }
    }
  }

  ////////////////////////
  // Partner track cut
  ////////////////////////
  if(!SingleTrackCut(trackPart)) return;

  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign==kOs)) {
      cuteffect[1] = 3.0;
      if(!fCuts[3]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[3] = kTRUE;
      }
    }
  }
  
  /////////////////////////
  // shared clusters cut
  /////////////////////////
  if(fSharedClusterCut && ShareCluster(track,trackPart)) return;

  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign==kOs)) {
      cuteffect[1] = 4.0;
      if(!fCuts[4]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[4] = kTRUE;
      }
    } 
  }

  ////////////////////////
  // PID Partner track 
  ////////////////////////
  if(!PIDTrackCut(trackPart)) return;

  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign==kOs)) {
      cuteffect[1] = 5.0;
      if(!fCuts[5]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[5] = kTRUE;
      }
    }
  }

  
  //////////////////////
  // DCA
  /////////////////////
  
  Double_t xthis,xp;
  Double_t dca = track->GetDCA(trackPart,fBz,xthis,xp);
  (dynamic_cast<TH1F *>(fList->At(kDatadca)))->Fill(dca);
  if(HasMCData()) {
    //printf("has MC data for DCA\n");
    //printf("IsPartner %d and isfrom %d\n",fIsPartner,fIsFrom);
    if(fIsFrom==kElectronFromBackground) (dynamic_cast<TH2F *>(fList->At(kMCdca)))->Fill(dca,fIsFrom);
    else {
      if(fIsPartner) (dynamic_cast<TH2F *>(fList->At(kMCdca)))->Fill(dca,fIsFrom);
    }
  }
  if(TMath::Abs(dca) > 3.0) return;

  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign==kOs)) {
      cuteffect[1] = 6.0;
      if(!fCuts[6]) {
	(dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	fCuts[6] = kTRUE;
      }
    }
  }

  ///////////////////////////////////
  // Calcul mother variables
  ///////////////////////////////////
  Double_t results[5];
  Double_t resultsr[5];

  
  if(!fUseAliKFCode) {
    
    /////////////////////////////
    // Propagate
    ////////////////////////////
    
    Double_t norradius = TMath::Sqrt(fkVertex->GetX()*fkVertex->GetX()+fkVertex->GetY()*fkVertex->GetY());
    
    AliESDtrack *trackCopy = new AliESDtrack(*track);
    AliESDtrack *trackPartCopy = new AliESDtrack(*trackPart);
    Bool_t propagateok = kTRUE;
    if((!(trackPartCopy->PropagateTo(norradius,fBz))) || (!(trackCopy->PropagateTo(norradius,fBz)))) propagateok = kFALSE;
    if(!propagateok) {
      if(trackCopy) delete trackCopy;
      if(trackPartCopy) delete trackPartCopy;
      return;
    }  
  
    CalculateMotherVariable(trackCopy,trackPartCopy,&results[0]);
    CalculateMotherVariableR(trackCopy,trackPartCopy,&resultsr[0]);
    
    if(trackCopy) delete trackCopy;
    if(trackPartCopy) delete trackPartCopy;
    
  }
  else {
    
    if(!CalculateMotherVariable(track,trackPart,&results[0])) return;
    if(fDebugLevel > 0) {     
      if(HasMCData() && fIsPartner && (sign==kOs)) {
	cuteffect[1] = 7.0;
	if(!fCuts[7]) {
	  (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	  fCuts[7] = kTRUE;
	}
      }
    }
  
  }
  
  /////////////////////////////////////
  // Fill
  /////////////////////////////////////
   
  FillOutput(results, resultsr, sign);

  if(fDebugLevel > 0) {  
    if(HasMCData() && fIsPartner && (sign==kOs)) {
      
      if(TMath::Abs(results[4]) < fOpeningAngleCut) {
	
	cuteffect[1] = 8.0;
	if(!fCuts[8]) {
	  (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	  fCuts[8] = kTRUE;
	}
	if(TMath::Abs(results[1]) < fInvMassCut) {
	  cuteffect[1] = 9.0;
	  if(!fCuts[9]) {
	    (dynamic_cast<THnSparseF *>(fList->At(kMCe)))->Fill(cuteffect);
	    fCuts[9] = kTRUE;
	  }
	}
      }
    }    
  }
  
 
}
//_____________________________________________________________________________________
Bool_t AliHFEelecbackground::CalculateMotherVariable(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
{
  //
  // variables mother and take the pt of the track
  //
  // results contain: ptmother, invmass, etamother, phimother, openingangle
  //
  // with a chi2Ndf cut for AliKF code
  //
  
  if(!fUseAliKFCode) {
    
    TVector3 v3Dtagged;
    TVector3 v3Dpart;
    
    Double_t *pxyz = new Double_t[3];
    track->PxPyPz(pxyz);
    v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
    
    Double_t *pxyzpart = new Double_t[3];
    trackpart->PxPyPz(pxyzpart);
    v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
    
    
    TVector3 motherrec = v3Dtagged + v3Dpart;
    
    Double_t etaESDmother = motherrec.Eta();
    Double_t ptESDmother  = motherrec.Pt();
    Double_t phiESDmother = motherrec.Phi();
    if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
    
    
    // openinganglepropagated
    Double_t openingangle = v3Dtagged.Angle(v3Dpart);
    
    // invmass
    Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
    Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
    
    // e propagate
    Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
    Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
    
    Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
    
    if(!results) return kFALSE;

    results[0] = ptESDmother;
    results[1] = etaESDmother;
    results[2] = phiESDmother;
    results[3] = invmass;
    results[4] = openingangle;
    
    return kTRUE;

  }
  else {
    
    AliKFParticle pair;
    pair.Initialize();
    
    // pid
    Int_t pid1 = -11;
    if(track->Charge() > 0.0) pid1 = 11;
    Int_t pid2 = -11;
    if(trackpart->Charge() > 0.0) pid2 = 11;
    
    
    // daughters
    AliKFParticle kf(*track,pid1);
    AliKFParticle kfpart(*trackpart,pid2);
    
    pair.AddDaughter(kf);
    pair.AddDaughter(kfpart);
    
    // variables
    Double_t openingangle = kf.GetAngle(kfpart);
    Double_t chi2ndf = pair.GetChi2()/pair.GetNDF();
    //Double_t decayLength = pair.GetDecayLength();
    Double_t radius = pair.GetR();
    //Double_t masserror = pair.GetErrMass()>0?pair.GetErrMass()/pair.GetMass():1000000;
    Double_t ptpair = pair.GetPt();
    Double_t etapair = pair.GetEta();
    Double_t phipair = pair.GetPhi();
    Double_t masspair = pair.GetMass();
    
    // Put them
    if(!results) return kFALSE;

    results[0] = ptpair;
    results[1] = etapair;
    results[2] = phipair;
    results[3] = masspair;
    results[4] = openingangle;

    // chi2Ndf cut
    (dynamic_cast<TH1F *>(fList->At(kDatachi2Ndf)))->Fill(chi2ndf);
    if(HasMCData()){
      if(fIsFrom==kElectronFromBackground) (dynamic_cast<TH2F *>(fList->At(kMCchi2Ndf)))->Fill(chi2ndf,fIsFrom);
      else {
	if(fIsPartner) (dynamic_cast<TH2F *>(fList->At(kMCchi2Ndf)))->Fill(chi2ndf,fIsFrom);
      }
    }
    if(chi2ndf > fChi2NdfCut) return kFALSE;
    else {
      (dynamic_cast<TH1F *>(fList->At(kDatar)))->Fill(radius);
      if(HasMCData()) {
	if(fIsFrom==kElectronFromBackground) (dynamic_cast<TH2F *>(fList->At(kMCr)))->Fill(radius,fIsFrom);
	else {
	  if(fIsPartner) (dynamic_cast<TH2F *>(fList->At(kMCr)))->Fill(radius,fIsFrom);
	}
      }
      return kTRUE;
    }
    
  }
  
}
//_____________________________________________________________________________________
void AliHFEelecbackground::CalculateMotherVariableR(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
{
  //
  // variables mother
  //
  // results contain: ptmother, invmass, etamother, phimother, openingangle
  // Implemented only for no AliKF
  //

  if(!fUseAliKFCode) {
    
    TVector3 v3Dtagged;
    TVector3 v3Dpart;
    
    Double_t *pxyz = new Double_t[3];
    track->PxPyPz(pxyz);
    v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
    Double_t *pxyzpart = new Double_t[3];
    trackpart->PxPyPz(pxyzpart);
    v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
    
    // rotate the partner
    v3Dpart.RotateZ(TMath::Pi());
    v3Dpart.GetXYZ(pxyzpart);
    
    
    TVector3 motherrec = v3Dtagged + v3Dpart;
    
    Double_t etaESDmother = motherrec.Eta();
    Double_t ptESDmother  = motherrec.Pt();
    Double_t phiESDmother = motherrec.Phi();
    if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
    
    
    // openinganglepropagated
    Double_t openingangle = v3Dtagged.Angle(v3Dpart);
    //printf("Openingangle %f\n",openingangle);
    
    // invmass
    Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
    Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
    // e propagate
    Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
    Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
    
    Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
    
    if(!results) return;
    
    results[0] = ptESDmother;
    results[1] = etaESDmother;
    results[2] = phiESDmother;
    results[3] = invmass;
    results[4] = openingangle;
    
  }
  
}
//_________________________________________________________________________________
void AliHFEelecbackground::FillOutput(Double_t *results, Double_t *resultsr, Int_t sign) 
{
  //
  // Fill the Data and MC THnSparseF 
  //

  if((!results) || (!resultsr)) return;
  
  Double_t co[6];
  co[0] = fPtESD;
  co[1] = results[0];
  co[2] = TMath::Abs(results[4]);
  co[3] = results[3];
  co[4] = sign;
  co[5] = 0.0;

  (dynamic_cast<THnSparseF *>(fList->At(kDatai)))->Fill(co);

  if((sign==kOs) && (!fUseAliKFCode)) {
    
    co[1] = resultsr[0];
    co[2] = TMath::Abs(resultsr[4]);
    co[3] = resultsr[3];
    co[4] = kR;
    co[5] = 0.0;

    (dynamic_cast<THnSparseF *>(fList->At(kDatai)))->Fill(co);
    
  }
  
  if(HasMCData()){

    // Reset
    co[1] = results[0];
    co[2] = TMath::Abs(results[4]);
    co[3] = results[3];

    // Origin
    co[4] = kElectronFromBackground;
    if((sign==kOs) && fIsPartner) co[4] = fIsFrom;
    
    // Splitted tracks
    co[5] = kNotSplitted;
    if(fIsSplittedTrack) {
      if(sign==kOs){
	co[5] = kSplittedOs;
      }
      else {
	co[5] = kSplittedSs;
      }
    }

    (dynamic_cast<THnSparseF *>(fList->At(kMCo)))->Fill(co);

  }
 
}    
//_______________________________________________________________________________________________
Bool_t AliHFEelecbackground::SingleTrackCut(AliESDtrack* const trackPart) const
{
  //
  // Return minimum quality for the partner
  //
  
  if(trackPart->GetKinkIndex(0)>0) return kFALSE;


  UInt_t status = trackPart->GetStatus();

  if(fRequireITSStandalone > 0) {

    /////////////////////
    // ITS Standalone
    ////////////////////
    
    if(fRequireITSStandalone==1) {
      if(((status & AliESDtrack::kITSin) == 0 || (trackPart->IsPureITSStandalone()) || ((status&AliESDtrack::kITSrefit) ==0))) return kFALSE;
    }
    
    if(fRequireITSStandalone==2) {
      if(!trackPart->IsPureITSStandalone() || ((status&AliESDtrack::kITSrefit)==0)) return kFALSE;
    }
    
    // Chi2
    Double_t chi2 = trackPart->GetITSchi2();
    if(chi2 > fMinITSChi2) return kFALSE;

    // Min Nb of clusters
    Int_t nbcl = trackPart->GetITSclusters(0);
    if(nbcl < fMinNbCls)  return kFALSE;  

    // Min Nb of points in SDD and SPD
    Int_t nbSDDSPD = 0;
    for(Int_t layer = 0; layer < 4; layer++){
      if(trackPart->HasPointOnITSLayer(layer)) nbSDDSPD++;
    }
    if(nbSDDSPD < fMinNbClsSDDSPD) return kFALSE;
    
    
  }
  else {

    /////////
    // TPC
    /////////
    
    if((status&AliESDtrack::kTPCrefit)==0) return kFALSE;

    // Min Nb of clusters
    Int_t nbcl = trackPart->GetTPCclusters(0);
    if(nbcl < fMinNbCls) return kFALSE;   
    
  }

  return kTRUE;

}  
//_______________________________________________________________________________________________
Bool_t AliHFEelecbackground::PIDTrackCut(AliESDtrack* const trackPart)
{
  //
  // PID for the partner using TPC or ITS
  //
  
  if(fRequireITSStandalone > 0) {

    /////////////////////
    // ITS Standalone
    ////////////////////
    
    // signal
    Double_t itsSignal = trackPart->GetITSsignal();
    Double_t p = trackPart->P();
    
    if(fDebugLevel > 1) {        
      (dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))->Fill(p,itsSignal);
    }

    ///////////
    // PID
    //////////
    if(fPIDPartner) {
      
      // Take signal trackPart
      Double_t dEdxSamplesPart[4];
      trackPart->GetITSdEdxSamples(dEdxSamplesPart);

      // Cut at 2 sigma
      if(!fPIDMethodPartnerITS) fPIDMethodPartnerITS = new AliESDpid;
      Float_t nsigma = fPIDMethodPartnerITS->NumberOfSigmasITS(trackPart, AliPID::kElectron);
      if(TMath::Abs(nsigma) > 2.0) return kFALSE;
      
      // fill signal
      if(fDebugLevel > 1) {        
	
	Double_t entries[5];
	entries[0] = p;
	entries[1] = dEdxSamplesPart[0];
	entries[2] = dEdxSamplesPart[1];
	entries[3] = dEdxSamplesPart[2];
	entries[4] = dEdxSamplesPart[3];

	(dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))->Fill(p,itsSignal);
	(dynamic_cast<THnSparseF *>(fList->At(kMCcutPart2)))->Fill(entries);
	
      }
      
    }
    
  }
  else {

    /////////
    // TPC
    /////////
    
    Double_t tpcSignal = trackPart->GetTPCsignal();
    Double_t p = trackPart->GetInnerParam() ? trackPart->GetInnerParam()->P() : trackPart->P();

    if(fDebugLevel > 1) {        
      //printf("tpcSignal %f\n",tpcSignal);
      (dynamic_cast<TH2F *>(fList->At(kMCcutPart0)))->Fill(p,tpcSignal);
    }

    // PID
    if(fPIDPartner) {
      if(!fPIDMethodPartner) return kFALSE;
      AliHFEpidObject hfetrack;
      hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
      hfetrack.fRecTrack = trackPart;
      //if(HasMCData()) hfetrack.fMCtrack = mctrack;
      if(!fPIDMethodPartner->IsSelected(&hfetrack)) return kFALSE;
      
      if(fDebugLevel > 1) {  
	(dynamic_cast<TH2F *>(fList->At(kMCcutPart1)))->Fill(p,tpcSignal);
      }
    }   
    
  }

  return kTRUE;

}
//__________________________________________________________________________________________
Bool_t AliHFEelecbackground::ShareCluster(AliESDtrack * const track1,AliESDtrack * const track2) const
{
  //
  // Look if the two tracks shared clusters in the TPC or in the ITS depending on the method
  //
  // For TPC:
  // hsfval: number of shared clusters 
  // hsmval: quality of the tracks
  //
  // For ITS:
  // compare the dEdx in the ITS
  //

  if(!fRequireITSStandalone) {

    //////////
    // TPC
    //////////

    
    Int_t nh = 0;
    Int_t an = 0;
    Int_t ns = 0;
    Float_t hsmval = 0.0;
    Float_t hsfval = 0.0;
   
    for(unsigned int imap=0;imap<track1->GetTPCClusterMap().GetNbits(); imap++) {
      if(track1->GetTPCClusterMap().TestBitNumber(imap) &&
         track2->GetTPCClusterMap().TestBitNumber(imap)) {
        // Do they share it ?
        if (track1->GetTPCSharedMap().TestBitNumber(imap) &&
            track2->GetTPCSharedMap().TestBitNumber(imap))
          {
            an++;
            nh+=2;
            ns+=2;
          }
	else {
          an--;
          nh+=2;
        }
      }
      else if (track1->GetTPCClusterMap().TestBitNumber(imap) ||
               track2->GetTPCClusterMap().TestBitNumber(imap)) {
	an++;
        nh++;
      }
    }

    if (nh >0) {
      hsmval = an*1.0/nh;
      hsfval = ns*1.0/nh;
    }
    
    
    if((hsfval > 0.15) || (hsmval > -0.15)) return kTRUE; //they share cluster 
    else return kFALSE;   
 
  
  }
  else {

    
    //////////
    // ITS
    /////////

    // Take signals 
    Double_t dEdxSamples1[4];
    track1->GetITSdEdxSamples(dEdxSamples1);
    Double_t dEdxSamples2[4];
    track2->GetITSdEdxSamples(dEdxSamples2);
    
    // If there are matching
    Int_t nbClusters = 0;
    Bool_t match[4] = {kTRUE,kTRUE,kTRUE,kTRUE};
    Double_t limit[4] = {1.5,1.5,1.5,1.5};
    for(Int_t layer = 0; layer < 4; layer++) {
      if(track1->HasPointOnITSLayer(layer+2) && track2->HasPointOnITSLayer(layer+2)) {
	if(TMath::Abs(dEdxSamples1[layer]-dEdxSamples2[layer])>limit[layer]) match[layer] = kFALSE;
	nbClusters++;
      }
    }
    //printf("nbClusters %d\n",nbClusters);
    
    // fill signal   
    if(fDebugLevel > 1) {        
      Double_t entriesSplit[5];
      entriesSplit[0] = 0.0;
      if(fIsSplittedTrack) entriesSplit[0] = 1.0; 
      
      for(Int_t layer = 0; layer < 4; layer++) {
	if(track1->HasPointOnITSLayer(layer+2) && track2->HasPointOnITSLayer(layer+2)) {
	  entriesSplit[layer+1] = dEdxSamples1[layer]-dEdxSamples2[layer];
	}
	else entriesSplit[layer+1] = -100.0;
      }
      (dynamic_cast<THnSparseF *>(fList->At(kMCcutPart3)))->Fill(entriesSplit);
    }

    // Return
    Int_t nbClustersNotClose = 0;
    for(Int_t layer = 0; layer < 4; layer++) {
      if(!match[layer]) nbClustersNotClose++;
    }
    if((nbClusters > 1) && (nbClustersNotClose > 0.75*nbClusters)) return kFALSE;
    else return kTRUE;
    
  }   
  
}
//____________________________________________________________________________________________________________
void AliHFEelecbackground::SetPIDPartner() {

  //
  // Init the stuff for PID on the partner track
  //

  fPIDPartner = kTRUE;

  if(fRequireITSStandalone == 0) {
    
    if(!fPIDMethodPartner) {
      fPIDMethodPartner = new AliHFEpid();
      fPIDMethodPartner->InitializePID("Strategy1");     // 3 sigma cut in TPC
    }

  }
  else {
    
    if(!fPIDMethodPartnerITS) fPIDMethodPartnerITS = new AliESDpid;
    
  }
  
}
//______________________________________________________________________________________________
void AliHFEelecbackground::SetEvent(AliESDEvent* const ESD)
{
  //
  // Set the AliESD Event, the magnetic field and the primary vertex
  //
  
  fESD1=ESD;
  fBz=fESD1->GetMagneticField();
  fkVertex = fESD1->GetPrimaryVertex();

}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::IsMotherGamma(Int_t tr) {

  //
  // Return the lab of gamma mother or -1 if not gamma
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // Take mother
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother   = particle->GetFirstMother(); 
  if((imother < 0) || (imother >= stack->GetNtrack())) return -1;  
  TParticle * mother = stack->Particle(imother);
  if(!mother) return -1;

  // Check gamma    
  Int_t pdg = mother->GetPdgCode();
  if(TMath::Abs(pdg) == 22) return imother;
  if(TMath::Abs(pdg) == 11) {
    return IsMotherGamma(imother);
  }
  return -1;
 
}
//
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::IsMotherPi0(Int_t tr) {

  //
  // Return the lab of pi0 mother or -1 if not pi0
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // Take mother
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother   = particle->GetFirstMother(); 
  if((imother < 0) || (imother >= stack->GetNtrack())) return -1;  
  TParticle * mother = stack->Particle(imother);
  if(!mother) return -1;

  // Check gamma    
  Int_t pdg = mother->GetPdgCode();
  if(TMath::Abs(pdg) == 111) return imother;
  if(TMath::Abs(pdg) == 11) {
    return IsMotherPi0(imother);
  }
  return -1;
 
}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::IsMotherEta(Int_t tr) {

  //
  // Return the lab of pi0 mother or -1 if not pi0
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // Take mother
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother   = particle->GetFirstMother(); 
  if((imother < 0) || (imother >= stack->GetNtrack())) return -1;  
  TParticle * mother = stack->Particle(imother);
  if(!mother) return -1;

  // Check gamma    
  Int_t pdg = mother->GetPdgCode();
  if(TMath::Abs(pdg) == 221) return imother;
  if(TMath::Abs(pdg) == 11) {
    return IsMotherEta(imother);
  }
  return -1;
 
}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::IsMotherC(Int_t tr) {

  //
  // Return the lab of signal mother or -1 if not signal
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // Take mother
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother   = particle->GetFirstMother(); 
  if((imother < 0) || (imother >= stack->GetNtrack())) return -1;  
  TParticle * mother = stack->Particle(imother);
  if(!mother) return -1;

  // Check gamma    
  Int_t pdg = mother->GetPdgCode();
  if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
  if(TMath::Abs(pdg) == 11) {
    return IsMotherC(imother);
  }
  return -1;
 
}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::IsMotherB(Int_t tr) {

  //
  // Return the lab of signal mother or -1 if not signal
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // Take mother
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother   = particle->GetFirstMother(); 
  if((imother < 0) || (imother >= stack->GetNtrack())) return -1;  
  TParticle * mother = stack->Particle(imother);
  if(!mother) return -1;

  // Check gamma    
  Int_t pdg = mother->GetPdgCode();
  if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
  if(TMath::Abs(pdg) == 11) {
    return IsMotherB(imother);
  }
  return -1;
 
}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::GetPdg(Int_t tr) {

  //
  // Simply pdg code
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // MC Information
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t pdg = particle->GetPdgCode();

  return pdg;
 
}
//____________________________________________________________________________________________________________
Int_t AliHFEelecbackground::GetLabMother(Int_t tr) {

  //
  // Simply lab mother
  //

  AliStack* stack = fMCEvent->Stack();
  if((tr < 0) || (tr >= stack->GetNtrack())) return -1;  

  // MC Information
  TParticle * particle = stack->Particle(tr);
  if(!particle) return -1;
  Int_t imother = particle->GetFirstMother(); 

  return imother;
 
}
//_______________________________________________________________________________________________
void AliHFEelecbackground::PostProcess()
{
  //
  // Post process the histos and extract the background pt spectra
  //

  if(!fList) return;

  gStyle->SetPalette(1);
  gStyle->SetOptStat(1111);
  gStyle->SetPadBorderMode(0);
  gStyle->SetCanvasColor(10);
  gStyle->SetPadLeftMargin(0.13);
  gStyle->SetPadRightMargin(0.13);

  /////////////////////////
  // Take the THnSparseF
  /////////////////////////
  THnSparseF *hsSparseData = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassData")); 
  THnSparseF *hsSparseMC = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassMC")); 
  THnSparseF *hsSparseCutPassedMC = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC")); 

  /////////////////////////////////
  // Cuts on the opening angle
  ////////////////////////////////
  TAxis *axisOpeningAngleData = hsSparseData->GetAxis(2);
  Int_t binCutData = axisOpeningAngleData->FindBin(fOpeningAngleCut);
  hsSparseData->GetAxis(2)->SetRange(1,binCutData);

  if(hsSparseMC) {
    TAxis *axisOpeningAngleMC = hsSparseMC->GetAxis(2);
    Int_t binCutMC = axisOpeningAngleMC->FindBin(fOpeningAngleCut);
    hsSparseMC->GetAxis(2)->SetRange(1,binCutMC);
  }
  
  /////////////////////////
  // Prepare the histos
  ////////////////////////  

  TAxis *ptaxisinvmass = hsSparseData->GetAxis(3);
  Int_t  nbinsptinvmass = ptaxisinvmass->GetNbins();  
  
  TH1D **invmassosptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassssptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassrptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassdiffptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassgammaptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmasspi0ptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassetaptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassCptproj = new TH1D*[nbinsptinvmass];
  TH1D **invmassBptproj = new TH1D*[nbinsptinvmass];

  TH1D *yieldPtFound = (TH1D *) hsSparseData->Projection(0);
  yieldPtFound->SetName("Found yield");
  yieldPtFound->Reset();

  TH1D *yieldPtSourcesMC = 0x0;
  TH1D *yieldPtSignalCutMC = 0x0;
  if(hsSparseMC) {
    yieldPtSourcesMC = (TH1D *) hsSparseMC->Projection(0);
    yieldPtSourcesMC->SetName("Found yield");
    yieldPtSourcesMC->Reset();
    
    yieldPtSignalCutMC = (TH1D *) hsSparseMC->Projection(0);
    yieldPtSignalCutMC->SetName("Found yield");
    yieldPtSignalCutMC->Reset();
  }
  
  ////////////
  // canvas
  ///////////
  Int_t nbrow = (Int_t) (nbinsptinvmass/5);
  TString namecanvas("InvMassSpectra");
  TCanvas * canvas =new TCanvas(namecanvas,namecanvas,800,800);
  canvas->Divide(5,nbrow+1);

  /////////////////////////////
  // Loop on pt bins
  /////////////////////////////

  for(Int_t k=1; k <= nbinsptinvmass; k++){

    Double_t lowedge = ptaxisinvmass->GetBinLowEdge(k);
    Double_t upedge  = ptaxisinvmass->GetBinUpEdge(k);

    // Pt bin
    hsSparseData->GetAxis(0)->SetRange(k,k);
    if(hsSparseMC) hsSparseMC->GetAxis(0)->SetRange(k,k);
    
    // 
    hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);
    invmassosptproj[k-1] = hsSparseData->Projection(3);
    hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);
    invmassssptproj[k-1] = hsSparseData->Projection(3);
    hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);
    invmassrptproj[k-1]  = hsSparseData->Projection(3);
    hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins()); 
    invmassgammaptproj[k-1] = 0x0;
    invmasspi0ptproj[k-1] = 0x0;
    invmassetaptproj[k-1] = 0x0;
    invmassCptproj[k-1] = 0x0;
    invmassBptproj[k-1] = 0x0;
    if(hsSparseMC) {
      hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);
      invmassgammaptproj[k-1] = hsSparseMC->Projection(3);
      hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);
      invmasspi0ptproj[k-1] = hsSparseMC->Projection(3);
      hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);
      invmassetaptproj[k-1] = hsSparseMC->Projection(3);
      hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);
      invmassCptproj[k-1] = hsSparseMC->Projection(3);
      hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);
      invmassBptproj[k-1] = hsSparseMC->Projection(3);
      hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins()); 
    }      

    invmassdiffptproj[k-1] = (TH1D *) invmassosptproj[k-1]->Clone();
    TString name("Invmassdiffptbin");
    name += k;
    invmassdiffptproj[k-1]->SetName(name);
    invmassdiffptproj[k-1]->Add(invmassssptproj[k-1],-1.0);

    TString namee("p_{T} tagged from ");
    namee += lowedge;
    namee += " GeV/c to ";
    namee += upedge;
    namee += " GeV/c";

    invmassosptproj[k-1]->SetTitle((const char*)namee);
    invmassssptproj[k-1]->SetTitle((const char*)namee);
    invmassrptproj[k-1]->SetTitle((const char*)namee);
    invmassdiffptproj[k-1]->SetTitle((const char*)namee);
    if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetTitle((const char*)namee);
    if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetTitle((const char*)namee);
    if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetTitle((const char*)namee);
    if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetTitle((const char*)namee);
    if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetTitle((const char*)namee);
    
    invmassosptproj[k-1]->SetStats(0);
    invmassssptproj[k-1]->SetStats(0);
    invmassrptproj[k-1]->SetStats(0);
    invmassdiffptproj[k-1]->SetStats(0);
    if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetStats(0);
    if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetStats(0);
    if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetStats(0);
    if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetStats(0);
    if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetStats(0);
        
    Double_t yieldf = invmassdiffptproj[k-1]->Integral();
    if(invmassetaptproj[k-1] && invmasspi0ptproj[k-1] && invmassgammaptproj[k-1] && invmassCptproj[k-1] && invmassBptproj[k-1]) {
      Double_t yieldg = invmassetaptproj[k-1]->Integral() + invmasspi0ptproj[k-1]->Integral() + invmassgammaptproj[k-1]->Integral();
      yieldPtSourcesMC->SetBinContent(k,yieldg);
      
      Double_t yieldsignal = invmassCptproj[k-1]->Integral() + invmassBptproj[k-1]->Integral();
      yieldPtSignalCutMC->SetBinContent(k,yieldsignal);
    }

    yieldPtFound->SetBinContent(k,yieldf);
    
    canvas->cd(k);
    invmassosptproj[k-1]->Draw();
    invmassssptproj[k-1]->Draw("same");
    invmassdiffptproj[k-1]->Draw("same");
    invmassrptproj[k-1]->Draw("same");
    TLegend *legiv = new TLegend(0.4,0.6,0.89,0.89);
    legiv->AddEntry(invmassosptproj[k-1],"Opposite signs","p"); 
    legiv->AddEntry(invmassssptproj[k-1],"Same signs","p"); 
    legiv->AddEntry(invmassdiffptproj[k-1],"(Opposite - Same) signs","p"); 
    legiv->AddEntry(invmassrptproj[k-1],"rotated","p"); 
    if(invmassgammaptproj[k-1]) legiv->AddEntry(invmassgammaptproj[k-1],"e^{+}e^{-} from #gamma","p"); 
    if(invmasspi0ptproj[k-1]) legiv->AddEntry(invmasspi0ptproj[k-1],"e^{+}e^{-} from #pi^{0}","p"); 
    if(invmassetaptproj[k-1]) legiv->AddEntry(invmassetaptproj[k-1],"e^{+}e^{-} from #eta","p"); 
    legiv->Draw("same");
  
    hsSparseData->GetAxis(0)->SetRange(1,hsSparseData->GetAxis(0)->GetNbins()); 
    if(hsSparseMC) hsSparseMC->GetAxis(0)->SetRange(1,hsSparseMC->GetAxis(0)->GetNbins()); 

  }

  ////////////////////////////////////////////////////
  // End of plotting: do subtraction of background
  ///////////////////////////////////////////////////

  yieldPtFound->SetStats(0);
  if(yieldPtSourcesMC) yieldPtSourcesMC->SetStats(0); 
  if(yieldPtSignalCutMC) yieldPtSignalCutMC->SetStats(0);  

  TCanvas * canvasfin =new TCanvas("ResultsElecBackGround","ResultsElecBackGround",800,800);
  canvasfin->cd(1);
  yieldPtFound->Draw();
  if(yieldPtSourcesMC && yieldPtSignalCutMC) {
    yieldPtSourcesMC->Draw("same");
    yieldPtSignalCutMC->Draw("same");
    TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
    lega->AddEntry(yieldPtFound,"Contributions found","l"); 
    lega->AddEntry(yieldPtSourcesMC,"Contributions of e^{+}e^{-} from #gamma, #pi^{0} and #eta","l"); 
    lega->AddEntry(yieldPtSignalCutMC,"Contributions of e^{+}e^{-} from C and B","l"); 
    lega->Draw("same");
  }
  
  if(hsSparseCutPassedMC){
    hsSparseCutPassedMC->GetAxis(1)->SetRange(1,1); 
    hsSparseCutPassedMC->GetAxis(2)->SetRange(1,4); 
    TH1D *hsSparseCutPassedMCproj = hsSparseCutPassedMC->Projection(0);

    TH1D *cYieldPtFound = (TH1D*)yieldPtFound->Clone("RatioEfficiency");
    if(hsSparseCutPassedMCproj->Integral() > 0.0) cYieldPtFound->Divide(hsSparseCutPassedMCproj);

    TCanvas * canvasfratio =new TCanvas("RatioEfficiency","RatioEfficiency",800,800);
    canvasfratio->cd(1);
    cYieldPtFound->Draw();
  }

  //////////////////////////
  // fListPostProcess
  /////////////////////////
  
  if(!fListPostProcess) fListPostProcess = new TList();
  fListPostProcess->SetName("ListPostProcess");
  
  for(Int_t k=0; k < nbinsptinvmass; k++){
    fListPostProcess->AddAt(invmassosptproj[k],kOos+kNOutput*k);
    fListPostProcess->AddAt(invmassssptproj[k],kOss+kNOutput*k);
    fListPostProcess->AddAt(invmassrptproj[k],kOr+kNOutput*k);
    fListPostProcess->AddAt(invmassdiffptproj[k],kOdiff+kNOutput*k);
    if(invmassgammaptproj[k]) fListPostProcess->AddAt(invmassgammaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromGamma);
    if(invmasspi0ptproj[k]) fListPostProcess->AddAt(invmasspi0ptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromPi0);
    if(invmassetaptproj[k]) fListPostProcess->AddAt(invmassetaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromEta);
    if(invmassCptproj[k]) fListPostProcess->AddAt(invmassCptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromC);
    if(invmassBptproj[k]) fListPostProcess->AddAt(invmassBptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromB);
  }

  fListPostProcess->AddAt(yieldPtFound,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass);
  if(yieldPtSourcesMC) fListPostProcess->AddAt(yieldPtSourcesMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+1);
  if(yieldPtSignalCutMC) fListPostProcess->AddAt(yieldPtSignalCutMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+2);
  
  // delete dynamic array
  delete[] invmassosptproj;
  delete[] invmassssptproj;
  delete[] invmassrptproj;
  delete[] invmassdiffptproj;
  delete[] invmassgammaptproj;
  delete[] invmasspi0ptproj;
  delete[] invmassetaptproj;
  delete[] invmassCptproj;
  delete[] invmassBptproj;

}
//_______________________________________________________________________________________________
void AliHFEelecbackground::Plot() const
{
  //
  // Plot the output
  //
  
  if(!fList) return;

  gStyle->SetPalette(1);
  gStyle->SetOptStat(1111);
  gStyle->SetPadBorderMode(0);
  gStyle->SetCanvasColor(10);
  gStyle->SetPadLeftMargin(0.13);
  gStyle->SetPadRightMargin(0.13);


  /////////////////////////
  // Take the THnSparseF
  /////////////////////////
  THnSparseF *hsSparseData = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassData")); 
  THnSparseF *hsSparseMC = dynamic_cast<THnSparseF *>(fList->FindObject("OpeningangleinvmassMC")); 
  
  ////////////////////
  // Opening angle
  ////////////////////

  // Opening angle one direction
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);  
  TH1D *openingangleppproj = hsSparseData->Projection(2);
  hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);  
  TH1D *openinganglennproj = hsSparseData->Projection(2);
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);  
  TH1D *openinganglessproj = hsSparseData->Projection(2);
  hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);  
  TH1D *openinganglerproj  = hsSparseData->Projection(2);
  hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);  
  TH1D *openingangleosproj = hsSparseData->Projection(2);
  hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins()); 

  TH1D *openinganglegammaproj = 0x0;
  TH1D *openinganglepi0proj = 0x0;
  TH1D *openingangleCproj = 0x0;
  TH1D *openingangleBproj = 0x0;
  TH1D *openingangleetaproj = 0x0;
  TH1D *openingangleSplittedTrackssproj = 0x0;
  TH1D *openingangleSplittedTrackosproj = 0x0;
  if(hsSparseMC) {
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);  
    openinganglegammaproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);  
    openinganglepi0proj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);  
    openingangleetaproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);  
    openingangleCproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);  
    openingangleBproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins());  
    hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);  
    openingangleSplittedTrackssproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);  
    openingangleSplittedTrackosproj = hsSparseMC->Projection(2);
    hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins()); 
  }

  // Projection pt-opening angle
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1); 
  TH2D *openingangleppproj2D = hsSparseData->Projection(0,2);
  hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1); 
  TH2D *openinganglennproj2D = hsSparseData->Projection(0,2);
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1); 
  TH2D *openinganglessproj2D = hsSparseData->Projection(0,2);
  hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1); 
  TH2D *openinganglerproj2D  = hsSparseData->Projection(0,2);
  hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1); 
  TH2D *openingangleosproj2D = hsSparseData->Projection(0,2);
  hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins()); 

  TH2D *openinganglegammaproj2D = 0x0;
  TH2D *openinganglepi0proj2D = 0x0;
  TH2D *openingangleCproj2D = 0x0;
  TH2D *openingangleBproj2D = 0x0;
  TH2D *openingangleetaproj2D = 0x0;
  TH2D *openingangleSplittedTrackssproj2D = 0x0;
  TH2D *openingangleSplittedTrackosproj2D = 0x0;
  if(hsSparseMC) {
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1); 
    openinganglegammaproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1); 
    openinganglepi0proj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1); 
    openingangleetaproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1); 
    openingangleCproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1); 
    openingangleBproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(4)->SetRange(1, hsSparseMC->GetAxis(4)->GetNbins());  
    hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1); 
    openingangleSplittedTrackssproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1); 
    openingangleSplittedTrackosproj2D = hsSparseMC->Projection(0,2);
    hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
  }

  openingangleppproj2D->SetStats(0);
  openinganglennproj2D->SetStats(0);
  openinganglessproj2D->SetStats(0);
  openinganglerproj2D->SetStats(0);
  openingangleosproj2D->SetStats(0);
  if(openinganglegammaproj2D) openinganglegammaproj2D->SetStats(0);
  if(openinganglepi0proj2D) openinganglepi0proj2D->SetStats(0);
  if(openingangleCproj2D) openingangleCproj2D->SetStats(0);
  if(openingangleBproj2D) openingangleBproj2D->SetStats(0);
  if(openingangleetaproj2D) openingangleetaproj2D->SetStats(0);
  if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->SetStats(0);
  if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->SetStats(0);

  openingangleppproj2D->SetTitle("openingangleppproj2D");
  openinganglennproj2D->SetTitle("openinganglennproj2D");
  openinganglessproj2D->SetTitle("openinganglessproj2D");
  openinganglerproj2D->SetTitle("openinganglerproj2D");
  openingangleosproj2D->SetTitle("openingangleosproj2D");
  if(openinganglegammaproj2D) openinganglegammaproj2D->SetTitle("openinganglegammaproj2D");
  if(openinganglepi0proj2D) openinganglepi0proj2D->SetTitle("openinganglepi0proj2D");
  if(openingangleCproj2D) openingangleCproj2D->SetTitle("openingangleCproj2D");
  if(openingangleBproj2D) openingangleBproj2D->SetTitle("openingangleBproj2D");
  if(openingangleetaproj2D) openingangleetaproj2D->SetTitle("openingangleetaproj2D");
  if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->SetTitle("openingangleSplittedTrackssproj2D");
  if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->SetTitle("openingangleSplittedTrackosproj2D");  

  openingangleppproj->SetStats(0);
  openinganglennproj->SetStats(0);
  openinganglessproj->SetStats(0);
  openinganglerproj->SetStats(0);
  openingangleosproj->SetStats(0);
  if(openinganglegammaproj) openinganglegammaproj->SetStats(0);
  if(openinganglepi0proj) openinganglepi0proj->SetStats(0);
  if(openingangleCproj) openingangleCproj->SetStats(0);
  if(openingangleBproj) openingangleBproj->SetStats(0);
  if(openingangleetaproj) openingangleetaproj->SetStats(0);
  if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->SetStats(0);
  if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->SetStats(0);

  openingangleppproj->SetTitle("");
  openinganglennproj->SetTitle("");
  openinganglessproj->SetTitle("");
  openinganglerproj->SetTitle("");
  openingangleosproj->SetTitle("");
  if(openinganglegammaproj) openinganglegammaproj->SetTitle("");
  if(openinganglepi0proj) openinganglepi0proj->SetTitle("");
  if(openingangleCproj) openingangleCproj->SetTitle("");
  if(openingangleBproj) openingangleBproj->SetTitle("");
  if(openingangleetaproj) openingangleetaproj->SetTitle("");
  if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->SetTitle("");
  if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->SetTitle("");

  ////////////////////////////
  // Invariant mass
  ///////////////////////////

  // Cuts on the opening angle
  TAxis *axisOpeningAngleData = hsSparseData->GetAxis(2);
  Int_t binCutData = axisOpeningAngleData->FindBin(fOpeningAngleCut);
  hsSparseData->GetAxis(2)->SetRange(1,binCutData);
  
  // Debug
  //printf("Get Bin low edge %f, Get Bin Up edge %f for hsSparseData\n",axisOpeningAngleData->GetBinLowEdge(binCutData),axisOpeningAngleData->GetBinUpEdge(binCutData));

  // Invariant mass
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1);  
  TH1D *invmassppproj = hsSparseData->Projection(3);
  hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1);  
  TH1D *invmassnnproj = hsSparseData->Projection(3);
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1);  
  TH1D *invmassssproj = hsSparseData->Projection(3);
  hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1);  
  TH1D *invmassrproj  = hsSparseData->Projection(3);
  hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1);  
  TH1D *invmassosproj = hsSparseData->Projection(3);
  hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins()); 

  TH1D *invmassgammaproj = 0x0;
  TH1D *invmasspi0proj = 0x0;
  TH1D *invmassCproj = 0x0;
  TH1D *invmassBproj = 0x0;
  TH1D *invmassetaproj = 0x0;
  TH1D *invmassSplittedTrackssproj = 0x0;
  TH1D *invmassSplittedTrackosproj = 0x0;
  if(hsSparseMC) {
    TAxis *axisOpeningAngleMC = hsSparseMC->GetAxis(2);
    Int_t binCutMC = axisOpeningAngleMC->FindBin(fOpeningAngleCut);
    hsSparseMC->GetAxis(2)->SetRange(1,binCutMC);
    
    // Debug
    //printf("Get Bin low edge %f, Get Bin Up edge %f for hsSparseMC\n",axisOpeningAngleMC->GetBinLowEdge(binCutMC),axisOpeningAngleMC->GetBinUpEdge(binCutMC));
    
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1);  
    invmassgammaproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1);  
    invmasspi0proj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1);  
    invmassetaproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1);  
    invmassCproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1);  
    invmassBproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins()); 
    hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1);  
    invmassSplittedTrackssproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1);  
    invmassSplittedTrackosproj = hsSparseMC->Projection(3);
    hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins());
  }
  
  invmassppproj->SetStats(0);
  invmassnnproj->SetStats(0);
  invmassssproj->SetStats(0);
  invmassrproj->SetStats(0);
  invmassosproj->SetStats(0);
  if(invmassgammaproj) invmassgammaproj->SetStats(0);
  if(invmasspi0proj) invmasspi0proj->SetStats(0);
  if(invmassCproj) invmassCproj->SetStats(0);
  if(invmassBproj) invmassBproj->SetStats(0);
  if(invmassetaproj) invmassetaproj->SetStats(0);
  if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->SetStats(0);
  if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->SetStats(0);

  invmassppproj->SetTitle("");
  invmassnnproj->SetTitle("");
  invmassssproj->SetTitle("");
  invmassrproj->SetTitle("");
  invmassosproj->SetTitle("");
  if(invmassgammaproj) invmassgammaproj->SetTitle("");
  if(invmasspi0proj) invmasspi0proj->SetTitle("");
  if(invmassCproj) invmassCproj->SetTitle("");
  if(invmassBproj) invmassBproj->SetTitle("");
  if(invmassetaproj) invmassetaproj->SetTitle("");
  if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->SetTitle("");
  if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->SetTitle("");

  // Projection pt-invariant mass angle
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kPp+1); 
  TH2D *invmassppproj2D = hsSparseData->Projection(0,3);
  hsSparseData->GetAxis(4)->SetRange(kNn+1,kNn+1); 
  TH2D *invmassnnproj2D = hsSparseData->Projection(0,3);
  hsSparseData->GetAxis(4)->SetRange(kPp+1,kNn+1); 
  TH2D *invmassssproj2D = hsSparseData->Projection(0,3);
  hsSparseData->GetAxis(4)->SetRange(kR+1,kR+1); 
  TH2D *invmassrproj2D  = hsSparseData->Projection(0,3);
  hsSparseData->GetAxis(4)->SetRange(kOs+1,kOs+1); 
  TH2D *invmassosproj2D = hsSparseData->Projection(0,3);
  hsSparseData->GetAxis(4)->SetRange(1,hsSparseData->GetAxis(4)->GetNbins()); 

  TH2D *invmassgammaproj2D = 0x0;
  TH2D *invmasspi0proj2D = 0x0;
  TH2D *invmassCproj2D = 0x0;
  TH2D *invmassBproj2D = 0x0;
  TH2D *invmassetaproj2D = 0x0;
  TH2D *invmassSplittedTrackssproj2D = 0x0;
  TH2D *invmassSplittedTrackosproj2D = 0x0;
  if(hsSparseMC) {
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromGamma+1,kElectronFromGamma+1); 
    invmassgammaproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromPi0+1,kElectronFromPi0+1); 
    invmasspi0proj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromEta+1,kElectronFromEta+1); 
    invmassetaproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromC+1,kElectronFromC+1); 
    invmassCproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(4)->SetRange(kElectronFromB+1,kElectronFromB+1); 
    invmassBproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(4)->SetRange(1,hsSparseMC->GetAxis(4)->GetNbins()); 
    hsSparseMC->GetAxis(5)->SetRange(kSplittedSs+1,kSplittedSs+1); 
    invmassSplittedTrackssproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(5)->SetRange(kSplittedOs+1,kSplittedOs+1); 
    invmassSplittedTrackosproj2D = hsSparseMC->Projection(0,3);
    hsSparseMC->GetAxis(5)->SetRange(1,hsSparseMC->GetAxis(5)->GetNbins()); 
  }

  
  invmassppproj2D->SetStats(0);
  invmassnnproj2D->SetStats(0);
  invmassssproj2D->SetStats(0);
  invmassrproj2D->SetStats(0);
  invmassosproj2D->SetStats(0);
  if(invmassgammaproj2D) invmassgammaproj2D->SetStats(0);
  if(invmasspi0proj2D) invmasspi0proj2D->SetStats(0);
  if(invmassCproj2D) invmassCproj2D->SetStats(0);
  if(invmassBproj2D) invmassBproj2D->SetStats(0);
  if(invmassetaproj2D) invmassetaproj2D->SetStats(0);
  if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->SetStats(0);
  if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->SetStats(0);
  
  invmassppproj2D->SetTitle("invmassppproj2D");
  invmassnnproj2D->SetTitle("invmassnnproj2D");
  invmassssproj2D->SetTitle("invmassssproj2D");
  invmassrproj2D->SetTitle("invmassrproj2D");
  invmassosproj2D->SetTitle("invmassosproj2D");
  if(invmassgammaproj2D) invmassgammaproj2D->SetTitle("invmassgammaproj2D");
  if(invmasspi0proj2D) invmasspi0proj2D->SetTitle("invmasspi0proj2D");
  if(invmassCproj2D) invmassCproj2D->SetTitle("invmassCproj2D");
  if(invmassBproj2D) invmassBproj2D->SetTitle("invmassBproj2D");
  if(invmassetaproj2D) invmassetaproj2D->SetTitle("invmassetaproj2D");
  if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->SetTitle("invmassSplittedTrackssproj2D");
  if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->SetTitle("invmassSplittedTrackosproj2D");


  /////////////
  // Plot 
  ////////////

  // Draw histograms for opening angle
  TCanvas * copeningangle =new TCanvas("openingangle","Openingangle",800,800);
  copeningangle->cd();
  //openingangleppproj->Draw();
  //openinganglennproj->Draw("same");
  openinganglessproj->Draw();
  //openinganglerproj->Draw("same");
  openingangleosproj->Draw("same");
  if(openinganglegammaproj) openinganglegammaproj->Draw("same");
  if(openinganglepi0proj) openinganglepi0proj->Draw("same");
  //if(openingangleCproj) openingangleCproj->Draw("same");
  //if(openingangleBproj) openingangleBproj->Draw("same");
  if(openingangleetaproj) openingangleetaproj->Draw("same");
  if(openingangleSplittedTrackssproj) openingangleSplittedTrackssproj->Draw("same");
  if(openingangleSplittedTrackosproj) openingangleSplittedTrackosproj->Draw("same");
  TLegend *lego = new TLegend(0.4,0.6,0.89,0.89);
  //lego->AddEntry(openingangleppproj,"positive-positive","p");
  //lego->AddEntry(openinganglennproj,"negative-negative","p");
  lego->AddEntry(openinganglessproj,"same-sign","p");
  //lego->AddEntry(openinganglerproj,"rotated","p");
  lego->AddEntry(openingangleosproj,"positive-negative","p");
  if(openinganglegammaproj) lego->AddEntry(openinganglegammaproj,"e^{+}e^{-} from #gamma","p");
  if(openinganglepi0proj) lego->AddEntry(openinganglepi0proj,"e^{+}e^{-} from #pi^{0}","p");
  //if(openingangleCproj) lego->AddEntry(openingangleCproj,"e^{+}e^{-} from c","p");
  //if(openingangleBproj) lego->AddEntry(openingangleBproj,"e^{+}e^{-} from b","p");
  if(openingangleetaproj) lego->AddEntry(openingangleetaproj,"e^{+}e^{-} from #eta","p");
  if(openingangleSplittedTrackssproj) lego->AddEntry(openingangleSplittedTrackssproj,"Splitted tracks same sign","p");
  if(openingangleSplittedTrackosproj) lego->AddEntry(openingangleSplittedTrackosproj,"Splitted tracks opposite sign","p");
  lego->Draw("same");
  
  // Draw histograms for invariant mass
  TCanvas * cinvmass =new TCanvas("invmass","Invmass",800,800);
  cinvmass->cd();
  //invmassppproj->Draw();
  //invmassnnproj->Draw("same");
  invmassssproj->Draw();
  //invmassrproj->Draw("same");
  invmassosproj->Draw("same");
  if(invmassgammaproj) invmassgammaproj->Draw("same");
  if(invmasspi0proj) invmasspi0proj->Draw("same");
  //if(invmassCproj) invmassCproj->Draw("same");
  //if(invmassBproj) invmassBproj->Draw("same");
  if(invmassetaproj) invmassetaproj->Draw("same");
  if(invmassSplittedTrackssproj) invmassSplittedTrackssproj->Draw("same");
  if(invmassSplittedTrackosproj) invmassSplittedTrackosproj->Draw("same");
  TLegend *legi = new TLegend(0.4,0.6,0.89,0.89);
  //legi->AddEntry(invmassppproj,"positive-positive","p");
  //legi->AddEntry(invmassnnproj,"negative-negative","p");
  legi->AddEntry(invmassssproj,"same-sign","p");
  //legi->AddEntry(invmassrproj,"rotated","p");
  legi->AddEntry(invmassosproj,"positive-negative","p");
  if(invmassgammaproj) legi->AddEntry(invmassgammaproj,"e^{+}e^{-} from #gamma","p");
  if(invmasspi0proj) legi->AddEntry(invmasspi0proj,"e^{+}e^{-} from #pi^{0}","p");
  //if(invmassCproj) legi->AddEntry(invmassCproj,"e^{+}e^{-} from c","p");
  //if(invmassBproj) legi->AddEntry(invmassBproj,"e^{+}e^{-} from b","p");
  if(invmassetaproj) legi->AddEntry(invmassetaproj,"e^{+}e^{-} from #eta","p");
  if(invmassSplittedTrackssproj) legi->AddEntry(invmassSplittedTrackssproj,"Splitted tracks same sign","p");
  if(invmassSplittedTrackosproj) legi->AddEntry(invmassSplittedTrackosproj,"Splitted tracks opposite sign","p");
  legi->Draw("same");

  

  // Draw histograms for opening angle 2D
  TCanvas * copeningangle2D =new TCanvas("openingangle2D","Openingangle2D",800,800);
  copeningangle2D->Divide(6,2);
  copeningangle2D->cd(1);
  openingangleppproj2D->Draw("lego");
  copeningangle2D->cd(2);
  openinganglennproj2D->Draw("lego");
  copeningangle2D->cd(3);
  openinganglessproj2D->Draw("lego");
  copeningangle2D->cd(4);
  openinganglerproj2D->Draw("lego");
  copeningangle2D->cd(5);
  openingangleosproj2D->Draw("lego");
  copeningangle2D->cd(6);
  if(openinganglegammaproj2D) openinganglegammaproj2D->Draw("lego");
  copeningangle2D->cd(7);
  if(openinganglepi0proj2D) openinganglepi0proj2D->Draw("lego");
  copeningangle2D->cd(8);
  if(openingangleCproj2D) openingangleCproj2D->Draw("lego");
  copeningangle2D->cd(9);
  if(openingangleBproj2D) openingangleBproj2D->Draw("lego");
  copeningangle2D->cd(10);
  if(openingangleetaproj2D) openingangleetaproj2D->Draw("lego");
  copeningangle2D->cd(11);
  if(openingangleSplittedTrackssproj2D) openingangleSplittedTrackssproj2D->Draw("lego");
  copeningangle2D->cd(12);
  if(openingangleSplittedTrackosproj2D) openingangleSplittedTrackosproj2D->Draw("lego");
  
  // Draw histograms for invariant mass 2D
  TCanvas * cinvmass2D =new TCanvas("invmass2D","Invmass2D",800,800);
  cinvmass2D->Divide(6,2);
  cinvmass2D->cd(1);
  invmassppproj2D->Draw("lego");
  cinvmass2D->cd(2);
  invmassnnproj2D->Draw("lego");
  cinvmass2D->cd(3);
  invmassssproj2D->Draw("lego");
  cinvmass2D->cd(4);
  invmassrproj2D->Draw("lego");
  cinvmass2D->cd(5);
  invmassosproj2D->Draw("lego");
  cinvmass2D->cd(6);
  if(invmassgammaproj2D) invmassgammaproj2D->Draw("lego");
  cinvmass2D->cd(7);
  if(invmasspi0proj2D) invmasspi0proj2D->Draw("lego");
  cinvmass2D->cd(8);
  if(invmassCproj2D) invmassCproj2D->Draw("lego");
  cinvmass2D->cd(9);
  if(invmassBproj2D) invmassBproj2D->Draw("lego");
  cinvmass2D->cd(10);
  if(invmassetaproj2D) invmassetaproj2D->Draw("lego");
  cinvmass2D->cd(11);
  if(invmassSplittedTrackssproj2D) invmassSplittedTrackssproj2D->Draw("lego");
  cinvmass2D->cd(12);
  if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->Draw("lego");


  ////////////////////////
  // Cut efficiencies
  ////////////////////////

  THnSparseF *hsSparseMCe = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
  
  if(hsSparseMCe) {

    // init histos
    TAxis *axissources = hsSparseMCe->GetAxis(2);
    Int_t  nbsources = axissources->GetNbins(); 
    TAxis *axiscuts = hsSparseMCe->GetAxis(1);
    Int_t  nbcuts = axiscuts->GetNbins(); 
    TH1D **histopassedcuts = new TH1D*[nbsources*nbcuts];  
    Double_t *nbEntriesCuts = new Double_t[nbsources*nbcuts]; 

    //printf("Number of cuts %d\n",nbcuts);

    // canvas
    TCanvas * chsSparseMCeeff =new TCanvas("hsSparseMCeeffDebug","hsSparseMCeeffDebug",800,800);
    chsSparseMCeeff->Divide(3,1);
 
    // histos
    for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
      hsSparseMCe->GetAxis(2)->SetRange(sourceid+1,sourceid+1);  
      for(Int_t cut = 0; cut < nbcuts; cut++){
	hsSparseMCe->GetAxis(1)->SetRange(cut+1,cut+1); 
	histopassedcuts[sourceid*nbcuts+cut] = hsSparseMCe->Projection(0);
     	hsSparseMCe->GetAxis(1)->SetRange(1,hsSparseMCe->GetAxis(1)->GetNbins()); 
      }
      hsSparseMCe->GetAxis(2)->SetRange(1,hsSparseMCe->GetAxis(2)->GetNbins());  
    }

    // calcul efficiencies
    ///////////////////////
    // histos
    for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
      // Next is compared to the partner tracked
      for(Int_t cut = 2; cut < nbcuts; cut++){
	nbEntriesCuts[sourceid*nbcuts+cut] = histopassedcuts[sourceid*nbcuts+cut]->GetEntries();
	if(histopassedcuts[sourceid*nbcuts+1]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+cut]->Divide(histopassedcuts[sourceid*nbcuts+1]);
      }
      // First one is if the partner is tracked.
      nbEntriesCuts[sourceid*nbcuts+1] = histopassedcuts[sourceid*nbcuts+1]->GetEntries();
      if(histopassedcuts[sourceid*nbcuts]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+1]->Divide(histopassedcuts[sourceid*nbcuts]);
      // First one is input
      nbEntriesCuts[sourceid*nbcuts] = histopassedcuts[sourceid*nbcuts]->GetEntries();
    }
    
    /////////////
    // ratios
    ////////////
    for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
      for(Int_t cut = 1; cut < nbcuts; cut++){
	if(nbEntriesCuts[sourceid*nbcuts] > 0.0) nbEntriesCuts[sourceid*nbcuts+cut] = nbEntriesCuts[sourceid*nbcuts+cut]/nbEntriesCuts[sourceid*nbcuts]; 
      }
    }
    TH1F *ratioHistoEntriesGamma = new TH1F("ratioHistoEntriesGamma","", nbcuts-1, 0.0, nbcuts-1.0);
    TH1F *ratioHistoEntriesPi0 = new TH1F("ratioHistoEntriesPi0","", nbcuts-1, 0.0, nbcuts-1.0);
    TH1F *ratioHistoEntriesC = new TH1F("ratioHistoEntriesC","", nbcuts-1, 0.0, nbcuts-1.0);
    for(Int_t k = 1; k < nbcuts; k++){
      ratioHistoEntriesGamma->SetBinContent(k,nbEntriesCuts[nbcuts+k]);
      ratioHistoEntriesPi0->SetBinContent(k,nbEntriesCuts[2*nbcuts+k]);
      ratioHistoEntriesC->SetBinContent(k,nbEntriesCuts[4*nbcuts+k]);
    }     
    //
    TAxis *xAxisGamma = ratioHistoEntriesGamma->GetXaxis();
    xAxisGamma->SetBinLabel(1,"Partner tracked");
    xAxisGamma->SetBinLabel(2,"Opposite sign");
    xAxisGamma->SetBinLabel(3,"Single Track Cut");
    xAxisGamma->SetBinLabel(4,"Shared Clusters");
    xAxisGamma->SetBinLabel(5,"PID");
    xAxisGamma->SetBinLabel(6,"DCA");
    xAxisGamma->SetBinLabel(7,"Chi^{2}/Ndf");
    xAxisGamma->SetBinLabel(8,"Opening angle");
    xAxisGamma->SetBinLabel(9,"Invariant mass");
    //
    TAxis *xAxisPi0 = ratioHistoEntriesPi0->GetXaxis();
    xAxisPi0->SetBinLabel(1,"Partner tracked");
    xAxisPi0->SetBinLabel(2,"Opposite sign");
    xAxisPi0->SetBinLabel(3,"Single Track Cut");
    xAxisPi0->SetBinLabel(4,"Shared Clusters");
    xAxisPi0->SetBinLabel(5,"PID");
    xAxisPi0->SetBinLabel(6,"DCA");
    xAxisPi0->SetBinLabel(7,"Chi^{2}/Ndf");
    xAxisPi0->SetBinLabel(8,"Opening angle");
    xAxisPi0->SetBinLabel(9,"Invariant mass");
    //
    TAxis *xAxisC = ratioHistoEntriesC->GetXaxis();
    xAxisC->SetBinLabel(1,"Partner tracked");
    xAxisC->SetBinLabel(2,"Opposite sign");
    xAxisC->SetBinLabel(3,"Single Track Cut");
    xAxisC->SetBinLabel(4,"Shared Clusters");
    xAxisC->SetBinLabel(5,"PID");
    xAxisC->SetBinLabel(6,"DCA");
    xAxisC->SetBinLabel(7,"Chi^{2}/Ndf");
    xAxisC->SetBinLabel(8,"Opening angle");
    xAxisC->SetBinLabel(9,"Invariant mass");
    //
    TCanvas * cRatioHistoEntries =new TCanvas("cRatioHistoEntries","cRatioHistoEntries",800,800);
    cRatioHistoEntries->cd(1);
    ratioHistoEntriesGamma->SetStats(0);
    ratioHistoEntriesGamma->Draw();
    ratioHistoEntriesPi0->SetStats(0);
    ratioHistoEntriesPi0->Draw("same");
    ratioHistoEntriesC->SetStats(0);
    //ratioHistoEntriesC->Draw("same");
    TLegend *legEntries = new TLegend(0.4,0.6,0.89,0.89);
    legEntries->AddEntry(ratioHistoEntriesGamma,"#gamma","l");
    legEntries->AddEntry(ratioHistoEntriesPi0,"#pi^{0}","l");
    //legEntries->AddEntry(ratioHistoEntriesC,"c","p");
    legEntries->Draw("same"); 

    ////////////////////
    // plot Debug
    ///////////////////
    Int_t source = 1;
    chsSparseMCeeff->cd(1);
    histopassedcuts[source*nbcuts+0]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+1]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+2]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+3]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+4]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+5]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+6]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+7]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+8]->SetTitle("#gamma");
    histopassedcuts[source*nbcuts+9]->SetTitle("#gamma");
    //histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+1]->SetStats(0);
    histopassedcuts[source*nbcuts+2]->SetStats(0);
    histopassedcuts[source*nbcuts+3]->SetStats(0);
    histopassedcuts[source*nbcuts+4]->SetStats(0);
    histopassedcuts[source*nbcuts+5]->SetStats(0);
    histopassedcuts[source*nbcuts+6]->SetStats(0);
    histopassedcuts[source*nbcuts+7]->SetStats(0);
    histopassedcuts[source*nbcuts+8]->SetStats(0);
    histopassedcuts[source*nbcuts+9]->SetStats(0);
    //histopassedcuts[source*nbcuts+0]->Draw();
    //histopassedcuts[source*nbcuts+1]->Draw("");
    histopassedcuts[source*nbcuts+2]->Draw();
    histopassedcuts[source*nbcuts+3]->Draw("same");
    //histopassedcuts[source*nbcuts+4]->Draw("same");
    histopassedcuts[source*nbcuts+5]->Draw("same");
    histopassedcuts[source*nbcuts+6]->Draw("same");
    //histopassedcuts[source*nbcuts+7]->Draw("same");
    histopassedcuts[source*nbcuts+8]->Draw("same");
    histopassedcuts[source*nbcuts+9]->Draw("same");
    TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
    //legb->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
    //legb->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
    //legb->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
    //legb->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
    legb->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
    legb->Draw("same");

    source = 2;
    chsSparseMCeeff->cd(2);
    histopassedcuts[source*nbcuts+0]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+1]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+2]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+3]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+4]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+5]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+6]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+7]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+8]->SetTitle("#pi^{0}");
    histopassedcuts[source*nbcuts+9]->SetTitle("#pi^{0}");
    //histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+1]->SetStats(0);
    histopassedcuts[source*nbcuts+2]->SetStats(0);
    histopassedcuts[source*nbcuts+3]->SetStats(0);
    histopassedcuts[source*nbcuts+4]->SetStats(0);
    histopassedcuts[source*nbcuts+5]->SetStats(0);
    histopassedcuts[source*nbcuts+6]->SetStats(0);
    histopassedcuts[source*nbcuts+7]->SetStats(0);
    histopassedcuts[source*nbcuts+8]->SetStats(0);
    histopassedcuts[source*nbcuts+9]->SetStats(0);
    //histopassedcuts[source*nbcuts+0]->Draw();
    //histopassedcuts[source*nbcuts+1]->Draw();
    histopassedcuts[source*nbcuts+2]->Draw();
    histopassedcuts[source*nbcuts+3]->Draw("same");
    //histopassedcuts[source*nbcuts+4]->Draw("same");
    histopassedcuts[source*nbcuts+5]->Draw("same");
    histopassedcuts[source*nbcuts+6]->Draw("same");
    //histopassedcuts[source*nbcuts+7]->Draw("same");
    histopassedcuts[source*nbcuts+8]->Draw("same");
    histopassedcuts[source*nbcuts+9]->Draw("same");
    TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
    //legc->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
    //legc->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
    //legc->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
    //legc->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
    legc->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
    legc->Draw("same");

    source = 4;
    chsSparseMCeeff->cd(3);
    histopassedcuts[source*nbcuts+0]->SetTitle("C");
    histopassedcuts[source*nbcuts+1]->SetTitle("C");
    histopassedcuts[source*nbcuts+2]->SetTitle("C");
    histopassedcuts[source*nbcuts+3]->SetTitle("C");
    histopassedcuts[source*nbcuts+4]->SetTitle("C");
    histopassedcuts[source*nbcuts+5]->SetTitle("C");
    histopassedcuts[source*nbcuts+6]->SetTitle("C");
    histopassedcuts[source*nbcuts+7]->SetTitle("C");
    histopassedcuts[source*nbcuts+8]->SetTitle("C");
    histopassedcuts[source*nbcuts+9]->SetTitle("C");
    //histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+1]->SetStats(0);
    histopassedcuts[source*nbcuts+2]->SetStats(0);
    histopassedcuts[source*nbcuts+3]->SetStats(0);
    histopassedcuts[source*nbcuts+4]->SetStats(0);
    histopassedcuts[source*nbcuts+5]->SetStats(0);
    histopassedcuts[source*nbcuts+6]->SetStats(0);
    histopassedcuts[source*nbcuts+7]->SetStats(0);
    histopassedcuts[source*nbcuts+8]->SetStats(0);
    histopassedcuts[source*nbcuts+9]->SetStats(0);
    //histopassedcuts[source*nbcuts+0]->Draw();
    //histopassedcuts[source*nbcuts+1]->Draw();
    histopassedcuts[source*nbcuts+2]->Draw();
    histopassedcuts[source*nbcuts+3]->Draw("same");
    //histopassedcuts[source*nbcuts+4]->Draw("same");
    histopassedcuts[source*nbcuts+5]->Draw("same");
    histopassedcuts[source*nbcuts+6]->Draw("same");
    //histopassedcuts[source*nbcuts+7]->Draw("same");
    histopassedcuts[source*nbcuts+8]->Draw("same");
    histopassedcuts[source*nbcuts+9]->Draw("same");
    TLegend *lege = new TLegend(0.4,0.6,0.89,0.89);
    //lege->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
    //lege->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
    //lege->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
    //lege->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
    lege->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
    lege->Draw("same");

    //////////////////////
    // Input
    //////////////////////

    TCanvas * chsSparseMCein =new TCanvas("hsSparseMCeinput","hsSparseMCeinput",800,800);
    chsSparseMCein->cd(1);
    Double_t nbGamma = 0.0;
    source = 1;
    nbGamma = histopassedcuts[source*nbcuts+0]->GetEntries();
    histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+0]->Draw();
    TLegend *leginput = new TLegend(0.4,0.6,0.89,0.89);
    leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#gamma","p");
    Double_t nbPi0 = 0.0;
    source = 2;
    nbPi0 = histopassedcuts[source*nbcuts+0]->GetEntries();
    histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+0]->Draw("same");
    leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#pi^{0}","p");
    Double_t nbEta = 0.0;
    source = 3;
    nbEta = histopassedcuts[source*nbcuts+0]->GetEntries();
    histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+0]->Draw("same");
    leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#eta","p");
    Double_t nbC = 0.0;
    source = 4;
    nbC = histopassedcuts[source*nbcuts+0]->GetEntries();
    histopassedcuts[source*nbcuts+0]->SetStats(0);
    histopassedcuts[source*nbcuts+0]->Draw("same");
    leginput->AddEntry(histopassedcuts[source*nbcuts+0],"c","p");
    leginput->Draw("same");

    //printf("Gamma %f, pi^{0} %f and #eta %f, c %f\n",nbGamma,nbPi0,nbEta,nbC);

    //////////////////////
    // Tracked
    //////////////////////

    TCanvas * cTracked = new TCanvas("cTracked","cTracked",800,800);
    cTracked->cd(1);
    source = 1;
    histopassedcuts[source*nbcuts+1]->Draw();
    TLegend *legTracked = new TLegend(0.4,0.6,0.89,0.89);
    legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#gamma","p");
    source = 2;
    histopassedcuts[source*nbcuts+1]->Draw("same");
    legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#pi^{0}","p");
    legTracked->Draw("same");
  
  }

  /////////////////////////////////////
  // Data Radius and chi2Ndf if AliKF
  ////////////////////////////////////
  
  TH1F *hDataRadius = dynamic_cast<TH1F *>(fList->FindObject("DataRadius")); 
  TH1F *hDataChi2Ndf = dynamic_cast<TH1F *>(fList->FindObject("DataChi2Ndf")); 

  if(hDataRadius || hDataChi2Ndf) {
    TCanvas * cDataRadiusChi2Ndf =new TCanvas("CanvasDataRadiusChi2Ndf","CanvasDataRadiusChi2Ndf",800,800);
    cDataRadiusChi2Ndf->Divide(2,1);
    cDataRadiusChi2Ndf->cd(1);
    if(hDataRadius) hDataRadius->Draw();
    cDataRadiusChi2Ndf->cd(2);
    if(hDataChi2Ndf) hDataChi2Ndf->Draw();
  }

  ///////////////////////
  // Data DCA
  //////////////////////
 
  TH1F *hDataDCA = dynamic_cast<TH1F *>(fList->FindObject("DataDCA")); 
  
  if(hDataDCA) {
    TCanvas * cDataDCA =new TCanvas("CanvasDataDCA","CanvasDataDCA",800,800);
    cDataDCA->cd(1);
    hDataDCA->Draw();
  }

  /////////////////////////////////////
  // MC Radius and chi2Ndf if AliKF
  ////////////////////////////////////
  
  TH2F *hMCRadius = dynamic_cast<TH2F *>(fList->FindObject("MCRadius")); 
  TH2F *hMCChi2Ndf = dynamic_cast<TH2F *>(fList->FindObject("MCChi2Ndf")); 

  if(hMCRadius || hMCChi2Ndf) {
    TCanvas * cMCRadiusChi2Ndf =new TCanvas("CanvasMCRadiusChi2Ndf","CanvasMCRadiusChi2Ndf",800,800);
    cMCRadiusChi2Ndf->Divide(2,1);
    cMCRadiusChi2Ndf->cd(1);
    //TH1D *hMCRadiusBackground = hMCRadius->ProjectionX("MCRadiusBackGround",1,1,"e");
    TH1D *hMCRadiusGamma = hMCRadius->ProjectionX("MCRadiusGamma",2,2,"e");
    TH1D *hMCRadiusPi0 = hMCRadius->ProjectionX("MCRadiusPi0",3,3,"e");
    TH1D *hMCRadiusEta = hMCRadius->ProjectionX("MCRadiusEta",4,4,"e");
    TH1D *hMCRadiusC = hMCRadius->ProjectionX("MCRadiusC",5,5,"e");
    TH1D *hMCRadiusB = hMCRadius->ProjectionX("MCRadiusB",6,6,"e");
    //hMCRadiusBackground->Draw();
    hMCRadiusGamma->Draw();
    hMCRadiusPi0->Draw("same");
    hMCRadiusEta->Draw("same");
    hMCRadiusC->Draw("same");
    hMCRadiusB->Draw("same");
    TLegend *legRadius = new TLegend(0.4,0.6,0.89,0.89);
    //legRadius->AddEntry(hMCRadiusBackground,"Background","p");
    legRadius->AddEntry(hMCRadiusGamma,"#gamma","p");
    legRadius->AddEntry(hMCRadiusPi0,"#pi^{0}","p");
    legRadius->AddEntry(hMCRadiusEta,"#eta","p");
    legRadius->AddEntry(hMCRadiusC,"c","p");
    legRadius->AddEntry(hMCRadiusB,"b","p");
    legRadius->Draw("same");
    cMCRadiusChi2Ndf->cd(2);
    //TH1D *hMCChi2NdfBackground = hMCChi2Ndf->ProjectionX("MCChi2NdfBackGround",1,1,"e");
    TH1D *hMCChi2NdfGamma = hMCChi2Ndf->ProjectionX("MCChi2NdfGamma",2,2,"e");
    TH1D *hMCChi2NdfPi0 = hMCChi2Ndf->ProjectionX("MCChi2NdfPi0",3,3,"e");
    TH1D *hMCChi2NdfEta = hMCChi2Ndf->ProjectionX("MCChi2NdfEta",4,4,"e");
    TH1D *hMCChi2NdfC = hMCChi2Ndf->ProjectionX("MCChi2NdfC",5,5,"e");
    TH1D *hMCChi2NdfB = hMCChi2Ndf->ProjectionX("MCChi2NdfB",6,6,"e");
    //hMCChi2NdfBackground->Draw();
    hMCChi2NdfGamma->Draw();
    hMCChi2NdfPi0->Draw("same");
    hMCChi2NdfEta->Draw("same");
    hMCChi2NdfC->Draw("same");
    hMCChi2NdfB->Draw("same");
    TLegend *legChi2Ndf = new TLegend(0.4,0.6,0.89,0.89);
    //legChi2Ndf->AddEntry(hMCChi2NdfBackground,"Background","p");
    legChi2Ndf->AddEntry(hMCChi2NdfGamma,"#gamma","p");
    legChi2Ndf->AddEntry(hMCChi2NdfPi0,"#pi^{0}","p");
    legChi2Ndf->AddEntry(hMCChi2NdfEta,"#eta","p");
    legChi2Ndf->AddEntry(hMCChi2NdfC,"c","p");
    legChi2Ndf->AddEntry(hMCChi2NdfB,"b","p");
    legChi2Ndf->Draw("same");
  }

  ///////////////////////
  // MC DCA
  //////////////////////
  
  TH2F *hMCDCA = dynamic_cast<TH2F *>(fList->FindObject("MCDCA")); 
  
  if(hMCDCA) {
    TCanvas * cMCDCA =new TCanvas("CanvasMCDCA","CanvasMCDCA",800,800);
    cMCDCA->cd(1);
    //TH1D *hMCDCABackground = hMCDCA->ProjectionX("MCDCABackGround",1,1,"e");
    TH1D *hMCDCAGamma = hMCDCA->ProjectionX("MCDCAGamma",2,2,"e");
    TH1D *hMCDCAPi0 = hMCDCA->ProjectionX("MCDCAPi0",3,3,"e");
    TH1D *hMCDCAEta = hMCDCA->ProjectionX("MCDCAEta",4,4,"e");
    TH1D *hMCDCAC = hMCDCA->ProjectionX("MCDCAC",5,5,"e");
    TH1D *hMCDCAB = hMCDCA->ProjectionX("MCDCAB",6,6,"e");
    //hMCDCABackground->Draw();
    hMCDCAGamma->Draw();
    hMCDCAPi0->Draw("same");
    hMCDCAEta->Draw("same");
    hMCDCAC->Draw("same");
    hMCDCAB->Draw("same");
    TLegend *legDCA = new TLegend(0.4,0.6,0.89,0.89);
    //legDCA->AddEntry(hMCDCABackground,"Background","p");
    legDCA->AddEntry(hMCDCAGamma,"#gamma","p");
    legDCA->AddEntry(hMCDCAPi0,"#pi^{0}","p");
    legDCA->AddEntry(hMCDCAEta,"#eta","p");
    legDCA->AddEntry(hMCDCAC,"c","p");
    legDCA->AddEntry(hMCDCAB,"b","p");
    legDCA->Draw("same");
  }


  /////////////////////////
  // PID Partner Signal
  /////////////////////////
  TF1 *betheBlochElectron = 0x0;
  TF1 *betheBlochMuon = 0x0;
  TF1 *betheBlochPion = 0x0;
  TF1 *betheBlochKaon = 0x0;
  TF1 *betheBlochProton = 0x0;

  TH2F *hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner0")); 
  TH2F *hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("TPCPartner1"));
  if((!hsignalPidPartner0) && (!hsignalPidPartner1)) {
    hsignalPidPartner0 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner0")); 
    hsignalPidPartner1 = dynamic_cast<TH2F *>(fList->FindObject("ITSPartner1"));
    
    betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronITS,0.1,10.0,0);
    betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonITS,0.1,10.0,0);
    betheBlochPion = new TF1("betheBlochPion",BetheBlochPionITS,0.1,10.0,0);
    betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonITS,0.1,10.0,0);
    betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonITS,0.1,10.0,0);

  }
  else {

    betheBlochElectron = new TF1("betheBlochElectron",BetheBlochElectronTPC,0.1,10.0,0);
    betheBlochMuon = new TF1("betheBlochMuon",BetheBlochMuonTPC,0.1,10.0,0);
    betheBlochPion = new TF1("betheBlochPion",BetheBlochPionTPC,0.1,10.0,0);
    betheBlochKaon = new TF1("betheBlochKaon",BetheBlochKaonTPC,0.1,10.0,0);
    betheBlochProton = new TF1("betheBlochProton",BetheBlochProtonTPC,0.1,10.0,0);

  }


  if((hsignalPidPartner0) || (hsignalPidPartner1)) {
    TCanvas * cPidSignal =new TCanvas("cPidSignal","cPidSignal",800,800);
    cPidSignal->Divide(2,1);
    cPidSignal->cd(1);
    if(hsignalPidPartner0) hsignalPidPartner0->Draw("colz");
    if(betheBlochElectron) betheBlochElectron->Draw("same");
    if(betheBlochMuon) betheBlochMuon->Draw("same");
    if(betheBlochPion) betheBlochPion->Draw("same");
    if(betheBlochKaon) betheBlochKaon->Draw("same");
    if(betheBlochProton) betheBlochProton->Draw("same");
    cPidSignal->cd(2);
    if(hsignalPidPartner1) hsignalPidPartner1->Draw("colz");
    if(betheBlochElectron) betheBlochElectron->Draw("same");
    if(betheBlochMuon) betheBlochMuon->Draw("same");
    if(betheBlochPion) betheBlochPion->Draw("same");
    if(betheBlochKaon) betheBlochKaon->Draw("same");
    if(betheBlochProton) betheBlochProton->Draw("same");
  }
  
  THnSparseF *hsSparseITSsignal = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignal"));
  if(hsSparseITSsignal) {

   
    TH2D *sddsdd = hsSparseITSsignal->Projection(1,2);
    TH2D *ssdssd = hsSparseITSsignal->Projection(3,4);
    TH2D *sddssda = hsSparseITSsignal->Projection(1,3);    
    TH2D *sddssdb = hsSparseITSsignal->Projection(2,4);    
    TH2D *sddssdc = hsSparseITSsignal->Projection(1,4);   
    TH2D *sddssdd = hsSparseITSsignal->Projection(2,3);    

    TCanvas * cITSSignal =new TCanvas("cITSSignal","cITSSignal",800,800);
    cITSSignal->Divide(2,3);
    cITSSignal->cd(1);
    sddsdd->Draw("colz");
    cITSSignal->cd(2);
    ssdssd->Draw("colz");
    cITSSignal->cd(3);
    sddssda->Draw("colz");
    cITSSignal->cd(4);
    sddssdb->Draw("colz");
    cITSSignal->cd(5);
    sddssdc->Draw("colz");
    cITSSignal->cd(6);
    sddssdd->Draw("colz");
  
  } 

  THnSparseF *hsSparseITSsignalSplit = dynamic_cast<THnSparseF *>(fList->FindObject("SparseITSsignalSplit"));
  if(hsSparseITSsignalSplit) {

    // no splitted
    hsSparseITSsignalSplit->GetAxis(0)->SetRange(1,1);
    
    TH1D *layerITS2 = hsSparseITSsignalSplit->Projection(1);
    TH1D *layerITS3 = hsSparseITSsignalSplit->Projection(2);
    TH1D *layerITS4 = hsSparseITSsignalSplit->Projection(3);
    TH1D *layerITS5 = hsSparseITSsignalSplit->Projection(4);

    // splitted
    hsSparseITSsignalSplit->GetAxis(0)->SetRange(2,2);
    
    TH1D *layerITS2s = hsSparseITSsignalSplit->Projection(1);
    TH1D *layerITS3s = hsSparseITSsignalSplit->Projection(2);
    TH1D *layerITS4s = hsSparseITSsignalSplit->Projection(3);
    TH1D *layerITS5s = hsSparseITSsignalSplit->Projection(4);
    
    TCanvas * cITSSignalSplit =new TCanvas("cITSSignalSplit","cITSSignalSplit",800,800);
    cITSSignalSplit->Divide(2,2);
    cITSSignalSplit->cd(1);
    layerITS2->Draw();
    layerITS2s->Draw("same");
    TLegend *legITS2 = new TLegend(0.4,0.6,0.89,0.89);
    legITS2->AddEntry(layerITS2,"No splitted","p");
    legITS2->AddEntry(layerITS2s,"Splitted","p");
    legITS2->Draw("same");
    cITSSignalSplit->cd(2);
    layerITS3->Draw();
    layerITS3s->Draw("same");
    TLegend *legITS3 = new TLegend(0.4,0.6,0.89,0.89);
    legITS3->AddEntry(layerITS3,"No splitted","p");
    legITS3->AddEntry(layerITS3s,"Splitted","p");
    legITS3->Draw("same");
    cITSSignalSplit->cd(3);
    layerITS4->Draw();
    layerITS4s->Draw("same");
    TLegend *legITS4 = new TLegend(0.4,0.6,0.89,0.89);
    legITS4->AddEntry(layerITS4,"No splitted","p");
    legITS4->AddEntry(layerITS4s,"Splitted","p");
    legITS4->Draw("same");
    cITSSignalSplit->cd(4);
    layerITS5->Draw();
    layerITS5s->Draw("same");
    TLegend *legITS5 = new TLegend(0.4,0.6,0.89,0.89);
    legITS5->AddEntry(layerITS5,"No splitted","p");
    legITS5->AddEntry(layerITS5s,"Splitted","p");
    legITS5->Draw("same");

  
  }  

}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochElectronITS(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for ITS
  //
  static AliITSPIDResponse itsPidResponse;
  return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(0));
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochMuonITS(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for ITS
  //
  static AliITSPIDResponse itsPidResponse;
  return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(1));
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochPionITS(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for ITS
  //
  static AliITSPIDResponse itsPidResponse;
  return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(2));
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochKaonITS(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for ITS
  //
  static AliITSPIDResponse itsPidResponse;
  return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(3));
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochProtonITS(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for ITS
  //
  static AliITSPIDResponse itsPidResponse;
  return itsPidResponse.Bethe(x[0],AliPID::ParticleMass(4));
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochElectronTPC(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for TPC
  //
  static AliTPCPIDResponse tpcPidResponse;
  return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kElectron);
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochMuonTPC(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for TPC
  //
  static AliTPCPIDResponse tpcPidResponse;
  return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kMuon);
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochPionTPC(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for TPC
  //
  static AliTPCPIDResponse tpcPidResponse;
  return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kPion);
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochKaonTPC(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for TPC
  //
  static AliTPCPIDResponse tpcPidResponse;
  return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kKaon);
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochProtonTPC(const Double_t *x, const Double_t * /*par*/) 
{
  //
  // Bethe Bloch for TPC
  //
  static AliTPCPIDResponse tpcPidResponse;
  return tpcPidResponse.GetExpectedSignal(x[0],AliPID::kProton);
}