]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx
New class for normalization studies (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
index 8acfa3a57af1a9679909d72a8e273cd510b56dbf..7fabe861a6fd0b291a072a32c38dc72d2c12a09c 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/////////////////////////////////////////////////////////////
-//
-// AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
-// decay candidates with the MC truth.
+//*************************************************************************
+// 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 <TH1F.h>
 #include <TDatabasePDG.h>
 
 #include "AliAnalysisManager.h"
+#include "AliRDHFCutsDplustoKpipi.h"
 #include "AliAODHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODVertex.h"
@@ -40,6 +43,7 @@
 #include "AliAnalysisVertexingHF.h"
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisTaskSEDplus.h"
+#include "AliNormalizationCounter.h"
 
 ClassImp(AliAnalysisTaskSEDplus)
 
@@ -47,84 +51,386 @@ ClassImp(AliAnalysisTaskSEDplus)
 //________________________________________________________________________
 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
 AliAnalysisTaskSE(),
-fOutput(0), 
-fHistNEvents(0),
-fNtupleDplus(0),
-fHistLS(0),
-fHistLStrip(0),
-fHistLStripcut(0),
-fHistOS(0),
-fHistOSbkg(0),
-fHistLSnoweight(0),
-fHistLSsinglecut(0),
-fFillNtuple(kFALSE),
-fReadMC(kFALSE),
-fDoLS(kFALSE),
-fVHF(0)
+  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(0),
+  fRDCutsAnalysis(0),
+  fCounter(0),
+  fFillNtuple(kFALSE),
+  fReadMC(kFALSE),
+  fUseStrangeness(kFALSE),
+  fDoLS(kFALSE)
 {
-  // Default constructor
+   // Default constructor
 }
 
 //________________________________________________________________________
-AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
+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),
-fHistLS(0),
-fHistLStrip(0),
-fHistLStripcut(0),
-fHistOS(0),
-fHistOSbkg(0),
-fHistLSnoweight(0),
-fHistLSsinglecut(0),
+fUpmasslimit(1.965),
+fLowmasslimit(1.765),
+fNPtBins(0),
+fBinWidth(0.002),
+fListCuts(0),
+fRDCutsProduction(dpluscutsprod),
+fRDCutsAnalysis(dpluscutsana),
+fCounter(0),
 fFillNtuple(fillNtuple),
 fReadMC(kFALSE),
-fDoLS(kFALSE),
-fVHF(0)
+fUseStrangeness(kFALSE),
+fDoLS(kFALSE)
 {
+  // 
+  // Standrd constructor
+  //
+  //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
+   //SetPtBinLimit(5, ptlim);
+  SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
   // 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
+  DefineOutput(3,AliNormalizationCounter::Class());
 
   if(fFillNtuple){
-    // Output slot #2 writes into a TNtuple container
-    DefineOutput(2,TNtuple::Class());  //My private output
+    // Output slot #4 writes into a TNtuple container
+    DefineOutput(4,TNtuple::Class());  //My private output
   }
 }
 
 //________________________________________________________________________
 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
 {
+  //
   // Destructor
+  //
   if (fOutput) {
     delete fOutput;
     fOutput = 0;
   }
+  if(fHistNEvents){
+    delete fHistNEvents;
+    fHistNEvents=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;}
+    if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
+    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(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(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;}
+    if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
+    if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
+    if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
+    if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
+    if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
+  }
+  if(fPtVsMass){
+    delete fPtVsMass;
+    fPtVsMass=0;
+  }
+  if(fPtVsMassTC){
+    delete fPtVsMassTC;
+    fPtVsMassTC=0;
+  }
+  if(fYVsPt){
+    delete fYVsPt;
+    fYVsPt=0;
+  }
+  if(fYVsPtTC){
+    delete fYVsPtTC;
+    fYVsPtTC=0;
+  }
+  if(fYVsPtSig){
+    delete fYVsPtSig;
+    fYVsPtSig=0;
+  }
+  if(fYVsPtSigTC){
+    delete fYVsPtSigTC;
+    fYVsPtSigTC=0;
+  }
   if(fNtupleDplus){
     delete fNtupleDplus;
     fNtupleDplus=0;
   }
-  if (fVHF) {
-    delete fVHF;
-    fVHF = 0;
+  if (fListCuts) {
+    delete fListCuts;
+    fListCuts = 0;
+  }
+  if(fRDCutsProduction){
+    delete fRDCutsProduction;
+    fRDCutsProduction = 0;
+  }
+  if(fRDCutsAnalysis){
+    delete fRDCutsAnalysis;
+    fRDCutsAnalysis = 0;
+  }
+
+  if(fCounter){
+    delete fCounter;
+    fCounter = 0;
   }
 
-}  
 
+}  
+//_________________________________________________________________
+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::Init()
-{
-  // Initialization
+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;
+  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;
+}
+//_________________________________________________________________
+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);
+}
+//_________________________________________________________________
+void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
+  //
+  //
+  // Fill the Like Sign histograms
+  //
 
-  if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
+  //count pos/neg tracks
+  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();
 
-  gROOT->LoadMacro("ConfigVertexingHF.C");
+  Int_t nDplusLS=0;
+  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
+  Int_t index; 
+
+  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;
+    }
+    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate))nDplusLS++;
+    if(unsetvtx) d->UnsetOwnPrimaryVtx();
+  }
 
-  fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
-  fVHF->PrintStatus();
+ 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;
+    }
+    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.);
+
+      }
+
+      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(unsetvtx) d->UnsetOwnPrimaryVtx();
+  }
+  
+  //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");
+
+}
+
+
+//__________________________________________
+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();
+  AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi();
+  production=fRDCutsProduction;
+  AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi();
+  analysis=fRDCutsAnalysis;
+  
+  fListCuts->Add(production);
+  fListCuts->Add(analysis);
+  PostData(2,fListCuts);
+  
   return;
 }
 
@@ -140,78 +446,287 @@ void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
   fOutput->SetOwner();
   fOutput->SetName("OutputHistos");
 
-  Int_t nPtBins=4;
   TString hisname;
-  for(Int_t i=0;i<nPtBins;i++){
+  Int_t index=0;
+  Int_t indexLS=0;
+  Int_t nbins=GetNBinsHistos();
+  for(Int_t i=0;i<fNPtBins;i++){
+
+    index=GetHistoIndex(i);
+    indexLS=GetLSHistoIndex(i);
+
     hisname.Form("hMassPt%d",i);
-    TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hm->Sumw2();
-    hisname.Form("hSigPt%d",i);
-    TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hs->Sumw2();
-    hisname.Form("hBkgPt%d",i);
-    TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hb->Sumw2();
+    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(),nbins,0.5,1.);
+    fCosPHist[index]->Sumw2();
+    hisname.Form("hDLenAllPt%d",i);
+    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
+    fDLenHist[index]->Sumw2();
+    hisname.Form("hSumd02AllPt%d",i);
+    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
+    fSumd02Hist[index]->Sumw2();
+    hisname.Form("hSigVertAllPt%d",i);
+    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fSigVertHist[index]->Sumw2();
+    hisname.Form("hPtMaxAllPt%d",i);
+    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtMaxHist[index]->Sumw2();
+
+    hisname.Form("hDCAAllPt%d",i);
+    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fDCAHist[index]->Sumw2();
+
+
 
     hisname.Form("hMassPt%dTC",i);
-    TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hmtc->Sumw2();
+    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();
+
+
+
+    
+    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();
+    hisname.Form("hCosPSigPt%d",i);
+    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
+    fCosPHist[index]->Sumw2();
+    hisname.Form("hDLenSigPt%d",i);
+    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
+    fDLenHist[index]->Sumw2();
+    hisname.Form("hSumd02SigPt%d",i);
+    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
+    fSumd02Hist[index]->Sumw2();
+    hisname.Form("hSigVertSigPt%d",i);
+    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fSigVertHist[index]->Sumw2();
+    hisname.Form("hPtMaxSigPt%d",i);
+    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtMaxHist[index]->Sumw2();    
+
+    hisname.Form("hDCASigPt%d",i);
+    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fDCAHist[index]->Sumw2();    
+
+
     hisname.Form("hSigPt%dTC",i);
-    TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hstc->Sumw2();
+    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();
+
+    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();
+    hisname.Form("hCosPBkgPt%d",i);
+    fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
+    fCosPHist[index]->Sumw2();
+    hisname.Form("hDLenBkgPt%d",i);
+    fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
+    fDLenHist[index]->Sumw2();
+    hisname.Form("hSumd02BkgPt%d",i);
+    fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
+    fSumd02Hist[index]->Sumw2();
+    hisname.Form("hSigVertBkgPt%d",i);
+    fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fSigVertHist[index]->Sumw2();
+    hisname.Form("hPtMaxBkgPt%d",i);
+    fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
+    fPtMaxHist[index]->Sumw2();
+
+    hisname.Form("hDCABkgPt%d",i);
+    fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
+    fDCAHist[index]->Sumw2();
+
+
     hisname.Form("hBkgPt%dTC",i);
-    TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
-    hbtc->Sumw2();
+    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();
+
+    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();
 
-    fOutput->Add(hm);
-    fOutput->Add(hs);
-    fOutput->Add(hb);
-    fOutput->Add(hmtc);
-    fOutput->Add(hstc);
-    fOutput->Add(hbtc);
+    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]);
+    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]);
+  }
+  
+
   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
   fHistNEvents->Sumw2();
   fHistNEvents->SetMinimum(0);
   fOutput->Add(fHistNEvents);
-  if(fDoLS){
-    Double_t massDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
-  
-    fHistLS = new TH1F("fHistLS","fHistLS",100,massDplus-0.2,massDplus+0.2);
-    fHistLS->Sumw2();
-    fOutput->Add(fHistLS);
 
-    fHistLStrip = new TH1F("fHistLStrip","fHistLStrip",100,massDplus-0.2,massDplus+0.2);
-    fHistLStrip->Sumw2();
-    fOutput->Add(fHistLStrip);
+  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.);  
+  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.);
+  fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 
-    fHistLStripcut = new TH1F("fHistLStripcut","fHistLStripcut",100,massDplus-0.2,massDplus+0.2);
-    fHistLStripcut->Sumw2();
-    fOutput->Add(fHistLStripcut);
+  fOutput->Add(fPtVsMass);
+  fOutput->Add(fPtVsMassTC);
+  fOutput->Add(fYVsPt);
+  fOutput->Add(fYVsPtTC);
+  fOutput->Add(fYVsPtSig);
+  fOutput->Add(fYVsPtSigTC);
 
-    fHistOS = new TH1F("fHistOS","fHistOS",100,massDplus-0.2,massDplus+0.2);
-    fHistOS->Sumw2();
-    fOutput->Add(fHistOS);
 
-    fHistOSbkg = new TH1F("fHistOSbkg","fHistOSbkg",100,massDplus-0.2,massDplus+0.2);
-    fHistOSbkg->Sumw2();
-    fOutput->Add(fHistOSbkg);
-
-    fHistLSnoweight = new TH1F("fHistLSnoweight","fHistLSnoweight",100,massDplus-0.2,massDplus+0.2);
-    fHistLSnoweight->Sumw2();
-    fOutput->Add(fHistLSnoweight);
-
-    fHistLSsinglecut = new TH1F("fHistLSsinglecut","fHistLSsinglecut",100,massDplus-0.2,massDplus+0.2);
-    fHistLSsinglecut->Sumw2();
-    fOutput->Add(fHistLSsinglecut);
-  }
+  //Counter for Normalization
+  fCounter = new AliNormalizationCounter("NormalizationCounter");//new line
 
   if(fFillNtuple){
     OpenFile(2); // 2 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");
+   
+    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");  
+    
   }
-
+  
   return;
 }
 
@@ -221,13 +736,9 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   // Execute analysis for current event:
   // heavy flavor candidates association to MC truth
 
-  
-  
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
-  fHistNEvents->Fill(0); // count event
-  // Post the data already here
-  PostData(1,fOutput);
   
 
   TClonesArray *array3Prong = 0;
   TClonesArray *arrayLikeSign =0;
@@ -235,8 +746,8 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
     // 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
+     // 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()) {
@@ -256,9 +767,17 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   }
   if(!arrayLikeSign) {
     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
-    // do in any case the OS analysis.
+    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;
+  fCounter->StoreEvent(aod,fReadMC);
+  fHistNEvents->Fill(0); // count event
+  // Post the data already here
+  PostData(1,fOutput);
  
   TClonesArray *arrayMC=0;
   AliAODMCHeader *mcHeader=0;
@@ -273,14 +792,14 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
     if(!arrayMC) {
       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
-      return;
+      //    return;
     }
     
   // load MC header
     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
     if(!mcHeader) {
-      printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
-      return;
+    printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
+    return;
     }
   }
   
@@ -289,8 +808,11 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
   
   
   Int_t nOS=0;
+  Int_t index;
   Int_t pdgDgDplustoKpipi[3]={321,211,211};
-  Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
+  // 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);
     
@@ -301,36 +823,20 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
       unsetvtx=kTRUE;
     }
 
-    if(d->SelectDplus(fVHF->GetDplusCuts())) {
-      Int_t iPtBin=0;
+    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) {
+
+  
+
+      Int_t iPtBin = -1;
       Double_t ptCand = d->Pt();
-      if(ptCand<2.){
-       iPtBin=0;
-       cutsDplus[7]=0.08;
-       cutsDplus[8]=0.5;
-       cutsDplus[9]=0.979;
-       cutsDplus[10]=0.0000055;
-      }
-      else if(ptCand>2. && ptCand<3){ 
-       iPtBin=1;
-       cutsDplus[7]=0.08;
-       cutsDplus[8]=0.5;
-       cutsDplus[9]=0.991;
-       cutsDplus[10]=0.000005;
-      }else if(ptCand>3. && ptCand<5){ 
-       iPtBin=2;
-       cutsDplus[7]=0.1;
-       cutsDplus[8]=0.5;
-       cutsDplus[9]=0.9955;
-       cutsDplus[10]=0.0000035;
-      }else{
-       iPtBin=3;
-       cutsDplus[7]=0.1;
-       cutsDplus[8]=0.5;
-       cutsDplus[9]=0.997;
-       cutsDplus[10]=0.000001;
+
+      for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
+       if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
       }
-      Bool_t passTightCuts=d->SelectDplus(cutsDplus);
+      
+      Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
+     
+
       Int_t labDp=-1;
       Float_t deltaPx=0.;
       Float_t deltaPy=0.;
@@ -358,20 +864,15 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        }
       }
       Double_t invMass=d->InvMassDplus();
-
-      TString hisNameA(Form("hMassPt%d",iPtBin));
-      TString hisNameS(Form("hSigPt%d",iPtBin));
-      TString hisNameB(Form("hBkgPt%d",iPtBin));
-      TString hisNameATC(Form("hMassPt%dTC",iPtBin));
-      TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
-      TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
-      
-      ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
-      if(passTightCuts){
-       ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
+      Double_t rapid=d->YDplus();
+      fYVsPt->Fill(ptCand,rapid);
+      if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid);
+      Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
+      if(isFidAcc){
+       fPtVsMass->Fill(invMass,ptCand);
+       if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
       }
-
-      Float_t tmp[22];
+      Float_t tmp[24];
       if(fFillNtuple){           
        tmp[0]=pdgCode;
        tmp[1]=deltaPx;
@@ -394,159 +895,339 @@ void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
        tmp[18]=d->GetSigmaVert();
        tmp[19]=d->Getd0Prong(0);
        tmp[20]=d->Getd0Prong(1);
-       tmp[21]=d->Getd0Prong(2);         
+       tmp[21]=d->Getd0Prong(2);
+       tmp[22]=d->GetDCA();
+       tmp[23]=d->Prodd0d0(); 
        fNtupleDplus->Fill(tmp);
-       PostData(2,fNtupleDplus);
+       PostData(3,fNtupleDplus);
       }
-
-      if(fReadMC){
-       if(labDp>=0) {
-         ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
-         if(passTightCuts){
-           ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
-         }
+      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);
+      }
+      if(iPtBin>=0){
+      
+       index=GetHistoIndex(iPtBin);
+       if(isFidAcc){
+         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);
          
-       }else{
-         ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
          if(passTightCuts){
-           ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
+           nSelectedtight++;
+           fMassHistTC[index]->Fill(invMass);
+           if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
+           else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
          }
-         fHistOSbkg->Fill(invMass);
        }
-      }
-
-      fHistOS->Fill(invMass);
-      nOS++;
+       
+       if(fReadMC){
+         if(labDp>=0) {
+           index=GetSignalHistoIndex(iPtBin);
+           if(isFidAcc){
+             Float_t factor[3]={1.,1.,1.};
+             if(fUseStrangeness){
+               for(Int_t iprong=0;iprong<3;iprong++){
+                 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];
+             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(passTightCuts){
+               fMassHistTC[index]->Fill(invMass);            
+               if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
+               else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
+             }
+           }       
+           fYVsPtSig->Fill(ptCand,rapid);
+           if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
+         }else{
+           index=GetBackgroundHistoIndex(iPtBin);
+           if(isFidAcc){
+             Float_t factor[3]={1.,1.,1.};
+             if(fUseStrangeness){
+               for(Int_t iprong=0;iprong<3;iprong++){
+                 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];
+             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(passTightCuts){
+               fMassHistTC[index]->Fill(invMass);
+               if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
+               else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
+             }
+           }   
+         }
+       }
+      }  
     }
     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(3,fCounter);    
   return;
 }
 
 
 
-//_________________________________________________________________
-void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
-
-/*
- * Fill the Like Sign histograms
- */
-
-
-
-
-  //count pos/neg tracks
-  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++;
-      }
-    }else{      
-      nNegTrks++;
-      if(track->Pt()>=0.4){
-       nspcminus++;
-      }
-    }
-  }
-
-  Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
-
-  Int_t nDplusLS=0;
-  Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
-  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;
-    }
-    if(d->SelectDplus(fVHF->GetDplusCuts()))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  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;
-    }
-    if(d->SelectDplus(fVHF->GetDplusCuts())){
-      Int_t sign= d->GetCharge();
-      Float_t wei1=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
-       wei1=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
-       wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
-      }
-      if(sign<0&&nNegTrks>2&&nspcminus>2){
-       wei1=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
-       wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
-      }
-      Double_t invMass=d->InvMassDplus();
-      fHistLS->Fill(invMass,wei1);
-      fHistLSnoweight->Fill(invMass);
-      fHistLStrip->Fill(invMass,wei2);
-      fHistLStripcut->Fill(invMass,wei3);
-      fHistLSsinglecut->Fill(invMass,wei4);
-     }
-    if(unsetvtx) d->UnsetOwnPrimaryVtx();
-  }
-  
-  //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");
-
-}
-
-//_________________________________________________________________
-
-
+//________________________________________________________________________
 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(fFillNtuple){
-    fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
-  }
-  fHistOS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOS"));
-  fHistOSbkg =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOSbkg"));
-  if(fDoLS){
-    fHistLS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLS"));
-    fHistLStrip =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStrip"));
-    fHistLStripcut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStripcut"));
-    fHistLSnoweight =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSnoweight"));
-    fHistLSsinglecut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSsinglecut"));
+  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);
+    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");
+  hMassPt->SetLineColor(kBlue);
+  hMassPt->SetXTitle("M[GeV/c^{2}]"); 
+  hMassPt->Draw();
   return;
 }