]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
Fix (AndreaR)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSED0Correlations.cxx
index 86232993345e83cb539546d9bcd78c7964e2c6d3..409b6d76f4f9862faedd03b2aa3145a4c3b9d5ce 100644 (file)
@@ -52,6 +52,7 @@
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisTaskSED0Correlations.h"
 #include "AliNormalizationCounter.h"
+#include "AliVertexingHFUtils.h"
 
 using std::cout;
 using std::endl;
@@ -74,14 +75,26 @@ AliAnalysisTaskSE(),
   fNentries(0), 
   fCutsD0(0),
   fCutsTracks(0),
+  fCorrelatorTr(0),
+  fCorrelatorKc(0),
+  fCorrelatorK0(0),
   fReadMC(0),
+  fRecoTr(kTRUE),
+  fRecoD0(kTRUE),
+  fSelEvType(kFALSE),
+  fMixing(kFALSE),
   fCounter(0),
   fNPtBins(1),
   fFillOnlyD0D0bar(0),
   fIsSelectedCandidate(0),
   fSys(0),
+  fEtaForCorrel(0),
   fIsRejectSDDClusters(0),
-  fFillGlobal(kTRUE)
+  fFillGlobal(kFALSE),
+  fMultEv(0.),
+  fSoftPiCut(kTRUE),
+  fMEAxisThresh(kFALSE),
+  fKaonCorr(kFALSE)   
 {
   // Default constructor
 
@@ -102,14 +115,26 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *nam
   fNentries(0),
   fCutsD0(0),
   fCutsTracks(cutsTrk),
+  fCorrelatorTr(0),
+  fCorrelatorKc(0),
+  fCorrelatorK0(0),
   fReadMC(0),
+  fRecoTr(kTRUE),
+  fRecoD0(kTRUE),
+  fSelEvType(kFALSE),
+  fMixing(kFALSE),
   fCounter(0),
   fNPtBins(1),
   fFillOnlyD0D0bar(0),
   fIsSelectedCandidate(0),
   fSys(0),
+  fEtaForCorrel(0),
   fIsRejectSDDClusters(0),
-  fFillGlobal(kTRUE)
+  fFillGlobal(kFALSE),
+  fMultEv(0.),
+  fSoftPiCut(kTRUE),
+  fMEAxisThresh(kFALSE),
+  fKaonCorr(kFALSE)
 {
   // Default constructor
 
@@ -148,14 +173,26 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalys
   fNentries(source.fNentries), 
   fCutsD0(source.fCutsD0),
   fCutsTracks(source.fCutsTracks),
+  fCorrelatorTr(source.fCorrelatorTr),
+  fCorrelatorKc(source.fCorrelatorKc),
+  fCorrelatorK0(source.fCorrelatorK0),
   fReadMC(source.fReadMC),
+  fRecoTr(source.fRecoTr),
+  fRecoD0(source.fRecoD0),
+  fSelEvType(source.fSelEvType),
+  fMixing(source.fMixing),
   fCounter(source.fCounter),
   fNPtBins(source.fNPtBins),
   fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
   fIsSelectedCandidate(source.fIsSelectedCandidate),
   fSys(source.fSys),
+  fEtaForCorrel(source.fEtaForCorrel),
   fIsRejectSDDClusters(source.fIsRejectSDDClusters),
-  fFillGlobal(source.fFillGlobal)
+  fFillGlobal(source.fFillGlobal),
+  fMultEv(source.fMultEv),
+  fSoftPiCut(source.fSoftPiCut),
+  fMEAxisThresh(source.fMEAxisThresh),
+  fKaonCorr(source.fKaonCorr)
 {
   // Copy constructor
 }
@@ -183,7 +220,19 @@ AliAnalysisTaskSED0Correlations::~AliAnalysisTaskSED0Correlations()
     delete fNentries;
     fNentries = 0;
   }
-  if(fCounter){
+  if (fCorrelatorTr) {
+    delete fCorrelatorTr;
+    fCorrelatorTr = 0;
+  }
+  if (fCorrelatorKc) {
+    delete fCorrelatorKc;
+    fCorrelatorKc = 0;
+  }
+  if (fCorrelatorK0) {
+    delete fCorrelatorK0;
+    fCorrelatorK0 = 0;
+  }
+  if (fCounter){
     delete fCounter;
     fCounter=0;
   }
@@ -208,14 +257,26 @@ AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(cons
   fNentries = orig.fNentries; 
   fCutsD0 = orig.fCutsD0;
   fCutsTracks = orig.fCutsTracks;
+  fCorrelatorTr = orig.fCorrelatorTr;
+  fCorrelatorKc = orig.fCorrelatorKc;
+  fCorrelatorK0 = orig.fCorrelatorK0;
   fReadMC = orig.fReadMC;
+  fRecoTr = orig.fRecoTr;
+  fRecoD0 = orig.fRecoD0;
+  fSelEvType = orig.fSelEvType;
+  fMixing = orig.fMixing;
   fCounter = orig.fCounter;
   fNPtBins = orig.fNPtBins;
   fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
   fIsSelectedCandidate = orig.fIsSelectedCandidate;
   fSys = orig.fSys;
+  fEtaForCorrel = orig.fEtaForCorrel;
   fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
   fFillGlobal = orig.fFillGlobal;
+  fMultEv = orig.fMultEv;
+  fSoftPiCut = orig.fSoftPiCut;
+  fMEAxisThresh = orig.fMEAxisThresh;
+  fKaonCorr = orig.fKaonCorr;
 
   return *this; //returns pointer of the class
 }
@@ -232,6 +293,7 @@ void AliAnalysisTaskSED0Correlations::Init()
   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
   copyfCutsD0->SetName(nameoutput);
 
+  //needer to clear completely the objects inside with Clear() method
   // Post the data
   PostData(3,copyfCutsD0);
   PostData(7,fCutsTracks);
@@ -247,6 +309,36 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
   //
   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
 
+  //HFCorrelator creation and definition
+  fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys,fCutsD0);//fSys=0 use multiplicity, =1 use centrality
+  fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys,fCutsD0);
+  fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys,fCutsD0);
+  fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);// set the Delta Phi Interval you want (in this case -0.5Pi to 1.5 Pi)
+  fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
+  fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
+  fCorrelatorTr->SetEventMixing(fMixing);// sets the analysis on a single event (kFALSE) or mixed events (kTRUE)
+  fCorrelatorKc->SetEventMixing(fMixing);
+  fCorrelatorK0->SetEventMixing(fMixing);
+  fCorrelatorTr->SetAssociatedParticleType(1);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
+  fCorrelatorKc->SetAssociatedParticleType(2);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
+  fCorrelatorK0->SetAssociatedParticleType(3);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
+  fCorrelatorTr->SetApplyDisplacementCut(2); //0: don't calculate d0; 1: return d0; 2: return d0/d0err
+  fCorrelatorKc->SetApplyDisplacementCut(2);
+  fCorrelatorK0->SetApplyDisplacementCut(0);
+  fCorrelatorTr->SetUseMC(fReadMC);// sets Montecarlo flag
+  fCorrelatorKc->SetUseMC(fReadMC);
+  fCorrelatorK0->SetUseMC(fReadMC);
+  fCorrelatorTr->SetUseReco(fRecoTr);// sets (if MC analysis) wheter to analyze Reco or Kinem tracks
+  fCorrelatorKc->SetUseReco(fRecoTr);
+  fCorrelatorK0->SetUseReco(fRecoTr);
+  fCorrelatorKc->SetPIDmode(2); //switch for K+/- PID option
+  Bool_t pooldefTr = fCorrelatorTr->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
+  Bool_t pooldefKc = fCorrelatorKc->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
+  Bool_t pooldefK0 = fCorrelatorK0->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
+  if(!pooldefTr) AliInfo("Warning:: Event pool not defined properly");
+  if(!pooldefKc) AliInfo("Warning:: Event pool not defined properly");
+  if(!pooldefK0) AliInfo("Warning:: Event pool not defined properly");
+
   // Several histograms are more conveniently managed in a TList
   fOutputMass = new TList();
   fOutputMass->SetOwner();
@@ -260,45 +352,100 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
   fOutputStudy->SetOwner();
   fOutputStudy->SetName("MCstudyplots");
 
-  TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ";
+  TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassWg=" ",nameSgnWg=" ", nameBkgWg=" ", nameRflWg=" ";
 
+//for origin c case (or for data)
   for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
 
-    nameMass="histMass_";
+    nameMass="histMass_"; if(fReadMC) nameMass+="c_";
     nameMass+=i;
-    nameSgn="histSgn_";
+    nameMassWg="histMass_WeigD0Eff_"; if(fReadMC) nameMassWg+="c_";
+    nameMassWg+=i;
+    nameSgn="histSgn_"; if(fReadMC) nameSgn+="c_";
     nameSgn+=i;
-    nameBkg="histBkg_";
+    nameSgnWg="histSgn_WeigD0Eff_"; if(fReadMC) nameSgnWg+="c_";
+    nameSgnWg+=i;
+    nameBkg="histBkg_"; if(fReadMC) nameBkg+="c_";
     nameBkg+=i;
-    nameRfl="histRfl_";
+    nameBkgWg="histBkg_WeigD0Eff_"; if(fReadMC) nameBkgWg+="c_";
+    nameBkgWg+=i;
+    nameRfl="histRfl_"; if(fReadMC) nameRfl+="c_";
     nameRfl+=i;
+    nameRflWg="histRfl_WeigD0Eff_"; if(fReadMC) nameRflWg+="c_";
+    nameRflWg+=i;
 
     //histograms of invariant mass distributions
 
     //MC signal
     if(fReadMC){
-      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass c - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
       tmpSt->Sumw2();
+      tmpStWg->Sumw2();
 
       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
-      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
-      TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass c - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpBtWg = new TH1F(nameBkgWg.Data(), "Background invariant mass c - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
       tmpBt->Sumw2();
+      tmpBtWg->Sumw2();
       tmpRt->Sumw2();
+      tmpRtWg->Sumw2();
       fOutputMass->Add(tmpSt);
+      fOutputMass->Add(tmpStWg);
       fOutputMass->Add(tmpRt);
+      fOutputMass->Add(tmpRtWg);
       fOutputMass->Add(tmpBt);
+      fOutputMass->Add(tmpBtWg);
     }
 
     //mass
-    TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",120,1.5648,2.1648);
+    TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass c; M [GeV]; Entries",120,1.5648,2.1648);
     tmpMt->Sumw2();
     fOutputMass->Add(tmpMt);
+    //mass weighted by 1/D0eff
+    TH1F* tmpMtwg = new TH1F(nameMassWg.Data(),"D^{0} invariant mass c - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+    tmpMtwg->Sumw2();
+    fOutputMass->Add(tmpMtwg);
+  }
+
+//for origin b case (no Bkg and Mass histos, here for weights you should use c+b efficiencies, while on data (on MC they're useless))
+  for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
+
+    nameSgn="histSgn_b_";
+    nameSgn+=i;
+    nameSgnWg="histSgn_WeigD0Eff_b_";
+    nameSgnWg+=i;
+    nameRfl="histRfl_b_";
+    nameRfl+=i;
+    nameRflWg="histRfl_WeigD0Eff_b_";
+    nameRflWg+=i;
+
+    //histograms of invariant mass distributions
+
+    //MC signal
+    if(fReadMC){
+      TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpStWg = new TH1F(nameSgnWg.Data(), "D^{0} invariant mass b - MC; M [GeV] - weight 1/D0eff; Entries",120,1.5648,2.1648);
+      tmpSt->Sumw2();
+      tmpStWg->Sumw2();
+
+      //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
+      TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass b - MC; M [GeV]; Entries",120,1.5648,2.1648);
+      TH1F* tmpRtWg = new TH1F(nameRflWg.Data(), "Reflected signal invariant mass b - MC - weight 1/D0eff; M [GeV]; Entries",120,1.5648,2.1648);
+      tmpRt->Sumw2();
+      tmpRtWg->Sumw2();
+      fOutputMass->Add(tmpSt);
+      fOutputMass->Add(tmpStWg);
+      fOutputMass->Add(tmpRt);
+      fOutputMass->Add(tmpRtWg);
+    }
   }
 
   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 
-  fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 18,-0.5,17.5);
+  fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 20,-0.5,19.5);
 
   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
@@ -316,6 +463,8 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
+  fNentries->GetXaxis()->SetBinLabel(19,"nEventsSelected");
+  if(fReadMC) fNentries->GetXaxis()->SetBinLabel(20,"nEvsWithProdMech");
   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 
   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
@@ -352,7 +501,6 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
   
-
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
   fEvents++;
 
@@ -360,6 +508,8 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
 
   TClonesArray *inputArray=0;
 
+  fMultEv = 0.; //reset event multiplicity
+
   if(!aod && AODEvent() && IsStandardAOD()) {
     // In case there is an AOD handler writing a standard AOD, use the AOD 
     // event in memory rather than the input (ESD) event.    
@@ -405,7 +555,7 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
       return;
     }
   }
-  
+
   //histogram filled with 1 for every AOD
   fNentries->Fill(0);
   fCounter->StoreEvent(aod,fCutsD0,fReadMC); 
@@ -422,6 +572,37 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
     return;
   }
 
+  fNentries->Fill(18); //event selected after selection
+
+  //Setting PIDResponse for associated tracks
+  fCorrelatorTr->SetPidAssociated();
+  fCorrelatorKc->SetPidAssociated();
+  fCorrelatorK0->SetPidAssociated();
+
+  //Selection on production type (MC)
+  if(fReadMC && fSelEvType){ 
+
+    Bool_t isMCeventgood = kFALSE;
+            
+    Int_t eventType = mcHeader->GetEventType();
+    Int_t NMCevents = fCutsTracks->GetNofMCEventType();
+               
+    for(Int_t k=0; k<NMCevents; k++){
+      Int_t * MCEventType = fCutsTracks->GetMCEventType();
+          
+      if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
+      ((TH1D*)fOutputStudy->FindObject("EventTypeMC"))->Fill(eventType);
+    }
+                
+    if(NMCevents && !isMCeventgood){
+      if(fDebug>2)std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
+      return; 
+    }
+    fNentries->Fill(19); //event with particular production type                
+  
+  } //end of selection
+
+
   // Check the Nb of SDD clusters
   if (fIsRejectSDDClusters) { 
     Bool_t skipEvent = kFALSE;
@@ -439,22 +620,36 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
     if (skipEvent) return;
   }
   
+  //HFCorrelators initialization (for this event)
+  fCorrelatorTr->SetAODEvent(aod); // set the AOD event from which you are processing
+  fCorrelatorKc->SetAODEvent(aod);
+  fCorrelatorK0->SetAODEvent(aod);
+  Bool_t correlatorONTr = fCorrelatorTr->Initialize(); // initialize the pool for event mixing
+  Bool_t correlatorONKc = fCorrelatorKc->Initialize();
+  Bool_t correlatorONK0 = fCorrelatorK0->Initialize();
+  if(!correlatorONTr) {AliInfo("AliHFCorrelator (tracks) didn't initialize the pool correctly or processed a bad event"); return;}
+  if(!correlatorONKc) {AliInfo("AliHFCorrelator (charged K) didn't initialize the pool correctly or processed a bad event"); return;}
+  if(!correlatorONK0) {AliInfo("AliHFCorrelator (K0) didn't initialize the pool correctly or processed a bad event"); return;}
+
+  if(fReadMC) {
+    fCorrelatorTr->SetMCArray(mcArray); // set the TClonesArray *fmcArray for analysis on monte carlo
+    fCorrelatorKc->SetMCArray(mcArray);
+    fCorrelatorK0->SetMCArray(mcArray);
+  }
+
   // AOD primary vertex
   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 
-  Bool_t isGoodVtx=kFALSE;
-
   //vtx1->Print();
   TString primTitle = vtx1->GetTitle();
   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
-    isGoodVtx=kTRUE;
     fNentries->Fill(3);
   }
 
   //Reset flag for tracks distributions fill
   fAlreadyFilled=kFALSE;
 
-  // loop over candidates
+  //***** Loop over D0 candidates *****
   Int_t nInD0toKpi = inputArray->GetEntriesFast();
   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 
@@ -467,9 +662,9 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
       for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
       AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
       if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
-        if(fReadMC && (v0->MatchToMC(311,mcArray,2,pdgCodes)<0)) continue;
+        if(fReadMC && fRecoTr && (v0->MatchToMC(310,mcArray,2,pdgCodes)<0)) continue; //310 = K0s, 311 = K0 generico!!
         ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
-        ((TH1F*)fOutputStudy->FindObject("hist_Pt_K_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
+        ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
       }
     }
 
@@ -478,7 +673,7 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
       //rejection of tracks
       if(track->GetID() < 0) continue; //discard negative ID tracks
       if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
-      if(!fCutsTracks->IsHadronSelected(track,vtx1,aod->GetMagneticField())) continue;
+      if(!fCutsTracks->IsHadronSelected(track) || !fCutsTracks->CheckHadronKinematic(track->Pt(),0.1)) continue; //0.1 = dummy (d0 value, no cut on it for me)
       //pT distribution (in all events), charg and hadr cases
       ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt()); 
       if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
@@ -487,37 +682,147 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
   } //end of loops for global plot fill
 
   Int_t nSelectedloose=0,nSelectedtight=0;  
-  for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
-    AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
+  
+  //Fill Event Multiplicity (needed only in Reco)
+  fMultEv = (Double_t)(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.));
+
+  //RecoD0 case ************************************************
+  if(fRecoD0) {
+
+    for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
+      AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
  
-    if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
-       fNentries->Fill(2);
-       continue; //skip the D0 from Dstar
+      if(d->Pt()<2.) continue; //to save time and merging memory...
+
+      if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
+       fNentries->Fill(2);
+       continue; //skip the D0 from Dstar  
       }
     
-    if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
-      nSelectedloose++;
-      nSelectedtight++;      
-      if(fSys==0){
-       if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);       
-      }
-      Int_t ptbin=fCutsD0->PtBin(d->Pt());
-      if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
-
-      fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
-      if(!fIsSelectedCandidate) continue;
+      if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
+        nSelectedloose++;
+        nSelectedtight++;      
+        if(fSys==0){
+         if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);       
+        }  
+        Int_t ptbin=fCutsD0->PtBin(d->Pt());
+        if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
+
+        fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
+        if(!fIsSelectedCandidate) continue;
+
+        //D0 infos
+        Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(d->Phi());
+                 phiD0 = fCorrelatorKc->SetCorrectPhiRange(d->Phi());  //bad usage, but returns a Double_t...
+                 phiD0 = fCorrelatorK0->SetCorrectPhiRange(d->Phi());
+        fCorrelatorTr->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta()); // sets the parameters of the trigger particles that are needed
+        fCorrelatorKc->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
+        fCorrelatorK0->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
+        fCorrelatorTr->SetD0Properties(d,fIsSelectedCandidate); //sets special properties for D0
+        fCorrelatorKc->SetD0Properties(d,fIsSelectedCandidate);
+        fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
+
+        if(!fReadMC) {
+          if (TMath::Abs(d->Eta())<fEtaForCorrel) {
+           CalculateCorrelations(d); //correlations on real data
+           if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv);
+         }
+        } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
+          if (TMath::Abs(d->Eta())<fEtaForCorrel) {
+            Int_t pdgDgD0toKpi[2]={321,211};
+           Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
+            if (labD0>-1) {
+             CalculateCorrelations(d,labD0,mcArray);
+             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
+           }
+          }
+        }
 
-      if(!fReadMC) CalculateCorrelations(aod,d); //correlations on real data
-      else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
-        Int_t pdgDgD0toKpi[2]={321,211};
-       Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
-        if (labD0>-1) CalculateCorrelations(aod,d,labD0,mcArray);
+        FillMassHists(d,mcArray,fCutsD0,fOutputMass);
       }
-      
-      FillMassHists(d,mcArray,fCutsD0,fOutputMass);
     }
+  }
+  //End RecoD0 case ************************************************
+  
+  //MCKineD0 case ************************************************
+  if(fReadMC && !fRecoD0) {
+
+    for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { //Loop over all the tracks of MCArray
+      AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
+      if (!mcPart) {
+        AliWarning("Particle not found in tree, skipping"); 
+        continue;
+      } 
+  
+      if(TMath::Abs(mcPart->GetPdgCode()) == 421){  // THIS IS A D0
+        if (fCutsD0->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y()) ) {
+          nSelectedloose++;
+          nSelectedtight++;      
+
+         //Removal of cases in which D0 decay is not in Kpi!
+         if(mcPart->GetNDaughters()!=2) continue;
+         AliAODMCParticle* mcDau1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
+         AliAODMCParticle* mcDau2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)));
+         if(!mcDau1 || !mcDau2) continue;
+         Int_t pdg1 = TMath::Abs(mcDau1->GetPdgCode());
+         Int_t pdg2 = TMath::Abs(mcDau2->GetPdgCode());
+          if(!((pdg1 == 211 && pdg2 == 321) || (pdg2 == 211 && pdg1 == 321))) continue;
+          if(TMath::Abs(mcDau1->Eta())>0.8||TMath::Abs(mcDau2->Eta())>0.8) continue;
+            //Check momentum conservation (to exclude 4-prong decays with tracks outside y=1.5)
+            Double_t p1[3]  = {mcDau1->Px(),mcDau1->Py(),mcDau1->Pz()};
+            Double_t p2[3]  = {mcDau2->Px(),mcDau2->Py(),mcDau2->Pz()};
+            Double_t pD0[3] = {mcPart->Px(),mcPart->Py(),mcPart->Pz()};
+            if(TMath::Abs( (p1[0]+p2[0]-pD0[0])*(p1[0]+p2[0]-pD0[0]) + (p1[1]+p2[1]-pD0[1])*(p1[1]+p2[1]-pD0[1]) + (p1[2]+p2[2]-pD0[2])*(p1[2]+p2[2]-pD0[2]) )>0.1) continue;
+
+          if(fSys==0) fNentries->Fill(6);
+          Int_t ptbin=fCutsD0->PtBin(mcPart->Pt());
+          if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds  
+  
+          //D0 infos
+          Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(mcPart->Phi());
+                   phiD0 = fCorrelatorKc->SetCorrectPhiRange(mcPart->Phi());  //bad usage, but returns a Double_t...
+                   phiD0 = fCorrelatorK0->SetCorrectPhiRange(mcPart->Phi());
+          fCorrelatorTr->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta()); // sets the parameters of the trigger particles that are needed
+          fCorrelatorKc->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
+          fCorrelatorK0->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
+          //fCorrelatorTr->SetD0Properties(mcPart,fIsSelectedCandidate); //needed for D* soft pions rejection, useless in MCKine
+          //fCorrelatorKc->SetD0Properties(mcPart,fIsSelectedCandidate);
+          //fCorrelatorK0->SetD0Properties(mcPart,fIsSelectedCandidate);
+  
+          if (TMath::Abs(mcPart->Eta())<fEtaForCorrel) {
+  
+            //Removal of D0 from D* feeddown! This solves also the problem of soft pions, now excluded
+         /*   Int_t mother = mcPart->GetMother();
+           AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
+            if(!mcMoth) continue;
+           if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;*/
+
+            if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
+           else fIsSelectedCandidate = 2;
+
+           TString fillthis="histSgn_"; 
+           if(CheckD0Origin(mcArray,mcPart)==4) fillthis+="c_";
+           else if(CheckD0Origin(mcArray,mcPart)==5) fillthis+="b_";
+            else continue;
+            fillthis+=ptbin;
+           ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
+         
+            CalculateCorrelationsMCKine(mcPart,mcArray);
+            if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
+          }
+        }
+      }
+    }
+  }
+  //End MCKineD0 case ************************************************
+
+  if(fMixing /* && fAlreadyFilled*/) { // update the pool for Event Mixing, if: enabled,  event is ok, at least a SelD0 found! (fAlreadyFilled's role!)
+    Bool_t updatedTr = fCorrelatorTr->PoolUpdate();
+    Bool_t updatedKc = fCorrelatorKc->PoolUpdate();
+    Bool_t updatedK0 = fCorrelatorK0->PoolUpdate();
+    if(!updatedTr || !updatedKc || !updatedK0) AliInfo("Pool was not updated");
+  }
 
-  } //end for prongs
   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);  
   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);  
 
@@ -556,56 +861,128 @@ void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *par
   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
 
     if(fReadMC){
-      if(labD0>=0) {
-       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+      if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
+       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
        Int_t pdgD0 = partD0->GetPdgCode();
        if (pdgD0==421){ //D0
-         fillthis="histSgn_";
+         fillthis="histSgn_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
-       } else{ //it was a D0bar
-         fillthis="histRfl_";
+          fillthis="histSgn_WeigD0Eff_c_";
+          fillthis+=ptbin;
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+       } else{ //it was a D0bar
+         fillthis="histRfl_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
-       }
+          fillthis="histRfl_WeigD0Eff_c_";
+          fillthis+=ptbin;
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+       }
+      } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
+         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+         Int_t pdgD0 = partD0->GetPdgCode();
+         if (pdgD0==421){ //D0
+           fillthis="histSgn_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+            fillthis="histSgn_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+         } else{ //it was a D0bar
+           fillthis="histRfl_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+            fillthis="histRfl_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
+         }
       } else {//background
-       fillthis="histBkg_";
+       fillthis="histBkg_c_";
        fillthis+=ptbin;
        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+        fillthis="histBkg_WeigD0Eff_c_";
+        fillthis+=ptbin;
+        Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
       }
     }else{
       fillthis="histMass_";
       fillthis+=ptbin;
       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+      fillthis="histMass_WeigD0Eff_";
+      fillthis+=ptbin;
+      Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
     }
      
   }
   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
 
     if(fReadMC){
-      if(labD0>=0) {
-       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+      if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
+       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
        Int_t pdgD0 = partD0->GetPdgCode();
-
-       if (pdgD0==-421){ //D0bar
-         fillthis="histSgn_";
+       if (pdgD0==-421){ //D0
+         fillthis="histSgn_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
-       } else{
-         fillthis="histRfl_";
+          fillthis="histSgn_WeigD0Eff_c_";
+          fillthis+=ptbin;
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+       } else{ //it was a D0bar
+         fillthis="histRfl_c_";
          fillthis+=ptbin;
          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
-       }
-      } else {//background or LS
-       fillthis="histBkg_";
+          fillthis="histRfl_WeigD0Eff_c_";
+          fillthis+=ptbin;
+          Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+          ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+       }
+      } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
+         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
+         Int_t pdgD0 = partD0->GetPdgCode();
+         if (pdgD0==-421){ //D0
+           fillthis="histSgn_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+            fillthis="histSgn_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+         } else{ //it was a D0bar
+           fillthis="histRfl_b_";
+           fillthis+=ptbin;
+           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+            fillthis="histRfl_WeigD0Eff_b_";
+            fillthis+=ptbin;
+            Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
+            ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
+         }
+      } else {//background
+       fillthis="histBkg_c_";
        fillthis+=ptbin;
        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+        fillthis="histBkg_WeigD0Eff_c_";
+        fillthis+=ptbin;
+        Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
       }
     }else{
       fillthis="histMass_";
       fillthis+=ptbin;
-      ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
+      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+      fillthis="histMass_WeigD0Eff_";
+      fillthis+=ptbin;
+      Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
+      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
     }
+
   }
 
   return;
@@ -666,7 +1043,7 @@ Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliA
   //
   // checking whether the mother of the particles come from a charm or a bottom quark
   //
-  printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
+  printf("AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
        
   Int_t pdgGranma = 0;
   Int_t mother = 0;
@@ -676,15 +1053,15 @@ Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliA
   Bool_t isQuarkFound=kFALSE;
 
   while (mother > 0){
-    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
-    if (mcGranma){
-      pdgGranma = mcGranma->GetPdgCode();
+    AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
+    if (mcMoth){
+      pdgGranma = mcMoth->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();
+      mother = mcMoth->GetMother();
     }else{
       AliError("Failed casting the mother particle!");
       break;
@@ -705,276 +1082,235 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
   TString namePlot = "";
 
   //These for limits in THnSparse (one per bin, same limits). 
-  //Vars: DeltaPhi, InvMass, PtTrack
-  Int_t nBinsPhi[4] = {32,150,30,3};
-  Double_t binMinPhi[4] = {-1.6,1.6,0.,0.};  //is the minimum for all the bins
-  Double_t binMaxPhi[4] = {4.8,2.2,3.0,3.};  //is the maximum for all the bins
-
-  for(Int_t i=0;i<fNPtBinsCorr;i++){
-
-    //THnSparse plots: correlations for various invariant mass (MC and data)
-    namePlot="hPhi_K_Bin";
-    namePlot+=i;
+  //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
+  Int_t nBinsPhi[5] = {32,150,6,3,16};
+  Double_t binMinPhi[5] = {-TMath::Pi()/2.,1.6,0.,0.,-1.6};  //is the minimum for all the bins
+  Double_t binMaxPhi[5] = {3.*TMath::Pi()/2.,2.2,3.0,3.,1.6};  //is the maximum for all the bins
 
-    THnSparseI *hPhiK = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-    hPhiK->Sumw2();
-    fOutputCorr->Add(hPhiK);
+  //Vars: DeltaPhi, InvMass, DeltaEta
+  Int_t nBinsMix[4] = {32,150,16,6};
+  Double_t binMinMix[4] = {-TMath::Pi()/2.,1.6,-1.6,0.};  //is the minimum for all the bins
+  Double_t binMaxMix[4] = {3.*TMath::Pi()/2.,2.2,1.6,3.};  //is the maximum for all the bins
 
-    namePlot="hPhi_Kcharg_Bin";
-    namePlot+=i;
-
-    THnSparseI *hPhiH = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-    hPhiH->Sumw2();
-    fOutputCorr->Add(hPhiH);
-
-    namePlot="hPhi_Charg_Bin";
-    namePlot+=i;
-
-    THnSparseI *hPhiC = new THnSparseI(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-    hPhiC->Sumw2();
-    fOutputCorr->Add(hPhiC);
-
-    //histos for c/b origin for D0 (MC only)
-    if (fReadMC) {
+  for(Int_t i=0;i<fNPtBinsCorr;i++){
 
-      //generic origin for tracks
-      namePlot="hPhi_K_From_c_Bin";
+    if(!fMixing) {
+      //THnSparse plots: correlations for various invariant mass (MC and data)
+      namePlot="hPhi_K0_Bin";
       namePlot+=i;
 
-      THnSparseI *hPhiK_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiK_c->Sumw2();
-      fOutputCorr->Add(hPhiK_c);
+      THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+      hPhiK->Sumw2();
+      fOutputCorr->Add(hPhiK);
 
-      namePlot="hPhi_Kcharg_From_c_Bin";
+      namePlot="hPhi_Kcharg_Bin";
       namePlot+=i;
 
-      THnSparseI *hPhiH_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiH_c->Sumw2();
-      fOutputCorr->Add(hPhiH_c);
+      THnSparseF *hPhiH = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+      hPhiH->Sumw2();
+      fOutputCorr->Add(hPhiH);
 
-      namePlot="hPhi_Charg_From_c_Bin";
+      namePlot="hPhi_Charg_Bin";
       namePlot+=i;
 
-      THnSparseI *hPhiC_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_c->Sumw2();
-      fOutputCorr->Add(hPhiC_c);
+      THnSparseF *hPhiC = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+      hPhiC->Sumw2();
+      fOutputCorr->Add(hPhiC);
   
-      namePlot="hPhi_K_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiK_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiK_b->Sumw2();
-      fOutputCorr->Add(hPhiK_b);
-
-      namePlot="hPhi_Kcharg_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiH_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiH_b->Sumw2();
-      fOutputCorr->Add(hPhiH_b);
-
-      namePlot="hPhi_Charg_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiC_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_b->Sumw2();
-      fOutputCorr->Add(hPhiC_b);
-
-      //HF-only tracks (c for c->D0, b for b->D0)
-      namePlot="hPhi_K_HF_From_c_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiK_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiK_HF_c->Sumw2();
-      fOutputCorr->Add(hPhiK_HF_c);
-
-      namePlot="hPhi_Kcharg_HF_From_c_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiH_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiH_HF_c->Sumw2();
-      fOutputCorr->Add(hPhiH_HF_c);
-
-      namePlot="hPhi_Charg_HF_From_c_Bin";
-      namePlot+=i;
-      THnSparseI *hPhiC_HF_c = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_HF_c->Sumw2();
-      fOutputCorr->Add(hPhiC_HF_c);
-
-      namePlot="hPhi_K_HF_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiK_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiK_HF_b->Sumw2();
-      fOutputCorr->Add(hPhiK_HF_b);
-
-      namePlot="hPhi_Kcharg_HF_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiH_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiH_HF_b->Sumw2();
-      fOutputCorr->Add(hPhiH_HF_b);
-
-      namePlot="hPhi_Charg_HF_From_b_Bin";
-      namePlot+=i;
-
-      THnSparseI *hPhiC_HF_b = new THnSparseI(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_HF_b->Sumw2();
-      fOutputCorr->Add(hPhiC_HF_b);
-    }
-
-    //leading hadron correlations
-    namePlot="hPhi_Lead_Bin";
-    namePlot+=i;
-
-    TH2F *hCorrLead = new TH2F(namePlot.Data(), "Leading particle correlation; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-    hCorrLead->Sumw2();
-    fOutputCorr->Add(hCorrLead);
-
-    if (fReadMC) {
-      namePlot="hPhi_Lead_From_c_Bin";
-      namePlot+=i;
-
-      TH2F *hCorrLead_c = new TH2F(namePlot.Data(), "Leading particle correlation - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      hCorrLead_c->Sumw2();
-      fOutputCorr->Add(hCorrLead_c);
+      //histos for c/b origin for D0 (MC only)
+      if (fReadMC) {
 
-      namePlot="hPhi_Lead_From_b_Bin";
-      namePlot+=i;
+        //generic origin for tracks
+        namePlot="hPhi_K0_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrLead_b = new TH2F(namePlot.Data(), "Leading particle correlation - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      hCorrLead_b->Sumw2();
-      fOutputCorr->Add(hCorrLead_b);
+        THnSparseF *hPhiK_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_c->Sumw2();
+        fOutputCorr->Add(hPhiK_c);
 
-      namePlot="hPhi_Lead_HF_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrLead_HF_c = new TH2F(namePlot.Data(), "Leading particle correlation HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      hCorrLead_HF_c->Sumw2();
-      fOutputCorr->Add(hCorrLead_HF_c);
+        THnSparseF *hPhiH_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_c->Sumw2();
+        fOutputCorr->Add(hPhiH_c);
 
-      namePlot="hPhi_Lead_HF_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrLead_HF_b = new TH2F(namePlot.Data(), "Leading particle correlation HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      hCorrLead_HF_b->Sumw2();
-      fOutputCorr->Add(hCorrLead_HF_b);
-    }
-    
-    //pT weighted correlations
-    namePlot="hPhi_Weig_Bin";
-    namePlot+=i;
+        THnSparseF *hPhiC_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_c->Sumw2();
+        fOutputCorr->Add(hPhiC_c);
+  
+        namePlot="hPhi_K0_From_b_Bin";
+        namePlot+=i;
 
-    TH2F *hCorrWeig = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-    fOutputCorr->Add(hCorrWeig);
+        THnSparseF *hPhiK_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_b->Sumw2();
+        fOutputCorr->Add(hPhiK_b);
 
-    if (fReadMC) {
-      namePlot="hPhi_Weig_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_From_b_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeig_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeig_c);
+        THnSparseF *hPhiH_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_b->Sumw2();
+        fOutputCorr->Add(hPhiH_b);
 
-      namePlot="hPhi_Weig_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_From_b_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeig_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeig_b);
+        THnSparseF *hPhiC_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_b->Sumw2();
+        fOutputCorr->Add(hPhiC_b);
 
-      namePlot="hPhi_Weig_HF_From_c_Bin";
-      namePlot+=i;
+        //HF-only tracks (c for c->D0, b for b->D0)
+        namePlot="hPhi_K0_HF_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeig_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeig_HF_c);
+        THnSparseF *hPhiK_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_HF_c->Sumw2();
+        fOutputCorr->Add(hPhiK_HF_c);
 
-      namePlot="hPhi_Weig_HF_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_HF_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeig_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeig_HF_b);
-    }
+        THnSparseF *hPhiH_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_HF_c->Sumw2();
+        fOutputCorr->Add(hPhiH_HF_c);
 
-    //pT weighted correlations (squared weights)
-    namePlot="hPhi_WeigSqr_Bin";
-    namePlot+=i;
+        namePlot="hPhi_Charg_HF_From_c_Bin";
+        namePlot+=i;
 
-    TH2F *hCorrWeigSqr = new TH2F(namePlot.Data(), "Charged particle correlation (pT Weighted - squared); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-    fOutputCorr->Add(hCorrWeigSqr);
+        THnSparseF *hPhiC_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_HF_c->Sumw2();
+        fOutputCorr->Add(hPhiC_HF_c);
 
-    if (fReadMC) {
-      namePlot="hPhi_WeigSqr_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_K0_HF_From_b_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigSqr_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigSqr_c);
+        THnSparseF *hPhiK_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_HF_b->Sumw2();
+        fOutputCorr->Add(hPhiK_HF_b);
 
-      namePlot="hPhi_WeigSqr_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_HF_From_b_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigSqr_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigSqr_b);
+        THnSparseF *hPhiH_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_HF_b->Sumw2();
+        fOutputCorr->Add(hPhiH_HF_b);
 
-      namePlot="hPhi_WeigSqr_HF_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_HF_From_b_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigSqr_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigSqr_HF_c);
+        THnSparseF *hPhiC_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_HF_b->Sumw2();
+        fOutputCorr->Add(hPhiC_HF_b);
 
-      namePlot="hPhi_WeigSqr_HF_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_K0_NonHF_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigSqr_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - squared) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigSqr_HF_b);
-    }
+        THnSparseF *hPhiK_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_Non->Sumw2();
+        fOutputCorr->Add(hPhiK_Non);
 
-    //pT weighted correlations (inverse of pT distribution weights)
-    namePlot="hPhi_WeigDist_Bin";
-    namePlot+=i;
+        namePlot="hPhi_Kcharg_NonHF_Bin";
+        namePlot+=i;
 
-    TH2F *hCorrWeigDist = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.); #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-    fOutputCorr->Add(hCorrWeigDist);
+        THnSparseF *hPhiH_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_Non->Sumw2();
+        fOutputCorr->Add(hPhiH_Non);
 
-    if (fReadMC) {
-      namePlot="hPhi_WeigDist_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_NonHF_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigDist_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigDist_c);
+        THnSparseF *hPhiC_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_Non->Sumw2();
+        fOutputCorr->Add(hPhiC_Non);
+      }
 
-      namePlot="hPhi_WeigDist_From_b_Bin";
+      //leading hadron correlations
+      namePlot="hPhi_Lead_Bin";
       namePlot+=i;
 
-      TH2F *hCorrWeigDist_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigDist_b);
+      THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+      hCorrLead->Sumw2();
+      fOutputCorr->Add(hCorrLead);
 
-      namePlot="hPhi_WeigDist_HF_From_c_Bin";
-      namePlot+=i;
+      if (fReadMC) {
+        namePlot="hPhi_Lead_From_c_Bin";
+        namePlot+=i;
 
-      TH2F *hCorrWeigDist_HF_c = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) HF - c origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigDist_HF_c);
+        THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        hCorrLead_c->Sumw2();
+        fOutputCorr->Add(hCorrLead_c);
+  
+        namePlot="hPhi_Lead_From_b_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        hCorrLead_b->Sumw2();
+        fOutputCorr->Add(hCorrLead_b);
+  
+        namePlot="hPhi_Lead_HF_From_c_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        hCorrLead_HF_c->Sumw2();
+        fOutputCorr->Add(hCorrLead_HF_c);
+  
+        namePlot="hPhi_Lead_HF_From_b_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        hCorrLead_HF_b->Sumw2();
+        fOutputCorr->Add(hCorrLead_HF_b);
 
-      namePlot="hPhi_WeigDist_HF_From_b_Bin";
+        namePlot="hPhi_Lead_NonHF_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        hCorrLead_Non->Sumw2();
+        fOutputCorr->Add(hCorrLead_Non);
+      }
+      
+      //pT weighted correlations
+      namePlot="hPhi_Weig_Bin";
       namePlot+=i;
+  
+      THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+      fOutputCorr->Add(hCorrWeig);
+  
+      if (fReadMC) {
+        namePlot="hPhi_Weig_From_c_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        fOutputCorr->Add(hCorrWeig_c);
+  
+        namePlot="hPhi_Weig_From_b_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        fOutputCorr->Add(hCorrWeig_b);
+  
+        namePlot="hPhi_Weig_HF_From_c_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        fOutputCorr->Add(hCorrWeig_HF_c);
+  
+        namePlot="hPhi_Weig_HF_From_b_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        fOutputCorr->Add(hCorrWeig_HF_b);
 
-      TH2F *hCorrWeigDist_HF_b = new TH2F(namePlot.Data(), "Charged particle correlation (pT weighted - inv.dis.) HF - b origin; #Delta#phi",32,-1.6,4.8,300,1.6,2.2);
-      fOutputCorr->Add(hCorrWeigDist_HF_b);
-    }
-
-    //counters
-    namePlot = "hist_Count_Charg_Bin"; namePlot+=i;
-    TH1F *hCountC = new TH1F(namePlot.Data(), "Charged track counter; # Tracks",100,0.,100.);
-    hCountC->SetMinimum(0);
-    fOutputStudy->Add(hCountC);
-
-    namePlot = "hist_Count_Kcharg_Bin"; namePlot+=i;
-    TH1F *hCountH = new TH1F(namePlot.Data(), "Hadrons counter; # Tracks",100,0.,100.);
-    hCountH->SetMinimum(0);
-    fOutputStudy->Add(hCountH);
-
-    namePlot = "hist_Count_K_Bin"; namePlot+=i;
-    TH1F *hCountK = new TH1F(namePlot.Data(), "Kaons counter; # Tracks",10,0.,10.);
-    hCountK->SetMinimum(0);
-    fOutputStudy->Add(hCountK);
+        namePlot="hPhi_Weig_NonHF_Bin";
+        namePlot+=i;
+  
+        THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+        fOutputCorr->Add(hCorrWeig_Non);
+      }
 
     //pT distribution histos
     namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
@@ -987,26 +1323,72 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     hPtH->SetMinimum(0);
     fOutputStudy->Add(hPtH);
 
-    namePlot = "hist_Pt_K_Bin"; namePlot+=i;
+    namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
     TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
     hPtK->SetMinimum(0);
     fOutputStudy->Add(hPtK);
 
     //D* feeddown pions rejection histos
     namePlot = "hDstarPions_Bin"; namePlot+=i;
-    TH1F *hDstarPions = new TH1F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.);
+    TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,150,1.6,2.2);
     hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
     hDstarPions->SetMinimum(0);
-    fOutputStudy->Add(hDstarPions);
+    fOutputStudy->Add(hDstarPions); 
+
+    //Events multiplicity
+    Double_t yAxisMult[13] = {0, 4, 8, 12, 16, 20, 28, 36, 44, 100}; 
+    namePlot = "hMultiplEvt_Bin"; namePlot+=i;
+    TH1F *hMultEv = new TH1F(namePlot.Data(), "Event multiplicity",9,yAxisMult);
+    hMultEv->SetMinimum(0);
+    fOutputStudy->Add(hMultEv);
+    }
+
+    if(fMixing) {
+      //THnSparse plots for event mixing!
+      namePlot="hPhi_K0_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+
+      THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+      hPhiK_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiK_EvMix);
+
+      namePlot="hPhi_Kcharg_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+  
+      THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+      hPhiH_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiH_EvMix);
+
+      namePlot="hPhi_Charg_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+
+      THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
+      hPhiC_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiC_EvMix);
+    }
   }
 
   //out of bin loop
-  TH1F *hRejTracks = new TH1F("hRejTracks", "Tracks accepted/rejected; # Tracks",2,0.,2.);
-  hRejTracks->SetMinimum(0);
-  hRejTracks->GetXaxis()->SetBinLabel(1,"Accepted");
-  hRejTracks->GetXaxis()->SetBinLabel(2,"Rejected");
-  fOutputStudy->Add(hRejTracks);
+  if(!fMixing) {
+    TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
+    hCountC->SetMinimum(0);
+    fOutputStudy->Add(hCountC);
+
+    TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
+    hCountH->SetMinimum(0);
+    fOutputStudy->Add(hCountH);
+
+    TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
+    hCountK->SetMinimum(0);
+    fOutputStudy->Add(hCountK);
+  }
+
+  if (fReadMC) {
+    TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
+    fOutputStudy->Add(hEventTypeMC); 
+  }
 
   if (fFillGlobal) { //all-events plots
     //pt distributions
@@ -1018,7 +1400,7 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     hPtHAll->SetMinimum(0);
     fOutputStudy->Add(hPtHAll);
 
-    TH1F *hPtKAll = new TH1F("hist_Pt_K_AllEv", "Kaons  pT (All); p_{T} (GeV/c)",240,0.,12.);
+    TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
     hPtKAll->SetMinimum(0);
     fOutputStudy->Add(hPtKAll);
 
@@ -1028,18 +1410,37 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     fOutputStudy->Add(hK0MassInv);
   }
 
+  if(!fMixing) {
+    //phi distributions
+    TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
+    hPhiDistCAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistCAll);
+
+    TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
+    hPhiDistHAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistHAll);
+
+    TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
+    hPhiDistKAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistKAll);
+
+    TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
+    hPhiDistDAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistDAll);
+  }
+
   //for MC analysis only
-  if (fReadMC) {
+  for(Int_t i=0;i<fNPtBinsCorr;i++) {
 
-    for(Int_t i=0;i<fNPtBinsCorr;i++){
+    if (fReadMC && !fMixing) {
 
       //displacement histos
-      namePlot="histDispl_K_Bin"; namePlot+=i;
+      namePlot="histDispl_K0_Bin"; namePlot+=i;
       TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
       hDisplK->SetMinimum(0);
       fOutputStudy->Add(hDisplK);
   
-      namePlot="histDispl_K_HF_Bin";  namePlot+=i;
+      namePlot="histDispl_K0_HF_Bin";  namePlot+=i;
       TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
       hDisplK_HF->SetMinimum(0);
       fOutputStudy->Add(hDisplK_HF);
@@ -1064,12 +1465,12 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hDisplCharg_HF->SetMinimum(0);
       fOutputStudy->Add(hDisplCharg_HF);
 
-      namePlot="histDispl_K_From_c_Bin"; namePlot+=i;
+      namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
       TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
       hDisplK_c->SetMinimum(0);
       fOutputStudy->Add(hDisplK_c);
   
-      namePlot="histDispl_K_HF_From_c_Bin";  namePlot+=i;
+      namePlot="histDispl_K0_HF_From_c_Bin";  namePlot+=i;
       TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
       hDisplK_HF_c->SetMinimum(0);
       fOutputStudy->Add(hDisplK_HF_c);
@@ -1095,12 +1496,12 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hDisplCharg_HF_c->SetMinimum(0);
       fOutputStudy->Add(hDisplCharg_HF_c);
 
-      namePlot="histDispl_K_From_b_Bin"; namePlot+=i;
+      namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
       TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
       hDisplK_b->SetMinimum(0);
       fOutputStudy->Add(hDisplK_b);
   
-      namePlot="histDispl_K_HF_From_b_Bin";  namePlot+=i;
+      namePlot="histDispl_K0_HF_From_b_Bin";  namePlot+=i;
       TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
       hDisplK_HF_b->SetMinimum(0);
       fOutputStudy->Add(hDisplK_HF_b);
@@ -1132,11 +1533,11 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
       hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
       hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
-      hOrigin_Charm->GetXaxis()->SetBinLabel(4,"B->#");
-      hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
-      hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->D->#");
-      hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->X->#");
-      hOrigin_Charm->GetXaxis()->SetBinLabel(8,"c hadr.");
+      hOrigin_Charm->GetXaxis()->SetBinLabel(4,"c hadr.");
+      hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->#");
+      hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
+      hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->#");
+      hOrigin_Charm->GetXaxis()->SetBinLabel(8,"B->D->X->#");
       hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
       fOutputStudy->Add(hOrigin_Charm);
 
@@ -1146,28 +1547,30 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
       hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
       hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
-      hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"B->#");
-      hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
-      hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->D->#");
-      hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->X->#");
-      hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"c hadr.");
+      hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"c hadr.");
+      hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->#");
+      hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
+      hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->#");
+      hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"B->D->X->#");
       hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
       fOutputStudy->Add(hOrigin_Kcharg);
 
-      namePlot="histOrig_K_Bin";  namePlot+=i;
+      namePlot="histOrig_K0_Bin";  namePlot+=i;
       TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
       hOrigin_K->SetMinimum(0);
       hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
       hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
       hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
-      hOrigin_K->GetXaxis()->SetBinLabel(4,"B->#");
-      hOrigin_K->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
-      hOrigin_K->GetXaxis()->SetBinLabel(6,"B->D->#");
-      hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->X->#");
-      hOrigin_K->GetXaxis()->SetBinLabel(8,"c hadr.");
+      hOrigin_K->GetXaxis()->SetBinLabel(4,"c hadr.");
+      hOrigin_K->GetXaxis()->SetBinLabel(5,"B->#");
+      hOrigin_K->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
+      hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->#");
+      hOrigin_K->GetXaxis()->SetBinLabel(8,"B->D->X->#");
       hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
       fOutputStudy->Add(hOrigin_K);
+    }
 
+    if (fReadMC) {
       //origin of D0 histos
       namePlot="histOrig_D0_Bin";  namePlot+=i;
       TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
@@ -1175,332 +1578,604 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
       fOutputStudy->Add(hOrigin_D0);
+
+      //primary tracks (Kine & Reco)
+      namePlot="hPhysPrim_Bin";  namePlot+=i;
+      TH1F *hPhysPrim = new TH1F(namePlot.Data(), "Origin of hadrons",2,0.,2.);
+      hPhysPrim->SetMinimum(0);
+      hPhysPrim->GetXaxis()->SetBinLabel(1,"OK");
+      hPhysPrim->GetXaxis()->SetBinLabel(2,"NO");
+      fOutputStudy->Add(hPhysPrim);
     }
   }
-
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODEvent* aod, AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
+void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
 //
 // Method for correlations D0-hadrons study
 //
-
-  Double_t mD0, mD0bar, deltaphi, d0, d0err;
+  Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
+  Double_t mD0, mD0bar;
+  Int_t origD0 = 0, PDGD0 = 0, ptbin = 0;
   d->InvMassD0(mD0,mD0bar);
-  Int_t ptbin = PtBinCorr(d->Pt());
+  Double_t mInv[2] = {mD0, mD0bar};
+  ptbin = PtBinCorr(d->Pt());
+
   if(ptbin < 0) return;
-  AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
 
-  Int_t N_Charg = 0, N_Kcharg = 0, N_Kaons = 0;
+  //Fill of D0 phi distribution
+  if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());  
 
-  if(fReadMC == 0) {
-    Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
-    Double_t highPt = 0; Double_t lead[2] = {0,0};  //infos for leading particle (pt,deltaphi)
+  //Origin of D0
+  TString orig="";
+  if(fReadMC) {
+    origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
+    PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
+    switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
+      case 4:
+        orig = "_From_c";
+        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
+        break;
+      case 5:
+        orig = "_From_b";
+        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
+        break;
+      default:
+        return;
+    }
+  }
 
-    for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
-      AliAODTrack * track = aod->GetTrack(itrack);
-     
-      if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
+  Double_t highPt = 0; Double_t lead[4] = {0,0,0,1};  //infos for leading particle (pt,deltaphi)
+
+  //loop over the tracks in the pool 
+  Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+  Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+  Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+               
+  Int_t NofEventsinPool = 1;
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorTr->GetNofEventsInPool(); 
+    if(!execPoolTr) {
+      AliInfo("Mixed event analysis: track pool is not ready");
+      NofEventsinPool = 0;
+    }
+  }
+
+  //Charged tracks
+  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
+    Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
+    if(!analyzetracksTr) {
+      AliInfo("AliHFCorrelator::Cannot process the track array");
+      continue;
+    }
+       
+    for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
+
+      Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
+
+      if(!fMixing) {      
+       if(fSoftPiCut && !track->CheckSoftPi()) { //removal of soft pions
+          if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
+          if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
+          continue;
+        } else { //not a soft pion
+          if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
+          if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
+        }
+        Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
+        if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
+      }
+      if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+
+      if(fReadMC) {
+        AliAODMCParticle* trkKine = (AliAODMCParticle*)mcArray->At(track->GetLabel());
+        if (!trkKine) continue;
+        if (!trkKine->IsPhysicalPrimary()) {
+         ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);  
+         continue; //reject the Reco track if correspondent Kine track is not primary
+        } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
+      }
 
-      GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
+      Double_t effTr = track->GetWeight(); //extract track efficiency
+      Double_t effD0 = 1.;
+      if(fReadMC) {
+        if(origD0==4) effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
+        if(origD0==5) effD0 = fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv);
+      } else effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
+      Double_t eff = effTr*effD0;
+
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks
+
+      if(!fMixing) N_Charg++;
+
+      //retrieving leading info...
+      if(track->Pt() > highPt) {
+        if(fReadMC && track->GetLabel()<1) continue;
+        if(fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) continue;
+        lead[0] = fCorrelatorTr->GetDeltaPhi();
+        lead[1] = fCorrelatorTr->GetDeltaEta();
+        lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
+        if(fReadMC) {
+         if(origD0==4) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv)); //weight is 1./efficiency
+         if(origD0==5) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv)); //weight is 1./efficiency
+       } else lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv));
+        highPt = track->Pt();
+      }
 
-      ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
-      deltaphi = d->Phi()-track->Phi();
-      if (deltaphi < -1.571) deltaphi+=6.283;
-      if (deltaphi > 4.712) deltaphi-=6.283;
-      Double_t pttrack = track->Pt();
+    } // end of tracks loop
+  } //end of event loop for fCorrelatorTr
 
-      if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi;}  //for leading particle
+ if(fKaonCorr) { //loops for Kcharg and K0
 
-      //Lines needed to include overflow into THnSparse projections!
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
-      if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
-      if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorKc->GetNofEventsInPool(); 
+    if(!execPoolKc) {
+      AliInfo("Mixed event analysis: K+/- pool is not ready");
+      NofEventsinPool = 0;
+    }
+  }
 
-      //variables for filling histos
-      Double_t fillSpPhiD0[4] = {deltaphi,mD0,pttrack,d0/d0err};
-      Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,pttrack,d0/d0err};
-   
-      //generic charged tracks (NO PID selection)
-      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
-      }
-      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
-      }
-      if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
-      N_Charg++;
-        //pT_weighted track correlations fill
-        if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
-        }
-        if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
-        }       
-
-      //hadron identification
-      if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) {
-        if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
-        }
-        if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
+  //Charged Kaons loop
+  for (Int_t jMix = 0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
+    Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
+    if(!analyzetracksKc) {
+      AliInfo("AliHFCorrelator::Cannot process the K+/- array");
+      continue;
+    }  
+
+    for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
+
+      Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
+
+      if(!fMixing) {      
+       if(fSoftPiCut && !kCharg->CheckSoftPi()) { //removal of soft pions
+          if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
+          if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
+          continue;
+        } else {
+          if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
+          if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
         }
-        if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
-        N_Kcharg++;
-      } //end hadron case
+        Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
+        if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
+      }
+      if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
 
-    } //end tracks loop
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
 
-    //kaon identification
-    TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
-    Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
-    for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
-      for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
+      if(!fMixing) N_KCharg++;
 
-      AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
-      if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
-      if(SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
+    } // end of charged kaons loop
+  } //end of event loop for fCorrelatorKc
 
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t ptV0=v0->Pt(), deltaphiV0=d->Phi()-v0->Phi();
-      if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
-      deltaphiV0 = d->Phi()-v0->Phi();
-      if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
-      if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorK0->GetNofEventsInPool(); 
+    if(!execPoolK0) {
+      AliInfo("Mixed event analysis: K0 pool is not ready");
+      NofEventsinPool = 0;
+    }
+  }
 
-      Double_t fillSpPhiD0K0[4] = {deltaphiV0,mD0,ptV0,0.};
-      Double_t fillSpPhiD0barK0[4] = {deltaphiV0,mD0bar,ptV0,0.};
+  //K0 loop
+  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
+    Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
+    if(!analyzetracksK0) {
+      AliInfo("AliHFCorrelator::Cannot process the K0 array");
+      continue;
+    }  
 
-      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
+    for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
+
+      Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
+
+      if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+  
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
+
+      if(!fMixing) N_Kaons++;
+
+    } // end of charged kaons loop
+  } //end of event loop for fCorrelatorK0
+
+ } //end of 'if(fKaonCorr)'
+
+  Double_t fillSpLeadD0[4] = {lead[0],mD0,lead[1],0.4}; //dummy value for threshold of leading!
+  Double_t fillSpLeadD0bar[4] = {lead[0],mD0bar,lead[1],0.4};
+
+  //leading track correlations fill
+  if(!fMixing) {
+    if(fReadMC) {
+      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
+        if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);  
+        if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);  
+        if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);  //non HF  
       }
-      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
+      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
+        if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  
+        if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); 
+        if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  //non HF  
       }
-      if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K_Bin%d",ptbin)))->Fill(v0->Pt());
-      N_Kaons++;
-    } // end kaon case
+    } else {
+        if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); 
+        if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
+    }
 
-    //Leading track correlations fill
-    if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-     ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0);
+    //Fill of count histograms
+    if (!fAlreadyFilled) { 
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
     }
-    if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-     ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
+  }
+
+  fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCParticle* d, TClonesArray* mcArray) {
+//
+// Method for correlations D0-hadrons study
+//
+  Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
+  Double_t mD0 = 1.864, mD0bar = 1.864;
+  Double_t mInv[2] = {mD0, mD0bar};
+  Int_t origD0 = 0, PDGD0 = 0;
+  Int_t ptbin = PtBinCorr(d->Pt());
+
+  if(ptbin < 0) return;
+
+  //Fill of D0 phi distribution
+  if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi()); 
+  
+  //Origin of D0
+  TString orig="";
+  origD0=CheckD0Origin(mcArray,d);
+  PDGD0 = d->GetPdgCode();
+  switch (CheckD0Origin(mcArray,d)) {
+    case 4:
+      orig = "_From_c";
+      ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
+      break;
+    case 5:
+      orig = "_From_b";
+      ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
+      break;
+    default:
+      return;
+  }
+
+  Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
+
+  //loop over the tracks in the pool 
+  Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+  Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+  Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+               
+  Int_t NofEventsinPool = 1;
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorTr->GetNofEventsInPool(); 
+    if(!execPoolTr) {
+      AliInfo("Mixed event analysis: track pool is not ready");
+      NofEventsinPool = 0;
     }
+  }
 
-  } //end fReadMC = 0
+  //Charged tracks
+  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
 
-  if(fReadMC == 1) {
-    Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
-    Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi,origtrack)
+    Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
+    if(!analyzetracksTr) {
+      AliInfo("AliHFCorrelator::Cannot process the track array");
+      continue;
+    }
+       
+    for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
 
-    //Origin of D0
-    TString orig=""; Int_t origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
-    switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0)))
-    {
-      case 4:
-        orig = "_From_c";
-        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
-        break;
-      case 5:
-        orig = "_From_b";
-        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
-        break;
-      default:
-        return;
+      Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
+      if(track->GetLabel()<0) continue;
+      if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+      if(track->Pt() < 0.3 || TMath::Abs(track->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
+      if(!fMixing) N_Charg++;
+
+      AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
+      if(!trkMC) continue;
+
+      if (!trkMC->IsPhysicalPrimary()) {  //reject material budget, or other fake tracks
+       ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);  
+       continue;
+      } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
+
+      if (IsDDaughter(d,trkMC,mcArray)) continue;
+      if (fSoftPiCut && IsSoftPion_MCKine(d,trkMC,mcArray)) continue; //remove soft pions (if requestes, e.g. for templates)
+
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
+
+      //retrieving leading info...
+      if(track->Pt() > highPt) {
+        lead[0] = fCorrelatorTr->GetDeltaPhi();
+        lead[1] = fCorrelatorTr->GetDeltaEta();
+        lead[2] = fReadMC ? CheckTrackOrigin(mcArray,trkMC) : 0;
+        highPt = track->Pt();
+      }
+
+    } // end of tracks loop
+  } //end of event loop for fCorrelatorTr
+
+ if(fKaonCorr) { //loops for Kcharg and K0
+
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorKc->GetNofEventsInPool(); 
+    if(!execPoolKc) {
+      AliInfo("Mixed event analysis: K+/- pool is not ready");
+      NofEventsinPool = 0;
     }
+  }
 
-    for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
-      AliAODTrack * track = aod->GetTrack(itrack);
+  //Charged Kaons loop
+  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
+    Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
+    if(!analyzetracksKc) {
+      AliInfo("AliHFCorrelator::Cannot process the K+/- array");
+      continue;
+    }  
 
-      if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
+    for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
 
-      Int_t label = track->GetLabel();
-      if (label<0) continue; //discard negative label tracks
-      Int_t PDGtrack = ((AliAODMCParticle*)mcArray->At(label))->GetPdgCode();
+      Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
+      if(kCharg->GetLabel()<1) continue;
+      if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+      if(TMath::Abs(kCharg->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
+      if(!fMixing) N_KCharg++;
+
+      AliAODMCParticle *kChargMC = (AliAODMCParticle*)mcArray->At(kCharg->GetLabel());
+      if(!kChargMC) continue;
+
+      if (IsDDaughter(d,kChargMC,mcArray)) continue;
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
+
+    } // end of charged kaons loop
+  } //end of event loop for fCorrelatorKc
+
+  if(fMixing) {
+    NofEventsinPool = fCorrelatorK0->GetNofEventsInPool(); 
+    if(!execPoolK0) {
+      AliInfo("Mixed event analysis: K0 pool is not ready");
+      NofEventsinPool = 0;
+    }
+  }
 
-      GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
+  //K0 loop
+  for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
+    Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
+    if(!analyzetracksK0) {
+      AliInfo("AliHFCorrelator::Cannot process the K0 array");
+      continue;
+    }  
 
-      //Infos on track (origin, phi, eta)
-      Int_t origTr = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(label));
-      deltaphi = d->Phi()-track->Phi();
-      if (deltaphi < -1.571) deltaphi+=6.283;
-      if (deltaphi > 4.712) deltaphi-=6.283;
-      Double_t pttrack = track->Pt();
+    for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
 
-      if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi; lead[2] = origTr;}  //for leading particle
+      Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
+      if(!runcorrelation) continue;
+      
+      AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
+      if(k0->GetLabel()<1) continue;
+      if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+      if(TMath::Abs(k0->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
+  
+      AliAODMCParticle *k0MC = (AliAODMCParticle*)mcArray->At(k0->GetLabel());
+      if(!k0MC) continue;
 
-      //Lines needed to include overflow into THnSparse projections!
-      Double_t d0orig = d0;
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
-      if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
-      if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
+      if (IsDDaughter(d,k0MC,mcArray)) continue;
+      FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
 
-      //variables for filling histos
-      Double_t fillSpPhiD0[4] = {deltaphi,mD0,pttrack,d0/d0err};
-      Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,pttrack,d0/d0err};
+      if(!fMixing) N_Kaons++;
 
-      //generic charged tracks case
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
-         if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//c tr
-         if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);//b tr
-      }
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
-         if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
-         if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar); 
-      }
-      if(!fAlreadyFilled) {
-       ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
-        if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF_Bin%d",ptbin)))->Fill(d0orig);
-        if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
-       ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Charg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
-        ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Charg_Bin%d",ptbin)))->Fill(track->Pt());
-        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Charg_Bin%d",ptbin)))->Fill(origTr);
-      }
-      N_Charg++;
-        //pT_weighted track correlations fill
-        if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pttrack);//b tr
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pow(pttrack,2.));//b tr
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,PtWeig(ptbin,pttrack));//b tr
-        }
-        if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
-           ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pttrack);//b tr
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigSqr_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pow(pttrack,2.));//b tr
-           if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));//c tr
-           if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_WeigDist_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,PtWeig(ptbin,pttrack));//b tr
-        } 
-
-      //hadron identification (K/pi/p from MCTruth)
-      if(TMath::Abs(PDGtrack) == 321) {  
-        if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
-           if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
-           if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0);
-        }
-        if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
-           ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
-           if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
-           if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0bar); 
-        }
-        if(!fAlreadyFilled) {
-         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_Bin%d",ptbin)))->Fill(d0orig); //Fills displacement histos
-         if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF_Bin%d",ptbin)))->Fill(d0orig);
-          if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig);
-         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_Kcharg%s_Bin%d",orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
-          ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_Kcharg_Bin%d",ptbin)))->Fill(track->Pt());
-          ((TH1F*)fOutputStudy->FindObject(Form("histOrig_Kcharg_Bin%d",ptbin)))->Fill(origTr);
-       }
-        N_Kcharg++;
+    } // end of charged kaons loop
+  } //end of event loop for fCorrelatorK0
 
-      } //end hadron case
+ } //end of 'if(fKaonCorr)'
 
-    } //end tracks loop
+  Double_t fillSpLeadMC[4] = {lead[0],mD0,lead[1],0.4}; //mD0 = mD0bar = 1.864
 
-    //Kaon identification
-    TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
-    Int_t idArrayV0[v0array->GetEntriesFast()][2]; //array for skipping K0 double-counting
-    for(int iV0=0; iV0<v0array->GetEntriesFast();iV0++) {
-      for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
-      AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
+  //leading track correlations fill
+  if(!fMixing) {
+    if(d->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
+      ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
+      ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
+      if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
+      if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
+      if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
+    }
+    if(d->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
+      ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
+      ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
+      if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);  
+      if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); 
+      if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
+    }
+
+    //Fill of count histograms
+    if (!fAlreadyFilled) { 
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
+      ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
+    }
+
+    fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
+  }
 
-      if(v0->Pt() < fPtThreshLow.at(ptbin) || v0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
-      if(SelectV0(v0,vtx,1,idArrayV0)) continue; //option 1 = for correlations
-      Int_t pdgCodes[2] = {211,211};
-      Int_t labV0 = v0->MatchToMC(311,mcArray,2,pdgCodes);
-      if(labV0<=0) continue;
-
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t ptV0=v0->Pt(), deltaphiV0=d->Phi()-v0->Phi();
-      deltaphiV0 = d->Phi()-v0->Phi();
-      if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
-      if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
-      if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
-
-      Double_t fillSpPhiD0K0[4] = {deltaphiV0,mD0,ptV0,0.};
-      Double_t fillSpPhiD0barK0[4] = {deltaphiV0,mD0bar,ptV0,0.};
-
-      Int_t origV0 = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(labV0));
-
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
-         if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);  
-         if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);  
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Double_t mInv[], Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t ptbin, Int_t type, Double_t wg) {
+  //
+  //fills the THnSparse for correlations, calculating the variables
+  //
+
+  //Initialization of variables
+  Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
+  mD0 = mInv[0];
+  mD0bar = mInv[1];
+
+  if (fReadMC && track->GetLabel()<1) return;
+  if (fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) return;
+  Double_t ptTrack = track->Pt();
+  Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
+  Double_t phiTr = track->Phi();
+  Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
+
+  TString part = "", orig = "";
+
+  switch (type) {
+    case(kTrack): {
+      part = "Charg";
+      deltaphi = fCorrelatorTr->GetDeltaPhi();
+      deltaeta = fCorrelatorTr->GetDeltaEta();
+      break;
+    }
+    case(kKCharg): {
+      part = "Kcharg";
+      deltaphi = fCorrelatorKc->GetDeltaPhi();
+      deltaeta = fCorrelatorKc->GetDeltaEta();
+      break;
+    }
+    case(kK0): {
+      part = "K0";
+      deltaphi = fCorrelatorK0->GetDeltaPhi();
+      deltaeta = fCorrelatorK0->GetDeltaEta();
+      break;
+    }
+  }
+  
+  if(fMixing == kSE) {
+
+    //Fixes limits; needed to include overflow into THnSparse projections!
+    Double_t pTorig = track->Pt();
+    Double_t d0orig = track->GetImpPar();
+    Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
+    Double_t displLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
+    Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
+    if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
+    if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
+    if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
+    if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
+  
+    //variables for filling histos
+    Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
+    Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
+    Double_t fillSpWeigD0[4] = {deltaphi,mD0,deltaeta,ptTrack};
+    Double_t fillSpWeigD0bar[4] = {deltaphi,mD0bar,deltaeta,ptTrack};
+
+    if(fReadMC == 0) {
+      //sparse fill for data (tracks, K+-, K0) + weighted
+      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
       }
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
-         if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);  
-         if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0); 
+      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+        ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
       }
       if(!fAlreadyFilled) {
-         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_Bin%d",ptbin)))->Fill(0.); //Fills displacement histos
-        if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_HF_Bin%d",ptbin)))->Fill(0.);
-         if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(0.);
-        ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K%s_Bin%d",orig.Data(),ptbin)))->Fill(0.); //Fills displacement histos
-         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K_Bin%d",ptbin)))->Fill(v0->Pt());
-         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_K_Bin%d",ptbin)))->Fill(origV0);
+       ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
+       ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
       }
-      N_Kaons++;
-    } // end kaon case
-
-    //leading track correlations fill
-    if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
-      ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0); //c and b D0
-      ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0); //c or b D0
-      if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0);  
-      if(origD0==5&&(int)lead[2]>=3&&(int)lead[2]<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0);  
-    }
-    if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-      ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
-      ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar); //c or b D0
-      if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar);  
-      if(origD0==5&&(int)lead[2]>=3&&(int)lead[2]<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[1],mD0bar); 
     }
 
-  } //end fReadMC = 1
+    if(fReadMC) {
+
+      if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
+
+      //sparse fill for data (tracks, K+-, K0) + weighted
+      if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0 (from MCTruth)
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+         if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+         if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+         if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
+         if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
+         if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
+         if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
+      }
+      if(PdgD0==-421 && fIsSelectedCandidate>1) { //D0bar
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+         if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+         if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg); 
+         if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
+         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
+         if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
+         if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
+         if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
+      } 
+      if(!fAlreadyFilled) {
+       ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
+        if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
+        if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
+       ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
+        ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
+        ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
+       ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
+      }
+    }//end MC case
 
-  if (!fAlreadyFilled) {
-    ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Charg_Bin%d",ptbin)))->Fill(N_Charg);
-    ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_Kcharg_Bin%d",ptbin)))->Fill(N_Kcharg);
-    ((TH1F*)fOutputStudy->FindObject(Form("hist_Count_K_Bin%d",ptbin)))->Fill(N_Kaons);
-  }
+  } //end of SE fill
+
+  if(fMixing == kME) {
+
+    //Fixes limits; needed to include overflow into THnSparse projections!
+    Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
+    if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
+    if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
+    Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(3)->GetXmax(); //all plots have same axes...
+    if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
+
+    //variables for filling histos
+    Double_t fillSpPhiD0[4] = {deltaphi,mD0,deltaeta,0.4}; //dummy for ME threshold! unless explicitly set by flag...
+    Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,deltaeta,0.4};
+    if(fMEAxisThresh) {
+      fillSpPhiD0[3] = ptTrack;
+      fillSpPhiD0bar[3] = ptTrack;
+    }
 
-  fAlreadyFilled=kTRUE; //distribution plots for tracks filled  
+    if(fReadMC == 0) {
+      //sparse fill for data (tracks, K+-, K0)
+      if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+      if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+    }
+    if(fReadMC == 1) {
+      //sparse fill for data (tracks, K+-, K0)
+      if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3))  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
+      if(PdgD0==-421 && fIsSelectedCandidate>1) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
+    }//end MC case
 
+  } //end of ME fill
+  
+  return;
 }
 
 //_________________________________________________________________________________________________
@@ -1510,14 +2185,14 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
   // 0) Not HF
   // 1) D->#
   // 2) D->X->#
-  // 3) B->#
-  // 4) B->X-># (X!=D)
-  // 5) B->D->#
-  // 6) B->D->X->#
-  // 7) c hadronization
+  // 3) c hadronization
+  // 4) B->#
+  // 5) B->X-># (X!=D)
+  // 6) B->D->#
+  // 7) B->D->X->#
   // 8) b hadronization
   //
-  printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
+  if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
        
   Int_t pdgGranma = 0;
   Int_t mother = 0;
@@ -1531,11 +2206,12 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
   Bool_t isBchaindaugh=kFALSE;
   Bool_t isQuarkFound=kFALSE;
 
-  while (mother > 0){
+  if (mother<0) return -1;
+  while (mother >= 0){
     istep++;
-    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
-    if (mcGranma){
-      pdgGranma = mcGranma->GetPdgCode();
+    AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
+    if (mcMoth){
+      pdgGranma = mcMoth->GetPdgCode();
       abspdgGranma = TMath::Abs(pdgGranma);
       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
        isBchaindaugh=kTRUE;
@@ -1546,10 +2222,10 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
         if(istep==1) isDdaugh=kTRUE;
       }
       if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
-      mother = mcGranma->GetMother();
+      mother = mcMoth->GetMother();
     }else{
       AliError("Failed casting the mother particle!");
-      break;
+      return -1;
     }
   }
 
@@ -1558,23 +2234,49 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
     if(!isFromB) {  //charm
       if(isDdaugh) return 1; //charm immediate
       else if(isDchaindaugh) return 2; //charm chain
-      else return 7; //charm hadronization
+      else return 3; //charm hadronization
     }
     else { //beauty
-      if(isBdaugh) return 3; //b immediate
+      if(isBdaugh) return 4; //b immediate
       else if(isBchaindaugh) { //b chain
         if(isDchaindaugh) {
-          if(isDdaugh) return 5; //d immediate
-          return 6; //d chain
+          if(isDdaugh) return 6; //d immediate
+          return 7; //d chain
           }
-        else return 4; //b, not d
+        else return 5; //b, not d
       }
       else return 8; //b hadronization
     }
   }
-  else return 0; //no HF quark
+  else if(!isDdaugh && !isDchaindaugh && !isBdaugh && !isBchaindaugh) return 0; //no HF decay 
+     //in this case, it's !isQuarkFound, but not in 100% cases it's a non HF particle!
+     //rarely you can find a D/B meson which comes from a -1! It isn't a Non-HF, in that case! And I'll return -1...
+
+  return -1; //some problem spotted
 }
+//________________________________________________________________________
+Bool_t AliAnalysisTaskSED0Correlations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
+
+  //Daughter removal in MCKine case
+  Bool_t isDaughter = kFALSE;
+  Int_t labelD0 = d->GetLabel();
+
+  Int_t mother = track->GetMother();
+
+  //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
+  while (mother > 0){
+    AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother)); //it's the mother of the track!
+    if (mcMoth){
+      if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
+      mother = mcMoth->GetMother(); //goes back by one
+    } else{
+      AliError("Failed casting the mother particle!");
+      break;
+    }
+  }
 
+  return isDaughter;
+}
 
 //________________________________________________________________________
 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
@@ -1590,53 +2292,6 @@ Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
   return ptbin;
 }
 
-//________________________________________________________________________
-Double_t AliAnalysisTaskSED0Correlations::PtWeig(Int_t ptbin, Double_t x) const {
-  //
-  //gives the weight for Weighted Corrs (inverse of pT distribution, from 3 D0 pt bins)
-  //
-  if(x>11.5) x=11.5; if(x<0.3) x=0.3; //bounds for fit of distributions!
-  if(ptbin >= 3 && ptbin <= 5)  return 1/(0.22958*(1.507e+03-5.035e+02*x+5.681e+01*x*x-2.186e+00*x*x*x+exp(1.336e+01-2.146e+00*x+1.206e-01*x*x)));
-  if(ptbin >= 6 && ptbin <= 8)  return 1/(1.71828*(1.984e+02-6.113e+01*x+6.444e+00*x*x-2.342e-01*x*x*x+exp(1.023e+01-2.059e+00*x+1.194e-01*x*x)));
-  if(ptbin >= 9 && ptbin <= 10) return 1/(1.25905*(1.990e+02-6.306e+01*x+6.995e+00*x*x-2.681e-01*x*x*x+exp(9.125e+00-2.053e+00*x+1.276e-01*x*x)));
-
-  return 0; //for other bins plot is disabled!
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskSED0Correlations::GetImpParameter(AliAODTrack *track, AliAODEvent* aod, Double_t &d0, Double_t &d0err) const {
-  //
-  //calculates d0 and error on d0 for the track
-  //      
-  Double_t  d0z0[2],covd0z0[3];
-  AliESDtrack esdTrack(track);  // AliESDTrack conversion of AOD track
-  esdTrack.PropagateToDCA(aod->GetPrimaryVertex(),aod->GetMagneticField(),10000.,d0z0,covd0z0); 
-  d0 = TMath::Abs(d0z0[0]); // impact parameter
-  d0err = TMath::Sqrt(covd0z0[0]); // resolution of impact parameter
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskSED0Correlations::TrackSelectedInLoop(AliAODTrack* track, AliAODRecoDecayHF2Prong *d, AliAODEvent *aod, Int_t ptbin, Int_t idDaughs[]) const {
-  //
-  //rejection of tracks in loop for correlations
-  //      
-  Bool_t output = kTRUE;
-  AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
-  Double_t bz = aod->GetMagneticField();
-  if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1] || track->GetID() < 0) output = kFALSE; //discards daughters of candidate
-  if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) output = kFALSE; //discard tracks outside pt range for hadrons/K
-  if(!fCutsTracks->IsHadronSelected(track,vtx,bz)) { //track discarded by quality cuts
-    ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(1);
-    output = kFALSE;
-  } else ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
-  if(!fCutsTracks->InvMassDstarRejection(d,track,fIsSelectedCandidate)) {
-    ((TH1F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1);
-    output = kFALSE; 
-  } else ((TH1F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0); //rejects possible pions from D* using inv.mass
-
-  return output;
-}
-
 //---------------------------------------------------------------------------
 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
 {
@@ -1651,9 +2306,7 @@ Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx
   AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
 
   if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
-    if(v0->Pt() < 3. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0040) return kFALSE;
-    if(v0->Pt() > 3. && v0->Pt() < 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0052) return kFALSE;
-    if(v0->Pt() > 6. && TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.0075) return kFALSE;
+    if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
   }
 
   //This part removes double counting for swapped tracks!
@@ -1669,6 +2322,40 @@ Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx
   return kTRUE;
 }
 
+//---------------------------------------------------------------------------
+Bool_t AliAnalysisTaskSED0Correlations::IsSoftPion_MCKine(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* arrayMC) const
+{
+  //
+  // Removes soft pions in Kine
+
+  //Daughter removal in MCKine case
+  Bool_t isSoftPi = kFALSE;
+  Int_t labelD0 = d->GetLabel();
+
+  Int_t mother = track->GetMother();
+  if(mother<0) return isSoftPi; //safety check
+
+  AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); //it's the mother of the track!
+  if(!mcMoth){
+    return isSoftPi;
+  }
+  if(TMath::Abs(mcMoth->GetPdgCode())==413 && mcMoth->GetNDaughters()==2) { //mother is D* with 2 daughs
+    Int_t labdau1 = mcMoth->GetDaughter(0);
+    Int_t labdau2 = mcMoth->GetDaughter(1);
+    AliAODMCParticle* dau1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau1));
+    AliAODMCParticle* dau2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau2));
+    if(!dau1 || !dau2) return isSoftPi; //safety check
+    if(dau1->GetLabel()==labelD0 || dau2->GetLabel()==labelD0) { //one of the daughs is the D0 trigger
+      if((TMath::Abs(dau1->GetPdgCode())==421 && TMath::Abs(dau2->GetPdgCode())==211)||(TMath::Abs(dau1->GetPdgCode())==211 && TMath::Abs(dau2->GetPdgCode())==421)) {
+       isSoftPi = kTRUE; //ok, soft pion was found
+       return isSoftPi;
+      }
+    }
+  } 
+
+  return isSoftPi;
+}
+
 //________________________________________________________________________
 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
 
@@ -1688,7 +2375,17 @@ void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
     cout << fPtThreshUp.at(i) << "; ";
   }
   cout << "\n--------------------------\n";
-  cout << "MC Truth = "<<fReadMC<<"\n";
+  cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
+  cout << "--------------------------\n";
+  cout << "MC Truth = "<<fReadMC<<" - MC Reco: Trk = "<<fRecoTr<<", D0 = "<<fRecoD0<<"\n";
+  cout << "--------------------------\n";
+  cout << "Sel of Event tpye (PP,GS,FE,...)= "<<fSelEvType<<"\n";
+  cout << "--------------------------\n";
+  cout << "Ev Mixing = "<<fMixing<<"\n";
+  cout << "--------------------------\n";
+  cout << "ME thresh axis = "<<fMEAxisThresh<<"\n";
+  cout << "--------------------------\n";
+  cout << "Soft Pi Cut = "<<fSoftPiCut<<"\n";
   cout << "--------------------------\n";
 }