/**************************************************************************
 * Copyright(c) 1998-2008, 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.                  *
 **************************************************************************/

/* $Id$ */

//*************************************************************************
// Class AliAnalysisTaskSEDplus
// AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and 
//comparison of heavy-flavour decay candidates
// to MC truth (kinematics stored in the AOD)
// Authors: Renu Bala, bala@to.infn.it
// F. Prino, prino@to.infn.it
// G. Ortona, ortona@to.infn.it
/////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TNtuple.h>
#include <TCanvas.h>
#include <TList.h>
#include <TString.h>
#include <TDatabasePDG.h>

#include "AliAnalysisManager.h"
#include "AliRDHFCutsDplustoKpipi.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliAODRecoDecayHF3Prong.h"
#include "AliAnalysisVertexingHF.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSEDplus.h"
#include "AliNormalizationCounter.h"
#include "AliVertexingHFUtils.h"
ClassImp(AliAnalysisTaskSEDplus)

//________________________________________________________________________
AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
AliAnalysisTaskSE(),
  fOutput(0), 
  fHistNEvents(0),
  fPtVsMass(0),
  fPtVsMassTC(0),
  fYVsPt(0),
  fYVsPtTC(0),
  fYVsPtSig(0),
  fYVsPtSigTC(0),
  fPhiEtaCand(0),
  fPhiEtaCandSigReg(0),
  fSPDMult(0),
  fNtupleDplus(0),
  fUpmasslimit(1.965),
  fLowmasslimit(1.765),
  fNPtBins(0),
  fBinWidth(0.002),
  fListCuts(0),
  fRDCutsAnalysis(0),
  fCounter(0),
  fFillNtuple(kFALSE),
  fReadMC(kFALSE),
  fUseStrangeness(kFALSE),
  fUseBit(kTRUE),
  fCutsDistr(kFALSE),
  fDoImpPar(kFALSE),
  fNImpParBins(400),
  fLowerImpPar(-1000.),
  fHigherImpPar(1000.),
  fDoLS(0),
  fEtaSelection(0),
  fSystem(0)
{
  // Default constructor

  for(Int_t i=0;i<3;i++){
    fHistCentrality[i]=0;
    fCorreld0Kd0pi[i]=0;
  }
  
  for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0;
    
  
  for(Int_t i=0;i<3*kMaxPtBins;i++){
    fMassHist[i]=0;
    fCosPHist[i]=0;
    fDLenHist[i]=0;
    fSumd02Hist[i]=0;
    fSigVertHist[i]=0;
    fPtMaxHist[i]=0;
    fPtKHist[i]=0;
    fPtpi1Hist[i]=0;
    fPtpi2Hist[i]=0;
    fDCAHist[i]=0;
    fMassHistTC[i]=0;
    fMassHistTCPlus[i]=0;
    fMassHistTCMinus[i]=0;
    
    fDLxy[i]=0;
    fDLxyTC[i]=0;
    fCosxy[i]=0;
    fCosxyTC[i]=0;
    fCosPHistLS[i]=0;
    fDLenHistLS[i]=0;
    fSumd02HistLS[i]=0;
    fSigVertHistLS[i]=0;
    fPtMaxHistLS[i]=0;
    fDCAHistLS[i]=0;
  }
  for(Int_t i=0;i<5*kMaxPtBins;i++){
    fMassHistLS[i]=0;
    fMassHistLSTC[i]=0;
  }
  for(Int_t i=0;i<kMaxPtBins+1;i++){
    fArrayBinLimits[i]=0;
  }

}

//________________________________________________________________________
AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,Bool_t fillNtuple):
  AliAnalysisTaskSE(name),
  fOutput(0),
  fHistNEvents(0),
  fPtVsMass(0),
  fPtVsMassTC(0),
  fYVsPt(0),
  fYVsPtTC(0),
  fYVsPtSig(0),
  fYVsPtSigTC(0),
  fPhiEtaCand(0),
  fPhiEtaCandSigReg(0),
  fSPDMult(0),
  fNtupleDplus(0),
  fUpmasslimit(1.965),
  fLowmasslimit(1.765),
  fNPtBins(0),
  fBinWidth(0.002),
  fListCuts(0),
  fRDCutsAnalysis(dpluscutsana),
  fCounter(0),
  fFillNtuple(fillNtuple),
  fReadMC(kFALSE),
  fUseStrangeness(kFALSE),
  fUseBit(kTRUE),
  fCutsDistr(kFALSE),
  fDoImpPar(kFALSE),
  fNImpParBins(400),
  fLowerImpPar(-1000.),
  fHigherImpPar(1000.),
  fDoLS(0),
  fEtaSelection(0),
  fSystem(0)
{
  // 
  // Standrd constructor
  //
  fNPtBins=fRDCutsAnalysis->GetNPtBins();
  
  for(Int_t i=0;i<3;i++){
    fHistCentrality[i]=0;
    fCorreld0Kd0pi[i]=0;
  }
  
  for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0;
    
  for(Int_t i=0;i<3*kMaxPtBins;i++){
    fMassHist[i]=0;
    fCosPHist[i]=0;
    fDLenHist[i]=0;
    fSumd02Hist[i]=0;
    fSigVertHist[i]=0;
    fPtMaxHist[i]=0;
    fPtKHist[i]=0;
    fPtpi1Hist[i]=0;
    fPtpi2Hist[i]=0;
    fDCAHist[i]=0;
    fMassHistTC[i]=0;
    fMassHistTCPlus[i]=0;
    fMassHistTCMinus[i]=0;
    
    fDLxy[i]=0;
    fDLxyTC[i]=0;
    fCosxy[i]=0;
    fCosxyTC[i]=0;
    fCosPHistLS[i]=0;
    fDLenHistLS[i]=0;
    fSumd02HistLS[i]=0;
    fSigVertHistLS[i]=0;
    fPtMaxHistLS[i]=0;
    fDCAHistLS[i]=0;
  }
  for(Int_t i=0;i<5*kMaxPtBins;i++){
    fMassHistLS[i]=0;
    fMassHistLSTC[i]=0;
  }
  for(Int_t i=0;i<kMaxPtBins+1;i++){
    fArrayBinLimits[i]=0;
  }
  
  
  // Default constructor
  // Output slot #1 writes into a TList container
  DefineOutput(1,TList::Class());  //My private output
  // Output slot #2 writes cut to private output
  //  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
  DefineOutput(2,TList::Class());
  // Output slot #3 writes cut to private output
  DefineOutput(3,AliNormalizationCounter::Class());
  
  if(fFillNtuple){
    // Output slot #4 writes into a TNtuple container
    DefineOutput(4,TNtuple::Class());  //My private output
  }
}

//________________________________________________________________________
AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
{
  //
  // Destructor
  //
  if(fOutput && !fOutput->IsOwner()){
    delete fHistNEvents;
    for(Int_t i=0;i<3*fNPtBins;i++){    
      delete fMassHist[i];
      delete fCosPHist[i];
      delete fDLenHist[i];
      delete fSumd02Hist[i];
      delete fSigVertHist[i];
      delete fPtMaxHist[i];
      delete fPtKHist[i];
      delete fPtpi1Hist[i];
      delete fPtpi2Hist[i];
      delete fDCAHist[i];
      delete fDLxy[i];
      delete fDLxyTC[i];
      delete fCosxy[i];
      delete fCosxyTC[i];
      delete fMassHistTC[i];
      delete fMassHistTCPlus[i];
      delete fMassHistTCMinus[i];
      delete fCosPHistLS[i];
      delete fDLenHistLS[i];
      delete fSumd02HistLS[i];
      delete fSigVertHistLS[i];
      delete fPtMaxHistLS[i];
      delete fDCAHistLS[i];
    }
    for(Int_t i=0;i<5*fNPtBins;i++){
      delete fMassHistLS[i];
      delete fMassHistLSTC[i];
    }
    for(Int_t i=0;i<3;i++){
      delete fCorreld0Kd0pi[i];
      delete fHistCentrality[i];
    }
    for(Int_t i=0;i<5;i++){
      delete fHistMassPtImpParTC[i];
    }
    delete fPtVsMass;
    delete fPtVsMassTC;
    delete fYVsPt;
    delete fYVsPtTC;
    delete fYVsPtSig;
    delete fYVsPtSigTC;
    delete fPhiEtaCand;
    delete fPhiEtaCandSigReg;
    delete fSPDMult;
  }

  delete fOutput;  
  delete fNtupleDplus;
  delete fListCuts;
  delete fRDCutsAnalysis;
  delete fCounter;

}  
//_________________________________________________________________
void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
  // set invariant mass limits
  Float_t bw=GetBinWidth();
  fUpmasslimit = 1.865+range;
  fLowmasslimit = 1.865-range;
  SetBinWidth(bw);
}
//_________________________________________________________________
void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
  // set invariant mass limits
  if(uplimit>lowlimit)
    {
      Float_t bw=GetBinWidth();
      fUpmasslimit = lowlimit;
      fLowmasslimit = uplimit;
      SetBinWidth(bw);
    }
}
//________________________________________________________________
void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
  Float_t width=w;
  Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
  Int_t missingbins=4-nbins%4;
  nbins=nbins+missingbins;
  width=(fUpmasslimit-fLowmasslimit)/nbins;
  if(missingbins!=0){
    printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
  }
  else{
    if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
  }
  fBinWidth=width;
}
//_________________________________________________________________
Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
  return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
}
//_________________________________________________________________
void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
  //
  //
  // Fill the Like Sign histograms
  //
  if(fDebug>1)printf("started LS\n");
  
  //histograms for like sign
  Int_t nbins=GetNBinsHistos();;
  TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
  TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
  
  Int_t nPosTrks=0,nNegTrks=0;
  Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
  Int_t nDplusLS=0;
  Int_t nDminusLS=0;
  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
  Int_t index=0; 
  
  // loop over like sign candidates
  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
    if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
    Bool_t unsetvtx=kFALSE;
    if(!d->GetOwnPrimaryVtx()) {
      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
      unsetvtx=kTRUE;
    }
    Double_t ptCand = d->Pt();
    Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
    if(iPtBin<0)continue;
    Int_t sign= d->GetCharge();
    if(sign>0){
      nPosTrks++;
    }else{
      nNegTrks++;
    }
    fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
    Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
    if(passTightCuts>0){
      Float_t invMass = d->InvMassDplus();
      index=GetLSHistoIndex(iPtBin);
      fMassHistLS[index+1]->Fill(invMass);
      if(sign>0){
	histLSPlus->Fill(invMass);
	nDplusLS++;
      }else{
	histLSMinus->Fill(invMass);
	nDminusLS++;
      }
      if(fCutsDistr){
	Double_t dlen=d->DecayLength();
	Double_t cosp=d->CosPointingAngle();
	Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
	Double_t dca=d->GetDCA();   
	Double_t sigvert=d->GetSigmaVert();   
	Double_t ptmax=0;
	for(Int_t i=0;i<3;i++){
	  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
	}
	fCosPHistLS[iPtBin]->Fill(cosp);
	fDLenHistLS[iPtBin]->Fill(dlen);
	fSumd02HistLS[iPtBin]->Fill(sumD02);
	fSigVertHistLS[iPtBin]->Fill(sigvert);
	fPtMaxHistLS[iPtBin]->Fill(ptmax);
	fDCAHistLS[iPtBin]->Fill(dca);
      }
    }
    if(unsetvtx) d->UnsetOwnPrimaryVtx();
  }
  //wei:normalized to the number of combinations (production)
  //wei2:normalized to the number of  LS/OS (production)
  //wei3:normalized to the number of  LS/OS (analysis)
  //wei4:normalized to the number of combinations (analysis)
  Float_t wei2=0;
  if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
  Float_t wei3=0;
  if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
  Float_t weiplus=1.,weiminus=1.;
  Float_t wei4plus=1.,wei4minus=1.;
  //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
  if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
  if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
  if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
  if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);

  fMassHistLS[index]->Add(histLSPlus,weiplus);
  fMassHistLS[index]->Add(histLSMinus,weiminus);
  fMassHistLS[index+2]->Add(histLSPlus,wei2);
  fMassHistLS[index+2]->Add(histLSMinus,wei2);
  fMassHistLS[index+3]->Add(histLSPlus,wei3);
  fMassHistLS[index+3]->Add(histLSMinus,wei3);
  fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
  fMassHistLS[index+4]->Add(histLSMinus,wei4minus);

  delete histLSPlus;histLSPlus=0;
  delete histLSMinus;histLSMinus=0;
  
  if(fDebug>1) printf("LS analysis terminated\n");  
}

//__________________________________________
void AliAnalysisTaskSEDplus::Init(){
  //
  // Initialization
  //
  if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
  
  //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
  fListCuts=new TList();
  fListCuts->SetOwner();
  AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
  analysis->SetName("AnalysisCuts");
  
  fListCuts->Add(analysis);
  PostData(2,fListCuts);
  
  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
{
  // Create the output container
  //
  if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList();
  fOutput->SetOwner();
  fOutput->SetName("OutputHistos");

  fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
  fHistNEvents->GetXaxis()->SetBinLabel(3,"Rejected due to trigger");
  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected pileup events");
  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to centrality"); 
  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vtxz");
  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to Physics Sel");
  fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
  fHistNEvents->GetXaxis()->SetBinLabel(10,"D+ after loose cuts");
  fHistNEvents->GetXaxis()->SetBinLabel(11,"D+ after tight cuts");
 
  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);  
  fHistNEvents->Sumw2();
  fHistNEvents->SetMinimum(0);
  fOutput->Add(fHistNEvents);

  TString hisname;
  Int_t index=0;
  Int_t nbins=GetNBinsHistos();
  fHistCentrality[0]=new TH2F("hCentrMult","centrality",100,0.5,30000.5,40,0.,100.);
  fHistCentrality[1]=new TH2F("hCentrMult(selectedCent)","centrality(selectedCent)",100,0.5,30000.5,40,0.,100.);
  fHistCentrality[2]=new TH2F("hCentrMult(OutofCent)","centrality(OutofCent)",100,0.5,30000.5,40,0.,100.);
  for(Int_t i=0;i<3;i++){
    fHistCentrality[i]->Sumw2();
    fOutput->Add(fHistCentrality[i]);
  }
  for(Int_t i=0;i<fNPtBins;i++){

    index=GetHistoIndex(i);

    hisname.Form("hMassPt%d",i);
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHist[index]->Sumw2();
    hisname.Form("hCosPAllPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenAllPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hSumd02AllPt%d",i);
    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
    fSumd02Hist[index]->Sumw2();
    hisname.Form("hSigVertAllPt%d",i);
    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fSigVertHist[index]->Sumw2();
    hisname.Form("hPtMaxAllPt%d",i);
    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
    fPtMaxHist[index]->Sumw2();
    hisname.Form("hPtKPt%d",i);
    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtKHist[index]->Sumw2();
    hisname.Form("hPtpi1Pt%d",i);
    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi1Hist[index]->Sumw2();
    hisname.Form("hPtpi2Pt%d",i);
    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi2Hist[index]->Sumw2();
    hisname.Form("hDCAAllPt%d",i);
    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fDCAHist[index]->Sumw2();

    hisname.Form("hDLxyPt%d",i);
    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxy[index]->Sumw2();
    hisname.Form("hCosxyPt%d",i);
    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxy[index]->Sumw2();
    hisname.Form("hDLxyPt%dTC",i);
    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxyTC[index]->Sumw2();
    hisname.Form("hCosxyPt%dTC",i);
    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxyTC[index]->Sumw2();
   
    hisname.Form("hMassPt%dTC",i);
    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTC[index]->Sumw2();
    hisname.Form("hMassPt%dTCPlus",i);
    fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCPlus[index]->Sumw2();
    hisname.Form("hMassPt%dTCMinus",i);
    fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCMinus[index]->Sumw2();


       
    index=GetSignalHistoIndex(i);    
    hisname.Form("hSigPt%d",i);
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHist[index]->Sumw2();
    hisname.Form("hCosPSigPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenSigPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hSumd02SigPt%d",i);
    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
    fSumd02Hist[index]->Sumw2();
    hisname.Form("hSigVertSigPt%d",i);
    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fSigVertHist[index]->Sumw2();
    hisname.Form("hPtMaxSigPt%d",i);
    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
    fPtMaxHist[index]->Sumw2();  
    hisname.Form("hPtKSigPt%d",i);  
    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtKHist[index]->Sumw2();
    hisname.Form("hPtpi1SigPt%d",i);
    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi1Hist[index]->Sumw2();
    hisname.Form("hPtpi2SigPt%d",i);
    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi2Hist[index]->Sumw2();

    hisname.Form("hDCASigPt%d",i);
    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fDCAHist[index]->Sumw2();    

    hisname.Form("hDLxySigPt%d",i);
    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxy[index]->Sumw2();
    hisname.Form("hCosxySigPt%d",i);
    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxy[index]->Sumw2();
    hisname.Form("hDLxySigPt%dTC",i);
    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxyTC[index]->Sumw2();
    hisname.Form("hCosxySigPt%dTC",i);
    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxyTC[index]->Sumw2();
    hisname.Form("hSigPt%dTC",i);
    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTC[index]->Sumw2();
    hisname.Form("hSigPt%dTCPlus",i);
    fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCPlus[index]->Sumw2();
    hisname.Form("hSigPt%dTCMinus",i);
    fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCMinus[index]->Sumw2();


    index=GetBackgroundHistoIndex(i); 
    hisname.Form("hBkgPt%d",i);
    fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHist[index]->Sumw2();
    hisname.Form("hCosPBkgPt%d",i);
    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
    fCosPHist[index]->Sumw2();
    hisname.Form("hDLenBkgPt%d",i);
    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
    fDLenHist[index]->Sumw2();
    hisname.Form("hSumd02BkgPt%d",i);
    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
    fSumd02Hist[index]->Sumw2();
    hisname.Form("hSigVertBkgPt%d",i);
    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fSigVertHist[index]->Sumw2();
    hisname.Form("hPtMaxBkgPt%d",i);
    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
    fPtMaxHist[index]->Sumw2();
    hisname.Form("hPtKBkgPt%d",i);  
    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtKHist[index]->Sumw2();
    hisname.Form("hPtpi1BkgPt%d",i);
    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi1Hist[index]->Sumw2();
    hisname.Form("hPtpi2BkgPt%d",i);
    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
    fPtpi2Hist[index]->Sumw2();
    hisname.Form("hDCABkgPt%d",i);
    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
    fDCAHist[index]->Sumw2();

    hisname.Form("hDLxyBkgPt%d",i);
    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxy[index]->Sumw2();
    hisname.Form("hCosxyBkgPt%d",i);
    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxy[index]->Sumw2();
    hisname.Form("hDLxyBkgPt%dTC",i);
    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
    fDLxyTC[index]->Sumw2();
    hisname.Form("hCosxyBkgPt%dTC",i);
    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
    fCosxyTC[index]->Sumw2();
  

    hisname.Form("hBkgPt%dTC",i);
    fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTC[index]->Sumw2();
    hisname.Form("hBkgPt%dTCPlus",i);
    fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCPlus[index]->Sumw2();
    hisname.Form("hBkgPt%dTCMinus",i);
    fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistTCMinus[index]->Sumw2();
  }
    

  for(Int_t i=0; i<3*fNPtBins; i++){
    fOutput->Add(fMassHist[i]);
    if(fCutsDistr){  
      fOutput->Add(fCosPHist[i]);
      fOutput->Add(fDLenHist[i]);
      fOutput->Add(fSumd02Hist[i]);
      fOutput->Add(fSigVertHist[i]);
      fOutput->Add(fPtMaxHist[i]);
      fOutput->Add(fPtKHist[i]);
      fOutput->Add(fPtpi1Hist[i]);
      fOutput->Add(fPtpi2Hist[i]);
      fOutput->Add(fDCAHist[i]);
      fOutput->Add(fDLxy[i]);
      fOutput->Add(fDLxyTC[i]);
      fOutput->Add(fCosxy[i]);
      fOutput->Add(fCosxyTC[i]);
    }  
    fOutput->Add(fMassHistTC[i]);
    fOutput->Add(fMassHistTCPlus[i]);
    fOutput->Add(fMassHistTCMinus[i]);

  }
  
  fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
  fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
  fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
  if(fCutsDistr){  
    for(Int_t i=0; i<3; i++){
      fCorreld0Kd0pi[i]->Sumw2();
      fOutput->Add(fCorreld0Kd0pi[i]);
    }
  }
  

  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
  fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);  
  fYVsPt=new TH3F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
  fYVsPtTC=new TH3F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
  fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
  fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
  fPhiEtaCand=new TH2F("hPhiEtaCand","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
  fPhiEtaCandSigReg=new TH2F("hPhiEtaCandSigReg","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
  fSPDMult = new TH1F("hSPDMult", "Tracklets multiplicity; Tracklets ; Entries",200,0.,200.); 
  fOutput->Add(fPtVsMass);
  fOutput->Add(fPtVsMassTC);
  fOutput->Add(fYVsPt);
  fOutput->Add(fYVsPtTC);
  fOutput->Add(fYVsPtSig);
  fOutput->Add(fYVsPtSigTC);
  fOutput->Add(fPhiEtaCand);
  fOutput->Add(fPhiEtaCandSigReg);
  fOutput->Add(fSPDMult);


  //Counter for Normalization
  TString normName="NormalizationCounter";
  AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
  if(cont)normName=(TString)cont->GetName();
  fCounter = new AliNormalizationCounter(normName.Data());
  fCounter->Init();

  if(fDoLS) CreateLikeSignHistos();
  if(fDoImpPar) CreateImpactParameterHistos();

  if(fFillNtuple){
    OpenFile(4); // 4 is the slot number of the ntuple
   
    
    fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Pt:charge:piddau0:piddau1:piddau2:Ptpi:PtK:Ptpi2:mompi:momK:mompi2:cosp:cospxy:DecLen:NormDecLen:DecLenXY:NormDecLenXY:InvMass:sigvert:d0Pi:d0K:d0Pi2:maxdca:ntracks:centr:RunNumber:BadDau");

  }
  
  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event:
  // heavy flavor candidates association to MC truth

  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
  
  TClonesArray *array3Prong = 0;
  TClonesArray *arrayLikeSign =0;
  if(!aod && AODEvent() && IsStandardAOD()) {
    // In case there is an AOD handler writing a standard AOD, use the AOD 
    // event in memory rather than the input (ESD) event.    
    aod = dynamic_cast<AliAODEvent*> (AODEvent());
    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
    // have to taken from the AOD event hold by the AliAODExtension
    AliAODHandler* aodHandler = (AliAODHandler*) 
      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
    if(aodHandler->GetExtensions()) {
      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
      AliAODEvent *aodFromExt = ext->GetAOD();
      array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
      arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
    }
  } else if(aod) {
    array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
    arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
  }

  if(!aod || !array3Prong) {
    printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
    return;
  }
  if(!arrayLikeSign && fDoLS) {
    printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
    return;
  }

  // fix for temporary bug in ESDfilter 
  // the AODs with null vertex pointer didn't pass the PhysSel
  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;

  //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
  fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);

  fHistNEvents->Fill(0); // count event
  Int_t runNumber=aod->GetRunNumber();

  //Event selection
  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
  Float_t ntracks=aod->GetNumberOfTracks();
  Float_t evCentr=0;
  if(fRDCutsAnalysis->GetUseCentrality()>0) evCentr=fRDCutsAnalysis->GetCentrality(aod);
  fHistCentrality[0]->Fill(ntracks,evCentr);
  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
  if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3); 
  if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
  if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
  if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);

  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
  //TString trigclass=aod->GetFiredTriggerClasses();
  // Post the data already here  
  PostData(1,fOutput);
  if(!isEvSel)return;
  Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
  // printf("ntracklet===%d\n",tracklets);
  fSPDMult->Fill(tracklets);

  fHistCentrality[1]->Fill(ntracks,evCentr);
  fHistNEvents->Fill(1);

  TClonesArray *arrayMC=0;
  AliAODMCHeader *mcHeader=0;

  // AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  //    vtx1->Print();
  //   TString primTitle = vtx1->GetTitle();
  //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 
  // load MC particles
  if(fReadMC){
    arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
    if(!arrayMC) {
      printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
      return;
    }
    
    // load MC header
    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
    if(!mcHeader) {
      printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
      return;
    }
  }
  
  Int_t n3Prong = array3Prong->GetEntriesFast();
  //  printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
  
  Int_t nOS=0;
  Int_t index;
  Int_t pdgDgDplustoKpipi[3]={321,211,211};

  if(fDoLS>1){//Normalizations for LS
    for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
      AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
      if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
	if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod))nOS++;
      }
    }
  }else{//Standard analysis
    // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE
    //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
    Int_t nSelectedloose=0,nSelectedtight=0;
    for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
      AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
      fHistNEvents->Fill(7);
      if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
	fHistNEvents->Fill(8);
	continue;
      }
    
      Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);

      if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;

      Double_t etaD=d->Eta();
      Double_t phiD=d->Phi();
      if(fEtaSelection!=0){
	if(fEtaSelection==1 && etaD<0) continue;
	if(fEtaSelection==-1 && etaD>0) continue;
      }

      Bool_t unsetvtx=kFALSE;
      if(!d->GetOwnPrimaryVtx()){
	d->SetOwnPrimaryVtx(vtx1);
	unsetvtx=kTRUE;
      }

      Double_t ptCand = d->Pt();
      Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);     
    
      Bool_t recVtx=kFALSE;
      AliAODVertex *origownvtx=0x0;
      if(fRDCutsAnalysis->GetIsPrimaryWithoutDaughters()){
	if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());	
	if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
	else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
      }
    
      Int_t labDp=-1;
      Bool_t isPrimary=kTRUE;
      Float_t pdgCode=-2;
      Float_t trueImpParXY=0.;
      if(fReadMC){
	labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
	if(labDp>=0){
	  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
	  if(AliVertexingHFUtils::CheckOrigin(arrayMC,partDp,kFALSE)==5) isPrimary=kFALSE;
	  pdgCode=TMath::Abs(partDp->GetPdgCode());
	  if(!isPrimary){
	    trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
	  }
	}else{
	  pdgCode=-1;
	}
      }

      Double_t invMass=d->InvMassDplus();
      Double_t rapid=d->YDplus();
      fYVsPt->Fill(ptCand,rapid,invMass);
      if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid,invMass);nOS++;}
      Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
      if(isFidAcc){
	fPtVsMass->Fill(invMass,ptCand);
	if(passTightCuts){
	  fPtVsMassTC->Fill(invMass,ptCand);
	  fPhiEtaCand->Fill(etaD,phiD);
	  if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
	}
      }


      Double_t  dlen=0,cosp=0,maxdca=0,sigvert=0,sumD02=0,ptmax=0,dlxy=0,cxy=0;
      if(fCutsDistr||fFillNtuple||fDoImpPar){
	dlen=d->DecayLength();
	cosp=d->CosPointingAngle();
	sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
	maxdca=-9999.;
	for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
	sigvert=d->GetSigmaVert();         
	ptmax=0;
	for(Int_t i=0;i<3;i++){
	  if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
	}
	dlxy=d->NormalizedDecayLengthXY();
	cxy=d->CosPointingAngleXY();
      }
      Double_t impparXY=d->ImpParXY()*10000.;
      Double_t arrayForSparse[6]={invMass,ptCand,impparXY,cosp,dlen,static_cast<Double_t>(tracklets)};
      Double_t arrayForSparseTrue[6]={invMass,ptCand,trueImpParXY,cosp,dlen,static_cast<Double_t>(tracklets)};

      //Ntuple
      Float_t tmp[31];
      if(fFillNtuple){
	tmp[0]=pdgCode;
	if(!isPrimary) tmp[0]+=5000.;
	tmp[1]=d->Px();
	tmp[2]=d->Py();
	tmp[3]=d->Pz();
	tmp[4]=d->Pt();
	tmp[5]=d->GetCharge();
	//	tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);	  
	tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
	tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
	tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
	tmp[9]=d->PtProng(0);
	tmp[10]=d->PtProng(1);
	tmp[11]=d->PtProng(2);
	tmp[12]=d->PProng(0);
	tmp[13]=d->PProng(1);
	tmp[14]=d->PProng(2);
	tmp[15]=cosp;
	tmp[16]=cxy;
	tmp[17]=dlen;
	tmp[18]=d->NormalizedDecayLength();
	tmp[19]=d->DecayLengthXY();
	tmp[20]=dlxy;
	tmp[21]=d->InvMassDplus();
	tmp[22]=sigvert;
	tmp[23]=d->Getd0Prong(0);
	tmp[24]=d->Getd0Prong(1);
	tmp[25]=d->Getd0Prong(2);
	tmp[26]=maxdca;
	tmp[27]=ntracks;
	tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
	tmp[29]=runNumber;
	tmp[30]=d->HasBadDaughters();
	fNtupleDplus->Fill(tmp);
	PostData(4,fNtupleDplus);
      }

      //Fill histos
      index=GetHistoIndex(iPtBin);
      if(isFidAcc){
	fHistNEvents->Fill(9);
	nSelectedloose++;
	fMassHist[index]->Fill(invMass);
	if(fCutsDistr){	  
	  fCosPHist[index]->Fill(cosp);
	  fDLenHist[index]->Fill(dlen);
	  fSumd02Hist[index]->Fill(sumD02);
	  fSigVertHist[index]->Fill(sigvert);
	  fPtMaxHist[index]->Fill(ptmax);
	  fPtKHist[index]->Fill(d->PtProng(1));
	  fPtpi1Hist[index]->Fill(d->PtProng(0));
	  fPtpi2Hist[index]->Fill(d->PtProng(2));
	  fDCAHist[index]->Fill(maxdca);
	  fDLxy[index]->Fill(dlxy);
	  fCosxy[index]->Fill(cxy);
	  fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
				  d->Getd0Prong(2)*d->Getd0Prong(1));
	}
	if(passTightCuts){ 
	  fHistNEvents->Fill(10);
	  nSelectedtight++;
	  fMassHistTC[index]->Fill(invMass);
	  if(fCutsDistr){  
	    fDLxyTC[index]->Fill(dlxy);
	    fCosxyTC[index]->Fill(cxy);
	  }
	  if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
	  else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
	  if(fDoImpPar){
	    fHistMassPtImpParTC[0]->Fill(arrayForSparse);
	  }
	}
      }
      
      if(fReadMC){
	if(isFidAcc){
	  if(labDp>=0) {
	    index=GetSignalHistoIndex(iPtBin);
	    if(passTightCuts&&fDoImpPar){
	      if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
	      else{
		fHistMassPtImpParTC[2]->Fill(arrayForSparse);
		fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
	      }
	    }
	  }else{
	    index=GetBackgroundHistoIndex(iPtBin);
	    if(passTightCuts&&fDoImpPar)fHistMassPtImpParTC[4]->Fill(arrayForSparse);
	  }
	  
	  fMassHist[index]->Fill(invMass);
	  if(fCutsDistr){
	    Float_t fact=1.;
	    Float_t factor[3]={1.,1.,1.};
	    if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
	    fCosPHist[index]->Fill(cosp,fact);
	    fDLenHist[index]->Fill(dlen,fact);
	    fDLxy[index]->Fill(dlxy);
	    fCosxy[index]->Fill(cxy);
	    
	    Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
	    fSumd02Hist[index]->Fill(sumd02s);
	    fSigVertHist[index]->Fill(sigvert,fact);
	    fPtMaxHist[index]->Fill(ptmax,fact);
	    fPtKHist[index]->Fill(d->PtProng(1),fact);
	    fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
	    fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
	    fDCAHist[index]->Fill(maxdca,fact);
	    fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
				    d->Getd0Prong(2)*d->Getd0Prong(1));
	  }
	  if(passTightCuts){
	    fMassHistTC[index]->Fill(invMass);	      
	    if(fCutsDistr){
	      fDLxyTC[index]->Fill(dlxy);
	      fCosxyTC[index]->Fill(cxy);
	    }	      
	    if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
	    else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
	  }
	}else{//outside fidAcc
	  if(labDp>=0){
	    fYVsPtSig->Fill(ptCand,rapid, invMass);
	    if(passTightCuts)fYVsPtSigTC->Fill(ptCand,rapid, invMass);
	  }
	}
      }//readmc
    
      if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
      
      if(unsetvtx) d->UnsetOwnPrimaryVtx();
    }
    fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
    fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
  }
  //start LS analysis
  if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
  
  PostData(1,fOutput); 
  PostData(2,fListCuts);
  PostData(3,fCounter);    
  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
  // Histos for Like Sign bckground

  TString hisname;
  Int_t indexLS=0;
  Int_t index=0;
  Int_t nbins=GetNBinsHistos();
  for(Int_t i=0;i<fNPtBins;i++){

    index=GetHistoIndex(i);
    indexLS=GetLSHistoIndex(i);

    hisname.Form("hLSPt%dLC",i);
    fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLS[indexLS]->Sumw2();
    hisname.Form("hLSPt%dTC",i);
    fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLSTC[indexLS]->Sumw2();

    hisname.Form("hCosPAllPt%dLS",i);
    fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
    fCosPHistLS[index]->Sumw2();
    hisname.Form("hDLenAllPt%dLS",i);
    fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
    fDLenHistLS[index]->Sumw2();
    hisname.Form("hSumd02AllPt%dLS",i);
    fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
    fSumd02HistLS[index]->Sumw2();
    hisname.Form("hSigVertAllPt%dLS",i);
    fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fSigVertHistLS[index]->Sumw2();
    hisname.Form("hPtMaxAllPt%dLS",i);
    fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
    fPtMaxHistLS[index]->Sumw2();
    hisname.Form("hDCAAllPt%dLS",i);
    fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fDCAHistLS[index]->Sumw2();    

    index=GetSignalHistoIndex(i);    
    indexLS++;
 
    hisname.Form("hLSPt%dLCnw",i);
    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLS[indexLS]->Sumw2();
    hisname.Form("hLSPt%dTCnw",i);
    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLSTC[indexLS]->Sumw2();

    hisname.Form("hCosPSigPt%dLS",i);
    fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
    fCosPHistLS[index]->Sumw2();
    hisname.Form("hDLenSigPt%dLS",i);
    fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
    fDLenHistLS[index]->Sumw2();
    hisname.Form("hSumd02SigPt%dLS",i);
    fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
    fSumd02HistLS[index]->Sumw2();
    hisname.Form("hSigVertSigPt%dLS",i);
    fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fSigVertHistLS[index]->Sumw2();
    hisname.Form("hPtMaxSigPt%dLS",i);
    fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
    fPtMaxHistLS[index]->Sumw2();
    hisname.Form("hDCASigPt%dLS",i);
    fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fDCAHistLS[index]->Sumw2();

    index=GetBackgroundHistoIndex(i); 
    indexLS++;

    hisname.Form("hLSPt%dLCntrip",i);
    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLS[indexLS]->Sumw2();
    hisname.Form("hLSPt%dTCntrip",i);
    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLSTC[indexLS]->Sumw2();

    hisname.Form("hCosPBkgPt%dLS",i);
    fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
    fCosPHistLS[index]->Sumw2();
    hisname.Form("hDLenBkgPt%dLS",i);
    fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
    fDLenHistLS[index]->Sumw2();
    hisname.Form("hSumd02BkgPt%dLS",i);
    fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
    fSumd02HistLS[index]->Sumw2();
    hisname.Form("hSigVertBkgPt%dLS",i);
    fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fSigVertHistLS[index]->Sumw2();
    hisname.Form("hPtMaxBkgPt%dLS",i);
    fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
    fPtMaxHistLS[index]->Sumw2();
    hisname.Form("hDCABkgPt%dLS",i);
    fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
    fDCAHistLS[index]->Sumw2();

    indexLS++;
    hisname.Form("hLSPt%dLCntripsinglecut",i);
    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLS[indexLS]->Sumw2();
    hisname.Form("hLSPt%dTCntripsinglecut",i);
    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLSTC[indexLS]->Sumw2();

    indexLS++;
    hisname.Form("hLSPt%dLCspc",i);
    fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLS[indexLS]->Sumw2();
    hisname.Form("hLSPt%dTCspc",i);
    fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
    fMassHistLSTC[indexLS]->Sumw2();
  }

  for(Int_t i=0; i<3*fNPtBins; i++){
    fOutput->Add(fCosPHistLS[i]);
    fOutput->Add(fDLenHistLS[i]);
    fOutput->Add(fSumd02HistLS[i]);
    fOutput->Add(fSigVertHistLS[i]);
    fOutput->Add(fPtMaxHistLS[i]);  
    fOutput->Add(fDCAHistLS[i]);  
  }
  for(Int_t i=0; i<5*fNPtBins; i++){
    fOutput->Add(fMassHistLS[i]);
    fOutput->Add(fMassHistLSTC[i]);
  }
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
  // Histos for impact paramter study

  Int_t nmassbins=GetNBinsHistos();
  Double_t maxmult;
  if(fSystem==1) maxmult=5000;
  else maxmult=200;
  Int_t nbins[6]={nmassbins,200,fNImpParBins,5,50,100};
  Double_t xmin[6]={fLowmasslimit,0.,fLowerImpPar,0.95,0.,-0.5};
  Double_t xmax[6]={fUpmasslimit,40.,fHigherImpPar,1.,1.,maxmult};

  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
					"Mass vs. pt vs.imppar - All",
					6,nbins,xmin,xmax);
  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
					"Mass vs. pt vs.imppar - promptD",
					6,nbins,xmin,xmax);
  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
					"Mass vs. pt vs.imppar - DfromB",
					6,nbins,xmin,xmax);
  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
					"Mass vs. pt vs.true imppar -DfromB",
					6,nbins,xmin,xmax);
  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
				        "Mass vs. pt vs.imppar - backgr.",
					6,nbins,xmin,xmax);
  for(Int_t i=0; i<5;i++){
    fOutput->Add(fHistMassPtImpParTC[i]);
  }
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");

  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    printf("ERROR: fOutput not available\n");
    return;
  }

  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
  if(fHistNEvents){
    printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
  }else{
    printf("ERROR: fHistNEvents not available\n");
    return;
  }

  return;
}
//_________________________________________________________________________________________________
Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
  // true impact parameter calculation

  Double_t vtxTrue[3];
  mcHeader->GetVertex(vtxTrue);
  Double_t origD[3];
  partDp->XvYvZv(origD);
  Short_t charge=partDp->Charge();
  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
  for(Int_t iDau=0; iDau<3; iDau++){
    pXdauTrue[iDau]=0.;
    pYdauTrue[iDau]=0.;
    pZdauTrue[iDau]=0.;
  }

  Int_t nDau=partDp->GetNDaughters();
  Int_t labelFirstDau = partDp->GetDaughter(0); 
  if(nDau==3){
    for(Int_t iDau=0; iDau<3; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return 99999.;
      } 
      pXdauTrue[iDau]=part->Px();
      pYdauTrue[iDau]=part->Py();
      pZdauTrue[iDau]=part->Pz();
    }
  }else if(nDau==2){
    Int_t theDau=0;
    for(Int_t iDau=0; iDau<2; iDau++){
      Int_t ind = labelFirstDau+iDau;
      AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
      if(!part){
	AliError("Daughter particle not found in MC array");
	return 99999.;
      } 
      Int_t pdgCode=TMath::Abs(part->GetPdgCode());
      if(pdgCode==211 || pdgCode==321){
	pXdauTrue[theDau]=part->Px();
	pYdauTrue[theDau]=part->Py();
	pZdauTrue[theDau]=part->Pz();
	++theDau;
      }else{
	Int_t nDauRes=part->GetNDaughters();
	if(nDauRes==2){
	  Int_t labelFirstDauRes = part->GetDaughter(0); 	
	  for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	    Int_t indDR = labelFirstDauRes+iDauRes;
	    AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
	    if(!partDR){
	      AliError("Daughter particle not found in MC array");
	      return 99999.;
	    } 
	    
	    Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 	    if(pdgCodeDR==211 || pdgCodeDR==321){
	      pXdauTrue[theDau]=partDR->Px();
	      pYdauTrue[theDau]=partDR->Py();
	      pZdauTrue[theDau]=partDR->Pz();
	      ++theDau;
	    }
	  }
	}
      }
    }
    if(theDau!=3){
      AliError("Wrong number of decay prongs");
      return 99999.;
    }
  }

  Double_t d0dummy[3]={0.,0.,0.};
  AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
  return aodDplusMC.ImpParXY();

}
//_________________________________________________________________________________________________
Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
  // Computes weights to adapt strangeness in MC to data

  for(Int_t iprong=0;iprong<3;iprong++){
    factor[iprong]=1;
    AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
    Int_t labd= trad->GetLabel();
    if(labd>=0){
      AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
      if(dau){
	Int_t labm = dau->GetMother();
	if(labm>=0){
	  AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
	  if(mot){
	    if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
	      if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
	      else factor[iprong]=1./.6;
	      //          fNentries->Fill(12);
	    }
	    if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
	      factor[iprong]=1./0.25;
	      //		  fNentries->Fill(13);
	    }//if 3122
	  }//if(mot)
	}//if labm>0
      }//if(dau)
    }//if labd>=0
  }//prong loop

  Float_t fact=1.;
  for(Int_t k=0;k<3;k++)fact=fact*factor[k];
  return fact;

}