New class for steering D-hadron correlation analysis (Sandro)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 10 Jul 2012 08:49:05 +0000 (08:49 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 10 Jul 2012 08:49:05 +0000 (08:49 +0000)
14 files changed:
PWGHF/CMakelibPWGHFcorrelationHF.pkg
PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx
PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h
PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx
PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.h
PWGHF/correlationHF/AliHFAssociatedTrackCuts.cxx
PWGHF/correlationHF/AliHFAssociatedTrackCuts.h
PWGHF/correlationHF/AliHFCorrelator.cxx [new file with mode: 0644]
PWGHF/correlationHF/AliHFCorrelator.h [new file with mode: 0644]
PWGHF/correlationHF/AliReducedParticle.cxx
PWGHF/correlationHF/AliReducedParticle.h
PWGHF/correlationHF/macros/AddTaskD0Correlations.C
PWGHF/correlationHF/macros/AddTaskDStarCorrelations.C
PWGHF/correlationHF/macros/makeTFileAssociatedTrackCuts.C

index 765c160..47d4c1c 100644 (file)
@@ -34,6 +34,7 @@ set ( CLASS_HDRS
     AliAnalysisTaskDxHFEParticleSelection.h
     AliAnalysisTaskDxHFECorrelation.h
     AliHFAssociatedTrackCuts.h
+    AliHFCorrelator.h
     AliReducedParticle.h
     AliAnalysisTaskDStarCorrelations.h
     AliAnalysisTaskSED0Correlations.h
index e9bc77e..a9886a2 100644 (file)
@@ -27,7 +27,7 @@
 
 /* $Id$ */
 
-#include <TDatabasePDG.h>
+//#include <TDatabasePDG.h>
 #include <TParticle.h>
 #include <TVector3.h>
 #include <TChain.h>
@@ -40,7 +40,7 @@
 #include "AliAODRecoCascadeHF.h"
 #include "AliAODRecoDecayHF2Prong.h"
 #include "AliAODPidHF.h"
-#include "AliEventPoolManager.h"
+//#include "AliEventPoolManager.h"
 #include "AliVParticle.h"
 #include "AliAnalysisManager.h"
 #include "AliAODInputHandler.h"
 #include "AliAODMCParticle.h"
 #include "AliNormalizationCounter.h"
 #include "AliReducedParticle.h"
+#include "AliHFCorrelator.h"
 
 
 
-using std::cout;
-using std::endl;
 
 ClassImp(AliAnalysisTaskDStarCorrelations)
 
@@ -62,14 +61,17 @@ ClassImp(AliAnalysisTaskDStarCorrelations)
 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
 AliAnalysisTaskSE(),
 fhandler(0x0),
-fPoolMgr(0x0),      
+//fPoolMgr(0x0),      
 fmcArray(0x0),
 fCounter(0x0),
+fCorrelator(0x0),
 fselect(0),
 fmontecarlo(kFALSE),
 fmixing(kFALSE),
+fSystem(kFALSE),
 fEvents(0),
 fDebug(0),
+fDisplacement(0),
 
 fOutput(0x0),            
 fCuts(0),
@@ -83,14 +85,17 @@ AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t*
 AliAnalysisTaskSE(name),
 
 fhandler(0x0),
-fPoolMgr(0x0),   
+//fPoolMgr(0x0),   
 fmcArray(0x0),
 fCounter(0x0),
+fCorrelator(0x0),
 fselect(0),
 fmontecarlo(kFALSE),
 fmixing(kFALSE),
+fSystem(kFALSE),
 fEvents(0),
 fDebug(0),
+fDisplacement(0),
 
 fOutput(0x0),                
 fCuts(0),
@@ -115,9 +120,10 @@ AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
        Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  
        
        if(fhandler) {delete fhandler; fhandler = 0;}
-       if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
+       //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
        if(fmcArray) {delete fmcArray; fmcArray = 0;}
        if(fCounter) {delete fCounter; fCounter = 0;}
+       if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
        if(fOutput) {delete fOutput; fOutput = 0;}
        if(fCuts) {delete fCuts; fCuts = 0;}
        if(fCuts2) {delete fCuts2; fCuts2=0;}
@@ -132,6 +138,10 @@ void AliAnalysisTaskDStarCorrelations::Init(){
        if(fDebug > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
        
        AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
+       
+       
+       
+       
        // Post the D* cuts
        PostData(2,copyfCuts);
        
@@ -153,39 +163,43 @@ void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
        
        // define histograms
        DefineHistoForAnalysis();
-       
        fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(3)->GetContainer()->GetName()));
        fCounter->Init();
        
-   
-       // definition of the Pool Manager for Event Mixing
+    Double_t Pi = TMath::Pi();
+       fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fSystem); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb
+       fCorrelator->SetDeltaPhiInterval((-0.5-1./32)*Pi,(1.5-1./32)*Pi); // set correct phi interval
+       //fCorrelator->SetDeltaPhiInterval(-1.57,4.71);
+       fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
+       fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
+       fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
+       fCorrelator->SetUseMC(fmontecarlo);
+       Bool_t pooldef = fCorrelator->DefineEventPool();
        
-       Int_t MaxNofEvents = 200;
-       Int_t MinNofTracks = 1000;
-       
-       Int_t NofMultiplicityBins = 5;
-       Double_t MBins[]={0,20,40,60,80,500};
-       Double_t * MultiplicityBins = MBins;
-       
-       Int_t NofZVrtxBins = 5;
-       Double_t ZBins[]={-10,-5,-2.5,2.5,5,10};
-       Double_t *ZVrtxBins = ZBins;
-       
-       
-       fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofMultiplicityBins, MultiplicityBins, NofZVrtxBins, ZVrtxBins);
+       if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
+
        
        
 }
 //_________________________________________________
 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
 
-       cout << " " << endl;
-       cout << "=================================================================================" << endl;
+       cout << "USER EXEC CHECKPOINT 1" << endl;
        
-       if(fselect==1) cout << "TASK::Correlation with hadrons "<< endl;
-       if(fselect==2) cout << "TASK::Correlation with kaons "<< endl;
-       if(fselect==3) cout << "TASK::Correlation with kzeros "<< endl;
+       //Double_t pi = TMath::Pi();
        
+       cout << " " << endl;
+       cout << "=================================================================================" << endl;
+       if(!fmixing){
+       if(fselect==1) cout << "TASK::Correlation with hadrons on SE "<< endl;
+       if(fselect==2) cout << "TASK::Correlation with kaons on SE "<< endl;
+       if(fselect==3) cout << "TASK::Correlation with kzeros on SE "<< endl;
+       }
+       if(fmixing){
+               if(fselect==1) cout << "TASK::Correlation with hadrons on ME "<< endl;
+               if(fselect==2) cout << "TASK::Correlation with kaons on ME "<< endl;
+               if(fselect==3) cout << "TASK::Correlation with kzeros on ME "<< endl;
+       }
        if (!fInputEvent) {
                Error("UserExec","NO EVENT FOUND!");
                return;
@@ -196,35 +210,30 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
          AliError("AOD event not found!");
          return;
        }
-       Double_t pi = TMath::Pi();
+       
+       fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
        
        fEvents++; // event counter
        ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
        fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+       
        if(fmontecarlo && !fmcArray){
          AliError("Array of MC particles not found");
          return;
        }
-       
-       // initialize the pool for event mixing
-       Int_t multiplicity = aodEvent->GetNTracks();
-       AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
-       Double_t zvertex = vtx->GetZ();
-       Double_t multip = multiplicity;
-       
-       if(TMath::Abs(zvertex)>=10 || multip>500 || multip == 0) {
-               AliInfo(Form("Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,multip));
+       Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
+       if(!isEvSel) return;
+       ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
+       //
+       Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
+       if(!correlatorON) {
+               AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
                return;
        }
        
-       
-       AliEventPool* pool = fPoolMgr->GetEventPool(multip, zvertex);
-       if (!pool) AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f cm", multip, zvertex));
-       
-       
-
+       if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);
        // D* reconstruction
-       
+       cout << "USER EXEC CHECKPOINT 2" << endl;
        TClonesArray *arrayDStartoD0pi=0;
 
        
@@ -257,20 +266,21 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
        Double_t deltainvMDStar;
 
        
-       Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-       Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
+       Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
+       Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
        
+       cout << "USER EXEC CHECKPOINT 3" << endl;
        
-       
-       if(fselect ==3){// check the K0 invariant mass
-       TObjArray * fillkzerospectra = AcceptAndReduceKZero(aodEvent, 0,0);
-       delete fillkzerospectra;
-       }
+       //if(fselect ==3){// check the K0 invariant mass
+       //TObjArray * fillkzerospectra = AcceptAndReduceKZero(aodEvent, 0,0);
+       //delete fillkzerospectra;
+       //}
        
        //MC tagging for DStar
        //D* and D0 prongs needed to MatchToMC method
        Int_t pdgDgDStartoD0pi[2]={421,211};
        Int_t pdgDgD0toKpi[2]={321,211};
+
        
        //loop on D* candidates
        for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
@@ -283,26 +293,24 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                invMassDZero = - 999;
                deltainvMDStar = -998;
                
-       
-
                AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
                if(!dstarD0pi->GetSecondaryVtx()) continue;
                AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
                if (!theD0particle) continue;
+
                
                // track quality cuts
                Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
                // region of interest + topological cuts + PID
                Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
-               
                //apply selections
                if(!isTkSelected) continue;
                if(!isSelected) continue;
                if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;    
-
+               Int_t mcLabelDStar = -999;
                if(fmontecarlo){
                        // find associated MC particle for D* ->D0toKpi
-                       Int_t mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
+                       mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray,kFALSE);
                        if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
                }
                
@@ -310,31 +318,45 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                phiDStar = dstarD0pi->Phi(); 
                etaDStar = dstarD0pi->Eta();
                
+               phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the phi of the D meson in the correct range
+               
                Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
                
                Double_t dmDStarWindow =0.0019;// 0.0019 = 3 sigma
                Double_t mD0Window=0.074;
                
-               if (ptbin==1)  mD0Window = 0.026; //0.5-1
-               if (ptbin==2)  mD0Window = 0.022; //1-2
-               if (ptbin==3)  mD0Window = 0.024; //2-3
-               if (ptbin==4)  mD0Window = 0.032;
-               if (ptbin==5)  mD0Window = 0.032;
-               if (ptbin==6)  mD0Window = 0.036;
-               if (ptbin==7)  mD0Window = 0.036;
-               if (ptbin==8)  mD0Window = 0.036;
-               if (ptbin==9)  mD0Window = 0.058;
-               if (ptbin==10) mD0Window = 0.058;
-               if (ptbin>10)  mD0Window = 0.074;
-
+               if (!fSystem){ // pp
+                       if (ptbin==1)  mD0Window = 0.026; //0.5-1
+                       if (ptbin==2)  mD0Window = 0.022; //1-2
+                       if (ptbin==3)  mD0Window = 0.024; //2-3
+                       if (ptbin==4)  mD0Window = 0.032;
+                       if (ptbin==5)  mD0Window = 0.032;
+                       if (ptbin==6)  mD0Window = 0.036;
+                       if (ptbin==7)  mD0Window = 0.036;
+                       if (ptbin==8)  mD0Window = 0.036;
+                       if (ptbin==9)  mD0Window = 0.058;
+                       if (ptbin==10) mD0Window = 0.058;
+                       if (ptbin>10)  mD0Window = 0.074;
+               }
+               if(fSystem){// PbPb
+                       if (ptbin==0)  mD0Window = 0.032; //1-1
+                       if (ptbin==1)  mD0Window = 0.032; //2-3
+                       if (ptbin==2)  mD0Window = 0.032; //3-4
+                       if (ptbin==3)  mD0Window = 0.032; //4-5
+                       if (ptbin==4)  mD0Window = 0.036; //5-6
+                       if (ptbin==5)  mD0Window = 0.036; //6-8
+                       if (ptbin==6)  mD0Window = 0.055; //8-12
+                       if (ptbin==7)  mD0Window = 0.074; //12-16
+                       if (ptbin==8)  mD0Window = 0.074; //16-24
+                       if (ptbin==9)  mD0Window = 0.074; //24-35 
+                }
                
                invMassDZero = dstarD0pi->InvMassD0();
                ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
                
                deltainvMDStar = dstarD0pi->DeltaInvMass();
                
-               
-               
+       
                //good candidates
                if (TMath::Abs(invMassDZero-mPDGD0)<mD0Window){
                        
@@ -343,6 +365,7 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                                
                                ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar);
                                isInPeak = kTRUE;
+                               ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);
                        }
                }// end if good candidates
                
@@ -354,155 +377,113 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                        if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
                                ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar);
                                isInSideBand = kTRUE;
+                               ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);
                        }
                        
                }//end if sidebands
-        
+               // getting the number of triggers in the MCtag D* case 
+               
 
+        if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);
+               
                
                if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
                
-       
+               
+           fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
+                               
+               
+               
                Int_t trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
                Int_t trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
                Int_t trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
                
-               ptDStar = dstarD0pi->Pt();
-               phiDStar = dstarD0pi->Phi(); 
-               etaDStar = dstarD0pi->Eta();
+               Bool_t execPool = fCorrelator->ProcessEventPool();
+               if(fmixing && !execPool) {
+                       AliInfo("Mixed event analysis: pool is not ready");
+                       continue;
+               }
+               
+               Int_t NofEventsinPool = 1;
+               if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
                                
-               if(!fmixing){ // analysis on Single Event
-                                               
-                       TObjArray* selectedtracks = new TObjArray();
-                       if(fselect==1 || fselect ==2)   selectedtracks = AcceptAndReduceTracks(aodEvent);
-                       if(fselect==3) {selectedtracks = AcceptAndReduceKZero(aodEvent,iDStartoD0pi,1);}        
-                       Int_t noftracks = selectedtracks->GetEntriesFast(); 
-                       Int_t counterPeak =0;
-                       Int_t counterSB = 0;
-
-                       for(Int_t i =0; i<noftracks; i++){ // loop on tracks/k0 candidates in aod event
-
-                               AliReducedParticle *redpart = (AliReducedParticle*)selectedtracks->At(i);
-                               Double_t phiHad=redpart->Phi();
-
-                               Double_t ptHad=redpart->Pt();
-                               Double_t etaHad=redpart->Eta();
-
-                               Int_t label = redpart->GetLabel();
-
-                               Int_t trackid = redpart->GetID();
-
-
-                               // check that the track is not a D* daughter
-                               if(trackid == trackiddaugh0) continue;
-                               if(trackid == trackiddaugh1) continue;
-                               if(trackid == trackidsoftPi) continue;
+               for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
+               
+                       Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
+                       
+                       if(!analyzetracks) {
+                               AliInfo("AliHFCorrelator::Cannot process the track array");
+                               continue;
+                       }
+               
+               
+                       Int_t NofTracks = fCorrelator->GetNofTracks();
+               
+                       for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
+                       
+                               Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
+                               if(!runcorrelation) continue;
+                       
+                               Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
+                               Double_t DeltaEta = fCorrelator->GetDeltaEta();
+                       
+                               AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
+                               
+                               Double_t ptHad = hadron->Pt();
+                               Double_t phiHad = hadron->Phi();
+                               Double_t etaHad = hadron->Eta(); 
+                               Double_t label = hadron->GetLabel(); 
+                               Double_t trackid = hadron->GetID();
+                               
+                               phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
+                               
+                               if(!fmixing){ // skip D* Daughetrs
+                                       if(trackid == trackiddaugh0) continue;
+                                       if(trackid == trackiddaugh1) continue;
+                                       if(trackid == trackidsoftPi) continue;
+                               }
+                       
+                               // from here on it is up to the user to decide what object to fill
                                
                                if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
-                                       
+                               
                                        Int_t PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
-
-                                       FillMCTagCorrelations(ptDStar,phiDStar,etaDStar,ptHad,phiHad,etaHad,PartSource);
+                               
+                                       FillMCTagCorrelations(ptDStar,DeltaPhi,DeltaEta,ptHad,PartSource);
+                               
                                }
+                       
                                if(isInPeak)  {
-                                       FillCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad);// function for correlations
-                                       counterPeak++;
-                                       if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
-                                       if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
-
-                                       ((TH1F*)fOutput->FindObject("PhiTrigger"))->Fill(phiDStar);
-
-                                       
-                                               
-                                               if (phiHad > 1.5*pi) phiHad = phiHad - 2*pi;
-                                               if (phiHad < -0.5*pi) phiHad = phiHad + 2*pi;
-                                               ((TH1F*)fOutput->FindObject("PhiPart"))->Fill(phiHad);
-                                       
-                               }
                                
-                               if(isInSideBand) {
-
-                                       FillSideBandCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad); // function for sidebands
-                                       if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
-                                       if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
-                                       ((TH1F*)fOutput->FindObject("PhiSideBand"))->Fill(phiDStar);
-
-                                       counterSB++;
-                               }
+                                       if(fselect==1) ((TH3D*)fOutput->FindObject("DPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                                       if(fselect==2) ((TH3D*)fOutput->FindObject("DPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                                       if(fselect==3) ((TH3D*)fOutput->FindObject("DPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                                   ((TH2F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad);
+                                       //counterPeak++; // count tracks per peak per event
                                
-                       } // end loop on tracks 
-                       
-                       if(counterPeak) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(counterPeak);
-                       if(counterSB) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(counterSB);
-                       
-                       
-                       
+                               }
                        
-               } // end if SE Analysis
-               
-               if(fmixing) { // analysis on Mixed Events
-                       if (pool->IsReady()|| pool->GetCurrentNEvents()>=5){ // check if the pool is ready
-               
-                               pool->PrintInfo();
-       
-                               Int_t multbinflag = pool->MultBinIndex();
-                               Int_t zvtxflag = pool->ZvtxBinIndex();
-                               if(isInPeak) ((TH2I*)fOutput->FindObject("EventMixingCheck"))->Fill(multbinflag,zvtxflag); 
+                               if(isInSideBand) {
+                               
+                                       if(fselect==1) ((TH3D*)fOutput->FindObject("bkgDPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                                       if(fselect==2) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                                       if(fselect==3) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
+                               
                                        
-
-                               TObjArray* mixedtracks = 0x0;
                                
-                               for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) {//loop over the events in the pool
-                                       mixedtracks = pool->GetEvent(jMix);
-                                       if(!mixedtracks) cout << "No Mixed tracks " << endl;
-                                       Int_t jMax = mixedtracks->GetEntriesFast();
-
-                                       for(Int_t iMix =0; iMix<jMax; iMix++){ //loop on the tracks of the event
-                                               AliVParticle *redpart = (AliVParticle*)mixedtracks->At(iMix);
-                                               Double_t phiHad=redpart->Phi();
-                                               Double_t etaHad=redpart->Eta();
-                                               Double_t ptHad=redpart->Pt();
-                                               Int_t label = redpart->GetLabel();
-                                               
-                                               
-                                               if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
-                                                       
-                                                       Int_t PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
-                                                       
-                                                       FillMCTagCorrelations(ptDStar,phiDStar,etaDStar,ptHad,phiHad,etaHad,PartSource);
-                                               }
-
-                                               if(isInPeak) {
-                                                       FillCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad);// function for correlations
-
-                                                               if (phiDStar > 1.5*pi) phiDStar = phiDStar - 2*pi;
-                                                               if (phiDStar < -0.5*pi) phiDStar = phiDStar + 2*pi;
-                                                               
-                                                               ((TH1F*)fOutput->FindObject("PhiTrigger"))->Fill(phiDStar);
-                                                               if (phiHad > 1.5*pi) phiHad = phiHad - 2*pi;
-                                                               if (phiHad < -0.5*pi) phiHad = phiHad + 2*pi;
-                                                               ((TH1F*)fOutput->FindObject("PhiPart"))->Fill(phiHad);
-                                                       
-
-                                                       
-                                               }
-                                               
-                                               if(isInSideBand) FillSideBandCorrelations(ptDStar,phiDStar,etaDStar,phiHad,etaHad); // function for sidebands
-
-                                               } // end loop on tracks
-                               }// end loop on events in pool
-                       } // end if pool is ready
+                                       //counterSB++;
+                               }
                        
-               } // end ME analysis
-               
+                       
+                       } // end loop on track candidates
+               } // end loop on events in the pool
+                               
        }// end loop on D* candidates           
+       cout << "USER EXEC CHECKPOINT 4" << endl;
+       Bool_t updated = fCorrelator->PoolUpdate();
+       if(!updated) AliInfo("Pool was not updated");
        
-       if(fmixing) { // update the pool for Event Mixing
-               if(fselect==1 || fselect==2)pool->UpdatePool(AcceptAndReduceTracks(aodEvent)); // updating the pool for hadrons
-               if(fselect==3) pool->UpdatePool(AcceptAndReduceKZero(aodEvent,0,0)); // updating the pool for K0s
-       }
-       //cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END OF THE EVENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
+               //cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END OF THE EVENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
        
        PostData(1,fOutput); // set the outputs
        PostData(3,fCounter); // set the outputs
@@ -529,93 +510,6 @@ void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
 
 
 //_____________________________________________________
-TObjArray*  AliAnalysisTaskDStarCorrelations::AcceptAndReduceTracks(AliAODEvent* inputEvent){
-       
-       Int_t nTracks = inputEvent->GetNTracks();
-       AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
-       Double_t Bz = inputEvent->GetMagneticField();
-
-       
-       TObjArray* tracksClone = new TObjArray;
-       tracksClone->SetOwner(kTRUE);
-       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
-               AliAODTrack* track = inputEvent->GetTrack(iTrack);
-               if (! track) continue;
-
-               if(!fCuts2->IsHadronSelected(track,vtx,Bz)) continue; // apply selection of hadrons
-                       
-               
-               if(fselect ==2){        
-                       if(!fCuts2->CheckKaonCompatibility(track,fmontecarlo,fmcArray)) continue; // check if it is a Kaon - data and MC
-                       }
-      
-               AliVParticle * part = (AliVParticle*)track;
-               tracksClone->Add(new AliReducedParticle(part->Eta(), part->Phi(), part->Pt(),track->GetLabel(),track->GetID()));
-               
-       }
-       return tracksClone;
-}
-
-//_____________________________________________________
-TObjArray*  AliAnalysisTaskDStarCorrelations::AcceptAndReduceKZero(AliAODEvent* inputEvent, Int_t loopindex, Int_t plotassociation){
-       
-       Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
-       TObjArray* KZeroClone = new TObjArray;
-       AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();
-       Int_t v0label = -1;
-       Int_t pdgDgK0toPipi[2] = {211,211};
-       Double_t mPDGK0=TDatabasePDG::Instance()->GetParticle(310)->Mass();
-       const Int_t size = inputEvent->GetNumberOfV0s();
-       Int_t prevnegID[size];
-       Int_t prevposID[size];
-       for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
-               AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);
-               if(!v0) {cout << "is not a v0 at step " << iv0 << endl; continue;}
-               
-               if(!fCuts2->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
-           
-               // checks to avoid double counting
-               Int_t negID = -999;
-               Int_t posID = -998;
-               //Int_t a = 0;
-               prevnegID[iv0] = -997;
-               prevposID[iv0]= -996;
-               negID = v0->GetNegID();
-               posID = v0->GetPosID();
-               Bool_t DoubleCounting = kFALSE;
-               
-               for(Int_t k=0; k<iv0; k++){
-                       if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){
-                               DoubleCounting = kTRUE;
-                               //a=k;
-                               break;
-                       }//end if
-               } // end for
-               
-               if(DoubleCounting) {cout << "SKIPPING DOUBLE COUNTING" << endl;continue;}
-               else{ prevposID[iv0] = posID; prevnegID[iv0] = negID;}
-               
-               if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
-               Double_t k0pt = v0->Pt();
-               Double_t k0eta = v0->Eta();
-               Double_t k0Phi = v0->Phi();
-               Double_t k0InvMass = v0->MassK0Short(); 
-
-               if (loopindex == 0) {
-                       if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
-                       if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*
-               }
-               // if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
-               
-               if(TMath::Abs(k0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
-               KZeroClone->Add(new AliReducedParticle(k0eta , k0Phi, k0pt,v0label,0));
-               
-       }
-       
-       return KZeroClone;
-}
-
-//_____________________________________________________
 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        
        Double_t Pi = TMath::Pi();
@@ -647,6 +541,9 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",30,0,30);
        fOutput->Add(RecoPtBkg);
        
+       TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",30,0,30);
+       fOutput->Add(MCtagPtDStar);
+       
        TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
        if(fselect==3) fOutput->Add(KZeroSpectra);
        
@@ -672,14 +569,14 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        MCSources->GetXaxis()->SetBinLabel(5," from b->B");
        if(fmontecarlo) fOutput->Add(MCSources);
        
-       TH1F * PhiTrigger = new TH1F("PhiTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi);
-       fOutput->Add(PhiTrigger);
+       TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
+       fOutput->Add(PhiEtaTrigger);
        
-       TH1F * PhiSideBand = new TH1F("PhiSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi);
-       fOutput->Add(PhiSideBand);
+       TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
+       fOutput->Add(PhiEtaSideBand);
        
-       TH1F * PhiPart = new TH1F("PhiPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi);
-       fOutput->Add(PhiPart);
+       TH2F * PhiEtaPart = new TH2F("PhiEtaPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
+       fOutput->Add(PhiEtaPart);
 
 
        //correlations histograms
@@ -745,68 +642,23 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
 }
 
 
-//____________________________  Function for correlations ___________________________________________________
-void AliAnalysisTaskDStarCorrelations::FillCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack){
-       Double_t pi = TMath::Pi();
-       Double_t deltaPhi, deltaEta;
-       deltaPhi = phiTrig - phiTrack;
-       deltaEta = etaTrig - etaTrack;
-       // set correct Delta Phi range
-       if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
-       if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
-       
-               if(fselect==1) ((TH3D*)fOutput->FindObject("DPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
-               if(fselect==2) ((TH3D*)fOutput->FindObject("DPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
-               if(fselect==3) ((TH3D*)fOutput->FindObject("DPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
-       
-       
-       return;
-}
 
-//____________________________  Function for sidebands ___________________________________________________
-void AliAnalysisTaskDStarCorrelations::FillSideBandCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack){
-       
-       Double_t pi = TMath::Pi();
-       Double_t deltaPhi, deltaEta;
-       deltaPhi = phiTrig - phiTrack;
-       deltaEta = etaTrig - etaTrack;
-       // set correct Delta Phi range
-       if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
-       if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
-       
-               if(fselect==1) ((TH3D*)fOutput->FindObject("bkgDPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
-               if(fselect==2) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
-               if(fselect==3) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
-       
-       
-       return;
-       
-       
-}
+//____________________________  Function for MC correlations ___________________________________________________
+void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Int_t mcSource){
 
-//____________________________  Function for sidebands ___________________________________________________
-void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t phiTrig,  Double_t etaTrig, Double_t ptTrack, Double_t phiTrack, Double_t etaTrack, Int_t mcSource){
-Double_t deltaPhi = phiTrig - phiTrack;
-Double_t deltaEta = etaTrig - etaTrack;
-// set correct Delta Phi range
-       
-       Double_t pi = TMath::Pi();
-       
-if (deltaPhi > 1.5*pi -pi/32) deltaPhi = deltaPhi - 2*pi;
-if (deltaPhi < -0.5*pi -pi/32) deltaPhi = deltaPhi + 2*pi;
 
-       if(fselect==1) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarHadron"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKaon"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKZero"))->Fill(deltaPhi,ptTrig,deltaEta);
+       if(fselect==1) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKaon"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==3) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKZero"))->Fill(DelPhi,ptTrig,DelEta);
 
 
 
 ((TH1F*)fOutput->FindObject("MCSources"))->Fill(0);
 
 if (mcSource==44){ // is from charm ->D
-       if(fselect==1) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
+       if(fselect==1) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==3) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
        
        
        ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
@@ -814,17 +666,17 @@ if (mcSource==44){ // is from charm ->D
        }
 
 if (mcSource==54){ // is from beauty -> D
-       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
+       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
        if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
        if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(3);
        }
 
 if (mcSource==55){ // is from beauty ->B
-       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarHadronMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarKaonMC"))->Fill(deltaPhi,ptTrig,deltaEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyOriginBDPhiDStarKZeroMC"))->Fill(deltaPhi,ptTrig,deltaEta);
+       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyOriginBDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
        if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
        if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(4);
        }
index 3930e95..9a4ec13 100644 (file)
@@ -37,6 +37,7 @@
 #include "AliAODRecoCascadeHF.h"
 #include "AliHFAssociatedTrackCuts.h"
 #include "AliNormalizationCounter.h"
+#include "AliHFCorrelator.h"
 
 class TParticle ;
 class TClonesArray ;
@@ -66,22 +67,17 @@ class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE
        virtual void UserExec(Option_t *option);
        virtual void Terminate(Option_t *option);
        
-       
-       // methods to get the tracks to correlate
-       TObjArray*  AcceptAndReduceTracks(AliAODEvent* inputEvent);
-       TObjArray*  AcceptAndReduceKZero(AliAODEvent* inputEvent, Int_t loopindex, Int_t plotassociation);
 
        void DefineHistoForAnalysis();
        
        // correlators
-       void FillCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack);
-       void FillSideBandCorrelations(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig, Double_t phiTrack, Double_t etaTrack);
-       void FillMCTagCorrelations(Double_t ptTrig, Double_t phiTrig,  Double_t etaTrig, Double_t ptTrack, Double_t phiTrack, Double_t etaTrack, Int_t mcSource);
-       
+       void FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Int_t mcSource);       
        // setters
        void SetMonteCarlo(Bool_t k) {fmontecarlo = k;}
        void SetUseMixing (Bool_t j) {fmixing = j;}
-       void SetCorrelator(Int_t l) {fselect = l;}
+       void SetCorrelator(Int_t l) {fselect = l;} // select 1 for hadrons, 2 for Kaons, 3 for Kzeros
+       void SetUseDisplacement(Int_t m) {fDisplacement=m;} // select 0 for no displ, 1 for abs displ, 2 for d0/sigma_d0
+       void SetRunPbPb(Bool_t system){fSystem=system;} // select between pp (kFALSE) or PbPb (kTRUE)
        
 
 
@@ -91,23 +87,27 @@ class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE
        AliAnalysisTaskDStarCorrelations& operator=(const AliAnalysisTaskDStarCorrelations& source); 
 
        TObject* fhandler; //! Analysis Handler
-       AliEventPoolManager*     fPoolMgr;         //! event pool manager
+       //      AliEventPoolManager*     fPoolMgr;         //! event pool manager
        TClonesArray* fmcArray; //mcarray
        AliNormalizationCounter *fCounter; // counter
+    AliHFCorrelator * fCorrelator; // object for correlations
+
        
        
        Int_t fselect; // select what to correlate with a D* 1-chargedtracks,2-chargedkaons,3-k0s
        Bool_t fmontecarlo;//switch for MC
        Bool_t fmixing;// switch for event mixing
+       Bool_t fSystem; // pp or PbPb
        Int_t fEvents; //! number of event
        Int_t fDebug; //! debug level
+       Int_t fDisplacement; // set 0 for no displacement cut, 1 for absolute d0, 2 for d0/sigma_d0
        
        
        TList *fOutput;                  //! user output data
        AliRDHFCutsDStartoKpipi *fCuts;  // Cuts D*
        AliHFAssociatedTrackCuts *fCuts2; // cuts for associated 
                                          
-       ClassDef(AliAnalysisTaskDStarCorrelations,1); // class for D meson correlations                           
+       ClassDef(AliAnalysisTaskDStarCorrelations,2); // class for D meson correlations                           
 };
 
 #endif
index bd33b81..3087627 100644 (file)
@@ -74,7 +74,11 @@ AliAnalysisTaskSE(),
   fNentries(0), 
   fCutsD0(0),
   fCutsTracks(0),
+  fCorrelatorTr(0),
+  fCorrelatorKc(0),
+  fCorrelatorK0(0),
   fReadMC(0),
+  fMixing(kFALSE),
   fCounter(0),
   fNPtBins(1),
   fFillOnlyD0D0bar(0),
@@ -102,7 +106,11 @@ AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *nam
   fNentries(0),
   fCutsD0(0),
   fCutsTracks(cutsTrk),
+  fCorrelatorTr(0),
+  fCorrelatorKc(0),
+  fCorrelatorK0(0),
   fReadMC(0),
+  fMixing(kFALSE),
   fCounter(0),
   fNPtBins(1),
   fFillOnlyD0D0bar(0),
@@ -148,7 +156,11 @@ 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),
+  fMixing(source.fMixing),
   fCounter(source.fCounter),
   fNPtBins(source.fNPtBins),
   fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
@@ -183,7 +195,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,7 +232,11 @@ 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;
+  fMixing = orig.fMixing;
   fCounter = orig.fCounter;
   fNPtBins = orig.fNPtBins;
   fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
@@ -232,6 +260,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 +276,33 @@ void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
   //
   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
 
+  //HFCorrelator creation and definition
+  fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys);
+  fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys);
+  fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys);
+  fCorrelatorTr->SetDeltaPhiInterval(-1.57,4.71);// set the Delta Phi Interval you want (in this case -0.5Pi to 1.5 Pi)
+  fCorrelatorKc->SetDeltaPhiInterval(-1.57,4.71);
+  fCorrelatorK0->SetDeltaPhiInterval(-1.57,4.71);
+  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);
+  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();
@@ -405,7 +461,7 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
       return;
     }
   }
-  
+
   //histogram filled with 1 for every AOD
   fNentries->Fill(0);
   fCounter->StoreEvent(aod,fCutsD0,fReadMC); 
@@ -439,9 +495,25 @@ 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();
@@ -454,7 +526,7 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
   //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,7 +539,7 @@ 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 && (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_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
       }
@@ -478,7 +550,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());
@@ -504,20 +576,42 @@ void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
       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
+      fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
       if(!fIsSelectedCandidate) continue;
 
-      if(!fReadMC) CalculateCorrelations(aod,d); //correlations on real data
-      else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
+      //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) {
+        CalculateCorrelations(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);
+        if (labD0>-1) {
+          CalculateCorrelations(d,labD0,mcArray);
+        }
       }
       
       FillMassHists(d,mcArray,fCutsD0,fOutputMass);
     }
 
   } //end for prongs
+
+  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");
+  }
+
   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);  
   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);  
 
@@ -666,7 +760,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;
@@ -705,276 +799,201 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
   TString namePlot = "";
 
   //These for limits in THnSparse (one per bin, same limits). 
-  //Vars: DeltaPhi, InvMass, PtTrack
+  //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
   Int_t nBinsPhi[5] = {32,150,30,3,16};
   Double_t binMinPhi[5] = {-1.6,1.6,0.,0.,-1.6};  //is the minimum for all the bins
   Double_t binMaxPhi[5] = {4.8,2.2,3.0,3.,1.6};  //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_K0_Bin";
-    namePlot+=i;
+  //Vars: DeltaPhi, InvMass, DeltaEta
+  Int_t nBinsMix[3] = {32,150,16};
+  Double_t binMinMix[3] = {-1.6,1.6,-1.6};  //is the minimum for all the bins
+  Double_t binMaxMix[3] = {4.8,2.2,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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-    hPhiK->Sumw2();
-    fOutputCorr->Add(hPhiK);
-
-    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)",5,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)",5,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_K0_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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiK_c->Sumw2();
-      fOutputCorr->Add(hPhiK_c);
+      THnSparseI *hPhiK = new THnSparseI(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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiH_c->Sumw2();
-      fOutputCorr->Add(hPhiH_c);
+      THnSparseI *hPhiH = new THnSparseI(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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_c->Sumw2();
-      fOutputCorr->Add(hPhiC_c);
+      THnSparseI *hPhiC = new THnSparseI(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_K0_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)",5,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)",5,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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_b->Sumw2();
-      fOutputCorr->Add(hPhiC_b);
-
-      //HF-only tracks (c for c->D0, b for b->D0)
-      namePlot="hPhi_K0_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)",5,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)",5,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)",5,nBinsPhi,binMinPhi,binMaxPhi);
-      hPhiC_HF_c->Sumw2();
-      fOutputCorr->Add(hPhiC_HF_c);
-
-      namePlot="hPhi_K0_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)",5,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)",5,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)",5,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);
-
-      namePlot="hPhi_Lead_From_b_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);
-
-      namePlot="hPhi_Lead_HF_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);
-
-      namePlot="hPhi_Lead_HF_From_b_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;
+      //histos for c/b origin for D0 (MC only)
+      if (fReadMC) {
 
-    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);
+        //generic origin for tracks
+        namePlot="hPhi_K0_From_c_Bin";
+        namePlot+=i;
 
-    if (fReadMC) {
-      namePlot="hPhi_Weig_From_c_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)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiK_c->Sumw2();
+        fOutputCorr->Add(hPhiK_c);
 
-      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);
+        namePlot="hPhi_Kcharg_From_c_Bin";
+        namePlot+=i;
 
-      namePlot="hPhi_Weig_From_b_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)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiH_c->Sumw2();
+        fOutputCorr->Add(hPhiH_c);
 
-      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);
+        namePlot="hPhi_Charg_From_c_Bin";
+        namePlot+=i;
 
-      namePlot="hPhi_Weig_HF_From_c_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)",5,nBinsPhi,binMinPhi,binMaxPhi);
+        hPhiC_c->Sumw2();
+        fOutputCorr->Add(hPhiC_c);
+  
+        namePlot="hPhi_K0_From_b_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);
+        THnSparseI *hPhiK_b = new THnSparseI(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);
 
-      namePlot="hPhi_Weig_HF_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_From_b_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);
-    }
+        THnSparseI *hPhiH_b = new THnSparseI(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);
 
-    //pT weighted correlations (squared weights)
-    namePlot="hPhi_WeigSqr_Bin";
-    namePlot+=i;
+        namePlot="hPhi_Charg_From_b_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);
+        THnSparseI *hPhiC_b = new THnSparseI(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);
 
-    if (fReadMC) {
-      namePlot="hPhi_WeigSqr_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 *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);
+        THnSparseI *hPhiK_HF_c = new THnSparseI(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_WeigSqr_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Kcharg_HF_From_c_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);
+        THnSparseI *hPhiH_HF_c = new THnSparseI(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);
 
-      namePlot="hPhi_WeigSqr_HF_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_HF_From_c_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);
+        THnSparseI *hPhiC_HF_c = new THnSparseI(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);
 
-      namePlot="hPhi_WeigSqr_HF_From_b_Bin";
-      namePlot+=i;
+        namePlot="hPhi_K0_HF_From_b_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);
-    }
+        THnSparseI *hPhiK_HF_b = new THnSparseI(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);
 
-    //pT weighted correlations (inverse of pT distribution weights)
-    namePlot="hPhi_WeigDist_Bin";
-    namePlot+=i;
+        namePlot="hPhi_Kcharg_HF_From_b_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);
+        THnSparseI *hPhiH_HF_b = new THnSparseI(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);
 
-    if (fReadMC) {
-      namePlot="hPhi_WeigDist_From_c_Bin";
-      namePlot+=i;
+        namePlot="hPhi_Charg_HF_From_b_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);
+        THnSparseI *hPhiC_HF_b = new THnSparseI(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_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);
+      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);
 
-      namePlot="hPhi_WeigDist_HF_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);
+      if (fReadMC) {
+        namePlot="hPhi_Lead_From_c_Bin";
+        namePlot+=i;
 
-      namePlot="hPhi_WeigDist_HF_From_b_Bin";
+        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);
+  
+        namePlot="hPhi_Lead_From_b_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);
+  
+        namePlot="hPhi_Lead_HF_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);
+  
+        namePlot="hPhi_Lead_HF_From_b_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;
-
-      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_K0_Bin"; namePlot+=i;
-    TH1F *hCountK = new TH1F(namePlot.Data(), "Kaons counter; # Tracks",10,0.,10.);
-    hCountK->SetMinimum(0);
-    fOutputStudy->Add(hCountK);
+  
+      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);
+  
+      if (fReadMC) {
+        namePlot="hPhi_Weig_From_c_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);
+  
+        namePlot="hPhi_Weig_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);
+  
+        namePlot="hPhi_Weig_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);
+  
+        namePlot="hPhi_Weig_HF_From_b_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);
+      }
 
     //pT distribution histos
     namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
@@ -998,7 +1017,33 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
     hDstarPions->SetMinimum(0);
-    fOutputStudy->Add(hDstarPions);
+    fOutputStudy->Add(hDstarPions); 
+    }
+
+    if(fMixing) {
+      //THnSparse plots for event mixing!
+      namePlot="hPhi_K0_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+
+      THnSparseI *hPhiK_EvMix = new THnSparseI(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      hPhiK_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiK_EvMix);
+
+      namePlot="hPhi_Kcharg_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+  
+      THnSparseI *hPhiH_EvMix = new THnSparseI(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      hPhiH_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiH_EvMix);
+
+      namePlot="hPhi_Charg_Bin";
+      namePlot+=i;namePlot+="_EvMix";
+
+      THnSparseI *hPhiC_EvMix = new THnSparseI(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
+      hPhiC_EvMix->Sumw2();
+      fOutputCorr->Add(hPhiC_EvMix);
+    }
   }
 
   //out of bin loop
@@ -1008,6 +1053,18 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
   hRejTracks->GetXaxis()->SetBinLabel(2,"Rejected");
   fOutputStudy->Add(hRejTracks);
 
+  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 (fFillGlobal) { //all-events plots
     //pt distributions
     TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
@@ -1018,10 +1075,27 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     hPtHAll->SetMinimum(0);
     fOutputStudy->Add(hPtHAll);
 
-    TH1F *hPtKAll = new TH1F("hist_Pt_K0_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);
 
+    //phi distributions
+    TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
+    hPhiDistCAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistCAll);
+
+    TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
+    hPhiDistHAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistHAll);
+
+    TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
+    hPhiDistKAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistKAll);
+
+    TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
+    hPhiDistDAll->SetMinimum(0);
+    fOutputStudy->Add(hPhiDistDAll);
+
     //K0 Invariant Mass plots
     TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
     hK0MassInv->SetMinimum(0);
@@ -1178,345 +1252,317 @@ void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
     }
   }
 
+  if(fMixing) { //check for event mixing (events in bins)
+    TH2I* EventMixingCheckTr = new TH2I("EventMixingCheckTr","EventMixingCheckTr",5,-0.5,4.5,7,-0.5,6.5);
+    fOutputStudy->Add(EventMixingCheckTr);
+
+    TH2I* EventMixingCheckK0 = new TH2I("EventMixingCheckK0","EventMixingCheckK0",5,-0.5,4.5,7,-0.5,6.5);
+    fOutputStudy->Add(EventMixingCheckK0);
+  }
+
 }
 
 //________________________________________________________________________
-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, deltaeta;
+  Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
+  Double_t mD0, mD0bar;
+  Int_t origD0 = 0, PDGD0 = 0;
   d->InvMassD0(mD0,mD0bar);
   Int_t ptbin = PtBinCorr(d->Pt());
   if(ptbin < 0) return;
-  AliAODVertex *vtx = (AliAODVertex*)aod->GetPrimaryVertex();
 
-  Int_t N_Charg = 0, N_Kcharg = 0, N_Kaons = 0;
+  //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;
+    }
+  }
+
+  Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
 
-  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)
+  //loop over the tracks in the pool 
+  Bool_t execPool = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
+               
+  if(fMixing && !execPool) {
+    AliInfo("Mixed event analysis: pool is not ready");
+    return;
+  }
+               
+  Int_t NofEventsinPool = 1;
+  if(fMixing) NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
+               
+  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
+    Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
+    Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
+    if(!analyzetracksTr || !analyzetracksKc || !analyzetracksK0) {
+      AliInfo("AliHFCorrelator::Cannot process the track array");
+      continue;
+    }
 
-    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
-
-      GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
-
-      ((TH1F*)fOutputStudy->FindObject("hRejTracks"))->Fill(0); //track accepted by quality cuts
-      deltaphi = d->Phi()-track->Phi();
-      deltaeta = d->Eta()-track->Eta();
-      if (deltaphi < -1.571) deltaphi+=6.283;
-      if (deltaphi > 4.712) deltaphi-=6.283;
-      Double_t pttrack = track->Pt();
-
-      if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi;}  //for leading particle
-
-      //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();
-      Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
-      if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
-      if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
-      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,d0/d0err,deltaeta};
-      Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,pttrack,d0/d0err,deltaeta};
-   
-      //generic charged tracks (NO PID selection)
-      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0);
+    //Charged tracks           
+    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(!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);
+        }
       }
-      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->Fill(fillSpPhiD0bar);
+      if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+
+      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
+   
+      FillSparsePlots(mcArray,d,origD0,PDGD0,track,kTrack); //fills for charged tracks
+
+      if(!fMixing) N_Charg++;
+
+      //retrieving leading info...
+      if(track->Pt() > highPt) {
+        lead[0] = fCorrelatorTr->GetDeltaPhi();
+        lead[1] = fCorrelatorTr->GetDeltaEta();
+        lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
+        highPt = track->Pt();
       }
-      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));
-        }       
-
-      //Charged Kaon 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);
+
+    } // end of tracks loop
+
+    //Charged Kaons loop
+    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(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
+      }
+      if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
 
-    } //end tracks loop
+      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
+   
+      FillSparsePlots(mcArray,d,origD0,PDGD0,kCharg,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
-
-      Double_t ptV0=v0->Pt(), deltaphiV0, deltaetaV0;
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
-      if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
-      deltaphiV0 = d->Phi()-v0->Phi();
-      deltaetaV0 = d->Eta()-v0->Eta();
-      if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
-      if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
-      if(deltaetaV0 > EtaLim_Sparse) deltaetaV0 = EtaLim_Sparse-0.01;
-      if(deltaetaV0 < -EtaLim_Sparse) deltaetaV0 = -EtaLim_Sparse+0.01;
-
-      Double_t fillSpPhiD0K0[5] = {deltaphiV0,mD0,ptV0,0.,deltaetaV0};
-      Double_t fillSpPhiD0barK0[5] = {deltaphiV0,mD0bar,ptV0,0.,deltaetaV0};
+    } // end of charged kaons loop
 
-      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
-      }
-      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
-      }
-      if(!fAlreadyFilled) ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K0_Bin%d",ptbin)))->Fill(v0->Pt());
-      N_Kaons++;
-    } // end kaon case
+    //K0 loop
+    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();
 
-    //Leading track correlations fill
-    if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
-     ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0);
+      if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
+
+      Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
+      if(k0->GetID() == idDaughs[0] || k0->GetID() == idDaughs[1]) continue; //discards daughters of candidate
+   
+      FillSparsePlots(mcArray,d,origD0,PDGD0,k0,kK0); //fills for charged tracks
+
+      if(!fMixing) N_Kaons++;
+
+    } // end of charged kaons loop
+
+  } //end of events loop
+
+  //leading track correlations fill
+  if(fReadMC) {
+    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[0],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[0],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[0],mD0);  
     }
-    if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
-     ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
+    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[0],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[0],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[0],mD0bar); 
     }
+  } else {
+      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[0],mD0); //c and b D0
+      if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[0],mD0bar);
+  }
 
-  } //end fReadMC = 0
+  //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(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)
+  fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
 
-    //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;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, AliAODRecoDecayHF2Prong *d, Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t type) {
+  //
+  //fills the THnSparse for correlations, calculating the variables
+  //
+
+  //Initialization of variables
+  Int_t ptbin = PtBinCorr(d->Pt());
+  if(ptbin < 0) return;
+  Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
+  d->InvMassD0(mD0,mD0bar);
+
+  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;
     }
+  }
 
-    for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tracks
-      AliAODTrack * track = aod->GetTrack(itrack);
+  //Fixes limits; needed to include overflow into THnSparse projections!
+  Double_t pTorig = track->Pt();
+  Double_t d0orig = track->GetImpPar();
+  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();
+  Double_t EtaLim_Sparse = ((THnSparseI*)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;
+  
+  if(fMixing == kSE) {
+
+  //variables for filling histos
+  Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
+  Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
 
-      if(!TrackSelectedInLoop(track,d,aod,ptbin,idDaughs)) continue; //rejection of track (daughter of D0/quality cuts not passed/soft pion/negative ID
-
-      Int_t label = track->GetLabel();
-      if (label<0) continue; //discard negative label tracks
-      Int_t PDGtrack = ((AliAODMCParticle*)mcArray->At(label))->GetPdgCode();
-
-      GetImpParameter(track,aod,d0,d0err); //evaluates d0 and sigma_d0
-
-      //Infos on track (origin, phi, eta)
-      Int_t origTr = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(label));
-      deltaphi = d->Phi()-track->Phi();
-      deltaeta = d->Eta()-track->Eta();
-      if (deltaphi < -1.571) deltaphi+=6.283;
-      if (deltaphi > 4.712) deltaphi-=6.283;
-      Double_t pttrack = track->Pt();
-
-      if (pttrack > highPt) {highPt = pttrack; lead[0] = pttrack; lead[1] = deltaphi; lead[2] = origTr;}  //for leading particle
-
-      //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();
-      Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
-      if(pttrack > ptLim_Sparse) pttrack = ptLim_Sparse-0.01;
-      if(d0/d0err > displLim_Sparse) d0 = (displLim_Sparse-0.001)*d0err;
-      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,d0/d0err,deltaeta};
-      Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,pttrack,d0/d0err,deltaeta};
-
-      //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(fReadMC == 0) {
+      //sparse fill for data (tracks, K+-, K0) + weighted
+      if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
+        ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0);
+        ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pTorig);
       }
-      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(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
+        ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+        ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pTorig);
       }
       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);
+       ((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_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
-        } 
-
-      //Charged Kaon identification (K 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 hadron case
+    }
 
-    } //end tracks loop
+    if(fReadMC) {
 
-    //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);
+      if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
 
-      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 ptV0=v0->Pt(), deltaphiV0, deltaetaV0;
-      Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
-      Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
-      deltaphiV0 = d->Phi()-v0->Phi();
-      deltaetaV0 = d->Eta()-v0->Eta();
-      if (deltaphiV0 < -1.571) deltaphiV0+=6.283;
-      if (deltaphiV0 > 4.712) deltaphiV0-=6.283;
-      if (ptV0 > ptLim_Sparse) ptV0 = ptLim_Sparse-0.01;
-      if(deltaetaV0 > EtaLim_Sparse) deltaetaV0 = EtaLim_Sparse-0.01;
-      if(deltaetaV0 < -EtaLim_Sparse) deltaetaV0 = -EtaLim_Sparse+0.01;
-
-      Double_t fillSpPhiD0K0[5] = {deltaphiV0,mD0,ptV0,0.,deltaetaV0};
-      Double_t fillSpPhiD0barK0[5] = {deltaphiV0,mD0bar,ptV0,0.,deltaetaV0};
-
-      Int_t origV0 = CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(labV0));
-
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0 (from MCTruth)
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0K0);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);
-         if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);  
-         if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0K0);  
-      }
-      if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_Bin%d",ptbin)))->Fill(fillSpPhiD0barK0);
-         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);
-         if(origD0==4&&origV0>=1&&origV0<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0);  
-         if(origD0==5&&origV0>=3&&origV0<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpPhiD0barK0); 
+      //sparse fill for data (tracks, K+-, K0) + weighted
+      if(PdgD0==421) { //D0 (from MCTruth)
+         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0);
+         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
+         if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
+         if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
+         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pTorig);
+         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
+         if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
+         if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
       }
+      if(PdgD0==-421) { //D0bar
+         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+         if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+         if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar); 
+         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pTorig);
+         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
+         if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
+         if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
+      } 
       if(!fAlreadyFilled) {
-         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_Bin%d",ptbin)))->Fill(0.); //Fills displacement histos
-        if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_HF_Bin%d",ptbin)))->Fill(0.);
-         if (origV0>=1&&origV0<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(0.);
-        ((TH1F*)fOutputStudy->FindObject(Form("histDispl_K0%s_Bin%d",orig.Data(),ptbin)))->Fill(0.); //Fills displacement histos
-         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_K0_Bin%d",ptbin)))->Fill(v0->Pt());
-         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_K0_Bin%d",ptbin)))->Fill(origV0);
+       ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
+        if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
+        if (origTr>=1&&origTr<=6) ((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);
       }
-      N_Kaons++;
-    } // end kaon case
+    }//end MC 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 of SE fill
 
-  } //end fReadMC = 1
+  if(fMixing == kME) {
 
-  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_K0_Bin%d",ptbin)))->Fill(N_Kaons);
-  }
+  //variables for filling histos
+  Double_t fillSpPhiD0[3] = {deltaphi,mD0,deltaeta};
+  Double_t fillSpPhiD0bar[3] = {deltaphi,mD0bar,deltaeta};
 
-  fAlreadyFilled=kTRUE; //distribution plots for tracks filled  
+    if(fReadMC == 0) {
+      //sparse fill for data (tracks, K+-, K0)
+      if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0);
+      if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+    }
+    if(fReadMC == 1) {
+      //sparse fill for data (tracks, K+-, K0)
+      if(PdgD0==421)  ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0);
+      if(PdgD0==-421) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
+    }//end MC case
 
+  } //end of ME fill
+  
+  return;
 }
 
 //_________________________________________________________________________________________________
@@ -1533,7 +1579,7 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
   // 7) c hadronization
   // 8) b hadronization
   //
-  printf(" AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
+  if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
        
   Int_t pdgGranma = 0;
   Int_t mother = 0;
@@ -1591,7 +1637,6 @@ Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, A
   else return 0; //no HF quark
 }
 
-
 //________________________________________________________________________
 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
   //
@@ -1606,62 +1651,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();
-  Double_t mD0, mD0bar;
-  d->InvMassD0(mD0,mD0bar);
-  
-  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)) {
-    if (fReadMC == 0 && (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3)) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
-    if (fReadMC == 0 && fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
-    if (fReadMC == 1) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,1.864);
-    output = kFALSE; 
-  } else {
-    if (fReadMC == 0 && (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3)) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0); 
-    if (fReadMC == 0 && fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar); 
-    if (fReadMC == 1) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,1.864);
-  }
-
-  return output;
-}
-
 //---------------------------------------------------------------------------
 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
 {
@@ -1676,9 +1665,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!
@@ -1715,5 +1702,7 @@ void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
   cout << "\n--------------------------\n";
   cout << "MC Truth = "<<fReadMC<<"\n";
   cout << "--------------------------\n";
+  cout << "Ev Mixing = "<<fMixing<<"\n";
+  cout << "--------------------------\n";
 }
 
index d034e12..5d2c878 100644 (file)
@@ -14,7 +14,6 @@
 // F.Colamaria, fabio.colamaria@ba.infn.it
 //*************************************************************************
 
-#include <vector>
 #include <TROOT.h>
 #include <TSystem.h>
 #include <TNtuple.h>
 #include "AliAnalysisTaskSE.h"
 #include "AliRDHFCutsD0toKpi.h"
 #include "AliHFAssociatedTrackCuts.h"
+#include "AliHFCorrelator.h"
 #include "AliNormalizationCounter.h"
 
-using std::vector;
-
 class AliAODEvent;
 
 class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
@@ -66,22 +64,23 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   void SetPtTreshUp(Double_t* pttreshup) {for(int i=0;i<fNPtBinsCorr;i++) {fPtThreshUp.push_back(pttreshup[i]);}}
   void PrintBinsAndLimits();
   Int_t PtBinCorr(Double_t pt) const;
+  void SetEvMixing(Bool_t mix) {fMixing=mix;}
+
+  enum PartType {kTrack,kKCharg,kK0};
+  enum FillType {kSE, kME}; //for single event or event mixing histos fill
 
  private:
 
   AliAnalysisTaskSED0Correlations(const AliAnalysisTaskSED0Correlations &source);
   AliAnalysisTaskSED0Correlations& operator=(const AliAnalysisTaskSED0Correlations& source); 
   void FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout);
-  void FillVarHists(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout);
   Int_t CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const;
   //correlation methods
   void CreateCorrelationsObjs();
-  void CalculateCorrelations(AliAODEvent* aod, AliAODRecoDecayHF2Prong* d, Int_t labD0=-1, TClonesArray* mcArray=0x0);
-  void GetImpParameter(AliAODTrack *track, AliAODEvent* aod, Double_t &d0, Double_t &d0err) const;
-  Bool_t TrackSelectedInLoop(AliAODTrack* track, AliAODRecoDecayHF2Prong *d, AliAODEvent *aod, Int_t ptbin, Int_t idDaughs[]) const;
+  void CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0=-1, TClonesArray* mcArray=0x0);
+  void FillSparsePlots(TClonesArray* arrayMC, AliAODRecoDecayHF2Prong *d, Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t type);
   Int_t CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const;
   Bool_t SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t option, Int_t idArrayV0[][2]) const;
-  Double_t PtWeig(Int_t ptbin, Double_t ptTrack) const;
 
   Int_t             fNPtBinsCorr;        // number of pt bins per correlations
   vector<Double_t>  fBinLimsCorr;        // limits of pt bins per correlations
@@ -96,7 +95,11 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   TH1F     *fNentries;                 //!histogram with number of events on output slot 2
   AliRDHFCutsD0toKpi *fCutsD0;         // Cuts for D0, output 3
   AliHFAssociatedTrackCuts *fCutsTracks;// Cuts for tracks and K0, output 7
+  AliHFCorrelator* fCorrelatorTr;      // Correlator for tracks
+  AliHFCorrelator* fCorrelatorKc;      // Correlator for charged K
+  AliHFCorrelator* fCorrelatorK0;      // Correlator for K0
   Bool_t    fReadMC;                   // flag for MC array: kTRUE = read it, kFALSE = do not read it
+  Bool_t    fMixing;                   // flag to enable also event mixing
   AliNormalizationCounter *fCounter;   //!AliNormalizationCounter on output slot 4
   Int_t     fNPtBins;                  // Number of pt bins
   Int_t     fFillOnlyD0D0bar;          // flag to fill mass histogram with D0/D0bar only (0 = fill with both, 1 = fill with D0 only, 2 = fill with D0bar only)
@@ -105,7 +108,7 @@ class AliAnalysisTaskSED0Correlations : public AliAnalysisTaskSE
   Bool_t    fIsRejectSDDClusters;      // flag to reject events with SDD clusters
   Bool_t    fFillGlobal;               // flag to fill global plots (in loops on tracks and V0 for each event)
 
-  ClassDef(AliAnalysisTaskSED0Correlations,1); // AliAnalysisTaskSE for D0->Kpi
+  ClassDef(AliAnalysisTaskSED0Correlations,2); // AliAnalysisTaskSE for D0->Kpi
 };
 
 #endif
index accb795..aef609c 100644 (file)
@@ -41,6 +41,17 @@ AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts():
 AliAnalysisCuts(),
 fESDTrackCuts(0),
 fPidObj(0),
+
+fPoolMaxNEvents(0), 
+fPoolMinNTracks(0), 
+fMinEventsToMix(0),
+fNzVtxBins(0), 
+fNzVtxBinsDim(0), 
+fZvtxBins(0), 
+fNCentBins(0), 
+fNCentBinsDim(0), 
+fCentBins(0),
+
 fNTrackCuts(0),
 fAODTrackCuts(0),
 fTrackCutsNames(0),
@@ -63,6 +74,17 @@ AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts(const char* name, const char*
 AliAnalysisCuts(name,title),
 fESDTrackCuts(0),
 fPidObj(0),
+
+fPoolMaxNEvents(0), 
+fPoolMinNTracks(0), 
+fMinEventsToMix(0),
+fNzVtxBins(0), 
+fNzVtxBinsDim(0), 
+fZvtxBins(0), 
+fNCentBins(0), 
+fNCentBinsDim(0), 
+fCentBins(0),
+
 fNTrackCuts(0),
 fAODTrackCuts(0),
 fTrackCutsNames(0),
@@ -81,6 +103,17 @@ AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts(const AliHFAssociatedTrackCut
 AliAnalysisCuts(source),
 fESDTrackCuts(source.fESDTrackCuts),
 fPidObj(source.fPidObj),
+
+fPoolMaxNEvents(source.fPoolMaxNEvents), 
+fPoolMinNTracks(source.fPoolMinNTracks), 
+fMinEventsToMix(source.fMinEventsToMix),
+fNzVtxBins(source.fNzVtxBins), 
+fNzVtxBinsDim(source.fNzVtxBinsDim), 
+fZvtxBins(source.fZvtxBins), 
+fNCentBins(source.fNCentBins), 
+fNCentBinsDim(source.fNCentBinsDim), 
+fCentBins(source.fCentBins),
+
 fNTrackCuts(source.fNTrackCuts),
 fAODTrackCuts(source.fAODTrackCuts),
 fTrackCutsNames(source.fTrackCutsNames),
@@ -127,6 +160,8 @@ AliHFAssociatedTrackCuts::~AliHFAssociatedTrackCuts()
 {
        if(fESDTrackCuts) {delete fESDTrackCuts; fESDTrackCuts = 0;}
        if(fPidObj) {delete fPidObj; fPidObj = 0;}
+       if(fZvtxBins) {delete[] fZvtxBins; fZvtxBins=0;} 
+       if(fCentBins) {delete[] fCentBins; fCentBins=0;}
        if(fAODTrackCuts) {delete[] fAODTrackCuts; fAODTrackCuts=0;}
        if(fTrackCutsNames) {delete[] fTrackCutsNames; fTrackCutsNames=0;}
        if(fAODvZeroCuts){delete[] fAODvZeroCuts; fAODvZeroCuts=0;}
@@ -140,29 +175,24 @@ Bool_t AliHFAssociatedTrackCuts::IsInAcceptance()
        return kFALSE;
 }
 //--------------------------------------------------------------------------
-Bool_t AliHFAssociatedTrackCuts::IsHadronSelected(AliAODTrack * track, AliAODVertex *vtx1, Double_t bz) 
+Bool_t AliHFAssociatedTrackCuts::IsHadronSelected(AliAODTrack * track)
 {
-       
-       
-       
-       Double_t d0 = -999.;
-       Double_t d0err = -998.;
-       Double_t VeryBig = 100;
-    
        AliESDtrack esdtrack(track);
-       Double_t d0z0[2],covd0z0[3];
+       if(!fESDTrackCuts->IsSelected(&esdtrack)) return kFALSE;
+       return kTRUE;
        
-       esdtrack.PropagateToDCA(vtx1,bz,VeryBig,d0z0,covd0z0);
+}
+
+//--------------------------------------------------------------------------
+Bool_t AliHFAssociatedTrackCuts::CheckHadronKinematic(Double_t pt, Double_t d0) 
+{
        
-       d0 = TMath::Abs(d0z0[0]);
-       d0err = TMath::Sqrt(covd0z0[0]);
        
-       if(!fESDTrackCuts->IsSelected(&esdtrack)) return kFALSE;
        
-       if(track->Pt() < fAODTrackCuts[0]) return kFALSE;
-       if(track->Pt() > fAODTrackCuts[1]) return kFALSE;
-       if(d0 < fAODTrackCuts[2])  return kFALSE;
-       if(d0 > fAODTrackCuts[3])  return kFALSE;
+       if(pt < fAODTrackCuts[0]) return kFALSE;
+       if(pt > fAODTrackCuts[1]) return kFALSE;
+       if(d0 < fAODTrackCuts[2]) return kFALSE;
+       if(d0 > fAODTrackCuts[3]) return kFALSE;
        
        return kTRUE;
 
@@ -273,17 +303,57 @@ Int_t AliHFAssociatedTrackCuts::IsMCpartFromHF(Int_t label, TClonesArray*mcArray
   if (isBeauty && isB) return 55;
   return 0;
 }
+
+//--------------------------------------------------------------------------
+Bool_t AliHFAssociatedTrackCuts::InvMassDstarRejection(AliAODRecoDecayHF2Prong* d, AliAODTrack *track, Int_t hypD0) const {
+       //
+       // Calculates invmass of track+D0 and rejects if compatible with D*
+       // (to remove pions from D*)
+       // 
+       Double_t nsigma = 3.;
+       
+       Double_t mD0, mD0bar;
+       d->InvMassD0(mD0,mD0bar);
+       
+       Double_t invmassDstar1 = 0, invmassDstar2 = 0; 
+       Double_t e1Pi = d->EProng(0,211), e2K = d->EProng(1,321); //hyp 1 (pi,K) - D0
+       Double_t e1K = d->EProng(0,321), e2Pi = d->EProng(1,211); //hyp 2 (K,pi) - D0bar
+       Double_t psum2 = (d->Px()+track->Px())*(d->Px()+track->Px())
+       +(d->Py()+track->Py())*(d->Py()+track->Py())
+       +(d->Pz()+track->Pz())*(d->Pz()+track->Pz());
+       
+       switch(hypD0) {
+               case 1:
+                       invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
+                       if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
+                       break;
+               case 2:
+                       invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
+                       if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
+                       break;
+               case 3:
+                       invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
+                       invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
+                       if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
+                       if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
+                       break;
+       }
+       
+       return kTRUE;
+}
 //________________________________________________________
 void AliHFAssociatedTrackCuts::SetAODTrackCuts(Float_t *cutsarray)
 {
-       
+       cout << "Check 1" << endl;
        if(!fAODTrackCuts) fAODTrackCuts = new Float_t[fNTrackCuts];
-       
+       cout << "Check 2" << endl;
        for(Int_t i =0; i<fNTrackCuts; i++){
-               
+               cout << "Check 2." << i << endl;
                fAODTrackCuts[i] = cutsarray[i];
        }
+       cout << "Check 3" << endl;
        SetTrackCutsNames();
+       cout << "Check 4" << endl;
        return;
 }
 //________________________________________________________
@@ -357,42 +427,30 @@ void AliHFAssociatedTrackCuts::PrintAll()
                cout << fvZeroCutsNames[k] <<  fAODvZeroCuts[k] << endl;
        }
        cout << " " << endl;
+       PrintPoolParameters();
+
 }
 
 //--------------------------------------------------------------------------
-Bool_t AliHFAssociatedTrackCuts::InvMassDstarRejection(AliAODRecoDecayHF2Prong* d, AliAODTrack *track, Int_t hypD0) const {
-  //
-  // Calculates invmass of track+D0 and rejects if compatible with D*
-  // (to remove pions from D*)
-  // 
-  Double_t nsigma = 3.;
-  Double_t mD0, mD0bar;
-  d->InvMassD0(mD0,mD0bar);
-
-  Double_t invmassDstar1 = 0, invmassDstar2 = 0; 
-  Double_t e1Pi = d->EProng(0,211), e2K = d->EProng(1,321); //hyp 1 (pi,K) - D0
-  Double_t e1K = d->EProng(0,321), e2Pi = d->EProng(1,211); //hyp 2 (K,pi) - D0bar
-  Double_t psum2 = (d->Px()+track->Px())*(d->Px()+track->Px())
-                  +(d->Py()+track->Py())*(d->Py()+track->Py())
-                  +(d->Pz()+track->Pz())*(d->Pz()+track->Pz());
-
-  switch(hypD0) {
-    case 1:
-      invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
-      if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
-      break;
-    case 2:
-      invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
-      if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
-      break;
-    case 3:
-      invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
-      invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
-      if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
-      if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
-      break;
-  }
-
-  return kTRUE;
+void AliHFAssociatedTrackCuts::PrintPoolParameters()
+{
+       printf("\nEvent Pool settings: \n \n");
+       
+       printf("Number of zVtx Bins: %d\n", fNzVtxBins);
+       printf("\nzVtx Bins:\n");
+       //Double_t zVtxbinLims[fNzVtxBins+1] = fNzVtxBins;
+       for(Int_t k=0; k<fNzVtxBins; k++){
+               printf("Bin %d..............................................: %.1f - %.1f cm\n ", k, fZvtxBins[k], fZvtxBins[k+1]);     
+       }
+       
+       printf("\nNumber of Centrality(multiplicity) Bins: %d\n", fNCentBins);
+       printf("\nCentrality(multiplicity) Bins:\n");
+       for(Int_t k=0; k<fNCentBins; k++){
+               printf("Bin %d..............................................: %.1f - %.1f\n ", k, fCentBins[k], fCentBins[k+1]);
+       }
+       
+       
+       
 }
+
+
index 76eddd9..edd5e1b 100644 (file)
@@ -51,25 +51,55 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
        Bool_t IsSelected(TList*  list) {if(list) return kTRUE; return kFALSE;};
        Bool_t IsSelected(TObject*  obj) {if(obj) return kTRUE; return kFALSE;};
        Bool_t IsInAcceptance();
-       Bool_t IsHadronSelected(AliAODTrack * track, AliAODVertex *vtx1, Double_t bz);
+       Bool_t IsHadronSelected(AliAODTrack * track);
+       Bool_t CheckHadronKinematic(Double_t pt, Double_t d0); 
        Bool_t CheckKaonCompatibility(AliAODTrack * track, Bool_t useMc, TClonesArray* mcArray, Int_t method=1);
        Bool_t IsKZeroSelected(AliAODv0 *vzero, AliAODVertex *vtx1);
        Int_t IsMCpartFromHF(Int_t label, TClonesArray*mcArray);
        Bool_t InvMassDstarRejection(AliAODRecoDecayHF2Prong* d, AliAODTrack *track, Int_t hypD0) const;        
        
+       //getters
+       Int_t GetMaxNEventsInPool() {return fPoolMaxNEvents;}
+       Int_t GetMinNTracksInPool() {return fPoolMinNTracks;}
+       Int_t GetMinEventsToMix(){return fMinEventsToMix;}
+       Int_t GetNZvtxPoolBins() {return fNzVtxBins;}
+       Double_t *GetZvtxPoolBins(){return fZvtxBins;}
+       Int_t GetNCentPoolBins() {return fNCentBins;}
+       Double_t *GetCentPoolBins(){return fCentBins;}
+       
+       
        
        void AddTrackCuts(const AliESDtrackCuts *cuts) {
                delete fESDTrackCuts; 
                fESDTrackCuts=new AliESDtrackCuts(*cuts); 
                return;
        }
+       //setters
+       //event pool settings
+       void SetMaxNEventsInPool(Int_t events){fPoolMaxNEvents=events;}
+       void SetMinNTracksInPool(Int_t tracks){fPoolMinNTracks=tracks;}
+       void SetMinEventsToMix(Int_t events){fMinEventsToMix=events;}
        
+       void SetNofPoolBins(Int_t Nzvtxbins, Int_t Ncentbins){
+               fNzVtxBins=Nzvtxbins;
+               fNzVtxBinsDim=Nzvtxbins+1;
+               
+           fNCentBins=Ncentbins;
+               fNCentBinsDim=Ncentbins+1;
+       }
+       
+       void SetPoolBins(Double_t *ZvtxBins, Double_t* CentBins){
+               fZvtxBins=ZvtxBins; 
+               fCentBins=CentBins;
+       }
+       //cut settings
        void SetAODTrackCuts(Float_t *cutsarray);
        void SetTrackCutsNames(/*TString *namearray*/);
        void SetAODvZeroCuts(Float_t *cutsarray);
        void SetvZeroCutsNames(/*TString *namearray*/);
        void SetPidHF(AliAODPidHF* pid) {fPidObj = pid; return;}
        virtual void PrintAll();
+       virtual void PrintPoolParameters();
        Int_t GetNVarsTrack(){return fNTrackCuts;}
        
        
@@ -81,6 +111,20 @@ private:
        //AliESDtrackCuts *fTrackCuts;
        AliESDtrackCuts *fESDTrackCuts; // track cut object
        AliAODPidHF * fPidObj;     /// PID object
+       
+       Int_t fPoolMaxNEvents; // set maximum number of events in the pool
+       Int_t fPoolMinNTracks; // se minimum number of tracks in the pool
+       Int_t fMinEventsToMix; // set the minimum number of events you wanna mix
+       
+       Int_t fNzVtxBins; // number of z vrtx bins
+       Int_t fNzVtxBinsDim; // number of z vrtx bins +1 : necessary to initialize correctly the array
+       Double_t* fZvtxBins; // [fNzVtxBinsDim]
+       
+       
+       Int_t fNCentBins; //number of centrality bins
+       Int_t fNCentBinsDim; //number of centrality bins bins +1 : necessary to initialize correctly the array
+       Double_t* fCentBins; // [fNCentBinsDim]
+       
        Int_t fNTrackCuts;     // array dimension
        Float_t* fAODTrackCuts;//[fNTrackCuts]
        TString * fTrackCutsNames;//[fNTrackCuts]
@@ -89,7 +133,7 @@ private:
        TString * fvZeroCutsNames;//[fNvZeroCuts]
        
        
-       ClassDef(AliHFAssociatedTrackCuts,1);
+       ClassDef(AliHFAssociatedTrackCuts,2);
 };
 
 
diff --git a/PWGHF/correlationHF/AliHFCorrelator.cxx b/PWGHF/correlationHF/AliHFCorrelator.cxx
new file mode 100644 (file)
index 0000000..d458a9c
--- /dev/null
@@ -0,0 +1,425 @@
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+//
+//
+//             Base class for Heavy Flavour Correlations Analysis
+//             Single Event and Mixed Event Analysis are implemented
+//
+//-----------------------------------------------------------------------
+//          
+//
+//                                                Author S.Bjelogrlic
+//                         Utrecht University 
+//                      sandro.bjelogrlic@cern.ch
+//
+//-----------------------------------------------------------------------
+
+/* $Id$ */
+
+#include <TParticle.h>
+#include <TVector3.h>
+#include <TChain.h>
+#include "TROOT.h"
+#include "AliHFCorrelator.h"
+#include "AliRDHFCutsDStartoKpipi.h"
+#include "AliHFAssociatedTrackCuts.h"
+#include "AliEventPoolManager.h"
+#include "AliReducedParticle.h"
+#include "AliCentrality.h"
+
+//_____________________________________________________
+AliHFCorrelator::AliHFCorrelator() :
+//
+// default constructor
+//
+TNamed(),
+fPoolMgr(0x0),         
+fPool(0x0),
+fhadcuts(0x0),
+fAODEvent(0x0),
+fAssociatedTracks(0x0),
+fmcArray(0x0),
+fReducedPart(0x0),
+fD0cand(0x0), 
+fhypD0(0), 
+
+fmixing(kFALSE),
+fmontecarlo(kFALSE),
+fsystem(kFALSE),
+fselect(0),
+fUseImpactParameter(0),
+fPIDmode(0),
+
+fNofTracks(0),
+fPoolContent(0),
+
+fPhiMin(0),
+fPhiMax(0),
+
+fPtTrigger(0),
+fPhiTrigger(0),
+fEtaTrigger(0),
+
+
+fDeltaPhi(0),
+fDeltaEta(0),
+fk0InvMass(0)
+
+{
+       // default constructor  
+}
+
+
+
+//_____________________________________________________
+AliHFCorrelator::AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb) :
+TNamed(name,"title"),
+fPoolMgr(0x0),         
+fPool(0x0),
+fhadcuts(0x0),
+fAODEvent(0x0),
+fAssociatedTracks(0x0),
+fmcArray(0x0),
+fReducedPart(0x0),
+fD0cand(0x0), 
+fhypD0(0),
+
+fmixing(kFALSE),
+fmontecarlo(kFALSE),
+fsystem(ppOrPbPb),
+fselect(0),
+fUseImpactParameter(0),
+fPIDmode(0),
+
+fNofTracks(0),
+fPoolContent(0),
+
+fPhiMin(0),
+fPhiMax(0),
+
+fPtTrigger(0),
+fPhiTrigger(0),
+fEtaTrigger(0),
+
+
+fDeltaPhi(0),
+fDeltaEta(0),
+fk0InvMass(0)
+{
+       fhadcuts = cuts; 
+}
+
+//_____________________________________________________
+AliHFCorrelator::~AliHFCorrelator() 
+{
+//
+// destructor
+//     
+       
+       if(fPoolMgr)  {delete fPoolMgr; fPoolMgr=0;}       
+       if(fPool) {delete fPool; fPool=0;}
+       if(fhadcuts) {delete fhadcuts; fhadcuts=0;}
+       if(fAODEvent) {delete fAODEvent; fAODEvent=0;}
+       if(fAssociatedTracks) {delete fAssociatedTracks; fAssociatedTracks=0;}
+       if(fmcArray) {delete fmcArray; fmcArray=0;}
+       if(fReducedPart) {delete fReducedPart; fReducedPart=0;}
+       if(fD0cand) {delete fD0cand; fD0cand=0;}
+       
+       
+       if(fNofTracks) fNofTracks = 0;
+       
+       if(fPhiMin) fPhiMin = 0;
+       if(fPhiMax) fPhiMax = 0;
+       
+       if(fPtTrigger) fPtTrigger=0;
+       if(fPhiTrigger) fPhiTrigger=0;
+       if(fEtaTrigger) fEtaTrigger=0;
+       
+       if(fDeltaPhi) fDeltaPhi=0;
+       if(fDeltaEta) fDeltaEta=0;
+       
+       if(fk0InvMass) fk0InvMass=0;
+       
+}
+//_____________________________________________________
+Bool_t AliHFCorrelator::DefineEventPool(){
+       // definition of the Pool Manager for Event Mixing
+       
+
+       Int_t MaxNofEvents = fhadcuts->GetMaxNEventsInPool();
+       Int_t MinNofTracks = fhadcuts->GetMinNTracksInPool();
+       Int_t NofCentBins = fhadcuts->GetNCentPoolBins();
+       Double_t * CentBins = fhadcuts->GetCentPoolBins();
+       Int_t NofZVrtxBins = fhadcuts->GetNZvtxPoolBins();
+       Double_t *ZVrtxBins = fhadcuts->GetZvtxPoolBins();
+               
+                       
+       fPoolMgr = new AliEventPoolManager(MaxNofEvents, MinNofTracks, NofCentBins, CentBins, NofZVrtxBins, ZVrtxBins);
+       if(!fPoolMgr) return kFALSE;
+       return kTRUE;
+}
+//_____________________________________________________
+Bool_t AliHFCorrelator::Initialize(){
+       
+       cout << "AliHFCorrelator::Initialize"<< endl;
+       if(!fAODEvent) cout << "No AOD event" << endl;
+       
+       AliCentrality *centralityObj = 0;
+       Int_t multiplicity = -1;
+       Double_t MultipOrCent = -1;
+       
+       // initialize the pool for event mixing
+       if(!fsystem){ // pp
+       multiplicity = fAODEvent->GetNTracks();
+               MultipOrCent = multiplicity; // convert from Int_t to Double_t
+       }
+       if(fsystem){ // PbPb
+               
+               centralityObj = fAODEvent->GetHeader()->GetCentralityP();
+               MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
+               AliInfo(Form("Centrality is %f", MultipOrCent));
+       }
+       
+       AliAODVertex *vtx = fAODEvent->GetPrimaryVertex();
+       Double_t zvertex = vtx->GetZ(); // zvertex
+       Double_t * CentBins = fhadcuts->GetCentPoolBins();
+       Double_t poolmin=CentBins[0];
+       Double_t poolmax=CentBins[fhadcuts->GetNCentPoolBins()];
+       cout << "pool limits in centr/multipl " << poolmin << " - " << poolmax << endl; 
+/*     for (Int_t iM=0; iM<nMultBins; iM++) {
+               std::vector<AliEventPool*> evp;
+               for (Int_t iZ=0; iZ<nZvtxBins; iZ++) {
+                       evp.push_back(new AliEventPool(depth, 
+                                                                                  multbin[iM], multbin[iM+1], 
+                                                                                  zvtxbin[iZ], zvtxbin[iZ+1] ));
+               }
+               fEvPool.push_back(evp);
+*/     
+//     if(TMath::Abs(zvertex)>=10 || multip>500 || multip == 0) {
+//             AliInfo(Form("Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,multip));
+//             return kFALSE;
+//     }
+       
+               if(TMath::Abs(zvertex)>=10 || MultipOrCent>poolmax || MultipOrCent < poolmin) {
+               if(!fsystem)AliInfo(Form("pp Event with Zvertex = %.2f cm and multiplicity = %.0f out of pool bounds, SKIPPING",zvertex,MultipOrCent));
+               else AliInfo(Form("PbPb Event with Zvertex = %.2f cm and centrality = %.1f  out of pool bounds, SKIPPING",zvertex,MultipOrCent));
+
+                       return kFALSE;
+               }
+       
+       fPool = fPoolMgr->GetEventPool(MultipOrCent, zvertex);
+       
+       if (!fPool){
+               AliInfo(Form("No pool found for multiplicity = %f, zVtx = %f cm", MultipOrCent, zvertex));
+           return kFALSE;
+       }
+       fPool->PrintInfo();
+       return kTRUE;
+}
+
+//_____________________________________________________
+Bool_t AliHFCorrelator::ProcessEventPool(){
+        // analysis on Mixed Events
+       cout << "AliHFCorrelator::ProcessEventPool"<< endl;
+               if(!fmixing) return kFALSE;
+               if(!fPool->IsReady()) return kFALSE;
+               if(fPool->GetCurrentNEvents()<fhadcuts->GetMinEventsToMix()) return kFALSE;
+               fPool->PrintInfo();
+               fPoolContent = fPool->GetCurrentNEvents();
+               
+               return kTRUE;
+       
+}
+
+//_____________________________________________________
+Bool_t AliHFCorrelator::ProcessAssociatedTracks(Int_t EventLoopIndex){
+
+       fAssociatedTracks = new TObjArray();
+       
+       if(!fmixing){ // analysis on Single Event
+               
+               
+               if(fselect==1 || fselect ==2)   fAssociatedTracks = AcceptAndReduceTracks(fAODEvent);
+               if(fselect==3) {fAssociatedTracks = AcceptAndReduceKZero(fAODEvent);}   
+               
+               
+       }
+       
+       if(fmixing) { // analysis on Mixed Events
+               
+                       
+               fAssociatedTracks = fPool->GetEvent(EventLoopIndex);
+                               
+                               
+                       
+
+       } // end if mixing
+       
+       if(!fAssociatedTracks) return kFALSE;
+       
+       fNofTracks = fAssociatedTracks->GetEntriesFast(); 
+               
+       return kTRUE;
+       
+}
+//_____________________________________________________
+Bool_t AliHFCorrelator::Correlate(Int_t loopindex){
+
+       if(loopindex >= fNofTracks) return kFALSE;
+       if(!fAssociatedTracks) return kFALSE;
+       
+       fReducedPart = (AliReducedParticle*)fAssociatedTracks->At(loopindex);
+       
+
+       fDeltaPhi = SetCorrectPhiRange(fPhiTrigger - fReducedPart->Phi());
+       
+       fDeltaEta = fEtaTrigger - fReducedPart->Eta();
+
+       return kTRUE;
+       
+}
+               
+//_____________________________________________________
+Bool_t AliHFCorrelator::PoolUpdate(){
+
+       if(!fmixing) return kFALSE;
+       if(!fPool) return kFALSE;
+       if(fmixing) { // update the pool for Event Mixing
+               if(fselect==1 || fselect==2)fPool->UpdatePool(AcceptAndReduceTracks(fAODEvent)); // updating the pool for hadrons
+               if(fselect==3) fPool->UpdatePool(AcceptAndReduceKZero(fAODEvent)); // updating the pool for K0s
+       }
+
+       return kTRUE;
+}
+               
+//_____________________________________________________
+Double_t AliHFCorrelator::SetCorrectPhiRange(Double_t phi){
+       Double_t pi = TMath::Pi();
+       
+       if(phi<fPhiMin) phi = phi + 2*pi;
+       if(phi>fPhiMax) phi = phi - 2*pi;
+       
+       return phi;
+}
+
+//_____________________________________________________
+TObjArray*  AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){
+       
+       Int_t nTracks = inputEvent->GetNTracks();
+       AliAODVertex * vtx = inputEvent->GetPrimaryVertex();
+       Double_t Bz = inputEvent->GetMagneticField();
+       
+       
+       TObjArray* tracksClone = new TObjArray;
+       tracksClone->SetOwner(kTRUE);
+       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
+               
+               AliAODTrack* track = inputEvent->GetTrack(iTrack);
+               if (!track) continue;
+               if(!fhadcuts->IsHadronSelected(track)) continue; // apply ESD level selections
+
+               Double_t pT = track->Pt();
+               
+               //compute impact parameter
+               Double_t d0z0[2],covd0z0[3];
+               Double_t d0=-999999.;
+               if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
+               else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
+               
+               if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
+               if(fUseImpactParameter==2) { // use impact parameter over resolution
+                 if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Abs(covd0z0[0]); 
+                 else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
+                       
+               }
+               
+
+               
+               
+               if(!fhadcuts->CheckHadronKinematic(pT,d0))continue; // apply kinematic cuts
+               Bool_t rejectsoftpi = kFALSE;
+               if(fD0cand) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0);
+               
+               
+               if(fselect ==2){        
+                       if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
+               }
+               
+               //AliVParticle * part = (AliVParticle*)track;
+               tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi));
+       }
+       
+       
+       return tracksClone;
+}
+
+//_____________________________________________________
+TObjArray*  AliHFCorrelator::AcceptAndReduceKZero(AliAODEvent* inputEvent){
+       
+       Int_t nOfVZeros = inputEvent->GetNumberOfV0s();
+       TObjArray* KZeroClone = new TObjArray;
+       AliAODVertex *vertex1 = (AliAODVertex*)inputEvent->GetPrimaryVertex();
+       Int_t v0label = -1;
+       Int_t pdgDgK0toPipi[2] = {211,211};
+       Double_t mPDGK0=0.497614;//TDatabasePDG::Instance()->GetParticle(310)->Mass();
+       const Int_t size = inputEvent->GetNumberOfV0s();
+       Int_t prevnegID[size];
+       Int_t prevposID[size];
+       for(Int_t iv0 =0; iv0< nOfVZeros; iv0++){// loop on all v0 candidates
+               AliAODv0 *v0 = (static_cast<AliAODEvent*> (inputEvent))->GetV0(iv0);
+               if(!v0) {cout << "is not a v0 at step " << iv0 << endl; continue;}
+               
+               if(!fhadcuts->IsKZeroSelected(v0,vertex1)) continue; // check if it is a k0
+           
+               // checks to avoid double counting
+               Int_t negID = -999;
+               Int_t posID = -998;
+               //Int_t a = 0;
+               prevnegID[iv0] = -997;
+               prevposID[iv0]= -996;
+               negID = v0->GetNegID();
+               posID = v0->GetPosID();
+               Bool_t DoubleCounting = kFALSE;
+               
+               for(Int_t k=0; k<iv0; k++){
+                       if(((negID==prevnegID[k])&&(posID==prevposID[k]))||((negID==prevposID[k])&&(posID==prevnegID[k]))){
+                               DoubleCounting = kTRUE;
+                               //a=k;
+                               break;
+                       }//end if
+               } // end for
+               
+               if(DoubleCounting) {cout << "SKIPPING DOUBLE COUNTING" << endl;continue;}
+               else{ prevposID[iv0] = posID; prevnegID[iv0] = negID;}
+               
+               if(fmontecarlo) v0label = v0->MatchToMC(310,fmcArray, 0, pdgDgK0toPipi); //match a K0 short
+               Double_t k0pt = v0->Pt();
+               Double_t k0eta = v0->Eta();
+               Double_t k0Phi = v0->Phi();
+           fk0InvMass = v0->MassK0Short();     
+               
+               //if (loopindex == 0) {
+               //      if(!plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectra"))->Fill(k0InvMass,k0pt); // spectra for all k0
+               //      if(plotassociation) ((TH2F*)fOutput->FindObject("KZeroSpectraifHF"))->Fill(k0InvMass,k0pt); // spectra for k0 in association with a D*
+               //}
+               // if there are more D* candidates per event, loopindex == 0 makes sure you fill the mass spectra only once!
+               
+               if(TMath::Abs(fk0InvMass-mPDGK0)>3*0.004) continue; // select candidates within 3 sigma
+               KZeroClone->Add(new AliReducedParticle(k0eta,k0Phi,k0pt,v0label));
+               
+       }
+       
+       return KZeroClone;
+}
diff --git a/PWGHF/correlationHF/AliHFCorrelator.h b/PWGHF/correlationHF/AliHFCorrelator.h
new file mode 100644 (file)
index 0000000..0744edb
--- /dev/null
@@ -0,0 +1,160 @@
+#ifndef AliHFCorrelator_H
+#define AliHFCorrelator_H
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//
+//             Base class for Heavy Flavour Correlations Analysis
+//             Single Event and Mixed Event Analysis are implemented
+//-----------------------------------------------------------------------
+//          
+//
+//                                                Author S.Bjelogrlic
+//                         Utrecht University 
+//                      sandro.bjelogrlic@cern.ch
+//
+//-----------------------------------------------------------------------
+
+/* $Id$ */
+
+#include "AliHFAssociatedTrackCuts.h"
+#include "AliEventPoolManager.h"
+#include "AliVParticle.h"
+#include "AliReducedParticle.h"
+
+
+class AliHFCorrelator : public TNamed
+{
+       
+ public:
+       
+       AliHFCorrelator();
+       AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t ppOrPbPb);
+       virtual ~AliHFCorrelator();
+       
+       //setters
+       void SetDeltaPhiInterval (Double_t min, Double_t max){
+               fPhiMin = min; fPhiMax = max;
+               if(TMath::Abs(fPhiMin-fPhiMax) != 2*TMath::Pi()) AliInfo("AliHFCorrelator::Warning: the delta phi interval is not set to 2 Pi");
+       }
+       void SetEventMixing(Bool_t mixON){fmixing=mixON;}
+       void SetTriggerParticleProperties(Double_t ptTrig, Double_t phiTrig, Double_t etaTrig)
+       {fPtTrigger = ptTrig; fPhiTrigger = phiTrig; fEtaTrigger = etaTrig;}
+       
+       
+       void SetAssociatedParticleType(Int_t type){fselect = type;}
+       void SetAODEvent(AliAODEvent* inputevent){fAODEvent = inputevent;}
+       void SetMCArray(TClonesArray* mcArray){fmcArray = mcArray;}
+       void SetUseMC(Bool_t useMC){fmontecarlo = useMC;}
+       void SetApplyDisplacementCut(Int_t applycut){fUseImpactParameter = applycut;}
+       void SetPIDmode(Int_t mode){fPIDmode = mode;}
+       
+       void SetD0Properties(AliAODRecoDecayHF2Prong* d, Int_t D0hyp)
+       {fD0cand = d; fhypD0 = D0hyp;}
+       
+       
+       
+    Bool_t DefineEventPool(); // Definition of the Event pool parameters
+       Bool_t Initialize(); // function that initlize everything for the analysis      
+       Bool_t ProcessEventPool(); // processes the event pool
+       Bool_t ProcessAssociatedTracks(Int_t EventLoopIndex); //
+       Bool_t Correlate(Int_t loopindex); // function that computes the correlations between the trigger particle and the track n. loopindex
+       Bool_t PoolUpdate();// updates the event pool
+       Double_t SetCorrectPhiRange(Double_t phi); // sets all the angles in the correct range
+       
+       //getters
+       AliEventPool* GetPool() {return fPool;}
+       TObjArray * GetTrackArray(){return fAssociatedTracks;}
+       AliHFAssociatedTrackCuts* GetSelectionCuts() {return fhadcuts;}
+       AliReducedParticle* GetAssociatedParticle() {return fReducedPart;}
+       
+       Int_t GetNofTracks(){return fNofTracks;}
+       Int_t GetNofEventsInPool(){return fPoolContent;}
+       
+       /*Double_t GetAssociatedTrackPt() {return fPtAssoc;} // pt of the associated track
+       Double_t GetAssociatedTrackPhi() {return fPhiAssoc;} // phi of the associated track
+    Double_t GetAssociatedTrackEta() {return fEtaAssoc;} // Eta of the associated track
+    Int_t GetAssociatedTrackLabel() {return fTrackLabel;} // tracklabel associated track
+       Int_t GetAssociatedTrackID() {return fTrackID;} // trackID associated track
+       Double_t GetAssociatedTrackImpPar() {return fTrackImpPar;} // impact parameter of the track
+       Bool_t GetAssociatedTrackSoftPiCompatibility() {return fIsSoftPiCand;}
+       */
+       Double_t GetDeltaPhi(){return fDeltaPhi;} // Delta Phi, needs to be called after the method correlate 
+       Double_t GetDeltaEta(){return fDeltaEta;} // Delta Eta
+       
+       Double_t GetAssociatedKZeroInvariantmass(){return fk0InvMass;}
+       
+       
+       
+       // methods to reduce the tracks to correlate with track selection cuts applied here
+       TObjArray*  AcceptAndReduceTracks(AliAODEvent* inputEvent); // selecting hadrons and kaons
+       TObjArray*  AcceptAndReduceKZero(AliAODEvent* inputEvent); // selecting kzeros
+       
+       
+ private:
+
+       AliHFCorrelator(const AliHFCorrelator& vtxr);
+       AliHFCorrelator& operator=(const AliHFCorrelator& vtxr );
+
+       AliEventPoolManager* fPoolMgr;         //! event pool manager
+       AliEventPool * fPool; //! Pool for event mixing
+       AliHFAssociatedTrackCuts* fhadcuts;//! hadron cuts
+       AliAODEvent * fAODEvent;//! AOD Event
+       TObjArray* fAssociatedTracks; // Array of associated tracks
+       TClonesArray* fmcArray; //mcarray
+       AliReducedParticle * fReducedPart; // reduced AOD particle;
+       AliAODRecoDecayHF2Prong* fD0cand; //D0 candidate
+       Int_t fhypD0; //hypothesis necessary for
+       
+       Bool_t fmixing;// switch for event mixing
+       Bool_t fmontecarlo; // switch for MonteCarlo
+       Bool_t fsystem; // select pp (kFALSE) or PbPb (kTRUE)
+       Int_t fselect; // 1 for hadrons, 2 for kaons, 3 for KZeros
+       Int_t fUseImpactParameter; // switch to use the impact parameter cut
+       Int_t fPIDmode; // set the PID mode for Kaon identification
+       
+       Int_t fNofTracks; // number of tracks in track array
+       Int_t fPoolContent; //  n of events in pool
+       
+       Double_t fPhiMin; // min for phi
+       Double_t fPhiMax; // max for phi
+       
+       Double_t fPtTrigger; // pt of the trigger D meson
+       Double_t fPhiTrigger; // phi of the trigger D meson
+       Double_t fEtaTrigger; // Eta of the trigger D meson
+       
+//  Double_t fPtAssoc; // pt of the trigger D meson
+//     Double_t fPhiAssoc; // phi of the trigger D meson
+//     Double_t fEtaAssoc; // Eta of the trigger D meson
+//     Int_t fTrackLabel; // tracklabel
+//     Int_t fTrackID; // trackID
+//     Double_t fTrackImpPar; // Impact parameter of the track in respect to the vertex
+//     Bool_t fIsSoftPiCand; // is kTRUE if the track  is compatible with a soft pion hypothesis (necessary for D0 analysis)
+       
+       Double_t fDeltaPhi; // delta phi between D meson and associated track
+       Double_t fDeltaEta; // delta eta between D meson and associated track
+       
+       Double_t fk0InvMass; // KZero invariant mass
+       
+       
+       ClassDef(AliHFCorrelator,1); // class for HF correlations       
+};
+
+
+
+
+
+#endif
index f420e4b..3df0170 100644 (file)
 #include "AliVParticle.h"
 #include "AliReducedParticle.h"
 
-AliReducedParticle::AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t McLabel, Int_t trackid) : 
+AliReducedParticle::AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t mcLabel, Int_t trackid, Double_t impPar, Bool_t checkSoftPi) : 
 fEta(eta), 
 fPhi(phi), 
 fpT(pt), 
-fMcLabel(McLabel), 
-fid(trackid)
+fMcLabel(mcLabel), 
+fid(trackid),
+fImpPar(impPar),
+fCheckSoftPi(checkSoftPi)
+{
+       //
+       // default constructor
+       //
+}
+
+AliReducedParticle::AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t McLabel) : 
+fEta(eta), 
+fPhi(phi), 
+fpT(pt), 
+fMcLabel(McLabel),
+fid(0),
+fImpPar(0.),
+fCheckSoftPi(kFALSE)
 {
        //
        // default constructor
@@ -48,5 +64,5 @@ AliReducedParticle::~AliReducedParticle() {
        // destructor
        //
        
-       Info("AliReducedParticle","Calling Destructor");
+       //Info("AliReducedParticle","Calling Destructor");
 }
index 9c525db..7edf23a 100644 (file)
 
 // class to get the reduced hadron candidate
 class AliReducedParticle : public AliVParticle
-//class AliReducedParticle : public AliAODTrack
 {
 public:
-    AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t McLabel, Int_t trackid);
+    AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t mcLabel, Int_t trackid, Double_t impPar, Bool_t checkSoftPi);
+    AliReducedParticle(Double_t eta, Double_t phi, Double_t pt, Int_t mcLabel);
 
     ~AliReducedParticle();
     
     // kinematics
        
     virtual Double_t Pt()         const { return fpT; }
-       virtual Double_t Phi()        const { return fPhi; }
-       virtual Double_t Eta()        const { return fEta; }
-       virtual Int_t   GetLabel()    const { return fMcLabel; }
-       virtual Int_t  GetID()          const{return fid;}
+    virtual Double_t Phi()        const { return fPhi; }
+    virtual Double_t Eta()        const { return fEta; }
+    virtual Int_t   GetLabel()    const { return fMcLabel; }
+    virtual Int_t  GetID()             const{return fid;}
+    virtual Double_t GetImpPar() const{return fImpPar;}
+    virtual Bool_t CheckSoftPi() const{return fCheckSoftPi;}
        
        // kinematics
     virtual Double_t Px() const { AliFatal("Not implemented"); return 0; }
@@ -89,10 +91,11 @@ private:
     Double_t fEta;      // eta
     Double_t fPhi;      // phi
     Double_t fpT;       // pT
-       Int_t fMcLabel; //mclabel
-       Int_t fid;
-
+    Int_t fMcLabel; //mclabel
+    Int_t fid; // track ID
+    Double_t fImpPar; // impact parameter
+    Bool_t fCheckSoftPi; // check if the track is compatible with a softpion from D*
     
-    ClassDef(AliReducedParticle, 1); // class which contains only quantities requires for this analysis to reduce memory consumption for event mixing
+    ClassDef(AliReducedParticle, 2); // class which contains only quantities requires for this analysis to reduce memory consumption for event mixing
 };
 #endif
index 39e68c9..cf55c31 100644 (file)
@@ -1,4 +1,4 @@
-AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Int_t system=0/*0=pp,1=PbPb*/, Int_t flagD0D0bar=0,
+AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Bool_t mixing=kFALSE, Int_t system=0/*0=pp,1=PbPb*/, Int_t flagD0D0bar=0,
                                       Float_t minC=0, Float_t maxC=0, TString finDirname="Output", TString cutsfilename="D0toKpiCuts.root", 
                                        TString cutsfilename2="AssocPartCuts.root", TString cutsD0name="D0toKpiCuts", TString cutsTrkname="AssociatedTrkCuts",                                          Bool_t flagAOD049=kFALSE, Int_t standardbins=1, Bool_t stdcuts=kFALSE)
 {
@@ -93,8 +93,8 @@ AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Int
   //Cuts for D0
   AliRDHFCutsD0toKpi* RDHFD0Corrs=new AliRDHFCutsD0toKpi();
 //RDHFD0Corrs->SetUsePhysicsSelection(kFALSE);
+  RDHFD0Corrs->SetLowPt(kFALSE); //low-pt special PID disabled
   if(stdcuts) {
-    RDHFD0Corrs->SetLowPt(kFALSE); //low-pt special PID disabled
     if(system==0) RDHFD0Corrs->SetStandardCutsPP2010();
     else {
       RDHFD0Corrs->SetStandardCutsPbPb2010();
@@ -145,6 +145,7 @@ AliAnalysisTaskSED0Correlations *AddTaskD0Correlations(Bool_t readMC=kFALSE, Int
   AliAnalysisTaskSED0Correlations *massD0Task = new AliAnalysisTaskSED0Correlations(taskname.Data(),RDHFD0Corrs,corrCuts);
   massD0Task->SetDebugLevel(2);
   massD0Task->SetReadMC(readMC);
+  massD0Task->SetEvMixing(mixing);
   massD0Task->SetFillOnlyD0D0bar(flagD0D0bar);
   massD0Task->SetSystem(system); //0=pp, 1=PbPb
 
index b657814..f9603b1 100644 (file)
- //DEFINITION OF A FEW CONSTANTS
-//----------------------------------------------------
-
-/* $Id$ */
-
-AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Int_t theMCon =5, Int_t mixing=4, Int_t trackselect =1)
-{
-
-  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
-  if (!mgr) {
-    ::Error("AliAnalysisTaskDStarCorrelations", "No analysis manager to connect to.");
-    return NULL;
-  } 
-
-  TFile* filecuts=new TFile("DStartoKpipiCuts_corr.root");
-  if(!filecuts->IsOpen()){
-    cout<<"Input file not found: exit"<<endl;
-    return;
-  }  
-         TFile* filecuts2=new TFile("AssocPartCuts.root");
-         if(!filecuts2->IsOpen()){
-                 cout<<"Input file2 not found: exit"<<endl;
-                 return;
-  }
-
-  AliRDHFCutsDStartoKpipi* RDHFDStartoKpipi=new AliRDHFCutsDStartoKpipi();
-  RDHFDStartoKpipi = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");
-  RDHFDStartoKpipi->SetName("DStartoKpipiCuts");
-       
-       
-       AliHFAssociatedTrackCuts* corrCuts=new AliHFAssociatedTrackCuts();
-       corrCuts = (AliHFAssociatedTrackCuts*)filecuts2->Get("AssociatedCuts");
-       corrCuts->SetName("AssociatedCuts");
-       corrCuts->PrintAll();
-
-  // mm let's see if everything is ok
-  if(!RDHFDStartoKpipi){
-    cout<<"Specific AliRDHFCuts not found"<<endl;
-    return;
-  } 
-
-  //CREATE THE TASK
-  printf("CREATE TASK \n");
-  // create the task
-  AliAnalysisTaskDStarCorrelations *task = new AliAnalysisTaskDStarCorrelations("AliAnalysisTaskDStarCorrelations",RDHFDStartoKpipi,corrCuts);
-       
-       // Setters
-
-       task->SetMonteCarlo(theMCon);
-       task->SetUseMixing(mixing);
-       task->SetCorrelator(trackselect) ;
-       //task->SetDebugLevel(0);
-       
-
-       if(trackselect == 1) Info("AliAnalysisTaskDStarCorrelations","Correlating D* with charged hadrons \n");
-       else if(trackselect == 2) Info("AliAnalysisTaskDStarCorrelations","Correlating D* with charged kaons \n");
-       else if(trackselect == 3) Info("AliAnalysisTaskDStarCorrelations","Correlating D* with reconstructed K0s \n");
-       else Fatal("AliAnalysisTaskDStarCorrelations","Nothing to correlate with!");
-       if(mixing) Info ("AliAnalysisTaskDStarCorrelations","Event Mixing Analysis\n");
-       if(!mixing) Info ("AliAnalysisTaskDStarCorrelations","Single Event Analysis \n");
-
-  // Create and connect containers for input/output
-       //TString dcavalue = " ";
-       if(!theMCon) TString contname = "Data";
-       if(theMCon) TString contname = "MonteCarlo";
-       if(trackselect ==1) TString particle = "Hadron";
-       if(trackselect ==2) TString particle = "Kaon";
-       if(trackselect ==3) TString particle = "KZero";
-       
-       TString cutname = "cuts" ;
-       TString cutname2 = "hadroncuts" ;
-   TString outputfile = AliAnalysisManager::GetCommonFileName();
-       TString counter = "NormCounter";
-   outputfile += ":PWGHF_D2H_";
-       if(!mixing) {
-               outputfile += "SE";
-               contname += "SE";
-               cutname += "SE";
-               cutname2 += "SE";
-               counter+= "SE";
-       }
-       if(mixing){
-               outputfile += "ME";
-               contname += "ME";
-               cutname += "ME";
-               cutname2 += "ME";
-               counter+= "ME";
-       }
-       outputfile += "Dphi_DStar";
-       outputfile += particle;
-       cutname += particle;
-       cutname2 += particle;
-       contname += particle;
-       counter+= particle;
-
-       
-       //cout << "Contname = " << contname << endl;
-  mgr->AddTask(task);
-  // ------ input data ------
-  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
-  
-  // ----- output data -----
-  
-  // output TH1I for event counting
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(cutname,AliRDHFCutsDStartoKpipi::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts
-  AliAnalysisDataContainer *coutputDstarNorm = mgr->CreateContainer(counter,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
-  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(cutname2,AliHFAssociatedTrackCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts
-  
-  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
-  mgr->ConnectOutput(task,1,coutput1);
-  mgr->ConnectOutput(task,2,coutput2);
-  mgr->ConnectOutput(task,3,coutputDstarNorm);
-  mgr->ConnectOutput(task,4,coutput4);
-
-  return task ;
-
-}
-
+ //DEFINITION OF A FEW CONSTANTS\r
+//----------------------------------------------------\r
+\r
+/* $Id$ */\r
+\r
+AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool_t theMCon, Bool_t mixing, Int_t trackselect =1, Int_t usedispl =0, TString DCutObjName = "DStartoKpipiCuts_corr.root", TString TrackCutObjName = "AssocPartCuts.root")\r
+{\r
+\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error(" AliAnalysisTaskDStarCorrelations", "No analysis manager to connect to.");\r
+    return NULL;\r
+  } \r
+\r
+  TFile* filecuts=new TFile(DCutObjName.Data());\r
+  if(!filecuts->IsOpen()){\r
+    cout<<"Input file not found: exit"<<endl;\r
+    return;\r
+  }  \r
+         TFile* filecuts2=new TFile(TrackCutObjName.Data());\r
+         if(!filecuts2->IsOpen()){\r
+                 cout<<"Input file2 not found: exit"<<endl;\r
+                 return;\r
+  }\r
+\r
+  AliRDHFCutsDStartoKpipi* RDHFDStartoKpipi=new AliRDHFCutsDStartoKpipi();\r
+  RDHFDStartoKpipi = (AliRDHFCutsDStartoKpipi*)filecuts->Get("DStartoKpipiCuts");\r
+  RDHFDStartoKpipi->SetName("DStartoKpipiCuts");\r
+       \r
+       \r
+       AliHFAssociatedTrackCuts* corrCuts=new AliHFAssociatedTrackCuts();\r
+       corrCuts = (AliHFAssociatedTrackCuts*)filecuts2->Get("AssociatedCuts");\r
+       corrCuts->SetName("AssociatedCuts");\r
+       corrCuts->PrintAll();\r
+\r
+  // mm let's see if everything is ok\r
+  if(!RDHFDStartoKpipi){\r
+    cout<<"Specific AliRDHFCuts not found"<<endl;\r
+    return;\r
+  } \r
+\r
+  //CREATE THE TASK\r
+  printf("CREATE TASK \n");\r
+  // create the task\r
+  AliAnalysisTaskDStarCorrelations *task = new AliAnalysisTaskDStarCorrelations(" AliAnalysisTaskDStarCorrelations",RDHFDStartoKpipi,corrCuts);\r
+       \r
+       // Setters\r
+\r
+       task->SetMonteCarlo(theMCon);\r
+       task->SetUseMixing(mixing);\r
+       task->SetCorrelator(trackselect) ;\r
+       task->SetUseDisplacement(usedispl);\r
+       task->SetRunPbPb(runOnPbPb);\r
+       //task->SetDebugLevel(0);\r
+       \r
+\r
+       if(trackselect == 1) Info(" AliAnalysisTaskDStarCorrelations","Correlating D* with charged hadrons \n");\r
+       else if(trackselect == 2) Info(" AliAnalysisTaskDStarCorrelations","Correlating D* with charged kaons \n");\r
+       else if(trackselect == 3) Info(" AliAnalysisTaskDStarCorrelations","Correlating D* with reconstructed K0s \n");\r
+       else Fatal(" AliAnalysisTaskDStarCorrelations","Nothing to correlate with!");\r
+       if(mixing) Info (" AliAnalysisTaskDStarCorrelations","Event Mixing Analysis\n");\r
+       if(!mixing) Info (" AliAnalysisTaskDStarCorrelations","Single Event Analysis \n");\r
+\r
+  // Create and connect containers for input/output\r
+       //TString dcavalue = " ";\r
+       if(!theMCon) TString contname = "Data";\r
+       if(theMCon) TString contname = "MonteCarlo";\r
+       if(trackselect ==1) TString particle = "Hadron";\r
+       if(trackselect ==2) TString particle = "Kaon";\r
+       if(trackselect ==3) TString particle = "KZero";\r
+       \r
+       TString cutname = "cuts" ;\r
+       TString cutname2 = "hadroncuts" ;\r
+   TString outputfile = AliAnalysisManager::GetCommonFileName();\r
+       TString counter = "NormCounter";\r
+   outputfile += ":PWGHF_D2H_";\r
+       if(!mixing) {\r
+               outputfile += "SE";\r
+               contname += "SE";\r
+               cutname += "SE";\r
+               cutname2 += "SE";\r
+               counter+= "SE";\r
+       }\r
+       if(mixing){\r
+               outputfile += "ME";\r
+               contname += "ME";\r
+               cutname += "ME";\r
+               cutname2 += "ME";\r
+               counter+= "ME";\r
+       }\r
+       outputfile += "Dphi_DStar";\r
+       outputfile += particle;\r
+       cutname += particle;\r
+       cutname2 += particle;\r
+       contname += particle;\r
+       counter+= particle;\r
+\r
+       \r
+       cout << "Contname = " << contname << endl;\r
+  mgr->AddTask(task);\r
+  // ------ input data ------\r
+  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();\r
+  \r
+  // ----- output data -----\r
+  \r
+  // output TH1I for event counting\r
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());\r
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(cutname,AliRDHFCutsDStartoKpipi::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts\r
+  AliAnalysisDataContainer *coutputDstarNorm = mgr->CreateContainer(counter,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());\r
+  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(cutname2,AliHFAssociatedTrackCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts\r
+  \r
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(task,1,coutput1);\r
+  mgr->ConnectOutput(task,2,coutput2);\r
+  mgr->ConnectOutput(task,3,coutputDstarNorm);\r
+  mgr->ConnectOutput(task,4,coutput4);\r
+\r
+  return task ;\r
+\r
+}\r
+\r
index 7e19461..f9c946f 100644 (file)
-#include <Riostream.h>
-#include <TFile.h>
-#include "AliHFAssociatedTrackCuts.h"
-#include <TClonesArray.h>
-#include <TParameter.h>
-//#include "AliAODPidHF.h"
-
-/* $Id$ */
-
-void makeInputHFCorrelation(){
-       
-       
-       AliHFAssociatedTrackCuts* HFCorrelationCuts=new AliHFAssociatedTrackCuts();
-       HFCorrelationCuts->SetName("AssociatedCuts");
-       HFCorrelationCuts->SetTitle("Cuts for associated track");
-       Float_t eta = 0.9;
-
-       //______________________________ set ESD track cuts
-       AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
-
-       
-       esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
-       //esdTrackCuts->SetDCAToVertex2D(kFALSE);
-       
-       esdTrackCuts->SetRequireTPCRefit(kTRUE);
-       esdTrackCuts->SetRequireITSRefit(kTRUE);
-       esdTrackCuts->SetMinNClustersTPC(80);
-       esdTrackCuts->SetMinNClustersITS(2); 
-       //esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); // set requirement on Pixels
-       // default is kBoth, otherwise kAny
-       esdTrackCuts->SetPtRange(0.3,1.e10);
-       esdTrackCuts->SetEtaRange(-eta,eta);
-       
-        
-       HFCorrelationCuts->AddTrackCuts(esdTrackCuts);
-       
-       //______________________________ set kinematics cuts for AOD track 
-       const int nofcuts = 4;
-       Float_t* trackcutsarray;
-       trackcutsarray=new Float_t[nofcuts];
-       trackcutsarray[0] = 0.3;//track min pt
-       trackcutsarray[1] = 100.;//track max pt
-       trackcutsarray[2] = 0.;//track min impact parameter
-       trackcutsarray[3] = 100.;//track max impact parameter
-       
-       HFCorrelationCuts->SetNVarsTrack(nofcuts);
-       HFCorrelationCuts->SetAODTrackCuts(trackcutsarray);
-       
-       
-       
-       //______________________________ set kinematics cuts for AOD v0 
-       
-       const int nofcuts2 = 7;
-       
-       Float_t* vzerocutsarray;
-       vzerocutsarray=new Float_t[nofcuts2];
-       vzerocutsarray[0] = 0.2; // max dca between two daugters (cm)
-       vzerocutsarray[1] = 2; //  max chi square
-       vzerocutsarray[2] = 2.; // min decay length (cm) 
-       vzerocutsarray[3] = 15; // max decay length (cm)
-       vzerocutsarray[4] = 0.2; // min opening angle between two daugters
-       vzerocutsarray[5] = 0; // min pt of k0 (GeV/c)
-       vzerocutsarray[6] = 0.9; // set eta acceptance
-
-       HFCorrelationCuts->SetNVarsVzero(nofcuts2);
-       HFCorrelationCuts->SetAODvZeroCuts(vzerocutsarray);
-       
-       //______________________________ set PID
-        
-       Int_t mode =1;
-       AliAODPidHF* pidObj=new AliAODPidHF();
-       pidObj->SetMatch(mode);
-       pidObj->SetSigma(0,2); // TPC
-       pidObj->SetSigma(3,3); // TOF
-       pidObj->SetTPC(kTRUE);
-       pidObj->SetTOF(kTRUE);
-       pidObj->SetCompat(kTRUE);
-       
-       HFCorrelationCuts->SetPidHF(pidObj);
-       
-       
-    //______________________________ save to *.root file
-       HFCorrelationCuts->PrintAll();
-       
-       TFile* fout=new TFile("AssocPartCuts.root","recreate");   //set this!! 
-       fout->cd();
-       HFCorrelationCuts->Write();
-       fout->Close();
-
-
-}
+#include <Riostream.h>\r
+#include <TFile.h>\r
+#include "AliHFAssociatedTrackCuts.h"\r
+#include <TClonesArray.h>\r
+#include <TParameter.h>\r
+#include "AliAODPidHF.h"\r
+\r
+/* $Id$ */\r
+\r
+void makeInputHFCorrelation(){\r
+       \r
+       AliHFAssociatedTrackCuts* HFCorrelationCuts=new AliHFAssociatedTrackCuts();\r
+       HFCorrelationCuts->SetName("AssociatedCuts");\r
+       HFCorrelationCuts->SetTitle("Cuts for associated track");\r
+       Float_t eta = 0.9;\r
+       //______________________________ set ESD track cuts\r
+       AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();\r
+\r
+       esdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
+       esdTrackCuts->SetDCAToVertex2D(kFALSE);\r
+       \r
+       esdTrackCuts->SetRequireTPCRefit(kTRUE);\r
+       esdTrackCuts->SetRequireITSRefit(kTRUE);\r
+       esdTrackCuts->SetMinNClustersTPC(80);\r
+       esdTrackCuts->SetMinNClustersITS(2); \r
+       //esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); // set requirement on Pixels\r
+       // default is kBoth, otherwise kAny\r
+       esdTrackCuts->SetPtRange(0.3,1.e10);\r
+       esdTrackCuts->SetEtaRange(-eta,eta);\r
+       HFCorrelationCuts->AddTrackCuts(esdTrackCuts);\r
+       //______________________________ set Pool for event mixing \r
+       \r
+       \r
+       \r
+       HFCorrelationCuts->SetMaxNEventsInPool(200);\r
+       HFCorrelationCuts->SetMinNTracksInPool(1000);\r
+       HFCorrelationCuts->SetMinEventsToMix(8);\r
+       \r
+       HFCorrelationCuts->SetNofPoolBins(3,5);\r
+       \r
+       //Double_t MBins[]={0,20,40,60,80,500};\r
+       Double_t MBins[]={0,1,2,3,4,500};\r
+       Double_t * MultiplicityBins = MBins;\r
+       \r
+       //Double_t ZBins[]={-10,-5,-2.5,2.5,5,10};\r
+       Double_t ZBins[]={-10,-2.5,2.5,10};\r
+       Double_t *ZVrtxBins = ZBins;\r
+       \r
+       HFCorrelationCuts->SetPoolBins(ZVrtxBins,MultiplicityBins);\r
+       cout << "Crash 1 " << endl;\r
+       //______________________________ set kinematics cuts for AOD track \r
+       const int nofcuts = 4;\r
+       Float_t* trackcutsarray;\r
+       trackcutsarray=new Float_t[nofcuts];\r
+       cout << "Crash 1.1 " << endl;\r
+       trackcutsarray[0] = 0.3;//track min pt\r
+       trackcutsarray[1] = 100.;//track max pt\r
+       trackcutsarray[2] = 0.;//track min impact parameter\r
+       trackcutsarray[3] = 100.;//track max impact parameter\r
+       cout << "Crash 1.2 " << endl;\r
+       HFCorrelationCuts->SetNVarsTrack(nofcuts);\r
+       cout << "Crash 1.3 " << endl;\r
+       HFCorrelationCuts->SetAODTrackCuts(trackcutsarray);\r
+       \r
+       cout << "Crash 2 " << endl;\r
+       //______________________________ set kinematics cuts for AOD v0 \r
+       \r
+       const int nofcuts2 = 7;\r
+       \r
+       Float_t* vzerocutsarray;\r
+       vzerocutsarray=new Float_t[nofcuts2];\r
+       vzerocutsarray[0] = 0.2; // max dca between two daugters (cm)\r
+       vzerocutsarray[1] = 2; //  max chi square\r
+       vzerocutsarray[2] = 2.; // min decay length (cm) \r
+       vzerocutsarray[3] = 15; // max decay length (cm)\r
+       vzerocutsarray[4] = 0.2; // min opening angle between two daugters\r
+       vzerocutsarray[5] = 0; // min pt of k0 (GeV/c)\r
+       vzerocutsarray[6] = 0.9; // set eta acceptance\r
+\r
+       HFCorrelationCuts->SetNVarsVzero(nofcuts2);\r
+       HFCorrelationCuts->SetAODvZeroCuts(vzerocutsarray);\r
+       //______________________________ set PID\r
+       cout << "Crash 3 " << endl;\r
+\r
+       Int_t mode =1;\r
+       AliAODPidHF* pidObj=new AliAODPidHF();\r
+       pidObj->SetMatch(mode);\r
+       pidObj->SetSigma(0,2); // TPC\r
+       pidObj->SetSigma(3,3); // TOF\r
+       pidObj->SetTPC(kTRUE);\r
+       pidObj->SetTOF(kTRUE);\r
+       pidObj->SetCompat(kTRUE);\r
+       \r
+       HFCorrelationCuts->SetPidHF(pidObj);\r
+       cout << "Crash 4 " << endl;\r
+\r
+    //______________________________ save to *.root file\r
+       HFCorrelationCuts->PrintAll();\r
+       \r
+       TFile* fout=new TFile("AssocPartCuts.root","recreate");   //set this!! \r
+       fout->cd();\r
+       HFCorrelationCuts->Write();\r
+       fout->Close();\r
+       cout << "Crash 5 " << endl;\r
+\r
+}\r