Updates in D-hadron correlation code (Sandro)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Sep 2012 08:38:36 +0000 (08:38 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Sep 2012 08:38:36 +0000 (08:38 +0000)
In AliHFAssociatedTrackCuts:
- possibility to select a specific process for MC analysis (gluon splitting, pair production etc)
- IsMCpartFromHF returns an array of 4 booleans containing infos on track origin (from charm, from beauty, from D meson, from b meson)
- a string is added to store comments in the cut object

In AliHFCorrelator
- added flag to select if perform the analysis on reconstructed tracks or generated particles
- in AcceptAndReduceTracks adde the possibility to save a TObjArray of charged particles in a given pt range/PDG code

In AliAnalysisTaskDStarCorrelations:
- added control histograms for event mixing
- updated methods that calls AliHFAssociatedTrackCuts::IsMCpartFromHF
- added tha analysis with correlation to leading particles
- added possibility to select on event type (for MC)
- bug fix with normalizationcounter

In AliAnalysisTaskSEDplusCorrelations:
- updated to make it compatible with new method AliHFAssociatedTrackCuts::IsMCpartFromHF

PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx
PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h
PWGHF/correlationHF/AliAnalysisTaskSEDplusCorrelations.cxx
PWGHF/correlationHF/AliAnalysisTaskSEDplusCorrelations.h
PWGHF/correlationHF/AliHFAssociatedTrackCuts.cxx
PWGHF/correlationHF/AliHFAssociatedTrackCuts.h
PWGHF/correlationHF/AliHFCorrelator.cxx
PWGHF/correlationHF/AliHFCorrelator.h
PWGHF/correlationHF/macros/AddTaskDStarCorrelations.C
PWGHF/correlationHF/macros/makeTFileAssociatedTrackCuts.C

index 1f924e4..bc983da 100644 (file)
@@ -49,6 +49,8 @@
 #include "AliNormalizationCounter.h"
 #include "AliReducedParticle.h"
 #include "AliHFCorrelator.h"
+#include "AliAODMCHeader.h"
+#include "AliEventPoolManager.h"
 
 using std::cout;
 using std::endl;
@@ -61,7 +63,6 @@ ClassImp(AliAnalysisTaskDStarCorrelations)
 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
 AliAnalysisTaskSE(),
 fhandler(0x0),
-//fPoolMgr(0x0),      
 fmcArray(0x0),
 fCounter(0x0),
 fCorrelator(0x0),
@@ -69,11 +70,13 @@ fselect(0),
 fmontecarlo(kFALSE),
 fmixing(kFALSE),
 fSystem(kFALSE),
+fReco(kTRUE),
 fEvents(0),
-fDebug(0),
+fDebugLevel(0),
 fDisplacement(0),
 
-fOutput(0x0),            
+fOutput(0x0), 
+fOutputMC(0x0),
 fCuts(0),
 fCuts2(0)
 {
@@ -92,21 +95,25 @@ fselect(0),
 fmontecarlo(kFALSE),
 fmixing(kFALSE),
 fSystem(kFALSE),
+fReco(kTRUE),
 fEvents(0),
-fDebug(0),
+fDebugLevel(0),
 fDisplacement(0),
 
-fOutput(0x0),                
+fOutput(0x0),   
+fOutputMC(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
+       
+       DefineOutput(1,TList::Class()); // histos from data and MC
+       DefineOutput(2,TList::Class()); // histos from MC
+       DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
+       DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
+       DefineOutput(5,AliNormalizationCounter::Class());   // normalization
 }
 
 //__________________________________________________________________________
@@ -124,6 +131,7 @@ AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
        if(fCounter) {delete fCounter; fCounter = 0;}
        if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
        if(fOutput) {delete fOutput; fOutput = 0;}
+       if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
        if(fCuts) {delete fCuts; fCuts = 0;}
        if(fCuts2) {delete fCuts2; fCuts2=0;}
 
@@ -134,7 +142,7 @@ void AliAnalysisTaskDStarCorrelations::Init(){
        //
        // Initialization
        //
-       if(fDebug > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
+       if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
        
        AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
        
@@ -142,7 +150,7 @@ void AliAnalysisTaskDStarCorrelations::Init(){
        
        
        // Post the D* cuts
-       PostData(2,copyfCuts);
+       PostData(3,copyfCuts);
        
        // Post the hadron cuts
        PostData(4,fCuts2);
@@ -162,9 +170,12 @@ void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
        fOutput = new TList();
        fOutput->SetOwner();
        
+       fOutputMC = new TList();
+       fOutputMC->SetOwner();
+       
        // define histograms
        DefineHistoForAnalysis();
-       fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(3)->GetContainer()->GetName()));
+       fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
        fCounter->Init();
        
     Double_t Pi = TMath::Pi();
@@ -174,31 +185,40 @@ void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
        fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
        fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
        fCorrelator->SetUseMC(fmontecarlo);
+       fCorrelator->SetUseReco(fReco);
        Bool_t pooldef = fCorrelator->DefineEventPool();
        
        if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
 
+
        
        PostData(1,fOutput); // set the outputs
-       PostData(3,fCounter); // set the outputs
+       PostData(2,fOutputMC); // set the outputs
+       PostData(5,fCounter); // set the outputs
 }
 //_________________________________________________
 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
 
        
-       
-       std::cout << " " << std::endl;
-       std::cout << "=================================================================================" << std::endl;
-       if(!fmixing){
-       if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;
-       if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;
-       if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;
-       }
-       if(fmixing){
-               if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;
-               if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;
-               if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;
-       }
+       if(fDebugLevel){
+               
+               if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;
+               if(!fReco) std::cout << "USING MC TRUTH" << std::endl;
+        std::cout << " " << std::endl;
+        std::cout << "=================================================================================" << std::endl;
+        if(!fmixing){
+            if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;
+            if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;
+            if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;
+        }
+        if(fmixing){
+            if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;
+            if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;
+            if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;
+        }
+        
+    }// end if debug
+    
        if (!fInputEvent) {
                Error("UserExec","NO EVENT FOUND!");
                return;
@@ -210,19 +230,28 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
          return;
        }
        
-       fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
+       
        
        fEvents++; // event counter
        ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
+    fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);
+  
+       // load MC array
        fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-       
        if(fmontecarlo && !fmcArray){
          AliError("Array of MC particles not found");
          return;
        }
+       
+
+       
+       
        Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
        if(!isEvSel) return;
-       ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
+       
+    fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
+    
+    ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
        //
        Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
        if(!correlatorON) {
@@ -232,10 +261,40 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
        ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);
        
        if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);
+       
+       
+       // check the event type
+       // load MC header
+       
+       if(fmontecarlo){
+               AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+               if (fmontecarlo && !mcHeader) {
+                       AliError("Could not find MC Header in AOD");
+                       return;
+               }
+       
+               Bool_t isMCeventgood = kFALSE;
+       
+               
+               Int_t eventType = mcHeader->GetEventType();
+               Int_t NMCevents = fCuts2->GetNofMCEventType();
+               
+               for(Int_t k=0; k<NMCevents; k++){
+                       Int_t * MCEventType = fCuts2->GetMCEventType();
+                       
+                       if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
+                       ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
+               }
+               
+               if(NMCevents && !isMCeventgood){
+               if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
+                       return; 
+               }
+               
+       } // end if montecarlo
+       
        // 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.    
@@ -274,7 +333,7 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
        Int_t pdgDgDStartoD0pi[2]={421,211};
        Int_t pdgDgD0toKpi[2]={321,211};
 
-       
+       Bool_t isDStarCand = kFALSE;
        //loop on D* candidates
        for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
                isInPeak = kFALSE;
@@ -382,6 +441,7 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                
                if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
                
+               isDStarCand = kTRUE;
                
            fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
                
@@ -411,9 +471,14 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                                continue;
                        }
                
+            //initialization of variables for correlations with leading particles
+            Double_t DeltaPhiLeading = -999.;
+                       Double_t DeltaEtaLeading = -999.;
+                       //Double_t ptleading = -999.;
+            Int_t labelleading = -999;
                
                        Int_t NofTracks = fCorrelator->GetNofTracks();
-               
+                       
                        for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
                        
                                Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
@@ -428,7 +493,7 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                                Double_t phiHad = hadron->Phi();
                                Double_t etaHad = hadron->Eta(); 
                                Double_t label = hadron->GetLabel(); 
-                               Double_t trackid = hadron->GetID();
+                               Int_t trackid = hadron->GetID();
                                
                                phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
                                
@@ -442,9 +507,17 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                                
                                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)
-                               
+                                       Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
                                        FillMCTagCorrelations(ptDStar,DeltaPhi,DeltaEta,ptHad,PartSource);
+                                       
+                                       
+                                       ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,0);
+                                       if(PartSource[0]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,1);
+                                       if(PartSource[1]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,2);
+                                       if(PartSource[2]&&PartSource[0]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,3);
+                                       if(PartSource[2]&&PartSource[1]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,4);
+                                       if(PartSource[3]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,5);
+                               
                                
                                }
                        
@@ -471,13 +544,28 @@ void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
                        
                        
                        } // end loop on track candidates
+            
+            // fill the leading particle histograms
+            
+            if(isInPeak) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
+                       if(isInSideBand) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
+            
+                       if(fmontecarlo && isDStarMCtag){
+                Bool_t* LeadPartSource = fCuts2->IsMCpartFromHF(labelleading,fmcArray);
+                FillMCTagLeadingCorrelations(ptDStar,DeltaPhiLeading,DeltaEtaLeading,LeadPartSource);
+                
+            }
+            
                } // end loop on events in the pool
                                
        }// end loop on D* candidates           
+       
+       
        Bool_t updated = fCorrelator->PoolUpdate();
+       
+       if(updated) EventMixingChecks(aodEvent);
        if(!updated) AliInfo("Pool was not updated");
        
-               //std::cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> END OF THE EVENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
        
        
                
@@ -513,9 +601,13 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        // ========================= histograms for both Data and MonteCarlo
        
        
-       TH1D * NofEvents = new TH1D("NofEvents","NofEvents",3,0,3);
+       TH1D * NofEvents = new TH1D("NofEvents","NofEvents",11,0,11);
        fOutput->Add(NofEvents);
-               
+       
+       
+       
+       
+       
        TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
        fOutput->Add(D0InvMass);
        
@@ -528,13 +620,13 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        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);
+       TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);
        fOutput->Add(RecoPtDStar);
        
-       TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",30,0,30);
+       TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);
        fOutput->Add(RecoPtBkg);
        
-       TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",30,0,30);
+       TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
        fOutput->Add(MCtagPtDStar);
        
        TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
@@ -552,15 +644,9 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        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);
        
        TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
        fOutput->Add(PhiEtaTrigger);
@@ -570,112 +656,380 @@ void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
        
        TH2F * PhiEtaPart = new TH2F("PhiEtaPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
        fOutput->Add(PhiEtaPart);
-
-
+       
+       
        //correlations histograms
        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);
+       TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.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);
+       TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
        
        
        fOutput->Add(DPhiDStar);
-
+       
        if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
        
        fOutput->Add(bkgDPhiDStar);
-
+       
        if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
-
-
-       // ========================= histos for analysis on MC
+       
+       
+       // leading particle
+       TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       fOutput->Add(leadingcand);
+       fOutput->Add(leadingsidebands);
+       
+       // ========================= histos for analysis on MC only
+       
+       TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
+       if(fmontecarlo) fOutputMC->Add(EventTypeMC);
+       
+       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 hadron Heavy flavour");
+       MCSources->GetXaxis()->SetBinLabel(3," from c->D");
+       MCSources->GetXaxis()->SetBinLabel(4," from b->D");
+       MCSources->GetXaxis()->SetBinLabel(5," from b->B");
+       MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
+       MCSources->GetXaxis()->SetBinLabel(7," from c");
+       MCSources->GetXaxis()->SetBinLabel(8," from b");
+       
+       if(fmontecarlo) fOutputMC->Add(MCSources);
+    
+    // leading particle from mc source
+    TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);
+       LeadingMCSources->GetXaxis()->SetBinLabel(1,"All ");
+       LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");
+       LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");
+       LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");
+       LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");
+       LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
+       LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");
+       LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");
+       
+       if(fmontecarlo) fOutputMC->Add(LeadingMCSources);
+       
+    // all hadrons
        TString histoname3 = "MCTag";
        histoname3 += histoname1;
-       TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,30,0,30,19,-0.95,0.95);
-        
+       TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.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);
+       TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.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);
+       TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.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);
+       TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TString histoname4 = "CharmQuarkOrigin";
+       histoname4 += histoname1;
+       histoname4 += "MC";
+       TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TString histoname5 = "BeautyQuarkOrigin";
+       histoname5 += histoname1;
+       histoname5 += "MC";
+       TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+    
+    if(fmontecarlo){
+        
+        fOutputMC->Add(MCTagDPhiDStar);
+        fOutputMC->Add(CharmDOriginDPhiDStar);
+        fOutputMC->Add(BeautyDOriginDPhiDStar);
+        fOutputMC->Add(BeautyBOriginDPhiDStar);
+        fOutputMC->Add(CharmQuarkOriginDPhiDStar);
+        fOutputMC->Add(BeautyQuarkOriginDPhiDStar);
+        
+       }
+    
+    // ========================= histos for analysis on MC
+    // all leading hadron
+       TString Leadinghistoname3 = "LeadingMCTag";
+       Leadinghistoname3 += histoname1;
+       TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+    
+       TString Leadinghistoname44 = "LeadingCharmDOrigin";
+       Leadinghistoname44 += histoname1;
+       Leadinghistoname44 += "MC";
+       
+       TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       
+       TString Leadinghistoname54 = "LeadingBeautyDOrigin";
+       Leadinghistoname54 += histoname1;
+       Leadinghistoname54 += "MC";
+       TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TString Leadinghistoname55 = "LeadingBeautyBOrigin";
+       Leadinghistoname55 += histoname1;
+       Leadinghistoname55 += "MC";
+       TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";
+       Leadinghistoname4 += histoname1;
+       Leadinghistoname4 += "MC";
+       TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+       
+       TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";
+       Leadinghistoname5 += histoname1;
+       Leadinghistoname5 += "MC";
+       TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
+    
+    
+       
        
        if(fmontecarlo){
+               
+               fOutputMC->Add(LeadingMCTagDPhiDStar);
+               fOutputMC->Add(LeadingCharmDOriginDPhiDStar);
+               fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);
+               fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);
+               fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);
+               fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);
+               
+       }
        
-       fOutput->Add(MCTagDPhiDStar);
-       fOutput->Add(CharmDOriginDPhiDStar);
-       fOutput->Add(BeautyDOriginDPhiDStar);
-       fOutput->Add(BeautyBOriginDPhiDStar);
+       TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi,50,-2.5,2.5,6,-0.5,6.5);
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");
+       MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");
+       if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);
+       
+       // ============================= EVENT MIXING CHECKS ======================================
+       
+       Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();
+       Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();
+       Int_t NofCentBins = fCuts2->GetNCentPoolBins();
+       Double_t * CentBins = fCuts2->GetCentPoolBins();
+       Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();
+       Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();
+       
+       Int_t k =0;
+       
+       if(fSystem) k = 100; // PbPb centrality
+       if(!fSystem) k = NofCentBins; // pp multiplicity
+       
+       
+       Double_t minvalue = CentBins[0];
+       Double_t maxvalue = CentBins[NofCentBins+1];
+       Double_t Zminvalue = ZVrtxBins[0];
+       Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];
        
-       }
        
 
+       Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};
+       Double_t * events = Nevents;
+       
+       TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,events);
+       EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
+       EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
+       EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");
+       if(fmixing) fOutput->Add(EventsPerPoolBin);
+       
+       Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;
+       Int_t Diff = MaxNofTracks-MinNofTracks;
+       
+       Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};
+       Double_t  * trackN = Ntracks;
+       
+       TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,trackN);
+       NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
+       NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
+       NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");
+       
+       if(fmixing) fOutput->Add(NofTracksPerPoolBin);
+       
+       TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);
+       NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");
+       NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");
+       if(fmixing) fOutput->Add(NofPoolBinCalls);
+       
+
+       
+       TH2D * EventProps = new TH2D("EventProps","Number of tracks in bin pool",k,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);
+       EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");
+       EventProps->GetYaxis()->SetTitle("Z vertex [cm]");
+       if(fmixing) fOutput->Add(EventProps);
        
 }
 
 
 
 //____________________________  Function for MC correlations ___________________________________________________
-void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Int_t mcSource){
-
-
-       if(fselect==1) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKaon"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("MCTagDPhiDStarKZero"))->Fill(DelPhi,ptTrig,DelEta);
+void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Bool_t *mcSource){
 
+       
+       
+       
+                       
+               if(fselect==1) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
+               if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarKaon"))->Fill(DelPhi,ptTrig,DelEta);
+               if(fselect==3) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarKZero"))->Fill(DelPhi,ptTrig,DelEta);
+               
+       
+               
+               ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(0);
+               
+               if(fDebugLevel){
+                       std::cout << "MC source " << mcSource[0] << " "  << mcSource[1] << " " << mcSource[2] << " " << mcSource[3] << std::endl;
+               
+                       if(mcSource[0]) std::cout << "mcSource 0 " << std::endl;
+                       if(mcSource[1]) std::cout << "mcSource 1 " << std::endl;
+                       if(mcSource[2]) std::cout << "mcSource 2 " << std::endl;
+                       if(mcSource[3]) std::cout << "mcSource 3 " << std::endl;
+               
+               }
+               if(mcSource[0]){ // is from charm quark
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(5); // all HF quarks
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(6); //  charm quarks
+                       if(fselect==1) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==3) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
+               }
+               
+               if(mcSource[1]){ // is from b quark
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(5); // all HF quarks
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(7); // beauty quarks
+                       if(fselect==1) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==3) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       
+               }
+               
+               if(mcSource[2]&&mcSource[0]){ // is from D meson and charm quark
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(1); // all HF mesons
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(2); //  charm + D
+                       if(fselect==1) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==3) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
+               }
+               
+               if(mcSource[2]&&mcSource[1]){ // is from D meson and b quark
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(1); // all HF mesons
+                       ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(3); //  beauty + D
+                       if(fselect==1) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
+                       if(fselect==3) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
+               }
+               
+       return;
+}
 
+//____________________________  Function for MC leading part correlations ___________________________________________________
+void AliAnalysisTaskDStarCorrelations::FillMCTagLeadingCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Bool_t *mcSource){
+    // correlations with leading hadron on MC
+    
+       if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingMCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
+       
+    
+    
+       ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(0);
+       
+       if(fDebugLevel){ std::cout << "MC source " << mcSource[0] << " "  << mcSource[1] << " " << mcSource[2] << " " << mcSource[3] << std::endl;
+    
+               if(mcSource[0]) std::cout << "mcSource 0 " << std::endl;
+               if(mcSource[1]) std::cout << "mcSource 1 " << std::endl;
+               if(mcSource[2]) std::cout << "mcSource 2 " << std::endl;
+               if(mcSource[3]) std::cout << "mcSource 3 " << std::endl;
+    }
+    
+       if(mcSource[0]){ // is from charm quark
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(5); // all HF quarks
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(6); //  charm quarks
+               if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingCharmQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+               
+       }
+    
+       if(mcSource[1]){ // is from b quaLeadingrk
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(5); // all HF quarks
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(7); // beauty quarks
+               if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingBeautyQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+               
+               
+       }
+    
+       if(mcSource[2]&&mcSource[0]){ // is from D meson and charm quark
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(1); // all HF mesons
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(2); //  charm + D
+               if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingCharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+               
+       }
+    
+       if(mcSource[2]&&mcSource[1]){ // is from D meson and b quark
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(1); // all HF mesons
+               ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(3); //  beauty + D
+               if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingBeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
+               
+       }
+       
+    
+       return;
+}
 
-((TH1F*)fOutput->FindObject("MCSources"))->Fill(0);
 
-if (mcSource==44){ // is from charm ->D
-       if(fselect==1) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
+//____________________________  Run checks on event mixing ___________________________________________________
+void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
        
+       AliCentrality *centralityObj = 0;
+       Int_t multiplicity = -1;
+       Double_t MultipOrCent = -1;
        
-       ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       ((TH1F*)fOutput->FindObject("MCSources"))->Fill(2);
+       // get the pool for event mixing
+       if(!fSystem){ // pp
+               multiplicity = AOD->GetNTracks();
+               MultipOrCent = multiplicity; // convert from Int_t to Double_t
        }
-
-if (mcSource==54){ // is from beauty -> D
-       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(3);
+       if(fSystem){ // PbPb
+               
+               centralityObj = AOD->GetHeader()->GetCentralityP();
+               MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
+               AliInfo(Form("Centrality is %f", MultipOrCent));
        }
+       
+       AliAODVertex *vtx = AOD->GetPrimaryVertex();
+       Double_t zvertex = vtx->GetZ(); // zvertex
+       
+       
+       
+       
+       AliEventPool * pool = fCorrelator->GetPool();
+       
 
-if (mcSource==55){ // is from beauty ->B
-       if(fselect==1) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutput->FindObject("BeautyBOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH3D*)fOutput->FindObject("BeautyOriginBDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(1);
-       if(fselect==3) ((TH1F*)fOutput->FindObject("MCSources"))->Fill(4);
-       }
-       return;
+       
+       
+       ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
+       ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
+       
+       ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of events in the pool
+       ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of calls of pool
 }
-
+       
 
 
 
index 9a4ec13..1bd325f 100644 (file)
@@ -71,13 +71,19 @@ class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE
        void DefineHistoForAnalysis();
        
        // correlators
-       void FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Int_t mcSource);       
+       void FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Bool_t *mcSource);     
+       void FillMCTagLeadingCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Bool_t *mcSource);
+       // checker for event mixing
+       void EventMixingChecks(AliAODEvent * AOD); 
        // setters
        void SetMonteCarlo(Bool_t k) {fmontecarlo = k;}
        void SetUseMixing (Bool_t j) {fmixing = j;}
        void SetCorrelator(Int_t l) {fselect = l;} // select 1 for hadrons, 2 for Kaons, 3 for Kzeros
        void SetUseDisplacement(Int_t m) {fDisplacement=m;} // select 0 for no displ, 1 for abs displ, 2 for d0/sigma_d0
        void SetRunPbPb(Bool_t system){fSystem=system;} // select between pp (kFALSE) or PbPb (kTRUE)
+    void SetLevelOfDebug(Int_t debug){fDebugLevel=debug;} // set debug level
+       void SetUseReconstruction(Bool_t reco){fReco = reco;}
+
        
 
 
@@ -87,7 +93,6 @@ class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE
        AliAnalysisTaskDStarCorrelations& operator=(const AliAnalysisTaskDStarCorrelations& source); 
 
        TObject* fhandler; //! Analysis Handler
-       //      AliEventPoolManager*     fPoolMgr;         //! event pool manager
        TClonesArray* fmcArray; //mcarray
        AliNormalizationCounter *fCounter; // counter
     AliHFCorrelator * fCorrelator; // object for correlations
@@ -98,16 +103,19 @@ class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE
        Bool_t fmontecarlo;//switch for MC
        Bool_t fmixing;// switch for event mixing
        Bool_t fSystem; // pp or PbPb
+       Bool_t fReco; // use reconstruction or MC truth
+       
        Int_t fEvents; //! number of event
-       Int_t fDebug; //! debug level
+       Int_t fDebugLevel; //! debug level
        Int_t fDisplacement; // set 0 for no displacement cut, 1 for absolute d0, 2 for d0/sigma_d0
        
        
        TList *fOutput;                  //! user output data
+       TList *fOutputMC;                //! outpu for MC
        AliRDHFCutsDStartoKpipi *fCuts;  // Cuts D*
        AliHFAssociatedTrackCuts *fCuts2; // cuts for associated 
                                          
-       ClassDef(AliAnalysisTaskDStarCorrelations,2); // class for D meson correlations                           
+       ClassDef(AliAnalysisTaskDStarCorrelations,3); // class for D meson correlations                           
 };
 
 #endif
index 25bff54..caf6294 100644 (file)
@@ -839,7 +839,7 @@ void AliAnalysisTaskSEDplusCorrelations::UserExec(Option_t */*option*/)
 
                index=GetSignalHistoIndex(iPtBin);
 
-               Int_t PartSource = fCuts->IsMCpartFromHF(label,arrayMC); // check source of associated particle (hadron/kaon/K0)
+               Bool_t* PartSource = fCuts->IsMCpartFromHF(label,arrayMC); // check source of associated particle (hadron/kaon/K0)
                FillMCCorrelations(d,DeltaPhi,DeltaEta,index,PartSource,fSelect);   
                
                         
@@ -950,7 +950,7 @@ void AliAnalysisTaskSEDplusCorrelations::FillCorrelations(AliAODRecoDecayHF3Pron
 
 
 //________________________________________________________________________
-void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind,Int_t mcSource, Int_t sel) const{
+void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind,Bool_t* mcSource, Int_t sel) const{
   // Filling histos with MC information
 
   Double_t invMass=d->InvMassDplus();
@@ -963,15 +963,15 @@ void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Pr
        
     fMCSources->Fill(0);
        
-    if (mcSource==44){ // is from charm ->D
+    if(mcSource[2]&&mcSource[0]){ // is from charm ->D
       fMCSources->Fill(1);
       fMCSources->Fill(2);
     }
-    if (mcSource==54){ // is from beauty -> D
+       if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
       fMCSources->Fill(1);
       fMCSources->Fill(3);
     }  
-    if (mcSource==55){ // is from beauty ->B
+    if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
       fMCSources->Fill(1);
       fMCSources->Fill(4);
     }
@@ -981,15 +981,15 @@ void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Pr
     fMassVsdPhiHistKaon[ind]->Fill(invMass,deltaPhi);
     fMassVsdEtaHistKaon[ind]->Fill(invMass,deltaEta);
     fKaonOrigin->Fill(0);
-    if (mcSource==44){ // is from charm ->D
+    if(mcSource[2]&&mcSource[0]){ // is from charm ->D
       fKaonOrigin->Fill(1);
       fKaonOrigin->Fill(2);
     }  
-    if (mcSource==54){ // is from beauty -> D
+    if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
       fKaonOrigin->Fill(1);
       fKaonOrigin->Fill(3);
     }  
-    if (mcSource==55){ // is from beauty ->B
+    if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
       fKaonOrigin->Fill(1);
       fKaonOrigin->Fill(4);
     }
@@ -998,15 +998,15 @@ void AliAnalysisTaskSEDplusCorrelations::FillMCCorrelations(AliAODRecoDecayHF3Pr
     fMassVsdPhiHistKshort[ind]->Fill(invMass,deltaPhi);
     fMassVsdEtaHistKshort[ind]->Fill(invMass,deltaEta);
     fK0Origin->Fill(0);
-    if (mcSource==44){ // is from charm ->D
+    if(mcSource[2]&&mcSource[0]){ // is from charm ->D
       fK0Origin->Fill(1);
       fK0Origin->Fill(2);
     }  
-    if (mcSource==54){ // is from beauty -> D
+    if(mcSource[2]&&mcSource[1]){ // is from beauty -> D
       fK0Origin->Fill(1);
       fK0Origin->Fill(3);
     }
-    if (mcSource==55){ // is from beauty ->B
+    if(mcSource[3]&&mcSource[1]){ // is from beauty ->B
       fK0Origin->Fill(1);
       fK0Origin->Fill(4);
     }
index e743503..c907d94 100644 (file)
@@ -64,7 +64,7 @@ class AliAnalysisTaskSEDplusCorrelations : public AliAnalysisTaskSE
 
   
   void FillCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t sel) const;
-  void FillMCCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Int_t mcSource,Int_t sel) const;
+  void FillMCCorrelations(AliAODRecoDecayHF3Prong* d, Double_t deltaPhi, Double_t deltaEta, Int_t ind, Bool_t* mcSource,Int_t sel) const;
   
   
 
index a2f5b28..adfd78e 100644 (file)
@@ -52,6 +52,9 @@ fNCentBins(0),
 fNCentBinsDim(0), 
 fCentBins(0),
 
+fNofMCEventType(0),
+fMCEventType(0),
+
 fNTrackCuts(0),
 fAODTrackCuts(0),
 fTrackCutsNames(0),
@@ -59,7 +62,8 @@ fNvZeroCuts(0),
 fAODvZeroCuts(0),
 fvZeroCutsNames(0),
 fBit(0),
-fCharge(0) 
+fCharge(0),
+fDescription("")
 
 {
        //
@@ -87,6 +91,9 @@ fNCentBins(0),
 fNCentBinsDim(0), 
 fCentBins(0),
 
+fNofMCEventType(0),
+fMCEventType(0),
+
 fNTrackCuts(0),
 fAODTrackCuts(0),
 fTrackCutsNames(0),
@@ -94,7 +101,8 @@ fNvZeroCuts(0),
 fAODvZeroCuts(0),
 fvZeroCutsNames(0),
 fBit(0),
-fCharge(0)
+fCharge(0),
+fDescription("")
 
 {
        //
@@ -118,6 +126,9 @@ fNCentBins(source.fNCentBins),
 fNCentBinsDim(source.fNCentBinsDim), 
 fCentBins(source.fCentBins),
 
+fNofMCEventType(source.fNofMCEventType),
+fMCEventType(source.fMCEventType),
+
 fNTrackCuts(source.fNTrackCuts),
 fAODTrackCuts(source.fAODTrackCuts),
 fTrackCutsNames(source.fTrackCutsNames),
@@ -125,7 +136,8 @@ fNvZeroCuts(source.fNvZeroCuts),
 fAODvZeroCuts(source.fAODvZeroCuts),
 fvZeroCutsNames(source.fvZeroCutsNames),
 fBit(source.fBit),
-fCharge(source.fCharge)
+fCharge(source.fCharge),
+fDescription(source.fDescription)
 {
        //
        // copy constructor
@@ -290,18 +302,27 @@ Bool_t AliHFAssociatedTrackCuts::IsKZeroSelected(AliAODv0 *vzero, AliAODVertex *
        return kTRUE;
 }
 //--------------------------------------------------------------------------
-Int_t AliHFAssociatedTrackCuts::IsMCpartFromHF(Int_t label, TClonesArray*mcArray){
+Bool_t *AliHFAssociatedTrackCuts::IsMCpartFromHF(Int_t label, TClonesArray*mcArray){
   // Check origin in MC
 
   AliAODMCParticle* mcParticle;
   Int_t pdgCode = -1;
-  if (label<0) return 0;
+       
   Bool_t isCharmy = kFALSE;
   Bool_t isBeauty = kFALSE;
   Bool_t isD = kFALSE;
   Bool_t isB = kFALSE;
-
-  while(pdgCode!=2212){ // loops back to the collision to check the particle source
+    
+     Bool_t *originvect = new Bool_t[4];
+    
+    originvect[0] = kFALSE;
+       originvect[1] = kFALSE;
+       originvect[2] = kFALSE;
+       originvect[3] = kFALSE;
+
+       if (label<0) return originvect;
+  
+       while(pdgCode!=2212){ // loops back to the collision to check the particle source
 
     mcParticle = dynamic_cast<AliAODMCParticle*>(mcArray->At(label));
     if(!mcParticle) {AliError("NO MC PARTICLE"); break;}
@@ -319,10 +340,15 @@ Int_t AliHFAssociatedTrackCuts::IsMCpartFromHF(Int_t label, TClonesArray*mcArray
     if(label<0) break;
 
   }
-  if (isCharmy && isD) return 44; // the first 4(5) indicates the charm/beauty quark, the second the charmy meson
-  if (isBeauty && isD) return 54;
-  if (isBeauty && isB) return 55;
-  return 0;
+
+       
+       originvect[0] = isCharmy;
+       originvect[1] = isBeauty;
+       originvect[2] = isD;
+       originvect[3] = isB;
+    
+  return originvect;
 }
 
 //--------------------------------------------------------------------------
@@ -363,18 +389,25 @@ Bool_t AliHFAssociatedTrackCuts::InvMassDstarRejection(AliAODRecoDecayHF2Prong*
        return kTRUE;
 }
 //________________________________________________________
+void AliHFAssociatedTrackCuts::SetMCEventTypes(Int_t *MCEventTypeArray)
+// set the array of event types you want to process in MonteCarlo (gluon splitting, pair production etc.)
+{
+       if(!fMCEventType) fMCEventType = new Int_t[fNofMCEventType];
+       
+       for(Int_t k=0; k<fNofMCEventType; k++){
+               fMCEventType[k] = MCEventTypeArray[k];
+       }
+       return; 
+}
+
+//________________________________________________________
 void AliHFAssociatedTrackCuts::SetAODTrackCuts(Float_t *cutsarray)
 {
-       cout << "Check 1" << endl;
        if(!fAODTrackCuts) fAODTrackCuts = new Float_t[fNTrackCuts];
-       cout << "Check 2" << endl;
        for(Int_t i =0; i<fNTrackCuts; i++){
-               cout << "Check 2." << i << endl;
                fAODTrackCuts[i] = cutsarray[i];
        }
-       cout << "Check 3" << endl;
        SetTrackCutsNames();
-       cout << "Check 4" << endl;
        return;
 }
 //________________________________________________________
@@ -421,6 +454,8 @@ void AliHFAssociatedTrackCuts::SetvZeroCutsNames(/*TString *namearray*/){
 //--------------------------------------------------------------------------
 void AliHFAssociatedTrackCuts::PrintAll()
 {
+       
+       printf("\n=================================================");
        printf("\nCuts for the associated track: \n \n");
               
        printf("ITS Refit........................................: %s\n",fESDTrackCuts->GetRequireITSRefit() ? "Yes" : "No");
@@ -430,34 +465,42 @@ void AliHFAssociatedTrackCuts::PrintAll()
        printf("Min number of ITS clusters.......................: %d\n",fESDTrackCuts->GetMinNClustersITS());
        printf("Min number of TPC clusters.......................: %d\n",fESDTrackCuts->GetMinNClusterTPC());
        Int_t spd = fESDTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD);
-       if(spd==0) cout <<  "SPD..............................................: kOff"  << endl;
-       if(spd==1) cout <<  "SPD..............................................: kNone"  << endl;
-       if(spd==2) cout <<  "SPD..............................................: kAny"  << endl;
-       if(spd==3) cout <<  "SPD..............................................: kFirst"  << endl;
-       if(spd==4) cout <<  "SPD..............................................: kOnlyFirst"  << endl;
-       if(spd==5) cout <<  "SPD..............................................: kSecond"  << endl;
-       if(spd==6) cout <<  "SPD..............................................: kOnlySecond"  << endl;
-       if(spd==7) cout <<  "SPD..............................................: kBoth"  << endl;
-       
-       cout <<  "Filter Bit.......................................: " << fBit  << endl;
-       cout <<  "Charge...........................................: " << fCharge  << endl;
+       if(spd==0) std::cout <<  "SPD..............................................: kOff"  << std::endl;
+       if(spd==1) std::cout <<  "SPD..............................................: kNone"  << std::endl;
+       if(spd==2) std::cout <<  "SPD..............................................: kAny"  << std::endl;
+       if(spd==3) std::cout <<  "SPD..............................................: kFirst"  << std::endl;
+       if(spd==4) std::cout <<  "SPD..............................................: kOnlyFirst"  << std::endl;
+       if(spd==5) std::cout <<  "SPD..............................................: kSecond"  << std::endl;
+       if(spd==6) std::cout <<  "SPD..............................................: kOnlySecond"  << std::endl;
+       if(spd==7) std::cout <<  "SPD..............................................: kBoth"  << std::endl;
+       
+       std::cout <<  "Filter Bit.......................................: " << fBit  << std::endl;
+       std::cout <<  "Charge...........................................: " << fCharge  << std::endl;
                
        for(Int_t j=0;j<fNTrackCuts;j++){
-               cout << fTrackCutsNames[j] << fAODTrackCuts[j] << endl;
+               std::cout << fTrackCutsNames[j] << fAODTrackCuts[j] << std::endl;
        }
-       
+       printf("\n");
+       printf("=================================================");
        printf("\nCuts for the K0 candidates: \n \n");
        for(Int_t k=0;k<fNvZeroCuts;k++){
-               cout << fvZeroCutsNames[k] <<  fAODvZeroCuts[k] << endl;
+               std::cout << fvZeroCutsNames[k] <<  fAODvZeroCuts[k] << std::endl;
        }
-       cout << " " << endl;
+       std::cout << " " << std::endl;
        PrintPoolParameters();
+       PrintSelectedMCevents();
+
+       printf("=================================================");
+       printf("\nAdditional description\n");
+       std::cout << fDescription << std::endl;
+       printf("\n");
 
 }
 
 //--------------------------------------------------------------------------
 void AliHFAssociatedTrackCuts::PrintPoolParameters()
-{
+{   
+       printf("=================================================");
        printf("\nEvent Pool settings: \n \n");
        
        printf("Number of zVtx Bins: %d\n", fNzVtxBins);
@@ -466,15 +509,41 @@ void AliHFAssociatedTrackCuts::PrintPoolParameters()
        for(Int_t k=0; k<fNzVtxBins; k++){
                printf("Bin %d..............................................: %.1f - %.1f cm\n", k, fZvtxBins[k], fZvtxBins[k+1]);      
        }
-       
+       printf("\n");
        printf("\nNumber of Centrality(multiplicity) Bins: %d\n", fNCentBins);
        printf("\nCentrality(multiplicity) Bins:\n");
        for(Int_t k=0; k<fNCentBins; k++){
                printf("Bin %d..............................................: %.1f - %.1f\n", k, fCentBins[k], fCentBins[k+1]);
        }
+
+       
+       
+}
+
+//--------------------------------------------------------------------------
+void AliHFAssociatedTrackCuts::PrintSelectedMCevents()
+{
+       printf("\n=================================================");
+       
+       printf("\nSelected MC events: \n \n");
+       printf("Number of selected events: %d\n",fNofMCEventType);
+       
+       for(Int_t k=0; k<fNofMCEventType; k++){
+       if(fMCEventType[k]==28) printf("=> Flavour excitation \n");     
+       if(fMCEventType[k]==53) printf("=> Pair creation \n");  
+       if(fMCEventType[k]==68) printf("=> Gluon splitting \n");        
+       }
        
+       printf("\n");
+       for(Int_t k=0; k<fNofMCEventType; k++){
+               printf("MC process code %d \n",fMCEventType[k]);                
+       }
+       
+       printf("\n");
        
        
+               
+       
 }
 
 
index 62496a6..f9cf850 100644 (file)
@@ -32,6 +32,8 @@
 #include "AliAODRecoDecayHF2Prong.h"
 #include <TClonesArray.h>
 
+
+
 class AliAODTrack;
 class AliAODEvent;
 
@@ -56,7 +58,7 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
        Bool_t Charge(Short_t charge, AliAODTrack* track);
        Bool_t CheckKaonCompatibility(AliAODTrack * track, Bool_t useMc, TClonesArray* mcArray, Int_t method=1);
        Bool_t IsKZeroSelected(AliAODv0 *vzero, AliAODVertex *vtx1);
-       Int_t IsMCpartFromHF(Int_t label, TClonesArray*mcArray);
+       Bool_t *IsMCpartFromHF(Int_t label, TClonesArray*mcArray);
        Bool_t InvMassDstarRejection(AliAODRecoDecayHF2Prong* d, AliAODTrack *track, Int_t hypD0) const;        
        
        //getters
@@ -68,6 +70,11 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
        Int_t GetNCentPoolBins() {return fNCentBins;}
        Double_t *GetCentPoolBins(){return fCentBins;}
        
+       Int_t GetNofMCEventType(){return fNofMCEventType;}      
+       Int_t *GetMCEventType(){return fMCEventType;}
+       TString GetDescription(){return fDescription;}
+       Int_t GetNVarsTrack(){return fNTrackCuts;}
+       
        
        
        void AddTrackCuts(const AliESDtrackCuts *cuts) {
@@ -75,6 +82,7 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
                fESDTrackCuts=new AliESDtrackCuts(*cuts); 
                return;
        }
+       
        //setters
        //event pool settings
        void SetMaxNEventsInPool(Int_t events){fPoolMaxNEvents=events;}
@@ -93,6 +101,12 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
                fZvtxBins=ZvtxBins; 
                fCentBins=CentBins;
        }
+       
+       // set MC events to process
+       
+       void SetNofMCEventTypes(Int_t k) {fNofMCEventType=k;}
+       void SetMCEventTypes(Int_t *MCEventTypeArray);
+       
        //cut settings
        void SetAODTrackCuts(Float_t *cutsarray);
        void SetTrackCutsNames(/*TString *namearray*/);
@@ -103,15 +117,17 @@ class AliHFAssociatedTrackCuts : public AliAnalysisCuts
        void SetFilterBit(Int_t bit) {fBit = bit;}
        virtual void PrintAll();
        virtual void PrintPoolParameters();
-       Int_t GetNVarsTrack(){return fNTrackCuts;}
+       virtual void PrintSelectedMCevents();
+
        
        
        
        void SetNVarsTrack(Int_t nVars){fNTrackCuts=nVars;}
        void SetNVarsVzero(Int_t nVars){fNvZeroCuts=nVars;}
        
+       
+       
 private:
-       //AliESDtrackCuts *fTrackCuts;
        AliESDtrackCuts *fESDTrackCuts; // track cut object
        AliAODPidHF * fPidObj;     /// PID object
        
@@ -128,6 +144,9 @@ private:
        Int_t fNCentBinsDim; //number of centrality bins bins +1 : necessary to initialize correctly the array
        Double_t* fCentBins; // [fNCentBinsDim]
        
+       Int_t fNofMCEventType;// number of event types to be selected in MC simultaneously;
+       Int_t *fMCEventType;//[fNofMCEventType]
+       
        Int_t fNTrackCuts;     // array dimension
        Float_t* fAODTrackCuts;//[fNTrackCuts]
        TString * fTrackCutsNames;//[fNTrackCuts]
@@ -136,9 +155,10 @@ private:
        TString * fvZeroCutsNames;//[fNvZeroCuts]
        Int_t fBit; // filterBit
        Short_t fCharge; // charge (+1 or -1)
+       TString fDescription; // additional description to the cuts
        
        
-       ClassDef(AliHFAssociatedTrackCuts,3);
+       ClassDef(AliHFAssociatedTrackCuts,4);
 };
 
 
index 4887272..5694de5 100644 (file)
@@ -38,6 +38,7 @@
 #include "AliEventPoolManager.h"
 #include "AliReducedParticle.h"
 #include "AliCentrality.h"
+#include "AliAODMCParticle.h"
 
 using std::cout;
 using std::endl;
@@ -62,7 +63,9 @@ fDCharge(0),
 fmixing(kFALSE),
 fmontecarlo(kFALSE),
 fsystem(kFALSE),
+fUseReco(kTRUE),
 fselect(0),
+
 fUseImpactParameter(0),
 fPIDmode(0),
 
@@ -104,6 +107,7 @@ fDCharge(0),
 fmixing(kFALSE),
 fmontecarlo(kFALSE),
 fsystem(ppOrPbPb),
+fUseReco(kTRUE),
 fselect(0),
 fUseImpactParameter(0),
 fPIDmode(0),
@@ -322,42 +326,69 @@ TObjArray*  AliHFCorrelator::AcceptAndReduceTracks(AliAODEvent* inputEvent){
        
        TObjArray* tracksClone = new TObjArray;
        tracksClone->SetOwner(kTRUE);
-       for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
-               
-               AliAODTrack* track = inputEvent->GetTrack(iTrack);
-               if (!track) continue;
-               if(!fhadcuts->IsHadronSelected(track)) continue; // apply ESD level selections
-        if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
+       
+       //*******************************************************
+    // use reconstruction
+       if(fUseReco){
+               for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) {
+                       AliAODTrack* track = inputEvent->GetTrack(iTrack);
+                       if (!track) continue;
+                       if(!fhadcuts->IsHadronSelected(track)) continue; // apply ESD level selections
+                       if(!fhadcuts->Charge(fDCharge,track)) continue; // apply selection on charge, if required
                
-               Double_t pT = track->Pt();
+                       Double_t pT = track->Pt();
                
                //compute impact parameter
-               Double_t d0z0[2],covd0z0[3];
-               Double_t d0=-999999.;
-               if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
-               else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
+                       Double_t d0z0[2],covd0z0[3];
+                       Double_t d0=-999999.;
+                       if(fUseImpactParameter) track->PropagateToDCA(vtx,Bz,100,d0z0,covd0z0);
+                       else d0z0[0] = 1. ; // random number - be careful with the cuts you applied
                
-               if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
-               if(fUseImpactParameter==2) { // use impact parameter over resolution
+                       if(fUseImpactParameter==1) d0 = TMath::Abs(d0z0[0]); // use impact parameter
+                       if(fUseImpactParameter==2) { // use impact parameter over resolution
                  if(TMath::Abs(covd0z0[0])>0.00000001) d0 = TMath::Abs(d0z0[0])/TMath::Sqrt(covd0z0[0]); 
                  else d0 = -1.; // if the resoultion is Zero, rejects the track - to be on the safe side
                        
                }
                
-        
+                       
+               if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
+                       Bool_t rejectsoftpi = kFALSE;
+                       if(fD0cand) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0);
                
+                       
+                       if(fselect ==2){        
+                               if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
+                       }
                
-               if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
-               Bool_t rejectsoftpi = kFALSE;
-               if(fD0cand) rejectsoftpi = fhadcuts->InvMassDstarRejection(fD0cand,track,fhypD0);
+               tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge()));
+               } // end loop on tracks
+       } // end if use reconstruction kTRUE
+       
+       //*******************************************************
+       
+       //use MC truth
+       if(fmontecarlo && !fUseReco){
                
+               for (Int_t iPart=0; iPart<fmcArray->GetEntriesFast(); iPart++) { 
+                       AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iPart));
+                       if (!mcPart) {
+                               AliWarning("MC Particle not found in tree, skipping"); 
+                               continue;
+                       }
+                       if(!mcPart->Charge()) continue; // consider only charged tracks
+                       
+                       Int_t PDG =TMath::Abs(mcPart->PdgCode()); 
+                       if(!((PDG==321)||(PDG==211)||(PDG==2212)||(PDG==13)||(PDG==11))) continue; // select only if kaon, pion, proton, muon or electron
                
-               if(fselect ==2){        
-                       if(!fhadcuts->CheckKaonCompatibility(track,fmontecarlo,fmcArray,fPIDmode)) continue; // check if it is a Kaon - data and MC
+                       Double_t pT = mcPart->Pt();
+            Double_t d0 =1; // set 1 fot the moment - no displacement calculation implemented yet
+                       if(!fhadcuts->CheckHadronKinematic(pT,d0)) continue; // apply kinematic cuts
+                       
+                       tracksClone->Add(new AliReducedParticle(mcPart->Eta(), mcPart->Phi(), pT,mcPart->GetLabel(),-1,d0,kFALSE,mcPart->Charge()));
                }
                
-               tracksClone->Add(new AliReducedParticle(track->Eta(), track->Phi(), pT,track->GetLabel(),track->GetID(),d0,rejectsoftpi,track->Charge()));
-       }
+       } // end if use  MC truth
        
        
        return tracksClone;
index 35c521b..e7c2476 100644 (file)
@@ -66,6 +66,8 @@ class AliHFCorrelator : public TNamed
        void SetD0Properties(AliAODRecoDecayHF2Prong* d, Int_t D0hyp)
        {fD0cand = d; fhypD0 = D0hyp;}
        
+       void SetUseReco(Bool_t useReco) {fUseReco = useReco;}
+       
        
        
     Bool_t DefineEventPool(); // Definition of the Event pool parameters
@@ -116,6 +118,8 @@ class AliHFCorrelator : public TNamed
        Bool_t fmixing;// switch for event mixing
        Bool_t fmontecarlo; // switch for MonteCarlo
        Bool_t fsystem; // select pp (kFALSE) or PbPb (kTRUE)
+       Bool_t fUseReco; // switch to use reconstruction (kTRUE) or MC truth (kFALSE)
+       
        Int_t fselect; // 1 for hadrons, 2 for kaons, 3 for KZeros
        Int_t fUseImpactParameter; // switch to use the impact parameter cut
        Int_t fPIDmode; // set the PID mode for Kaon identification
@@ -137,7 +141,7 @@ class AliHFCorrelator : public TNamed
        Double_t fk0InvMass; // KZero invariant mass
        
        
-       ClassDef(AliHFCorrelator,2); // class for HF correlations       
+       ClassDef(AliHFCorrelator,3); // class for HF correlations       
 };
 
 
index a0e3680..e3d37ca 100644 (file)
@@ -3,7 +3,7 @@
 
 /* $Id$ */
 
-AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool_t theMCon, Bool_t mixing, Int_t trackselect =1, Int_t usedispl =0, TString DCutObjName = "DStartoKpipiCuts_corr.root", TString TrackCutObjNamePrefix = "")
+AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool_t theMCon, Bool_t mixing, Bool_t UseReco = kTRUE, Int_t trackselect =1, Int_t usedispl =0, TString DCutObjNamePrefix = "_corr", TString TrackCutObjNamePrefix = "")
 {
 
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
@@ -12,19 +12,36 @@ AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool
     return NULL;
   } 
 
+       
+       TString DCutObjPath = "~/CorrelationAnalysis/CutObjects/DStar/";
+       
+       
+       
+       TString DCutObjName = "DStartoKpipiCuts";
+       DCutObjName += DCutObjNamePrefix;
+       DCutObjName += ".root";
+       
+       DCutObjName.Prepend(DCutObjPath.Data());
+       
+       cout << "D* cut object is " << DCutObjName << endl;
   TFile* filecuts=new TFile(DCutObjName.Data());
   if(!filecuts->IsOpen()){
-    cout<<"Input file not found: exit"<<endl;
+    cout<<"DStar cut object file not found: exit"<<endl;
     return;
   }  
        
+       TString TrackCutObjPath = "~/CorrelationAnalysis/CutObjects/AssocTracks/";
+       
        TString TrackCutObjName = "AssocPartCuts";
        TrackCutObjName += TrackCutObjNamePrefix;
        TrackCutObjName += ".root";
        
+       TrackCutObjName.Prepend(TrackCutObjPath.Data());
+       
+       cout << "tracks cut object is " << TrackCutObjName << endl;
          TFile* filecuts2=new TFile(TrackCutObjName.Data());
          if(!filecuts2->IsOpen()){
-                 cout<<"Input file2 not found: exit"<<endl;
+                 cout<<"Track cut object file not found: exit"<<endl;
                  return;
   }
 
@@ -37,26 +54,49 @@ AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool
        corrCuts = (AliHFAssociatedTrackCuts*)filecuts2->Get("AssociatedCuts");
        corrCuts->SetName("AssociatedCuts");
        corrCuts->PrintAll();
+       
+       // mm let's see if everything is ok
+       if(!RDHFDStartoKpipi){
+               cout<<"Specific AliRDHFCuts not found"<<endl;
+               return;
+       } 
+       
+       if(!corrCuts){
+               cout<<"Specific associated track cuts not found"<<endl;
+               return;
+       } 
+       
+       TString selectMCproc = "";
+       
+       Int_t NMCevents = corrCuts->GetNofMCEventType();
+       for(Int_t k=0; k<NMCevents; k++){
+               Int_t * MCEventType = corrCuts->GetMCEventType();
+               selectMCproc += Form("%d",MCEventType[k]);
+       }
+
+       cout << "Select process string = " << selectMCproc << endl;
+       
 
-  // mm let's see if everything is ok
-  if(!RDHFDStartoKpipi){
-    cout<<"Specific AliRDHFCuts not found"<<endl;
-    return;
-  } 
 
   //CREATE THE TASK
   printf("CREATE TASK \n");
   // create the task
-  AliAnalysisTaskDStarCorrelations *task = new AliAnalysisTaskDStarCorrelations(" AliAnalysisTaskDStarCorrelations",RDHFDStartoKpipi,corrCuts);
+  AliAnalysisTaskDStarCorrelations *task = new AliAnalysisTaskDStarCorrelations("AliAnalysisTaskDStarCorrelations",RDHFDStartoKpipi,corrCuts);
        
        // Setters
 
+       if(!theMCon) {
+               printf("Analysis on Data - reconstruction only!");
+               UseReco = kTRUE;
+       }
+       
        task->SetMonteCarlo(theMCon);
        task->SetUseMixing(mixing);
        task->SetCorrelator(trackselect) ;
        task->SetUseDisplacement(usedispl);
        task->SetRunPbPb(runOnPbPb);
-       //task->SetDebugLevel(0);
+       task->SetLevelOfDebug(2);
+       task->SetUseReconstruction(UseReco); // set kTRUE for Using Reconstruction, kFALSe for MC Truth
        
 
        if(trackselect == 1) Info(" AliAnalysisTaskDStarCorrelations","Correlating D* with charged hadrons \n");
@@ -70,45 +110,95 @@ AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool
        //TString dcavalue = " ";
        if(!theMCon) TString contname = "Data";
        if(theMCon) TString contname = "MonteCarlo";
+       TString contname2 = "MC";
        if(trackselect ==1) TString particle = "Hadron";
        if(trackselect ==2) TString particle = "Kaon";
        if(trackselect ==3) TString particle = "KZero";
        
        TString cutname = "cuts" ;
        TString cutname2 = "hadroncuts" ;
-   TString outputfile = AliAnalysisManager::GetCommonFileName();
+       TString outputfile = AliAnalysisManager::GetCommonFileName();
+       TString outputfileMC = AliAnalysisManager::GetCommonFileName();
        TString counter = "NormCounter";
-   outputfile += ":PWGHF_D2H_";
+       outputfile += ":PWGHF_D2H_";
+       outputfileMC += ":PWGHF_D2H_";
+       
        if(!mixing) {
                outputfile += "SE";
+               outputfileMC += "SE";
                contname += "SE";
+               contname2 += "SE";
                cutname += "SE";
                cutname2 += "SE";
                counter+= "SE";
        }
        if(mixing){
                outputfile += "ME";
+               outputfileMC += "ME";
                contname += "ME";
+               contname2 += "ME";
                cutname += "ME";
                cutname2 += "ME";
                counter+= "ME";
        }
        outputfile += "Dphi_DStar";
+       outputfileMC += "Dphi_DStar";
        outputfile += particle;
+       outputfileMC += particle;
        cutname += particle;
        cutname2 += particle;
        contname += particle;
+       contname2 += particle;
        counter+= particle;
        
+       
+       
+       outputfile += DCutObjNamePrefix;
+       outputfileMC += DCutObjNamePrefix;
+       cutname += DCutObjNamePrefix;
+       cutname2 += DCutObjNamePrefix;
+       contname += DCutObjNamePrefix;
+       contname2 += DCutObjNamePrefix;
+       counter+= DCutObjNamePrefix;
+       
        outputfile += TrackCutObjNamePrefix;
+       outputfileMC += TrackCutObjNamePrefix;
        cutname += TrackCutObjNamePrefix;
        cutname2 += TrackCutObjNamePrefix;
        contname += TrackCutObjNamePrefix;
+       contname2 += TrackCutObjNamePrefix;
        counter+= TrackCutObjNamePrefix;
        
-
+       outputfile += selectMCproc;
+       outputfileMC += selectMCproc;
+       cutname += selectMCproc;
+       cutname2 += selectMCproc;
+       contname += selectMCproc;
+       contname2 += selectMCproc;
+       counter+= selectMCproc;
+       
+       TString reco = "";
+       
+       if(UseReco) reco = "_reco";
+       if(!UseReco) reco = "_MCTruth";
+       
+       outputfile += reco;
+       outputfileMC += reco;
+       cutname += reco;
+       cutname2 += reco;
+       contname += reco;
+       contname2 += reco;
+       counter+= reco;
+       
+       
+       cout << contname << endl;
+       cout << contname2 << endl;
+       cout << cutname << endl;
+       cout << cutname2 << endl;
+       cout << counter << endl;
+       cout << outputfile << endl;
+       //return;
        
-       cout << "Contname = " << contname << endl;
   mgr->AddTask(task);
   // ------ input data ------
   AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
@@ -116,16 +206,23 @@ AliAnalysisTaskDStarCorrelations *AddTaskDStarCorrelations(Bool_t runOnPbPb,Bool
   // ----- output data -----
   
   // output TH1I for event counting
+       
+       //TLists
   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(cutname,AliRDHFCutsDStartoKpipi::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts
-  AliAnalysisDataContainer *coutputDstarNorm = mgr->CreateContainer(counter,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
-  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(cutname2,AliHFAssociatedTrackCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(contname2, TList::Class(),AliAnalysisManager::kOutputContainer,outputfileMC.Data());
+   // Cut Objects
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(cutname,AliRDHFCutsDStartoKpipi::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts D
+  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(cutname2,AliHFAssociatedTrackCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); //cuts tracks
+   // Normalization Counter
+  AliAnalysisDataContainer *coutput5 = mgr->CreateContainer(counter,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data());
   
+       
   mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
   mgr->ConnectOutput(task,1,coutput1);
   mgr->ConnectOutput(task,2,coutput2);
-  mgr->ConnectOutput(task,3,coutputDstarNorm);
+  mgr->ConnectOutput(task,3,coutput3);
   mgr->ConnectOutput(task,4,coutput4);
+  mgr->ConnectOutput(task,5,coutput5);
 
   return task ;
 
index 749aecc..e540bb9 100644 (file)
@@ -7,7 +7,7 @@
 
 /* $Id$ */
 
-void makeInputHFCorrelation(){
+void makeInputHFCorrelation(TString suffix = ""){
        
        AliHFAssociatedTrackCuts* HFCorrelationCuts=new AliHFAssociatedTrackCuts();
        HFCorrelationCuts->SetName("AssociatedCuts");
@@ -46,6 +46,13 @@ void makeInputHFCorrelation(){
        
        HFCorrelationCuts->SetPoolBins(ZVrtxBins,MultiplicityBins);
        
+       //______________________________ MC event type
+       
+       
+       HFCorrelationCuts->SetNofMCEventTypes(2);
+       Int_t MCEvType[]={53,68};
+       Int_t *MCEvTypeArray = MCEvType;
+       HFCorrelationCuts->SetMCEventTypes(MCEvTypeArray);
        
        //______________________________ set kinematics cuts for AOD track 
        const int nofcuts = 4;
@@ -91,12 +98,29 @@ void makeInputHFCorrelation(){
        
        HFCorrelationCuts->SetPidHF(pidObj);
 
+        //______________________________ Add additional info
+       
+       TString description = "Cuts for hadron selection - low pt selection cut";
+       
+       HFCorrelationCuts->AddDescription(description);
     //______________________________ save to *.root file
        HFCorrelationCuts->PrintAll();
        
-       TFile* fout=new TFile("AssocPartCuts.root","recreate");   //set this!! 
+       TString outputfilename = "AssocPartCuts";
+       outputfilename += suffix;
+       outputfilename += ".root";
+       
+       TString pathOutputFile = "~/CorrelationAnalysis/CutObjects/AssocTracks/";
+       
+       outputfilename.Prepend(pathOutputFile.Data());
+       
+       
+       TFile* fout=new TFile(outputfilename.Data(),"recreate");   //set this!! 
        fout->cd();
        HFCorrelationCuts->Write();
        fout->Close();
+       cout << "***************************************** " << endl;
+       cout << "File " << outputfilename << "  created " << endl;
+       cout << "***************************************** " << endl;
 
 }