]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
Coveriry
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
index 7fabe861a6fd0b291a072a32c38dc72d2c12a09c..01b0ef6a45223f5303741a4d1499ce1b2809d7fa 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //*************************************************************************
 // Class AliAnalysisTaskSEDplus
 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and 
@@ -28,7 +30,6 @@
 #include <TCanvas.h>
 #include <TList.h>
 #include <TString.h>
-#include <TH1F.h>
 #include <TDatabasePDG.h>
 
 #include "AliAnalysisManager.h"
@@ -37,8 +38,6 @@
 #include "AliAODEvent.h"
 #include "AliAODVertex.h"
 #include "AliAODTrack.h"
-#include "AliAODMCHeader.h"
-#include "AliAODMCParticle.h"
 #include "AliAODRecoDecayHF3Prong.h"
 #include "AliAnalysisVertexingHF.h"
 #include "AliAnalysisTaskSE.h"
@@ -71,51 +70,101 @@ AliAnalysisTaskSE(),
   fFillNtuple(kFALSE),
   fReadMC(kFALSE),
   fUseStrangeness(kFALSE),
-  fDoLS(kFALSE)
+  fUseBit(kTRUE),
+  fCutsDistr(kFALSE),
+  fDoImpPar(kFALSE),
+  fNImpParBins(400),
+  fLowerImpPar(-2000.),
+  fHigherImpPar(2000.),
+  fDoLS(0)
 {
    // Default constructor
 }
 
 //________________________________________________________________________
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
-AliAnalysisTaskSE(name),
-fOutput(0),
-fHistNEvents(0),
-fPtVsMass(0),
-fPtVsMassTC(0),
-fYVsPt(0),
-fYVsPtTC(0),
-fYVsPtSig(0),
-fYVsPtSigTC(0),
-fNtupleDplus(0),
-fUpmasslimit(1.965),
-fLowmasslimit(1.765),
-fNPtBins(0),
-fBinWidth(0.002),
-fListCuts(0),
-fRDCutsProduction(dpluscutsprod),
-fRDCutsAnalysis(dpluscutsana),
-fCounter(0),
-fFillNtuple(fillNtuple),
-fReadMC(kFALSE),
-fUseStrangeness(kFALSE),
-fDoLS(kFALSE)
+  AliAnalysisTaskSE(name),
+  fOutput(0),
+  fHistNEvents(0),
+  fPtVsMass(0),
+  fPtVsMassTC(0),
+  fYVsPt(0),
+  fYVsPtTC(0),
+  fYVsPtSig(0),
+  fYVsPtSigTC(0),
+  fNtupleDplus(0),
+  fUpmasslimit(1.965),
+  fLowmasslimit(1.765),
+  fNPtBins(0),
+  fBinWidth(0.002),
+  fListCuts(0),
+  fRDCutsProduction(dpluscutsprod),
+  fRDCutsAnalysis(dpluscutsana),
+  fCounter(0),
+  fFillNtuple(fillNtuple),
+  fReadMC(kFALSE),
+  fUseStrangeness(kFALSE),
+  fUseBit(kTRUE),
+  fCutsDistr(kFALSE),
+  fDoImpPar(kFALSE),
+  fNImpParBins(400),
+  fLowerImpPar(-2000.),
+  fHigherImpPar(2000.),
+  fDoLS(0)
 {
   // 
   // Standrd constructor
   //
-  //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
-   //SetPtBinLimit(5, ptlim);
-  SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
+  fNPtBins=fRDCutsAnalysis->GetNPtBins();
+  
+  for(Int_t i=0;i<3;i++){
+   if(fHistCentrality[i])fHistCentrality[i]=0;
+   if(fCorreld0Kd0pi[i])fCorreld0Kd0pi[i]=0;
+  }
+  
+  for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0;
+    
+  
+  for(Int_t i=0;i<3*fNPtBins;i++){
+    if(fMassHist[i])fMassHist[i]=0;
+    if(fCosPHist[i])fCosPHist[i]=0;
+    if(fDLenHist[i])fDLenHist[i]=0;
+    if(fSumd02Hist[i])fSumd02Hist[i]=0;
+    if(fSigVertHist[i])fSigVertHist[i]=0;
+    if(fPtMaxHist[i])fPtMaxHist[i]=0;
+    if(fPtKHist[i])fPtKHist[i]=0;
+    if(fPtpi1Hist[i])fPtpi1Hist[i]=0;
+    if(fPtpi2Hist[i])fPtpi2Hist[i]=0;
+    if(fDCAHist[i])fDCAHist[i]=0;
+    if(fMassHistTC[i])fMassHistTC[i]=0;
+    if(fMassHistTCPlus[i])fMassHistTCPlus[i]=0;
+    if(fMassHistTCMinus[i])fMassHistTCMinus[i]=0;
+    
+    if(fDLxy[i])fDLxy[i]=0;
+    if(fDLxyTC[i])fDLxyTC[i]=0;
+    if(fCosxy[i])fCosxy[i]=0;
+    if(fCosxyTC[i])fCosxyTC[i]=0;
+    if(fMassHistLS[i])fMassHistLS[i]=0;
+    if(fCosPHistLS[i])fCosPHistLS[i]=0;
+    if(fDLenHistLS[i])fDLenHistLS[i]=0;
+    if(fSumd02HistLS[i])fSumd02HistLS[i]=0;
+    if(fSigVertHistLS[i])fSigVertHistLS[i]=0;
+    if(fPtMaxHistLS[i])fPtMaxHistLS[i]=0;
+    if(fDCAHistLS[i])fDCAHistLS[i]=0;
+    if(fMassHistLSTC[i])fMassHistLSTC[i]=0;
+  }
+  
+  
+  
   // Default constructor
-   // Output slot #1 writes into a TList container
+  // 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
+  // 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
@@ -136,6 +185,10 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     delete fHistNEvents;
     fHistNEvents=0;
   }  
+  for(Int_t i=0;i<3;i++){
+    if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
+  }
+
   for(Int_t i=0;i<3*fNPtBins;i++){
     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
@@ -143,11 +196,18 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
+    if(fPtKHist[i]){ delete fPtKHist[i]; fPtKHist[i]=0;}
+    if(fPtpi1Hist[i]){ delete fPtpi1Hist[i]; fPtpi1Hist[i]=0;}
+    if(fPtpi2Hist[i]){ delete fPtpi2Hist[i]; fPtpi2Hist[i]=0;}
     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
     if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
     if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
     if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
 
+    if(fDLxy[i]){delete fDLxy[i]; fDLxy[i]=0;}
+    if(fDLxyTC[i]){delete fDLxyTC[i]; fDLxyTC[i]=0;}
+    if(fCosxy[i]){delete fCosxy[i]; fCosxy[i]=0;}
+    if(fCosxyTC[i]){delete fCosxyTC[i]; fCosxyTC[i]=0;}
     if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
     if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
     if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
@@ -157,6 +217,11 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
     if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
   }
+
+  for(Int_t i=0;i<3;i++){
+    if(fCorreld0Kd0pi[i]){ delete fCorreld0Kd0pi[i]; fCorreld0Kd0pi[i]=0;}
+  }
+
   if(fPtVsMass){
     delete fPtVsMass;
     fPtVsMass=0;
@@ -197,7 +262,9 @@ AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
     delete fRDCutsAnalysis;
     fRDCutsAnalysis = 0;
   }
-
+  for(Int_t i=0; i<5; i++){
+    delete fHistMassPtImpParTC[i];
+  }
   if(fCounter){
     delete fCounter;
     fCounter = 0;
@@ -224,34 +291,6 @@ void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
       SetBinWidth(bw);
     }
 }
-//________________________________________________________________________
-void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){
-  // define pt bins for analysis
-  if(n>kMaxPtBins){
-    printf("Max. number of Pt bins = %d\n",kMaxPtBins);
-    fNPtBins=kMaxPtBins;
-    fArrayBinLimits[0]=0.;
-    fArrayBinLimits[1]=2.;
-    fArrayBinLimits[2]=3.;
-    fArrayBinLimits[3]=5.;
-    for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
-  }else{
-    fNPtBins=n-1;
-    fArrayBinLimits[0]=lim[0];
-    for(Int_t i=1; i<fNPtBins+1; i++) 
-      if(lim[i]>fArrayBinLimits[i-1]){
-       fArrayBinLimits[i]=lim[i];
-      }
-      else {
-       fArrayBinLimits[i]=fArrayBinLimits[i-1];
-      }
-    for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
-  }
-  if(fDebug > 1){
-    printf("Number of Pt bins = %d\n",fNPtBins);
-    for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);    
-  }
-}
 //________________________________________________________________
 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
   Float_t width=w;
@@ -268,12 +307,6 @@ void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
   fBinWidth=width;
 }
 //_________________________________________________________________
-Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
-  // get pt bin limit
-  if(ibin>fNPtBins)return -1;
-  return fArrayBinLimits[ibin];
-} 
-//_________________________________________________________________
 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
 }
@@ -283,136 +316,103 @@ void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesA
   //
   // Fill the Like Sign histograms
   //
-
-  //count pos/neg tracks
+  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;
- //counter for particles passing single particle cuts
-  Int_t nspcplus=0;
-  Int_t nspcminus=0;
-
-  for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
-    AliAODTrack *track = aod->GetTrack(it);
-    if(track->Charge()>0){
-      nPosTrks++;
-      if(track->Pt()>=0.4){
-       nspcplus++;
-      }
-    }
-    if(track->Charge()<0)
-      {
-       nNegTrks++;
-       if(track->Pt()>=0.4){
-         nspcminus++;
-       }
-      }
-  }
-
   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
-
   Int_t nDplusLS=0;
+  Int_t nDminusLS=0;
   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
-  Int_t index; 
-
+  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;
     }
-    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate))nDplusLS++;
-    if(unsetvtx) d->UnsetOwnPrimaryVtx();
-  }
-
- 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;
-
- // loop over like sign candidates
-  for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
-    AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
-    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++;
     }
-    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)){
-
-      //set tight cuts values
-      Int_t iPtBin=-1;
-      Double_t ptCand = d->Pt();
-      for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
-       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
-      }
-      
-      if(iPtBin<0){
-       return;
-      }
-
-      Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
-
-      Int_t sign= d->GetCharge();
-      Float_t wei=1;
-      Float_t wei4=1;
-      if(sign>0&&nPosTrks>2&&nspcplus>2) {  //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
-       
-       wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
-       wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
-      }
-      
-      if(sign<0&&nNegTrks>2&&nspcminus>2){     
-       wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
-       wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
-
-      }
-
+    //    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
+    fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+    Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
+    if(passTightCuts>0){
       Float_t invMass = d->InvMassDplus();
-      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);
-     }
-     
       index=GetLSHistoIndex(iPtBin);
-      fMassHistLS[index]->Fill(invMass,wei);
       fMassHistLS[index+1]->Fill(invMass);
-      fMassHistLS[index+2]->Fill(invMass,wei2);
-      fMassHistLS[index+3]->Fill(invMass,wei3);
-      fMassHistLS[index+4]->Fill(invMass,wei4);
-      
-      Int_t indexcut=GetHistoIndex(iPtBin);
-      fCosPHistLS[indexcut]->Fill(cosp);
-      fDLenHistLS[indexcut]->Fill(dlen);
-      fSumd02HistLS[indexcut]->Fill(sumD02);
-      fSigVertHistLS[indexcut]->Fill(sigvert);
-      fPtMaxHistLS[indexcut]->Fill(ptmax);
-      fDCAHistLS[indexcut]->Fill(dca);
-      
-      if(passTightCuts){
-       fMassHistLSTC[index]->Fill(invMass,wei);
-       fMassHistLSTC[index+1]->Fill(invMass);
-       fMassHistLSTC[index+2]->Fill(invMass,wei2);
-       fMassHistLSTC[index+3]->Fill(invMass,wei3);
-       fMassHistLSTC[index+4]->Fill(invMass,wei4);
+      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;
   
-  //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
-  //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
-
-  //  printf("LS analysis...done\n");
-
+  if(fDebug>1) printf("LS analysis terminated\n");  
 }
 
-
 //__________________________________________
 void AliAnalysisTaskSEDplus::Init(){
   //
@@ -422,10 +422,10 @@ void AliAnalysisTaskSEDplus::Init(){
   
   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
   fListCuts=new TList();
-  AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi();
-  production=fRDCutsProduction;
-  AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi();
-  analysis=fRDCutsAnalysis;
+  AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
+  production->SetName("ProductionCuts");
+  AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
+  analysis->SetName("AnalysisCuts");
   
   fListCuts->Add(production);
   fListCuts->Add(analysis);
@@ -448,12 +448,17 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 
   TString hisname;
   Int_t index=0;
-  Int_t indexLS=0;
   Int_t nbins=GetNBinsHistos();
+  fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
+  fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
+  fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
+  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);
-    indexLS=GetLSHistoIndex(i);
 
     hisname.Form("hMassPt%d",i);
     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
@@ -473,13 +478,32 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     hisname.Form("hPtMaxAllPt%d",i);
     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
     fPtMaxHist[index]->Sumw2();
-
+    hisname.Form("hPtKPt%d",i);
+    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtKHist[index]->Sumw2();
+    hisname.Form("hPtpi1Pt%d",i);
+    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi1Hist[index]->Sumw2();
+    hisname.Form("hPtpi2Pt%d",i);
+    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi2Hist[index]->Sumw2();
     hisname.Form("hDCAAllPt%d",i);
     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
     fDCAHist[index]->Sumw2();
 
-
-
+    hisname.Form("hDLxyPt%d",i);
+    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxy[index]->Sumw2();
+    hisname.Form("hCosxyPt%d",i);
+    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxy[index]->Sumw2();
+    hisname.Form("hDLxyPt%dTC",i);
+    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxyTC[index]->Sumw2();
+    hisname.Form("hCosxyPt%dTC",i);
+    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxyTC[index]->Sumw2();
+   
     hisname.Form("hMassPt%dTC",i);
     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHistTC[index]->Sumw2();
@@ -491,40 +515,8 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     fMassHistTCMinus[index]->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();
-    
-    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();
-
-
-    
+       
     index=GetSignalHistoIndex(i);    
-    indexLS++;
     hisname.Form("hSigPt%d",i);
     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHist[index]->Sumw2();
@@ -542,13 +534,33 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     fSigVertHist[index]->Sumw2();
     hisname.Form("hPtMaxSigPt%d",i);
     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
-    fPtMaxHist[index]->Sumw2();    
+    fPtMaxHist[index]->Sumw2();  
+    hisname.Form("hPtKSigPt%d",i);  
+    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtKHist[index]->Sumw2();
+    hisname.Form("hPtpi1SigPt%d",i);
+    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi1Hist[index]->Sumw2();
+    hisname.Form("hPtpi2SigPt%d",i);
+    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi2Hist[index]->Sumw2();
 
     hisname.Form("hDCASigPt%d",i);
     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
     fDCAHist[index]->Sumw2();    
 
-
+    hisname.Form("hDLxySigPt%d",i);
+    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxy[index]->Sumw2();
+    hisname.Form("hCosxySigPt%d",i);
+    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxy[index]->Sumw2();
+    hisname.Form("hDLxySigPt%dTC",i);
+    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxyTC[index]->Sumw2();
+    hisname.Form("hCosxySigPt%dTC",i);
+    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxyTC[index]->Sumw2();
     hisname.Form("hSigPt%dTC",i);
     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHistTC[index]->Sumw2();
@@ -559,39 +571,8 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHistTCMinus[index]->Sumw2();
 
-    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("hBkgPt%d",i);
     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHist[index]->Sumw2();
@@ -610,11 +591,32 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     hisname.Form("hPtMaxBkgPt%d",i);
     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
     fPtMaxHist[index]->Sumw2();
-
+    hisname.Form("hPtKBkgPt%d",i);  
+    fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtKHist[index]->Sumw2();
+    hisname.Form("hPtpi1BkgPt%d",i);
+    fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi1Hist[index]->Sumw2();
+    hisname.Form("hPtpi2BkgPt%d",i);
+    fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtpi2Hist[index]->Sumw2();
     hisname.Form("hDCABkgPt%d",i);
     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
     fDCAHist[index]->Sumw2();
 
+    hisname.Form("hDLxyBkgPt%d",i);
+    fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxy[index]->Sumw2();
+    hisname.Form("hCosxyBkgPt%d",i);
+    fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxy[index]->Sumw2();
+    hisname.Form("hDLxyBkgPt%dTC",i);
+    fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
+    fDLxyTC[index]->Sumw2();
+    hisname.Form("hCosxyBkgPt%dTC",i);
+    fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
+    fCosxyTC[index]->Sumw2();
+  
 
     hisname.Form("hBkgPt%dTC",i);
     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
@@ -625,85 +627,59 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
     hisname.Form("hBkgPt%dTCMinus",i);
     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
     fMassHistTCMinus[index]->Sumw2();
-
-    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(fMassHist[i]);
-     fOutput->Add(fCosPHist[i]);
-    fOutput->Add(fDLenHist[i]);
-    fOutput->Add(fSumd02Hist[i]);
-    fOutput->Add(fSigVertHist[i]);
-    fOutput->Add(fPtMaxHist[i]);
-    fOutput->Add(fDCAHist[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]);
+
   }
-  for(Int_t i=0; i<3*fNPtBins&&fDoLS; 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&&fDoLS; i++){
-    fOutput->Add(fMassHistLS[i]);
-    fOutput->Add(fMassHistLSTC[i]);
+  
+  if(fCutsDistr){  
+    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);
+    for(Int_t i=0; i<3; i++){
+      fCorreld0Kd0pi[i]->Sumw2();
+      fOutput->Add(fCorreld0Kd0pi[i]);
+    }
   }
   
-
-  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
+  fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
+  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
+  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
+  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
+  fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
+  fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
+  fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
+  fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
+  fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
+  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
+  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);  
   fHistNEvents->Sumw2();
   fHistNEvents->SetMinimum(0);
   fOutput->Add(fHistNEvents);
 
-  fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.);
-  fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.);  
+  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 TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
@@ -718,10 +694,17 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 
 
   //Counter for Normalization
-  fCounter = new AliNormalizationCounter("NormalizationCounter");//new line
+  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(2); // 2 is the slot number of the ntuple
+    OpenFile(4); // 4 is the slot number of the ntuple
    
     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");  
     
@@ -756,16 +739,16 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
     }
-  } else {
+  } else if(aod) {
     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
   }
 
-  if(!array3Prong) {
+  if(!aod || !array3Prong) {
     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
     return;
   }
-  if(!arrayLikeSign) {
+  if(!arrayLikeSign && fDoLS) {
     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
     return;
   }
@@ -774,25 +757,44 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   // 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;
-  fCounter->StoreEvent(aod,fReadMC);
+  fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
   fHistNEvents->Fill(0); // count event
-  // Post the data already here
+
+  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
+  Bool_t isEvSelP=kTRUE;
+  isEvSelP=fRDCutsProduction->IsEventSelected(aod); // to have proper PID object settings
+
+  Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
+  fHistCentrality[0]->Fill(centrality);
+  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
+  TString trigclass=aod->GetFiredTriggerClasses();
+  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
+  if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3); 
+  if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
+
+  // Post the data already here  
   PostData(1,fOutput);
+  if(!isEvSel)return;
+
+  fHistCentrality[1]->Fill(centrality);
+  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;
+      return;
     }
     
   // load MC header
@@ -804,40 +806,53 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   }
   
   Int_t n3Prong = array3Prong->GetEntriesFast();
-  if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
+  //  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){
+    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::kCandidate,aod))nOS++;
+      }
+    }
+  }else{
   // 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(4);
+    if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
+      fHistNEvents->Fill(8);
+      continue;
+    }
     Bool_t unsetvtx=kFALSE;
     if(!d->GetOwnPrimaryVtx()){
       d->SetOwnPrimaryVtx(vtx1);
       unsetvtx=kTRUE;
     }
 
-    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) {
+    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
 
-  
-
-      Int_t iPtBin = -1;
       Double_t ptCand = d->Pt();
-
-      for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
-       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
+      Int_t iPtBin = fRDCutsProduction->PtBin(ptCand);
+      
+      Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+      Bool_t recVtx=kFALSE;
+      AliAODVertex *origownvtx=0x0;
+      if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){
+       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());  
+       if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
+       else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
       }
       
-      Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
-     
-
       Int_t labDp=-1;
+      Bool_t isPrimary=kTRUE;
       Float_t deltaPx=0.;
       Float_t deltaPy=0.;
       Float_t deltaPz=0.;
@@ -846,10 +861,12 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       Float_t yDecay=0.;
       Float_t zDecay=0.;
       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(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
          AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
          deltaPx=partDp->Px()-d->Px();
          deltaPy=partDp->Py()-d->Py();
@@ -859,21 +876,37 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
          yDecay=dg0->Yv();       
          zDecay=dg0->Zv();
          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);
-      if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid);
+      if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
       if(isFidAcc){
        fPtVsMass->Fill(invMass,ptCand);
        if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
       }
       Float_t tmp[24];
-      if(fFillNtuple){           
+      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);
+      }
+      Double_t impparXY=d->ImpParXY()*10000.;
+      Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
+      Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
+      if(fFillNtuple){
        tmp[0]=pdgCode;
        tmp[1]=deltaPx;
        tmp[2]=deltaPy;
@@ -886,52 +919,61 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        tmp[9]=d->PtProng(1);
        tmp[10]=d->PtProng(2);
        tmp[11]=d->Pt();
-       tmp[12]=d->CosPointingAngle();
-       tmp[13]=d->DecayLength();
+       tmp[12]=cosp;
+       tmp[13]=dlen;
        tmp[14]=d->Xv();
        tmp[15]=d->Yv();
        tmp[16]=d->Zv();
        tmp[17]=d->InvMassDplus();
-       tmp[18]=d->GetSigmaVert();
+       tmp[18]=sigvert;
        tmp[19]=d->Getd0Prong(0);
        tmp[20]=d->Getd0Prong(1);
        tmp[21]=d->Getd0Prong(2);
-       tmp[22]=d->GetDCA();
+       tmp[22]=dca;
        tmp[23]=d->Prodd0d0(); 
        fNtupleDplus->Fill(tmp);
-       PostData(3,fNtupleDplus);
-      }
-      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);
+       PostData(4,fNtupleDplus);
       }
       if(iPtBin>=0){
-      
+       Float_t dlxy=d->NormalizedDecayLengthXY();
+       Float_t cxy=d->CosPointingAngleXY();
        index=GetHistoIndex(iPtBin);
        if(isFidAcc){
+         fHistNEvents->Fill(5);
          nSelectedloose++;
          fMassHist[index]->Fill(invMass);
-         fCosPHist[index]->Fill(cosp);
-         fDLenHist[index]->Fill(dlen);
-         fSumd02Hist[index]->Fill(sumD02);
-         fSigVertHist[index]->Fill(sigvert);
-         fPtMaxHist[index]->Fill(ptmax);
-         fDCAHist[index]->Fill(dca);
-         
-         if(passTightCuts){
+         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(dca);
+           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(6);
            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(fCutsDistr){
          if(labDp>=0) {
            index=GetSignalHistoIndex(iPtBin);
            if(isFidAcc){
@@ -964,17 +1006,38 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
              }
              Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
              fMassHist[index]->Fill(invMass);
-             fCosPHist[index]->Fill(cosp,fact);
-             fDLenHist[index]->Fill(dlen,fact);
-             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);
-             fDCAHist[index]->Fill(dca,fact);
+             if(fCutsDistr){
+               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(dca,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);
+               if(fDoImpPar){
+                 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
+                 else{
+                   fHistMassPtImpParTC[2]->Fill(arrayForSparse);
+                   fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
+                 }
+               }
              }
            }       
            fYVsPtSig->Fill(ptCand,rapid);
@@ -1009,39 +1072,217 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
                  }//if labd>=0
                }//prong loop
              }
+           
              Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
              fMassHist[index]->Fill(invMass);
-             fCosPHist[index]->Fill(cosp,fact);
-             fDLenHist[index]->Fill(dlen,fact);
-             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);
-             fDCAHist[index]->Fill(dca,fact);
+             if(fCutsDistr){
+               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(dca,fact);
+               fCorreld0Kd0pi[2]->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);
+               if(fDoImpPar){
+                 fHistMassPtImpParTC[4]->Fill(arrayForSparse);
+               }
              }
-           }   
+           }
          }
+         
        }
-      }  
+      }
+    
+      if(recVtx)fRDCutsProduction->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(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();
+  Int_t nbins[3]={nmassbins,200,fNImpParBins};
+  Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
+  Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
+
+  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
+                                       "Mass vs. pt vs.imppar - All",
+                                       3,nbins,xmin,xmax);
+  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
+                                       "Mass vs. pt vs.imppar - promptD",
+                                       3,nbins,xmin,xmax);
+  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
+                                       "Mass vs. pt vs.imppar - DfromB",
+                                       3,nbins,xmin,xmax);
+  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
+                                       "Mass vs. pt vs.true imppar -DfromB",
+                                       3,nbins,xmin,xmax);
+  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
+                                       "Mass vs. pt vs.imppar - backgr.",
+                                       3,nbins,xmin,xmax);
+
+  for(Int_t i=0; i<5;i++){
+    fOutput->Add(fHistMassPtImpParTC[i]);
+  }
+}
 
 //________________________________________________________________________
 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
@@ -1056,172 +1297,17 @@ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
     return;
   }
   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
-  fYVsPt = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPt"));
-  fYVsPtTC = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtTC"));
-  fYVsPtSig = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtSig"));
-  fYVsPtSigTC = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtSigTC"));
-  fPtVsMass = dynamic_cast<TH2F*>(fOutput->FindObject("hPtVsMass"));
-  fPtVsMassTC = dynamic_cast<TH2F*>(fOutput->FindObject("hPtVsMassTC"));
 
   TString hisname;
   Int_t index=0;
 
-  Int_t indexLS=0;
   for(Int_t i=0;i<fNPtBins;i++){
     index=GetHistoIndex(i);
-    if(fDoLS)indexLS=GetLSHistoIndex(i);
-    hisname.Form("hMassPt%d",i);
-    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hCosPAllPt%d",i);
-    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDLenAllPt%d",i);
-    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSumd02AllPt%d",i);
-    fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSigVertAllPt%d",i);
-    fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hPtMaxAllPt%d",i);
-    fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDCAAllPt%d",i);
-    fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hMassPt%dTC",i);
-    fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hMassPt%dTCPlus",i);
-    fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hMassPt%dTCMinus",i);
-    fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    if(fDoLS){
-      hisname.Form("hLSPt%dLC",i);
-      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hCosPAllPt%dLS",i);
-      fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDLenAllPt%dLS",i);
-      fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSumd02AllPt%dLS",i);
-      fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSigVertAllPt%dLS",i);
-      fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hPtMaxAllPt%dLS",i);
-      fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDCAAllPt%dLS",i);
-      fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      hisname.Form("hLSPt%dTC",i);
-      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-    } 
-    
-    index=GetSignalHistoIndex(i);    
-    if(fDoLS)indexLS++;
-    hisname.Form("hSigPt%d",i);
-    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hCosPSigPt%d",i);
-    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDLenSigPt%d",i);
-    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSumd02SigPt%d",i);
-    fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSigVertSigPt%d",i);
-    fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hPtMaxSigPt%d",i);
-    fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDCASigPt%d",i);
-    fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
     
-    hisname.Form("hSigPt%dTC",i);
+    hisname.Form("hMassPt%dTC",i);
     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSigPt%dTCPlus",i);
-    fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSigPt%dTCMinus",i);
-    fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    if(fDoLS){
-      hisname.Form("hLSPt%dLCnw",i);
-      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hCosPSigPt%dLS",i);
-      fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDLenSigPt%dLS",i);
-      fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSumd02SigPt%dLS",i);
-      fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSigVertSigPt%dLS",i);
-      fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hPtMaxSigPt%dLS",i);
-      fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDCASigPt%dLS",i);
-      fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-
-      hisname.Form("hLSPt%dTCnw",i);
-      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    }
+  } 
     
-    index=GetBackgroundHistoIndex(i); 
-    if(fDoLS)indexLS++;
-    hisname.Form("hBkgPt%d",i);
-    fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hCosPBkgPt%d",i);
-    fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDLenBkgPt%d",i);
-    fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSumd02BkgPt%d",i);
-    fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hSigVertBkgPt%d",i);
-    fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hPtMaxBkgPt%d",i);
-    fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hDCABkgPt%d",i);
-    fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hBkgPt%dTC",i);
-    fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hBkgPt%dTCPlus",i);
-    fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    hisname.Form("hBkgPt%dTCMinus",i);
-    fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    if(fDoLS){
-      hisname.Form("hLSPt%dLCntrip",i);
-      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-
-      hisname.Form("hCosPBkgPt%dLS",i);
-      fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDLenBkgPt%dLS",i);
-      fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSumd02BkgPt%dLS",i);
-      fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hSigVertBkgPt%dLS",i);
-      fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hPtMaxBkgPt%dLS",i);
-      fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      hisname.Form("hDCABkgPt%dLS",i);
-      fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      hisname.Form("hLSPt%dTCntrip",i);
-      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      
-      indexLS++;
-      hisname.Form("hLSPt%dLCntripsinglecut",i);
-      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      hisname.Form("hLSPt%dTCntripsinglecut",i);
-      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      
-      indexLS++;
-      hisname.Form("hLSPt%dLCspc",i);
-      fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-      
-      hisname.Form("hLSPt%dTCspc",i);
-      fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
-    }
-  }
-  fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(3)); 
-
-  if(fFillNtuple){
-    fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(4));
-  }
-
   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
   c1->cd();
   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
@@ -1231,3 +1317,115 @@ void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
  
   return;
 }
+//_________________________________________________________________________________________________
+Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {            
+  //
+  // checking whether the mother of the particles come from a charm or a bottom quark
+  //
+       
+  Int_t pdgGranma = 0;
+  Int_t mother = 0;
+  mother = mcPartCandidate->GetMother();
+  Int_t istep = 0;
+  Int_t abspdgGranma =0;
+  Bool_t isFromB=kFALSE;
+  Bool_t isQuarkFound=kFALSE;
+  while (mother >0 ){
+    istep++;
+    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
+    if (mcGranma){
+      pdgGranma = mcGranma->GetPdgCode();
+      abspdgGranma = TMath::Abs(pdgGranma);
+      if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+       isFromB=kTRUE;
+      }
+      if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+      mother = mcGranma->GetMother();
+    }else{
+      AliError("Failed casting the mother particle!");
+      break;
+    }
+  }
+  
+  if(isFromB) return 5;
+  else return 4;
+}
+//_________________________________________________________________________________________________
+Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, 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();
+
+}