]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx
Corrected EINCLUDE
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskDStarCorrelations.cxx
index 899557f32d66419ead863dd9e9312cc8591f093b..37a755b51b81e7847e312dc9844ab39d8f8ae856 100644 (file)
-/**************************************************************************
- * 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 DStar - Hadron Correlations Analysis
-//
-//-----------------------------------------------------------------------
-//          
-//
-//                                                Author S.Bjelogrlic
-//                         Utrecht University 
-//                      sandro.bjelogrlic@cern.ch
-//
-//-----------------------------------------------------------------------
-
-/* $Id$ */
-
-#include <TDatabasePDG.h>
-#include <TParticle.h>
-#include <TVector3.h>
-#include <TChain.h>
-#include "TROOT.h"
-
-#include "AliAnalysisTaskDStarCorrelations.h"
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliHFAssociatedTrackCuts.h"
-#include "AliAODRecoDecay.h"
-#include "AliAODRecoCascadeHF.h"
-#include "AliAODRecoDecayHF2Prong.h"
-#include "AliAODPidHF.h"
-#include "AliEventPoolManager.h"
-#include "AliVParticle.h"
-#include "AliAnalysisManager.h"
-#include "AliAODInputHandler.h"
-#include "AliAODHandler.h"
-#include "AliESDtrack.h"
-#include "AliAODMCParticle.h"
-#include "AliNormalizationCounter.h"
-#include "AliReducedParticle.h"
-
-
-
-ClassImp(AliAnalysisTaskDStarCorrelations)
-
-
-//__________________________________________________________________________
-AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
-AliAnalysisTaskSE(),
-fhandler(0x0),
-fPoolMgr(0x0),      
-fmcArray(0x0),
-fCounter(0x0),
-fselect(0),
-fmontecarlo(kFALSE),
-fmixing(kFALSE),
-fEvents(0),
-fDebug(0),
-
-fOutput(0x0),            
-fCuts(0),
-fCuts2(0)
-{
-// default constructor 
-}
-
-//__________________________________________________________________________
-AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts) :
-AliAnalysisTaskSE(name),
-
-fhandler(0x0),
-fPoolMgr(0x0),   
-fmcArray(0x0),
-fCounter(0x0),
-fselect(0),
-fmontecarlo(kFALSE),
-fmixing(kFALSE),
-fEvents(0),
-fDebug(0),
-
-fOutput(0x0),                
-fCuts(0),
-fCuts2(AsscCuts)
-{
-       fCuts=cuts;
-       Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
-       DefineInput(0, TChain::Class());
-       DefineOutput(1,TList::Class()); // histos from data
-       DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
-       DefineOutput(3,AliNormalizationCounter::Class());   // normalization
-       DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my cuts
-}
-
-//__________________________________________________________________________
-
-AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
-       //
-       // destructor
-       //
-       
-       Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  
-       
-       if(fhandler) {delete fhandler; fhandler = 0;}
-       if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    
-       if(fmcArray) {delete fmcArray; fmcArray = 0;}
-       if(fCounter) {delete fCounter; fCounter = 0;}
-       if(fOutput) {delete fOutput; fOutput = 0;}
-       if(fCuts) {delete fCuts; fCuts = 0;}
-       if(fCuts2) {delete fCuts2; fCuts2=0;}
-
-}
-
-//___________________________________________________________
-void AliAnalysisTaskDStarCorrelations::Init(){
-       //
-       // Initialization
-       //
-       if(fDebug > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
-       
-       AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
-       // Post the D* cuts
-       PostData(2,copyfCuts);
-       
-       // Post the hadron cuts
-       PostData(4,fCuts2);
-       
-       return;
-}
-
-
-//_________________________________________________
-void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
-       Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
-       
-       //slot #1  
-       //OpenFile(0);
-       fOutput = new TList();
-       fOutput->SetOwner();
-       
-       // define histograms
-       DefineHistoForAnalysis();
-       
-       fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(3)->GetContainer()->GetName()));
-       fCounter->Init();
-       
-   
-       // definition of the Pool Manager for Event Mixing
-       
-       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);
-       
-       
-}
-//_________________________________________________
-void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
-
-       cout << " " << endl;
-       cout << "=================================================================================" << 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;
-       
-       if (!fInputEvent) {
-               Error("UserExec","NO EVENT FOUND!");
-               return;
-       }
-       
-        AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
-       Double_t pi = TMath::Pi();
-       
-       fEvents++; // event counter
-       ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
-               
-       fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-               
-       // 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));
-               return;
-       }
-       
-       
-       AliEventPool* pool = fPoolMgr->GetEventPool(multip, zvertex);
-       if (!pool) AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f cm", multip, zvertex));
-       
-       
-
-       // D* reconstruction
-       
-       TClonesArray *arrayDStartoD0pi=0;
-
-       
-       if(!aodEvent && AODEvent() && IsStandardAOD()) {
-               // In case there is an AOD handler writing a standard AOD, use the AOD 
-               // event in memory rather than the input (ESD) event.    
-               aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
-               // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
-               // have to taken from the AOD event hold by the AliAODExtension
-               AliAODHandler* aodHandler = (AliAODHandler*) 
-               ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
-               if(aodHandler->GetExtensions()) {
-                       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
-                       AliAODEvent *aodFromExt = ext->GetAOD();
-                       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
-               }
-       } else {
-               arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
-       }
-       
-       if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
-       
-       // initialize variables you will need for the D*
-       
-       Double_t ptDStar;//
-       Double_t phiDStar;//
-       Double_t etaDStar;//
-       Bool_t isInPeak, isInSideBand, isDStarMCtag;
-       Double_t invMassDZero;
-       Double_t deltainvMDStar;
-
-       
-       Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-       Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
-       
-       
-       
-       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++) {
-               isInPeak = kFALSE;
-               isInSideBand = kFALSE;
-               isDStarMCtag = kFALSE;
-               ptDStar = -123.4;
-               phiDStar = -999;
-               etaDStar = -56.;
-               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;    
-
-               if(fmontecarlo){
-                       // find associated MC particle for D* ->D0toKpi
-                       Int_t mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray);
-                       if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
-               }
-               
-               ptDStar = dstarD0pi->Pt();
-               phiDStar = dstarD0pi->Phi(); 
-               etaDStar = dstarD0pi->Eta();
-               
-               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;
-
-               
-               invMassDZero = dstarD0pi->InvMassD0();
-               ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
-               
-               deltainvMDStar = dstarD0pi->DeltaInvMass();
-               
-               
-               
-               //good candidates
-               if (TMath::Abs(invMassDZero-mPDGD0)<mD0Window){
-                       
-                       ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
-                       if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
-                               
-                               ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar);
-                               isInPeak = kTRUE;
-                       }
-               }// end if good candidates
-               
-               //sidebands
-               if (TMath::Abs(invMassDZero-mPDGD0)>1.3*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<4.*mD0Window ){
-                       ((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
-                       ((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero);
-                       
-                       if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
-                               ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar);
-                               isInSideBand = kTRUE;
-                       }
-                       
-               }//end if sidebands
-        
-
-               
-               if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
-               
-       
-               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();
-                               
-               if(!fmixing){ // analysis on Single Event
-                                               
-                       TObjArray* selectedtracks = new TObjArray();
-                       if(fselect==1 || fselect ==2)   selectedtracks = AcceptAndReduceTracks(aodEvent);
-                       if(fselect==3) {cout << " 2 "<< endl; 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;
-                               
-                               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
-                                       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);
-                                       
-                               }
-                               cout << "fill6" << endl;
-                               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++;
-                               }
-                               cout << "fill7" << endl;
-                       } // 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){ cout << "check 1" << endl;
-                                       ((TH2I*)fOutput->FindObject("EventMixingCheck"))->Fill(multbinflag,zvtxflag); 
-                                       cout << "filling" << endl;}
-
-                               TObjArray* mixedtracks = new TObjArray();
-                               
-                               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
-                       
-               } // end ME analysis
-               
-       }// end loop on D* candidates           
-       
-       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;
-       
-       PostData(1,fOutput); // set the outputs
-       PostData(3,fCounter); // set the outputs
-               
-} //end the exec
-
-//________________________________________ terminate ___________________________
-void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
-{    
-       // The Terminate() function is the last function to be called during
-       // a query. It always runs on the client, it can be used to present
-       // the results graphically or save the results to file.
-       
-       AliAnalysisTaskSE::Terminate();
-       
-       fOutput = dynamic_cast<TList*> (GetOutputData(1));
-       if (!fOutput) {     
-               printf("ERROR: fOutput not available\n");
-               return;
-       }
-
-       return;
-}
-
-
-//_____________________________________________________
-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();
-       Int_t nbinscorr = 32;
-       Double_t lowcorrbin = -0.5*Pi - Pi/32; // shift the bin by half the width so that at 0 is it the bin center
-       Double_t upcorrbin = 1.5*Pi - Pi/32;    
-       
-       // ========================= histograms for both Data and MonteCarlo
-       
-       
-       TH1D * NofEvents = new TH1D("NofEvents","NofEvents",1,0,1);
-       fOutput->Add(NofEvents);
-               
-       TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
-       fOutput->Add(D0InvMass);
-       
-       TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);
-       fOutput->Add(D0InvMassinSB);
-       
-       TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
-       fOutput->Add(DeltaInvMass);
-       
-       TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
-       fOutput->Add(bkgDeltaInvMass);
-       
-       TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",30,0,30);
-       fOutput->Add(RecoPtDStar);
-       
-       TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",30,0,30);
-       fOutput->Add(RecoPtBkg);
-       
-       TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
-       if(fselect==3) fOutput->Add(KZeroSpectra);
-       
-       TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);
-       if(fselect==3) fOutput->Add(KZeroSpectraifHF);
-       
-       TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","NofTracksInPeak",500,0.5,500.5);
-       fOutput->Add(NofTracksInPeak);
-       
-       TH1D * NofTracksInSB = new TH1D("NofTracksInSB","NofTracksInSB",500,0.5,500.5);
-       fOutput->Add(NofTracksInSB);
-       
-       TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
-       if(fmixing) fOutput->Add(EventMixingCheck);
-       
-
-       
-       TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
-       MCSources->GetXaxis()->SetBinLabel(1,"All ");
-       MCSources->GetXaxis()->SetBinLabel(2," from Heavy flavour");
-       MCSources->GetXaxis()->SetBinLabel(3," from c->D");
-       MCSources->GetXaxis()->SetBinLabel(4," from b->D");
-       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);
-       
-       TH1F * PhiSideBand = new TH1F("PhiSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi);
-       fOutput->Add(PhiSideBand);
-       
-       TH1F * PhiPart = new TH1F("PhiPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi);
-       fOutput->Add(PhiPart);
-
-
-       //correlations histograms
-       TString histoname1 = "DPhiDStar";
-       if(fselect==1) histoname1 += "Hadron";
-       if(fselect==2) histoname1 += "Kaon";
-       if(fselect==3) histoname1 += "KZero";
-
-       
-       TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       //side band background histograms
-       TString histoname2 = "bkg";
-       histoname2 += histoname1;
-       TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       
-       fOutput->Add(DPhiDStar);
-
-       if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
-       
-       fOutput->Add(bkgDPhiDStar);
-
-       if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
-
-
-       // ========================= histos for analysis on MC
-       TString histoname3 = "MCTag";
-       histoname3 += histoname1;
-       TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-        
-       TString histoname44 = "CharmDOrigin";
-       histoname44 += histoname1;
-       histoname44 += "MC";
-       
-       TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       
-       TString histoname54 = "BeautyDOrigin";
-       histoname54 += histoname1;
-       histoname54 += "MC";
-       TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       TString histoname55 = "BeautyBOrigin";
-       histoname55 += histoname1;
-       histoname55 += "MC";
-       TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-       
-       if(fmontecarlo){
-       
-       fOutput->Add(MCTagDPhiDStar);
-       fOutput->Add(CharmDOriginDPhiDStar);
-       fOutput->Add(BeautyDOriginDPhiDStar);
-       fOutput->Add(BeautyBOriginDPhiDStar);
-       
-       }
-       
-
-       
-}
-
-
-//____________________________  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;
-       cout << "CRASH CHECK 1 " << endl;
-               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);
-       cout << "CRASH CHECK 2 " << endl;
-       
-       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;
-       cout << "CRASH CHECK bkg 1 " << endl;
-               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);
-       cout << "CRASH CHECK bkg 2 " << endl;
-       
-       return;
-       
-       
-}
-
-//____________________________  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;
-cout << "fill1" << endl;
-       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);
-
-
-
-((TH1F*)fOutput->FindObject("MCSources"))->Fill(0);
-cout << "fill2" << endl;
-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);
-       
-       
-       ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       ((TH1F*)fOutput->FindObject("MCSources"))->Fill(2);
-       }
-cout << "fill3" << endl;
-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==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(3);
-       }
-cout << "fill4" << endl;
-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==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(4);
-       }
-       return;
-}
-
-
-
-
-
+/**************************************************************************\r
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+//\r
+//\r
+//             Base class for DStar - Hadron Correlations Analysis\r
+//\r
+//-----------------------------------------------------------------------\r
+//          \r
+//\r
+//                                                Author S.Bjelogrlic\r
+//                         Utrecht University \r
+//                      sandro.bjelogrlic@cern.ch\r
+//\r
+//-----------------------------------------------------------------------\r
+\r
+/* $Id$ */\r
+\r
+//#include <TDatabasePDG.h>\r
+#include <TParticle.h>\r
+#include <TVector3.h>\r
+#include <TChain.h>\r
+#include "TROOT.h"\r
+\r
+#include "AliAnalysisTaskDStarCorrelations.h"\r
+#include "AliRDHFCutsDStartoKpipi.h"\r
+#include "AliHFAssociatedTrackCuts.h"\r
+#include "AliAODRecoDecay.h"\r
+#include "AliAODRecoCascadeHF.h"\r
+#include "AliAODRecoDecayHF2Prong.h"\r
+#include "AliAODPidHF.h"\r
+#include "AliVParticle.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliAODHandler.h"\r
+#include "AliESDtrack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliNormalizationCounter.h"\r
+#include "AliReducedParticle.h"\r
+#include "AliHFCorrelator.h"\r
+#include "AliAODMCHeader.h"\r
+#include "AliEventPoolManager.h"\r
+#include "AliVertexingHFUtils.h"\r
+\r
+using std::cout;\r
+using std::endl;\r
+\r
+\r
+ClassImp(AliAnalysisTaskDStarCorrelations)\r
+\r
+\r
+//__________________________________________________________________________\r
+AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :\r
+AliAnalysisTaskSE(),\r
+fhandler(0x0),\r
+fmcArray(0x0),\r
+fCounter(0x0),\r
+fCorrelator(0x0),\r
+fselect(0),\r
+fmontecarlo(kFALSE),\r
+fmixing(kFALSE),\r
+fFullmode(kFALSE),\r
+fSystem(pp),\r
+fEfficiencyVariable(kNone),\r
+fReco(kTRUE),\r
+fUseEfficiencyCorrection(kFALSE),\r
+fUseDmesonEfficiencyCorrection(kFALSE),\r
+fUseCentrality(kFALSE),\r
+fUseHadronicChannelAtKineLevel(kFALSE),\r
+fPhiBins(32),\r
+fEvents(0),\r
+fDebugLevel(0),\r
+fDisplacement(0),\r
+fDim(0),\r
+fNofPtBins(0),\r
+fMaxEtaDStar(0.9),\r
+fDMesonSigmas(0),\r
+\r
+fD0Window(0),\r
+\r
+fOutput(0x0),\r
+fOutputMC(0x0),\r
+fCuts(0),\r
+fCuts2(0),\r
+fUtils(0),\r
+fTracklets(0),\r
+fDeffMapvsPt(0),\r
+fDeffMapvsPtvsMult(0),\r
+fDeffMapvsPtvsEta(0)\r
+{\r
+    SetDim();\r
+    // default constructor\r
+}\r
+\r
+//__________________________________________________________________________\r
+AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :\r
+AliAnalysisTaskSE(name),\r
+\r
+fhandler(0x0),\r
+fmcArray(0x0),\r
+fCounter(0x0),\r
+fCorrelator(0x0),\r
+fselect(0),\r
+fmontecarlo(kFALSE),\r
+fmixing(kFALSE),\r
+fFullmode(mode),\r
+fSystem(syst),\r
+fEfficiencyVariable(kNone),\r
+fReco(kTRUE),\r
+fUseEfficiencyCorrection(kFALSE),\r
+fUseDmesonEfficiencyCorrection(kFALSE),\r
+fUseCentrality(kFALSE),\r
+fUseHadronicChannelAtKineLevel(kFALSE),\r
+fPhiBins(32),\r
+fEvents(0),\r
+fDebugLevel(0),\r
+fDisplacement(0),\r
+fDim(0),\r
+fNofPtBins(0),\r
+fMaxEtaDStar(0.9),\r
+fDMesonSigmas(0),\r
+fD0Window(0),\r
+\r
+fOutput(0x0),\r
+fOutputMC(0x0),\r
+fCuts(0),\r
+fCuts2(AsscCuts),\r
+fUtils(0),\r
+fTracklets(0),\r
+fDeffMapvsPt(0),\r
+fDeffMapvsPtvsMult(0),\r
+fDeffMapvsPtvsEta(0)\r
+{\r
+     Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");\r
+  \r
+    SetDim();\r
+    if(fSystem == AA)  fUseCentrality = kTRUE; else fUseCentrality = kFALSE;\r
+  \r
+    fCuts=cuts;\r
+    fNofPtBins= fCuts->GetNPtBins();\r
+    //cout << "Enlarging the DZero window " << endl;\r
+    EnlargeDZeroMassWindow();\r
+   // cout << "Done" << endl;\r
+    \r
+   \r
+  DefineInput(0, TChain::Class());\r
+  \r
+  DefineOutput(1,TList::Class()); // histos from data and MC\r
+  DefineOutput(2,TList::Class()); // histos from MC\r
+  DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts\r
+  DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts\r
+  DefineOutput(5,AliNormalizationCounter::Class());   // normalization\r
+}\r
+\r
+//__________________________________________________________________________\r
+\r
+AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {\r
+  //\r
+       // destructor\r
+       //\r
+       \r
+       Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");  \r
+       \r
+       if(fhandler) {delete fhandler; fhandler = 0;}\r
+       //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}    \r
+       if(fmcArray) {delete fmcArray; fmcArray = 0;}\r
+       if(fCounter) {delete fCounter; fCounter = 0;}\r
+       if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}\r
+       if(fOutput) {delete fOutput; fOutput = 0;}\r
+       if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}\r
+       if(fCuts) {delete fCuts; fCuts = 0;}\r
+       if(fCuts2) {delete fCuts2; fCuts2=0;}\r
+    if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}\r
+    if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}\r
+    if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}\r
+\r
+}\r
+\r
+//___________________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::Init(){\r
+       //\r
+       // Initialization\r
+       //\r
+       if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");\r
+       \r
+       AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);\r
+       \r
+       \r
+    \r
+       // Post the D* cuts\r
+       PostData(3,copyfCuts);\r
+       \r
+       // Post the hadron cuts\r
+       PostData(4,fCuts2);\r
+    \r
+       return;\r
+}\r
+\r
+\r
+//_________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){\r
+       Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
+       \r
+       //slot #1\r
+       //OpenFile(0);\r
+       fOutput = new TList();\r
+       fOutput->SetOwner();\r
+       \r
+       fOutputMC = new TList();\r
+       fOutputMC->SetOwner();\r
+       \r
+       // define histograms\r
+       DefineHistoForAnalysis();\r
+    DefineThNSparseForAnalysis();\r
+    \r
+\r
+       fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));\r
+       fCounter->Init();\r
+       \r
+    Double_t Pi = TMath::Pi();\r
+       fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb\r
+       fCorrelator->SetDeltaPhiInterval(  -0.5*Pi, 1.5*Pi); // set correct phi interval\r
+       //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval\r
+       fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On\r
+       fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros\r
+       fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut\r
+       fCorrelator->SetUseMC(fmontecarlo);\r
+       fCorrelator->SetUseReco(fReco);\r
+    // fCorrelator->SetKinkRemoval(kTRUE);\r
+       Bool_t pooldef = fCorrelator->DefineEventPool();\r
+       \r
+       if(!pooldef) AliInfo("Warning:: Event pool not defined properly");\r
+    \r
+    fUtils = new AliAnalysisUtils();\r
+    \r
+    \r
+       \r
+       PostData(1,fOutput); // set the outputs\r
+       PostData(2,fOutputMC); // set the outputs\r
+       PostData(5,fCounter); // set the outputs\r
+}\r
+\r
+//_________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){\r
+  \r
+  \r
+  if(fDebugLevel){\r
+    \r
+    if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;\r
+    if(!fReco) std::cout << "USING MC TRUTH" << std::endl;\r
+    std::cout << " " << std::endl;\r
+    std::cout << "=================================================================================" << std::endl;\r
+    if(!fmixing){\r
+      if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;\r
+      if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;\r
+      if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;\r
+    }\r
+    if(fmixing){\r
+      if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;\r
+      if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;\r
+      if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;\r
+    }\r
+    \r
+  }// end if debug\r
+  \r
+  \r
+  \r
+  if (!fInputEvent) {\r
+    Error("UserExec","NO EVENT FOUND!");\r
+               return;\r
+  }\r
+  \r
+  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
+  if(!aodEvent){\r
+    AliError("AOD event not found!");\r
+    return;\r
+  }\r
+  \r
+    fTracklets = aodEvent->GetTracklets();\r
+    \r
+    fEvents++; // event counter\r
+    ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);\r
+    \r
+    fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);\r
+    \r
+    // load MC array\r
+    fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
+    if(fmontecarlo && !fmcArray){\r
+      AliError("Array of MC particles not found");\r
+      return;\r
+    }\r
+    \r
+    \r
+    \r
+    \r
+    // ********************************************** START EVENT SELECTION ****************************************************\r
+    \r
+    Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);\r
+    \r
+    if(!isEvSel) {\r
+      \r
+      if(fCuts->IsEventRejectedDueToPileupSPD()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);\r
+      if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);\r
+      if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);\r
+      if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);\r
+      if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);\r
+      if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);\r
+      if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);\r
+      \r
+      return;\r
+    }\r
+    \r
+    // added event selection for pA\r
+    \r
+    if(fSystem == pA){\r
+      \r
+      if(fUtils->IsFirstEventInChunk(aodEvent)) {\r
+       AliInfo("Rejecting the event - first in the chunk");\r
+       ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);\r
+       return;\r
+        }\r
+      if(!fUtils->IsVertexSelected2013pA(aodEvent)) {\r
+       ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);\r
+       AliInfo("Rejecting the event - bad vertex");\r
+       return;\r
+      }\r
+    }\r
+    // ******************************** END event selections **************************************************\r
+    \r
+    AliCentrality *centralityObj = 0;\r
+    Double_t MultipOrCent = -1;\r
+    \r
+    if(fUseCentrality){\r
+      /* if(fSystem == AA ){   */      centralityObj = aodEvent->GetHeader()->GetCentralityP();\r
+      MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
+      //AliInfo(Form("Centrality is %f", MultipOrCent));\r
+    }\r
+    \r
+    if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.);\r
+    \r
+        \r
+    fCorrelator->SetAODEvent(aodEvent); // set the event to be processed\r
+    \r
+    ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);\r
+    \r
+    Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings\r
+       if(!correlatorON) {\r
+         AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");\r
+         return;\r
+       }\r
+       \r
+       if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);\r
+               \r
+       // check the event type\r
+       // load MC header\r
+       if(fmontecarlo){\r
+         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));\r
+         if (fmontecarlo && !mcHeader) {\r
+           AliError("Could not find MC Header in AOD");\r
+           return;\r
+         }\r
+         \r
+         Bool_t isMCeventgood = kFALSE;\r
+         \r
+         \r
+         Int_t eventType = mcHeader->GetEventType();\r
+         Int_t NMCevents = fCuts2->GetNofMCEventType();\r
+         \r
+         for(Int_t k=0; k<NMCevents; k++){\r
+           Int_t * MCEventType = fCuts2->GetMCEventType();\r
+           \r
+           if(eventType == MCEventType[k]) isMCeventgood= kTRUE;\r
+           ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);\r
+         }\r
+         \r
+         if(NMCevents && !isMCeventgood){\r
+            if(fDebugLevel)    std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;\r
+           return;\r
+         }\r
+         \r
+       } // end if montecarlo\r
+       \r
+       \r
+       // checks on vertex and multiplicity of the event\r
+       AliAODVertex *vtx = aodEvent->GetPrimaryVertex();\r
+       Double_t zVtxPosition = vtx->GetZ(); // zvertex\r
+       \r
+       \r
+       if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);\r
+       \r
+       \r
+       \r
+       // D* reconstruction\r
+       TClonesArray *arrayDStartoD0pi=0;\r
+       if(!aodEvent && AODEvent() && IsStandardAOD()) {\r
+         // In case there is an AOD handler writing a standard AOD, use the AOD\r
+         // event in memory rather than the input (ESD) event.\r
+         aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r
+         // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
+         // have to taken from the AOD event hold by the AliAODExtension\r
+         AliAODHandler* aodHandler = (AliAODHandler*)\r
+           ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
+         if(aodHandler->GetExtensions()) {\r
+           AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
+           AliAODEvent *aodFromExt = ext->GetAOD();\r
+           arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r
+         }\r
+       } else {\r
+         arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");\r
+       }\r
+       \r
+       if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
+       \r
+    \r
+    // get the poolbin\r
+    \r
+    Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition);\r
+  \r
+    \r
+       // initialize variables you will need for the D*\r
+       \r
+       Double_t ptDStar;//\r
+       Double_t phiDStar;//\r
+       Double_t etaDStar;//\r
+       Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag;\r
+       Double_t invMassDZero;\r
+       Double_t deltainvMDStar;\r
+    \r
+       \r
+       Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+       Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();\r
+       \r
+       \r
+       //MC tagging for DStar\r
+       //D* and D0 prongs needed to MatchToMC method\r
+       Int_t pdgDgDStartoD0pi[2]={421,211};\r
+       Int_t pdgDgD0toKpi[2]={321,211};\r
+       \r
+       Bool_t isDStarCand = kFALSE;\r
+    Bool_t isDfromB = kFALSE;\r
+       Bool_t isEventMixingFilledPeak = kFALSE;\r
+       Bool_t isEventMixingFilledSB = kFALSE;\r
+    Bool_t EventHasDStarCandidate = kFALSE;\r
+    Bool_t EventHasDZeroSideBandCandidate = kFALSE;\r
+    Bool_t EventHasDStarSideBandCandidate = kFALSE;\r
+       //loop on D* candidates\r
+       \r
+       Int_t looponDCands = 0;\r
+       if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast();\r
+       if(!fReco) looponDCands = fmcArray->GetEntriesFast();\r
+       \r
+       Int_t nOfDStarCandidates = 0;\r
+       Int_t nOfSBCandidates = 0;\r
+       \r
+       Double_t DmesonEfficiency = 1.;\r
+       Double_t DmesonWeight = 1.;\r
+    Double_t efficiencyvariable = -999;\r
+    \r
+       \r
+       \r
+       for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {\r
+         isInPeak = kFALSE;\r
+         isInDStarSideBand = kFALSE;\r
+         isInDZeroSideBand = kFALSE;\r
+         isDStarMCtag = kFALSE;\r
+         isDfromB = kFALSE;\r
+         ptDStar = -123.4;\r
+         phiDStar = -999;\r
+         etaDStar = -56.;\r
+         invMassDZero = - 999;\r
+         deltainvMDStar = -998;\r
+         AliAODRecoCascadeHF* dstarD0pi;\r
+         AliAODRecoDecayHF2Prong* theD0particle;\r
+         AliAODMCParticle* DStarMC=0x0;\r
+      Short_t daughtercharge = -2;\r
+         Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info\r
+         Int_t trackiddaugh1 = -1;\r
+         Int_t trackidsoftPi = -1;\r
+         \r
+         // start the if reconstructed candidates from here ************************************************\r
+         \r
+         if(fReco){//// if reconstruction is applied\r
+           dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);\r
+           if(!dstarD0pi->GetSecondaryVtx()) continue;\r
+           theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();\r
+           if (!theD0particle) continue;\r
+            \r
+           \r
+           // track quality cuts\r
+           Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks\r
+           // region of interest + topological cuts + PID\r
+           Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected\r
+            \r
+            //apply track selections\r
+           if(!isTkSelected) continue;\r
+           if(!isSelected) continue;\r
+        if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;\r
+          \r
+          \r
+          ptDStar = dstarD0pi->Pt();\r
+          phiDStar = dstarD0pi->Phi();\r
+          etaDStar = dstarD0pi->Eta();\r
+          if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
+          if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr)  efficiencyvariable = MultipOrCent;\r
+          if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;\r
+          if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();\r
+          if(fEfficiencyVariable == kNone) efficiencyvariable = 0;\r
+          \r
+        // get the D meson efficiency\r
+        DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);\r
+            \r
+         \r
+\r
+           if(fUseDmesonEfficiencyCorrection){\r
+             if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;\r
+             else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI\r
+               if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value\r
+               else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low\r
+             }\r
+           }\r
+            else DmesonWeight = 1.; \r
+            \r
+            // continue;\r
+          \r
+           Int_t mcLabelDStar = -999;\r
+          if(fmontecarlo){\r
+             // find associated MC particle for D* ->D0toKpi\r
+             mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);\r
+             if(mcLabelDStar>=0) isDStarMCtag = kTRUE;\r
+              if(!isDStarMCtag) continue;\r
+            AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);\r
+            //check if DStar from B\r
+            Int_t labelMother = MCDStar->GetMother();\r
+            AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
+            if(!mother) continue;\r
+            Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
+            if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
+           }\r
+                        \r
+           \r
+           phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar);\r
+            \r
+            // set the phi of the D meson in the correct range\r
+            \r
+           Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());\r
+            \r
+           \r
+           Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma\r
+           //  Double_t mD0Window=0.074/3;\r
+            \r
+            Double_t mD0Window= fD0Window[ptbin]/3;\r
+            //cout << "Check with new window " << fD0Window[ptbin]/3 << endl;\r
+            \r
+           \r
+            \r
+           invMassDZero = dstarD0pi->InvMassD0();\r
+           if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);\r
+            \r
+           deltainvMDStar = dstarD0pi->DeltaInvMass();\r
+            \r
+           //good D0 candidates\r
+           if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){\r
+             \r
+             if(!fmixing)      ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
+             // good D*?\r
+             if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow){\r
+               \r
+               if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight);\r
+               if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);\r
+               isInPeak = kTRUE;\r
+        EventHasDStarCandidate = kTRUE;\r
+               nOfDStarCandidates++;\r
+             } // end Good D*\r
+              \r
+                //  D* sideband?\r
+             if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)<fDMesonSigmas[3]*dmDStarWindow)){\r
+               isInDStarSideBand = kTRUE;\r
+              EventHasDStarSideBandCandidate = kTRUE;\r
+             } // end D* sideband\r
+              \r
+            }// end good D0 candidates\r
+            \r
+            //D0 sidebands\r
+           if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){\r
+             if(!fmixing)((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
+             if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight);\r
+              \r
+             if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0] *dmDStarWindow){ // is in DStar peak region?\r
+               if(!fmixing)    ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight);\r
+               isInDZeroSideBand = kTRUE;\r
+        EventHasDZeroSideBandCandidate = kTRUE;\r
+               nOfSBCandidates++;\r
+               if(!fmixing)    ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);\r
+             }\r
+              \r
+           }//end if sidebands\r
+            \r
+          \r
+            \r
+            \r
+           if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME\r
+           \r
+           \r
+            // check properties of the events containing the D*\r
+\r
+     \r
+            \r
+            isDStarCand = kTRUE;\r
+            \r
+            // charge of the daughter od the\r
+           daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();\r
+            \r
+           \r
+           trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();\r
+           trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();\r
+           trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();\r
+           \r
+           // end here the reco\r
+            \r
+           \r
+         }// end of if for applied reconstruction to D*\r
+         \r
+         Int_t DStarLabel = -1;\r
+         \r
+         if(!fReco){ // use pure MC information\r
+          \r
+        // get the DStar Particle\r
+           DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));\r
+           if (!DStarMC) {\r
+             AliWarning("Careful: DStar MC Particle not found in tree, skipping");\r
+             continue;\r
+           }\r
+           DStarLabel = DStarMC->GetLabel();\r
+        if(DStarLabel>0)isDStarMCtag = kTRUE;\r
+           \r
+           Int_t PDG =TMath::Abs(DStarMC->PdgCode());\r
+           if(PDG !=413) continue; // skip if it is not a DStar\r
+           // check fiducial acceptance\r
+        if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;\r
+            \r
+            //check if DStar from B\r
+            Int_t labelMother = DStarMC->GetMother();\r
+            AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
+                if(!mother) continue;\r
+            Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
+            if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
+          \r
+          Bool_t isDZero = kFALSE;\r
+          Bool_t isSoftPi = kFALSE;\r
+          \r
+          if(fUseHadronicChannelAtKineLevel){\r
+          //check decay channel on MC ************************************************\r
+                Int_t NDaugh = DStarMC->GetNDaughters();\r
+                if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs\r
+                \r
+                for(Int_t i=0; i<NDaugh;i++){ // loop on daughters\r
+                        Int_t daugh_label = DStarMC->GetDaughter(i);\r
+                        AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));\r
+                        if(!mcDaughter) continue;\r
+                        Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());\r
+                        if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;\r
+             \r
+                        if(daugh_pdg == 421) {isDZero = kTRUE;\r
+                            Int_t NDaughD0 = mcDaughter->GetNDaughters();\r
+                            if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs\r
+                            trackiddaugh0 = mcDaughter->GetDaughter(0);\r
+                            trackiddaugh1 = mcDaughter->GetDaughter(1);\r
+                            Bool_t isKaon = kFALSE;\r
+                            Bool_t isPion = kFALSE;\r
+               \r
+                            for(Int_t k=0;k<NDaughD0;k++){\r
+                                Int_t labelD0daugh = mcDaughter->GetDaughter(k);\r
+                                AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));\r
+                                if(!mcGrandDaughter) continue;\r
+                                Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());\r
+                                if(granddaugh_pdg==321) isKaon = kTRUE;\r
+                                if(granddaugh_pdg==211) isPion = kTRUE;\r
+                            }\r
+                            if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0\r
+                        }\r
+             \r
+                        if(daugh_pdg == 211) {\r
+                            isSoftPi = kTRUE;\r
+                            daughtercharge = mcDaughter->Charge();\r
+                            trackidsoftPi = daugh_label;}\r
+                    }\r
+              if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel\r
+          } // end check decay channel\r
+           \r
+           ptDStar = DStarMC->Pt();\r
+           phiDStar = DStarMC->Phi();\r
+           etaDStar = DStarMC->Eta();\r
+          \r
+          if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
+           \r
+         } // end use pure MC information\r
+         \r
+    \r
+        // getting the number of triggers in the MCtag D* case\r
+        if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);\r
+         if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);\r
+         if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);\r
+         \r
+         \r
+         fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters\r
+         fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);\r
+         \r
+         \r
+         // ************************************************ CORRELATION ANALYSIS STARTS HERE\r
+         \r
+         \r
+         Bool_t execPool = fCorrelator->ProcessEventPool();\r
+         \r
+         \r
+         if(fmixing && !execPool) {\r
+           AliInfo("Mixed event analysis: pool is not ready");\r
+            if(!isEventMixingFilledPeak && isInPeak)  {\r
+             ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1);\r
+             isEventMixingFilledPeak = kTRUE;\r
+            }\r
+            if (!isEventMixingFilledSB && isInDZeroSideBand)  {\r
+             ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3);\r
+             isEventMixingFilledSB=kTRUE;\r
+            }\r
+            \r
+           continue;\r
+         }\r
+         \r
+         // check event topology\r
+         if(fmixing&&execPool){\r
+            // pool is ready - run checks on bins filling\r
+            if(!isEventMixingFilledPeak && isInPeak)  {\r
+             ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0);\r
+             if(fFullmode) EventMixingChecks(aodEvent);\r
+             isEventMixingFilledPeak = kTRUE;\r
+            }\r
+            \r
+            if(!isEventMixingFilledSB && isInDZeroSideBand) {\r
+             ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2);\r
+             isEventMixingFilledSB=kTRUE;\r
+            }\r
+         }\r
+         \r
+         Int_t NofEventsinPool = 1;\r
+         if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();\r
+         \r
+         \r
+         for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one\r
+           \r
+           Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);\r
+           if(!analyzetracks) {\r
+             AliInfo("AliHFCorrelator::Cannot process the track array");\r
+             continue;\r
+           }\r
+            \r
+            //initialization of variables for correlations with leading particles\r
+            Double_t DeltaPhiLeading = -999.;\r
+          Double_t DeltaEtaLeading = -999.;\r
+       \r
+           \r
+           Int_t NofTracks = fCorrelator->GetNofTracks();\r
+            \r
+           \r
+           if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks);\r
+           if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks);\r
+           \r
+           \r
+           \r
+            Double_t arraytofill[5];\r
+            Double_t MCarraytofill[7];\r
+            \r
+            \r
+            Double_t weight;\r
+            \r
+          for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates\r
+              Bool_t runcorrelation = fCorrelator->Correlate(iTrack);\r
+              if(!runcorrelation) continue;\r
+              \r
+              Double_t DeltaPhi = fCorrelator->GetDeltaPhi();\r
+              Double_t DeltaEta = fCorrelator->GetDeltaEta();\r
+             \r
+              AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();\r
+              if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}\r
+              \r
+              Double_t ptHad = hadron->Pt();\r
+              Double_t phiHad = hadron->Phi();\r
+              Double_t etaHad = hadron->Eta();\r
+              Int_t label = hadron->GetLabel();\r
+              Int_t trackid = hadron->GetID();\r
+              Double_t efficiency = hadron->GetWeight();\r
+              \r
+              weight = 1;\r
+              if(fUseEfficiencyCorrection && efficiency){\r
+                  weight = DmesonWeight * (1./efficiency);\r
+              }\r
+        \r
+              phiHad = fCorrelator->SetCorrectPhiRange(phiHad);\r
+              \r
+              \r
+              if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency);\r
+              \r
+              arraytofill[0] = DeltaPhi;\r
+              arraytofill[1] = DeltaEta;\r
+              arraytofill[2] = ptDStar;\r
+              arraytofill[3] = ptHad;\r
+              arraytofill[4] = poolbin;\r
+                \r
+              \r
+              MCarraytofill[0] = DeltaPhi;\r
+              MCarraytofill[1] = DeltaEta;\r
+              MCarraytofill[2] = ptDStar;\r
+              MCarraytofill[3] = ptHad;\r
+              MCarraytofill[4] = poolbin;\r
+             \r
+            \r
+             if(fmontecarlo){\r
+               if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks);\r
+               if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks);\r
+               if(label<0) continue; // skip track with wrong label\r
+             }\r
+              \r
+               Bool_t isDdaughter = kFALSE;\r
+            // skip the D daughters in the correlation\r
+              if(!fmixing && fReco){\r
+                  if(trackid == trackiddaugh0) continue;\r
+                  if(trackid == trackiddaugh1) continue;\r
+                  if(trackid == trackidsoftPi) continue;\r
+              }\r
+              \r
+              if(!fmixing && !fReco){\r
+                  AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);\r
+                  if(!part) continue;\r
+                  if(IsDDaughter(DStarMC, part)) continue;\r
+                  cout << "Skipping DStar  daugheter " << endl;\r
+              }\r
+              if(!fmixing && !fReco && fmontecarlo){  // skip D* Daughetrs if it is Pure MCDStar\r
+                    Int_t hadronlabel = label;\r
+                    for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers\r
+                        if(DStarLabel<0){ break;}\r
+                        if(hadronlabel<0) { break;}\r
+                        AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fmcArray->At(hadronlabel));\r
+                        if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;}\r
+                        hadronlabel = mcParticle->GetMother();\r
+                        if(hadronlabel == DStarLabel) isDdaughter = kTRUE;\r
+                    }\r
+               \r
+                  if(isDdaughter && fDebugLevel){\r
+                      std::cout << "It is the D* daughter with label " << label << std::endl;\r
+                      std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl;\r
+                      std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl;\r
+                      std::cout << "Soft pi label = " << trackidsoftPi << std::endl;\r
+                  }\r
+                    \r
+                                       if(isDdaughter) continue; // skip if track is from DStar\r
+                               }\r
+                \r
+                               // ================ FILL CORRELATION HISTOGRAMS ===============================\r
+                               \r
+                // monte carlo case (mc tagged D*)\r
+             if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo\r
+                \r
+               Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)\r
+                             \r
+               MCarraytofill[5] = 0;\r
+               if(PartSource[0]) MCarraytofill[5] = 1;\r
+               if(PartSource[1]) MCarraytofill[5] = 2;\r
+               if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3;\r
+                    if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4;\r
+                    if(PartSource[3]) MCarraytofill[5] = 5;\r
+                    if(!isDfromB) MCarraytofill[6] = 0;\r
+                    if(isDfromB) MCarraytofill[6] = 1;\r
+                   if(!fReco && TMath::Abs(etaHad)>0.8) {\r
+                     delete [] PartSource;\r
+                     continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
+                   }\r
+                    ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill);\r
+                    \r
+                    delete[] PartSource;\r
+                               }\r
+              \r
+                // Good DStar canidates\r
+                               if(isInPeak)  {\r
+                    \r
+                                       if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
+                    if(fselect==1)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
+                    if(fselect==2)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
+                    if(fselect==3)  ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
+                                       \r
+                                   ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent);\r
+                                       if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad);\r
+                                    \r
+                               }\r
+              \r
+                // Sidebands from D0 candidate\r
+                               if(isInDZeroSideBand) {\r
+                                       \r
+                                       if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
+                    if(fselect==1)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
+                    if(fselect==2)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
+                    if(fselect==3)  ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
+                    \r
+                                       if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad);\r
+                    \r
+                }\r
+              \r
+              // Sidebands from D* candidate\r
+                if(isInDStarSideBand) {\r
+                                       \r
+                                       if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC  truth only if particles are in acceptance\r
+                    if(fselect==1 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
+                    if(fselect==2 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
+                    if(fselect==3 && fFullmode)  ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
+\r
+                               }\r
+                \r
+                \r
+                       } // end loop on track candidates\r
+            \r
+                       \r
+                       \r
+            // fill the leading particle histograms\r
+            \r
+            if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
+                       if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
+            \r
+               } // end loop on events in the pool\r
+        \r
+       }// end loop on D* candidates\r
+       \r
+    \r
\r
+        // check events with D* or SB canidates\r
+    if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition);\r
+     if(fFullmode && EventHasDZeroSideBandCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition);\r
+    \r
+    if(fFullmode && EventHasDStarCandidate)  ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition);\r
+    \r
+    \r
+       if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent);\r
+    if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent);\r
+       \r
+    // update event pool\r
+    Bool_t updated = fCorrelator->PoolUpdate();\r
+       \r
+    // if(updated) EventMixingChecks(aodEvent);\r
+       if(!updated) AliInfo("Pool was not updated");\r
+       \r
+    \r
+} //end the exec\r
+\r
+//________________________________________ terminate ___________________________\r
+void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)\r
+{    \r
+       // The Terminate() function is the last function to be called during\r
+       // a query. It always runs on the client, it can be used to present\r
+       // the results graphically or save the results to file.\r
+       \r
+       AliAnalysisTaskSE::Terminate();\r
+       \r
+       fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
+       if (!fOutput) {     \r
+               printf("ERROR: fOutput not available\n");\r
+               return;\r
+       }\r
+\r
+       return;\r
+}\r
+//_____________________________________________________\r
+Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {\r
+    \r
+    //Daughter removal in MCKine case\r
+    Bool_t isDaughter = kFALSE;\r
+    Int_t labelD0 = d->GetLabel();\r
+    \r
+    Int_t mother = track->GetMother();\r
+    \r
+    //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)\r
+    while (mother > 0){\r
+        AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!\r
+        if (mcMoth){\r
+            if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;\r
+            mother = mcMoth->GetMother(); //goes back by one\r
+        } else{\r
+            AliError("Failed casting the mother particle!");\r
+            break;\r
+        }\r
+    }\r
+    \r
+    return isDaughter;\r
+}\r
+\r
+//_____________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){\r
+    \r
+    //cout << "DEFINING THNSPARSES "<< endl;\r
+    \r
+    Double_t Pi = TMath::Pi();\r
+       Int_t nbinscorr = fPhiBins;\r
+       Double_t lowcorrbin = -0.5*Pi;\r
+       Double_t upcorrbin = 1.5*Pi;\r
+    // define the THnSparseF\r
+    \r
+    //sparse bins\r
+    \r
+    //1 delta_phi\r
+    //2 delta_eta\r
+    //3 D* pt\r
+    //4 multiplicity\r
+    //5 track pt\r
+    //6 zVtx position\r
+    \r
+    Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins());\r
+    \r
+    \r
+    Int_t nbinsSparse[5]=         {nbinscorr,   32,30,250,nbinsPool};\r
+    Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0,  0,-0.5};\r
+    Double_t binUpLimitSparse[5]= {upcorrbin,  1.6,30,25,nbinsPool-0.5};\r
+  \r
+    Int_t MCnbinsSparse[7]=         {nbinscorr,   32,30,250,nbinsPool,10,2};    \r
+    Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0,  0,-0.5,-0.5,-0.5};      //  \r
+    Double_t MCbinUpLimitSparse[7]= {upcorrbin,  1.6,30,25,nbinsPool-0.5,9.5,1.5};\r
+    \r
+    TString sparsename = "CorrelationsDStar";\r
+    if(fselect==1) sparsename += "Hadron";\r
+       if(fselect==2) sparsename += "Kaon";\r
+       if(fselect==3) sparsename += "KZero";\r
+    \r
+    TString D0Bkgsparsename = "DZeroBkg";\r
+    D0Bkgsparsename += sparsename;\r
+    \r
+    TString DStarBkgsparsename = "DStarBkg";\r
+    DStarBkgsparsename += sparsename;\r
+    \r
+    TString MCSparseName = "MCDStar";\r
+    MCSparseName += sparsename;\r
+    // signal correlations\r
+    THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
+    \r
+    // bkg correlations from D0 sidebands\r
+    THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
+    \r
+    // bkg correlations from D* sidebands\r
+    THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
+    \r
+    \r
+    THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse);\r
+    \r
+    MCCorrelations->GetAxis(5)->SetBinLabel(1," All ");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(7," from c");\r
+       MCCorrelations->GetAxis(5)->SetBinLabel(8," from b");\r
+    \r
+    MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c");\r
+    MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b");\r
+    \r
+    Correlations->Sumw2();\r
+    DZeroBkgCorrelations->Sumw2();\r
+    DStarBkgCorrelations->Sumw2();\r
+    \r
+    fOutput->Add(Correlations);\r
+    fOutput->Add(DZeroBkgCorrelations);\r
+    if(fFullmode) fOutput->Add(DStarBkgCorrelations);\r
+    if(fmontecarlo) fOutputMC->Add(MCCorrelations);\r
+    \r
+    \r
+}\r
+//__________________________________________________________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){\r
+       \r
+       Double_t Pi = TMath::Pi();\r
+       Int_t nbinscorr = fPhiBins;\r
+       Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center\r
+       Double_t upcorrbin = 1.5*Pi ;\r
+       \r
+       // ========================= histograms for both Data and MonteCarlo\r
+       \r
+       \r
+       TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);\r
+    NofEvents->GetXaxis()->SetBinLabel(1," All events");\r
+       NofEvents->GetXaxis()->SetBinLabel(2," Selected events");\r
+       NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");\r
+       NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");\r
+       NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");\r
+       NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");\r
+       NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");\r
+       NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");\r
+    NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");\r
+    NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");\r
+    NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");\r
+    fOutput->Add(NofEvents);\r
+       \r
+       \r
+       \r
+       \r
+       TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);\r
+       if(!fmixing && fFullmode) fOutput->Add(D0InvMass);\r
+       \r
+       TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);\r
+       if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB);\r
+       \r
+       //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);\r
+       //if(!fmixing) fOutput->Add(DeltaInvMass);\r
+    TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution; D* p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);\r
+       if(!fmixing) fOutput->Add(DeltaInvMass);\r
+       \r
+       TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution; SB p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2);\r
+       if(!fmixing) fOutput->Add(bkgDeltaInvMass);\r
+    \r
+    DeltaInvMass->Sumw2();\r
+    bkgDeltaInvMass->Sumw2();\r
+       \r
+       TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);\r
+       if(!fmixing) fOutput->Add(RecoPtDStar);\r
+       \r
+       TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);\r
+       if(!fmixing) fOutput->Add(RecoPtBkg);\r
+    \r
+    TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);\r
+    if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);\r
+    \r
+    TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);\r
+    if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);\r
+       \r
+       TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);\r
+       if(!fmixing) fOutput->Add(MCtagPtDStar);\r
+       \r
+       TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);\r
+    if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra);\r
+       \r
+       TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);\r
+       if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF);\r
+       \r
+       TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5);\r
+       if(fFullmode) fOutput->Add(NofTracksInPeak);\r
+       \r
+       TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5);\r
+       if(fFullmode) fOutput->Add(NofTracksInSB);\r
+       \r
+       TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
+       if(fFullmode)fOutput->Add(TracksInPeakSpectra);\r
+       \r
+       TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
+       if(fFullmode)fOutput->Add(TracksInSBSpectra);\r
+       \r
+       \r
+       //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);\r
+       //if(fmixing) fOutput->Add(EventMixingCheck);\r
+    \r
+    \r
+    TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
+       if(fFullmode)fOutput->Add(EventPropsCheck);\r
+    \r
+    TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
+       if(fFullmode)fOutput->Add(EventPropsCheckifDStar);\r
+    \r
+    TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
+       if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB);\r
+    \r
+    TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
+       if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB);\r
+    \r
+       \r
+    TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005);\r
+       if(fFullmode)fOutput->Add(WeightChecks);\r
+    \r
+       \r
+       \r
+       TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
+       fOutput->Add(PhiEtaTrigger);\r
+       \r
+       TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
+       fOutput->Add(PhiEtaSideBand);\r
+       \r
+       TH3F * PhiEtaPart = new TH3F("PhiEtaPart","#phi distribution of the associated particle; #phi; #eta; multiplicity",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9,100,0,1000);\r
+       fOutput->Add(PhiEtaPart);\r
+    \r
+    TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
+       if(fFullmode)fOutput->Add(DStarCandidates);\r
+    \r
+    TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
+       if(fFullmode)fOutput->Add(SBCandidates);\r
+       \r
+       \r
+       //correlations histograms\r
+       TString histoname1 = "DPhiDStar";\r
+       if(fselect==1) histoname1 += "Hadron";\r
+       if(fselect==2) histoname1 += "Kaon";\r
+       if(fselect==3) histoname1 += "KZero";\r
+       \r
+       /*\r
+     TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+     \r
+     TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+     \r
+     //side band background histograms\r
+     TString histoname2 = "bkg";\r
+     histoname2 += histoname1;\r
+     TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+     TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+     \r
+     \r
+     fOutput->Add(DPhiDStar);\r
+     \r
+     if(fselect==3){fOutput->Add(DPhiDStarKZero1);}\r
+     \r
+     fOutput->Add(bkgDPhiDStar);\r
+     \r
+     if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}\r
+     \r
+     */\r
+       // leading particle\r
+       TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       if(fFullmode)fOutput->Add(leadingcand);\r
+       if(fFullmode)fOutput->Add(leadingsidebands);\r
+       \r
+       // ========================= histos for analysis on MC only\r
+       \r
+       TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);\r
+       if(fmontecarlo) fOutputMC->Add(EventTypeMC);\r
+       \r
+       TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);\r
+       MCSources->GetXaxis()->SetBinLabel(1," All ");\r
+       MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
+       MCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
+       MCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
+       MCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
+       MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
+       MCSources->GetXaxis()->SetBinLabel(7," from c");\r
+       MCSources->GetXaxis()->SetBinLabel(8," from b");\r
+       \r
+       if(fmontecarlo) fOutputMC->Add(MCSources);\r
+    \r
+    // leading particle from mc source\r
+    TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(1," All ");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");\r
+       LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");\r
+       \r
+       if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources);\r
+       \r
+    // all hadrons\r
+       TString histoname3 = "MCTag";\r
+       histoname3 += histoname1;\r
+       TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString histoname44 = "CharmDOrigin";\r
+       histoname44 += histoname1;\r
+       histoname44 += "MC";\r
+       \r
+       TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       \r
+       TString histoname54 = "BeautyDOrigin";\r
+       histoname54 += histoname1;\r
+       histoname54 += "MC";\r
+       TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString histoname55 = "BeautyBOrigin";\r
+       histoname55 += histoname1;\r
+       histoname55 += "MC";\r
+       TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString histoname4 = "CharmQuarkOrigin";\r
+       histoname4 += histoname1;\r
+       histoname4 += "MC";\r
+       TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString histoname5 = "BeautyQuarkOrigin";\r
+       histoname5 += histoname1;\r
+       histoname5 += "MC";\r
+       TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString histoname6 = "NonHFOrigin";\r
+       histoname6 += histoname1;\r
+       histoname6 += "MC";\r
+       TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+    \r
+    if(fmontecarlo && fFullmode){\r
+        \r
+        fOutputMC->Add(MCTagDPhiDStar);\r
+        fOutputMC->Add(CharmDOriginDPhiDStar);\r
+        fOutputMC->Add(BeautyDOriginDPhiDStar);\r
+        fOutputMC->Add(BeautyBOriginDPhiDStar);\r
+        fOutputMC->Add(CharmQuarkOriginDPhiDStar);\r
+        fOutputMC->Add(BeautyQuarkOriginDPhiDStar);\r
+               fOutputMC->Add(NonHFOriginDPhiDStar);\r
+        \r
+       }\r
+    \r
+    // ========================= histos for analysis on MC\r
+    // all leading hadron\r
+       TString Leadinghistoname3 = "LeadingMCTag";\r
+       Leadinghistoname3 += histoname1;\r
+       TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+    \r
+       TString Leadinghistoname44 = "LeadingCharmDOrigin";\r
+       Leadinghistoname44 += histoname1;\r
+       Leadinghistoname44 += "MC";\r
+       \r
+       TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       \r
+       TString Leadinghistoname54 = "LeadingBeautyDOrigin";\r
+       Leadinghistoname54 += histoname1;\r
+       Leadinghistoname54 += "MC";\r
+       TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString Leadinghistoname55 = "LeadingBeautyBOrigin";\r
+       Leadinghistoname55 += histoname1;\r
+       Leadinghistoname55 += "MC";\r
+       TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";\r
+       Leadinghistoname4 += histoname1;\r
+       Leadinghistoname4 += "MC";\r
+       TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+       \r
+       TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";\r
+       Leadinghistoname5 += histoname1;\r
+       Leadinghistoname5 += "MC";\r
+       TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
+    \r
+    \r
+       \r
+       \r
+       if(fmontecarlo && fFullmode){\r
+               \r
+               fOutputMC->Add(LeadingMCTagDPhiDStar);\r
+               fOutputMC->Add(LeadingCharmDOriginDPhiDStar);\r
+               fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);\r
+               fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);\r
+               fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);\r
+               fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);\r
+               \r
+       }\r
+       \r
+       TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5);\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");\r
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");\r
+       if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);\r
+       \r
+       TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5);\r
+       if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels);\r
+       \r
+       // ============================= EVENT MIXING CHECKS ======================================\r
+       \r
+       Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();\r
+       Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();\r
+       Int_t NofCentBins = fCuts2->GetNCentPoolBins();\r
+       Double_t * CentBins = fCuts2->GetCentPoolBins();\r
+       Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();\r
+       Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();\r
+    \r
+    \r
+    \r
+       Int_t k =0;\r
+       \r
+       if(fSystem == AA) k = 100; // PbPb centrality\r
+       if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity\r
+       \r
+       \r
+       //Double_t minvalue = CentBins[0];\r
+       //Double_t maxvalue = CentBins[NofCentBins+1];\r
+       //Double_t Zminvalue = ZVrtxBins[0];\r
+       //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];\r
+       \r
+       Double_t minvalue, maxvalue;\r
+    Double_t Zminvalue, Zmaxvalue;\r
+    \r
+    Zminvalue = -15.;\r
+    Zmaxvalue = 15;\r
+    if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb\r
+    if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity\r
+    \r
+       //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
+    // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
+    // Double_t * events = Nevents;\r
+    Double_t eventsv[] ={0,1000000};\r
+    //Double_t * events = new Double_t[2];\r
+   // events[0] = 0;\r
+//     events[1] = 1000000;\r
+    Double_t *events = eventsv;\r
+    Int_t Nevents = 1000000;\r
+    //  TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events);\r
+    \r
+    TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]);\r
+    \r
+       EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
+       EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
+       EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");\r
+       if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin);\r
+       \r
+       Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;\r
+       //Int_t Diff = MaxNofTracks-MinNofTracks;\r
+       \r
+    //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};\r
+    // Double_t  * trackN = Ntracks;\r
+       \r
+       TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks);\r
+       NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
+       NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
+       NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");\r
+       \r
+       if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin);\r
+       \r
+       TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);\r
+       NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
+       NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");\r
+       if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls);\r
+       \r
+    \r
+       \r
+       TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);\r
+       EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
+       EventProps->GetYaxis()->SetTitle("Z vertex [cm]");\r
+       if(fmixing && fFullmode) fOutput->Add(EventProps);\r
+    \r
+    TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);\r
+    CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");\r
+       CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");\r
+    CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");\r
+       CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");\r
+       \r
+    if(fmixing) fOutput->Add(CheckPoolReadiness);\r
+    \r
+       \r
+}\r
+\r
+//__________________________________________________________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){\r
+    \r
+\r
+  //Float_t* ptbins = fCuts->GetPtBinLimits();\r
+  if(fD0Window) delete fD0Window;\r
+  fD0Window = new Float_t[fNofPtBins];\r
+    \r
+  AliInfo("Enlarging the D0 mass windows from cut object\n"); \r
+  Int_t nvars = fCuts->GetNVars();\r
+\r
+  if(nvars<1){\r
+    AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");\r
+    return;\r
+  }\r
+  Float_t** rdcutsvalmine=new Float_t*[nvars];\r
+  for(Int_t iv=0;iv<nvars;iv++){\r
+    rdcutsvalmine[iv]=new Float_t[fNofPtBins];\r
+  }\r
+    \r
+    \r
+    for (Int_t k=0;k<nvars;k++){\r
+      for (Int_t j=0;j<fNofPtBins;j++){\r
+            \r
+       // enlarge D0 window\r
+       if(k==0)        {\r
+         fD0Window[j] =fCuts->GetCutValue(0,j);\r
+         rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);\r
+         cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;\r
+       }\r
+       else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);\r
+                       \r
+       // set same windows\r
+       //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);\r
+      }\r
+    }\r
+    \r
+    fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);\r
+    \r
+    AliInfo("\n New windows set\n");     \r
+    fCuts->PrintAll();\r
+    \r
+    \r
+    for(Int_t iv=0;iv<nvars;iv++){\r
+      delete rdcutsvalmine[iv];\r
+    }\r
+    delete [] rdcutsvalmine;\r
+    \r
+}\r
+\r
+\r
+//____________________________  Run checks on event mixing ___________________________________________________\r
+void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){\r
+       \r
+    \r
+       AliCentrality *centralityObj = 0;\r
+       Int_t multiplicity = -1;\r
+       Double_t MultipOrCent = -1;\r
+       \r
+       // get the pool for event mixing\r
+       if(fSystem != AA){ // pp\r
+               multiplicity = AOD->GetNTracks();\r
+               MultipOrCent = multiplicity; // convert from Int_t to Double_t\r
+       }\r
+       if(fSystem == AA){ // PbPb\r
+               \r
+               centralityObj = AOD->GetHeader()->GetCentralityP();\r
+               MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
+               AliInfo(Form("Centrality is %f", MultipOrCent));\r
+       }\r
+       \r
+       AliAODVertex *vtx = AOD->GetPrimaryVertex();\r
+       Double_t zvertex = vtx->GetZ(); // zvertex\r
+       \r
+       \r
+       \r
+       \r
+       AliEventPool * pool = fCorrelator->GetPool();\r
+       \r
+    \r
+       \r
+       \r
+       ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool\r
+       ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties\r
+       \r
+       ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool\r
+       ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool\r
+}\r
+       \r
+\r
+\r
+\r
+\r