PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //
  16 //
  17 //             Base class for DStar - Hadron Correlations Analysis
  18 //
  19 //-----------------------------------------------------------------------
  20 //
  21 //
  22 //                                                 Author S.Bjelogrlic
  23 //                         Utrecht University
  24 //                      sandro.bjelogrlic@cern.ch
  25 //
  26 //-----------------------------------------------------------------------
  27
  28 /* $Id$ */
  29
  30 //#include <TDatabasePDG.h>
  31 #include <TParticle.h>
  32 #include <TVector3.h>
  33 #include <TChain.h>
  34 #include "TROOT.h"
  35
  36 #include "AliAnalysisTaskDStarCorrelations.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliHFAssociatedTrackCuts.h"
  39 #include "AliAODRecoDecay.h"
  40 #include "AliAODRecoCascadeHF.h"
  41 #include "AliAODRecoDecayHF2Prong.h"
  42 #include "AliAODPidHF.h"
  43 #include "AliVParticle.h"
  44 #include "AliAnalysisManager.h"
  45 #include "AliAODInputHandler.h"
  46 #include "AliAODHandler.h"
  47 #include "AliESDtrack.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliNormalizationCounter.h"
  50 #include "AliReducedParticle.h"
  51 #include "AliHFCorrelator.h"
  52 #include "AliAODMCHeader.h"
  53 #include "AliEventPoolManager.h"
  54
  55 using std::cout;
  56 using std::endl;
  57
  58
  59 ClassImp(AliAnalysisTaskDStarCorrelations)
  60
  61
  62 //__________________________________________________________________________
  63 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :
  64 AliAnalysisTaskSE(),
  65 fhandler(0x0),
  66 fmcArray(0x0),
  67 fCounter(0x0),
  68 fCorrelator(0x0),
  69 fselect(0),
  70 fmontecarlo(kFALSE),
  71 fmixing(kFALSE),
  72 fSystem(kFALSE),
  73 fReco(kTRUE),
  74 fEvents(0),
  75 fDebugLevel(0),
  76 fDisplacement(0),
  77
  78 fOutput(0x0),
  79 fOutputMC(0x0),
  80 fCuts(0),
  81 fCuts2(0)
  82 {
  83 // default constructor
  84 }
  85
  86 //__________________________________________________________________________
  87 AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts) :
  88 AliAnalysisTaskSE(name),
  89
  90 fhandler(0x0),
  91 fmcArray(0x0),
  92 fCounter(0x0),
  93 fCorrelator(0x0),
  94 fselect(0),
  95 fmontecarlo(kFALSE),
  96 fmixing(kFALSE),
  97 fSystem(kFALSE),
  98 fReco(kTRUE),
  99 fEvents(0),
 100 fDebugLevel(0),
 101 fDisplacement(0),
 102
 103 fOutput(0x0),
 104 fOutputMC(0x0),
 105 fCuts(0),
 106 fCuts2(AsscCuts)
 107 {
 108         fCuts=cuts;
 109         Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
 110         DefineInput(0, TChain::Class());
 111
 112         DefineOutput(1,TList::Class()); // histos from data and MC
 113         DefineOutput(2,TList::Class()); // histos from MC
 114         DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
 115         DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
 116         DefineOutput(5,AliNormalizationCounter::Class());   // normalization
 117 }
 118
 119 //__________________________________________________________________________
 120
 121 AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {
 122         //
 123         // destructor
 124         //
 125
 126         Info("AliAnalysisTaskDStarCorrelations","Calling Destructor");
 127
 128         if(fhandler) {delete fhandler; fhandler = 0;}
 129         //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;}
 130         if(fmcArray) {delete fmcArray; fmcArray = 0;}
 131         if(fCounter) {delete fCounter; fCounter = 0;}
 132         if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}
 133         if(fOutput) {delete fOutput; fOutput = 0;}
 134         if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}
 135         if(fCuts) {delete fCuts; fCuts = 0;}
 136         if(fCuts2) {delete fCuts2; fCuts2=0;}
 137
 138 }
 139
 140 //___________________________________________________________
 141 void AliAnalysisTaskDStarCorrelations::Init(){
 142         //
 143         // Initialization
 144         //
 145         if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");
 146
 147         AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
 148
 149
 150
 151
 152         // Post the D* cuts
 153         PostData(3,copyfCuts);
 154
 155         // Post the hadron cuts
 156         PostData(4,fCuts2);
 157
 158
 159
 160         return;
 161 }
 162
 163
 164 //_________________________________________________
 165 void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){
 166         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 167
 168         //slot #1
 169         //OpenFile(0);
 170         fOutput = new TList();
 171         fOutput->SetOwner();
 172
 173         fOutputMC = new TList();
 174         fOutputMC->SetOwner();
 175
 176         // define histograms
 177         DefineHistoForAnalysis();
 178         fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
 179         fCounter->Init();
 180
 181     Double_t Pi = TMath::Pi();
 182         fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fSystem); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb
 183         fCorrelator->SetDeltaPhiInterval((-0.5-1./32)*Pi,(1.5-1./32)*Pi); // set correct phi interval
 184         fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On
 185         fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros
 186         fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut
 187         fCorrelator->SetUseMC(fmontecarlo);
 188         fCorrelator->SetUseReco(fReco);
 189         Bool_t pooldef = fCorrelator->DefineEventPool();
 190
 191         if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
 192
 193
 194
 195         PostData(1,fOutput); // set the outputs
 196         PostData(2,fOutputMC); // set the outputs
 197         PostData(5,fCounter); // set the outputs
 198 }
 199 //_________________________________________________
 200 void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
 201
 202
 203         if(fDebugLevel){
 204
 205                 if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;
 206                 if(!fReco) std::cout << "USING MC TRUTH" << std::endl;
 207         std::cout << " " << std::endl;
 208         std::cout << "=================================================================================" << std::endl;
 209         if(!fmixing){
 210             if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;
 211             if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;
 212             if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;
 213         }
 214         if(fmixing){
 215             if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;
 216             if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;
 217             if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;
 218         }
 219
 220     }// end if debug
 221
 222         if (!fInputEvent) {
 223                 Error("UserExec","NO EVENT FOUND!");
 224                 return;
 225         }
 226
 227         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 228         if(!aodEvent){
 229           AliError("AOD event not found!");
 230           return;
 231         }
 232
 233
 234
 235         fEvents++; // event counter
 236         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
 237     fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);
 238
 239         // load MC array
 240         fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 241         if(fmontecarlo && !fmcArray){
 242           AliError("Array of MC particles not found");
 243           return;
 244         }
 245
 246
 247
 248
 249         Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
 250         if(!isEvSel) return;
 251
 252     fCorrelator->SetAODEvent(aodEvent); // set the event to be processed
 253
 254     ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);
 255         //
 256         Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
 257         if(!correlatorON) {
 258                 AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
 259                 return;
 260         }
 261         ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);
 262
 263         if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);
 264
 265
 266         // check the event type
 267         // load MC header
 268
 269         if(fmontecarlo){
 270                 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 271                 if (fmontecarlo && !mcHeader) {
 272                         AliError("Could not find MC Header in AOD");
 273                         return;
 274                 }
 275
 276                 Bool_t isMCeventgood = kFALSE;
 277
 278
 279                 Int_t eventType = mcHeader->GetEventType();
 280                 Int_t NMCevents = fCuts2->GetNofMCEventType();
 281
 282                 for(Int_t k=0; k<NMCevents; k++){
 283                         Int_t * MCEventType = fCuts2->GetMCEventType();
 284
 285                         if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
 286                         ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);
 287                 }
 288
 289                 if(NMCevents && !isMCeventgood){
 290                 if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
 291                         return;
 292                 }
 293
 294         } // end if montecarlo
 295
 296         // D* reconstruction
 297         TClonesArray *arrayDStartoD0pi=0;
 298         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 299                 // In case there is an AOD handler writing a standard AOD, use the AOD
 300                 // event in memory rather than the input (ESD) event.
 301                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 302                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 303                 // have to taken from the AOD event hold by the AliAODExtension
 304                 AliAODHandler* aodHandler = (AliAODHandler*)
 305                 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 306                 if(aodHandler->GetExtensions()) {
 307                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 308                         AliAODEvent *aodFromExt = ext->GetAOD();
 309                         arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 310                 }
 311         } else {
 312                 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 313         }
 314
 315         if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 316
 317         // initialize variables you will need for the D*
 318
 319         Double_t ptDStar;//
 320         Double_t phiDStar;//
 321         Double_t etaDStar;//
 322         Bool_t isInPeak, isInSideBand, isDStarMCtag;
 323         Double_t invMassDZero;
 324         Double_t deltainvMDStar;
 325
 326
 327         Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();
 328         Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();
 329
 330
 331         //MC tagging for DStar
 332         //D* and D0 prongs needed to MatchToMC method
 333         Int_t pdgDgDStartoD0pi[2]={421,211};
 334         Int_t pdgDgD0toKpi[2]={321,211};
 335
 336         Bool_t isDStarCand = kFALSE;
 337         //loop on D* candidates
 338         for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
 339                 isInPeak = kFALSE;
 340                 isInSideBand = kFALSE;
 341                 isDStarMCtag = kFALSE;
 342                 ptDStar = -123.4;
 343                 phiDStar = -999;
 344                 etaDStar = -56.;
 345                 invMassDZero = - 999;
 346                 deltainvMDStar = -998;
 347
 348                 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 349                 if(!dstarD0pi->GetSecondaryVtx()) continue;
 350                 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 351                 if (!theD0particle) continue;
 352
 353
 354                 // track quality cuts
 355                 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
 356                 // region of interest + topological cuts + PID
 357                 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
 358                 //apply selections
 359                 if(!isTkSelected) continue;
 360                 if(!isSelected) continue;
 361                 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
 362                 Int_t mcLabelDStar = -999;
 363                 if(fmontecarlo){
 364                         // find associated MC particle for D* ->D0toKpi
 365                         mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray,kFALSE);
 366                         if(mcLabelDStar>=0) isDStarMCtag = kTRUE;
 367                 }
 368
 369                 ptDStar = dstarD0pi->Pt();
 370                 phiDStar = dstarD0pi->Phi();
 371                 etaDStar = dstarD0pi->Eta();
 372
 373                 phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the phi of the D meson in the correct range
 374
 375                 Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
 376
 377                 Double_t dmDStarWindow =0.0019;// 0.0019 = 3 sigma
 378                 Double_t mD0Window=0.074;
 379
 380                 if (!fSystem){ // pp
 381                         if (ptbin==1)  mD0Window = 0.026; //0.5-1
 382                         if (ptbin==2)  mD0Window = 0.022; //1-2
 383                         if (ptbin==3)  mD0Window = 0.024; //2-3
 384                         if (ptbin==4)  mD0Window = 0.032;
 385                         if (ptbin==5)  mD0Window = 0.032;
 386                         if (ptbin==6)  mD0Window = 0.036;
 387                         if (ptbin==7)  mD0Window = 0.036;
 388                         if (ptbin==8)  mD0Window = 0.036;
 389                         if (ptbin==9)  mD0Window = 0.058;
 390                         if (ptbin==10) mD0Window = 0.058;
 391                         if (ptbin>10)  mD0Window = 0.074;
 392                 }
 393                 if(fSystem){// PbPb
 394                         if (ptbin==0)  mD0Window = 0.032; //1-1
 395                         if (ptbin==1)  mD0Window = 0.032; //2-3
 396                         if (ptbin==2)  mD0Window = 0.032; //3-4
 397                         if (ptbin==3)  mD0Window = 0.032; //4-5
 398                         if (ptbin==4)  mD0Window = 0.036; //5-6
 399                         if (ptbin==5)  mD0Window = 0.036; //6-8
 400                         if (ptbin==6)  mD0Window = 0.055; //8-12
 401                         if (ptbin==7)  mD0Window = 0.074; //12-16
 402                         if (ptbin==8)  mD0Window = 0.074; //16-24
 403                         if (ptbin==9)  mD0Window = 0.074; //24-35
 404                  }
 405
 406                 invMassDZero = dstarD0pi->InvMassD0();
 407                 ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);
 408
 409                 deltainvMDStar = dstarD0pi->DeltaInvMass();
 410
 411
 412                 //good candidates
 413                 if (TMath::Abs(invMassDZero-mPDGD0)<mD0Window){
 414
 415                         ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
 416                         if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
 417
 418                                 ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar);
 419                                 isInPeak = kTRUE;
 420                                 ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);
 421                         }
 422                 }// end if good candidates
 423
 424                 //sidebands
 425                 if (TMath::Abs(invMassDZero-mPDGD0)>1.3*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<4.*mD0Window ){
 426                         ((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar);
 427                         ((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero);
 428
 429                         if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<dmDStarWindow){ // is in DStar peak region?
 430                                 ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar);
 431                                 isInSideBand = kTRUE;
 432                                 ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);
 433                         }
 434
 435                 }//end if sidebands
 436                 // getting the number of triggers in the MCtag D* case
 437
 438
 439         if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);
 440
 441
 442                 if(!isInPeak && !isInSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME
 443
 444                 isDStarCand = kTRUE;
 445
 446             fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
 447
 448                 Short_t daughtercharge = ((AliAODTrack*)theD0particle->GetDaughter(0))->Charge();
 449                 fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
 450
 451
 452                 Int_t trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
 453                 Int_t trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
 454                 Int_t trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
 455
 456                 Bool_t execPool = fCorrelator->ProcessEventPool();
 457                 if(fmixing && !execPool) {
 458                         AliInfo("Mixed event analysis: pool is not ready");
 459                         continue;
 460                 }
 461
 462                 Int_t NofEventsinPool = 1;
 463                 if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
 464
 465                 for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one
 466
 467                         Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);
 468
 469                         if(!analyzetracks) {
 470                                 AliInfo("AliHFCorrelator::Cannot process the track array");
 471                                 continue;
 472                         }
 473
 474             //initialization of variables for correlations with leading particles
 475             Double_t DeltaPhiLeading = -999.;
 476                         Double_t DeltaEtaLeading = -999.;
 477                         //Double_t ptleading = -999.;
 478             Int_t labelleading = -999;
 479
 480                         Int_t NofTracks = fCorrelator->GetNofTracks();
 481
 482                         for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
 483
 484                                 Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
 485                                 if(!runcorrelation) continue;
 486
 487                                 Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
 488                                 Double_t DeltaEta = fCorrelator->GetDeltaEta();
 489
 490                                 AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
 491
 492                                 Double_t ptHad = hadron->Pt();
 493                                 Double_t phiHad = hadron->Phi();
 494                                 Double_t etaHad = hadron->Eta();
 495                                 Double_t label = hadron->GetLabel();
 496                                 Int_t trackid = hadron->GetID();
 497
 498                                 phiHad = fCorrelator->SetCorrectPhiRange(phiHad);
 499
 500                                 if(!fmixing){ // skip D* Daughetrs
 501                                         if(trackid == trackiddaugh0) continue;
 502                                         if(trackid == trackiddaugh1) continue;
 503                                         if(trackid == trackidsoftPi) continue;
 504                                 }
 505
 506                                 // from here on it is up to the user to decide what object to fill
 507
 508                                 if(fmontecarlo && isDStarMCtag){ // check correlations of MC tagged DStars in MonteCarlo
 509
 510                                         Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)
 511                                         FillMCTagCorrelations(ptDStar,DeltaPhi,DeltaEta,ptHad,PartSource);
 512
 513
 514                                         ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,0);
 515                                         if(PartSource[0]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,1);
 516                                         if(PartSource[1]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,2);
 517                                         if(PartSource[2]&&PartSource[0]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,3);
 518                                         if(PartSource[2]&&PartSource[1]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,4);
 519                                         if(PartSource[3]) ((TH3F*)fOutputMC->FindObject("MCPhiEtaPart"))->Fill(phiHad,etaHad,5);
 520
 521
 522                                 }
 523
 524                                 if(isInPeak)  {
 525
 526                                         if(fselect==1) ((TH3D*)fOutput->FindObject("DPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 527                                         if(fselect==2) ((TH3D*)fOutput->FindObject("DPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 528                                         if(fselect==3) ((TH3D*)fOutput->FindObject("DPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 529                                     ((TH2F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad);
 530                                         //counterPeak++; // count tracks per peak per event
 531
 532                                 }
 533
 534                                 if(isInSideBand) {
 535
 536                                         if(fselect==1) ((TH3D*)fOutput->FindObject("bkgDPhiDStarHadron"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 537                                         if(fselect==2) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKaon"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 538                                         if(fselect==3) ((TH3D*)fOutput->FindObject("bkgDPhiDStarKZero"))->Fill(DeltaPhi,ptDStar,DeltaEta);
 539
 540
 541
 542                                         //counterSB++;
 543                                 }
 544
 545
 546                         } // end loop on track candidates
 547
 548             // fill the leading particle histograms
 549
 550             if(isInPeak) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
 551                         if(isInSideBand) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);
 552
 553                         if(fmontecarlo && isDStarMCtag){
 554                 Bool_t* LeadPartSource = fCuts2->IsMCpartFromHF(labelleading,fmcArray);
 555                 FillMCTagLeadingCorrelations(ptDStar,DeltaPhiLeading,DeltaEtaLeading,LeadPartSource);
 556
 557             }
 558
 559                 } // end loop on events in the pool
 560
 561         }// end loop on D* candidates
 562
 563
 564         Bool_t updated = fCorrelator->PoolUpdate();
 565
 566         if(updated) EventMixingChecks(aodEvent);
 567         if(!updated) AliInfo("Pool was not updated");
 568
 569
 570
 571
 572 } //end the exec
 573
 574 //________________________________________ terminate ___________________________
 575 void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)
 576 {
 577         // The Terminate() function is the last function to be called during
 578         // a query. It always runs on the client, it can be used to present
 579         // the results graphically or save the results to file.
 580
 581         AliAnalysisTaskSE::Terminate();
 582
 583         fOutput = dynamic_cast<TList*> (GetOutputData(1));
 584         if (!fOutput) {
 585                 printf("ERROR: fOutput not available\n");
 586                 return;
 587         }
 588
 589         return;
 590 }
 591
 592
 593 //_____________________________________________________
 594 void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
 595
 596         Double_t Pi = TMath::Pi();
 597         Int_t nbinscorr = 32;
 598         Double_t lowcorrbin = -0.5*Pi - Pi/32; // shift the bin by half the width so that at 0 is it the bin center
 599         Double_t upcorrbin = 1.5*Pi - Pi/32;
 600
 601         // ========================= histograms for both Data and MonteCarlo
 602
 603
 604         TH1D * NofEvents = new TH1D("NofEvents","NofEvents",11,0,11);
 605         fOutput->Add(NofEvents);
 606
 607
 608
 609
 610
 611         TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);
 612         fOutput->Add(D0InvMass);
 613
 614         TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);
 615         fOutput->Add(D0InvMassinSB);
 616
 617         TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
 618         fOutput->Add(DeltaInvMass);
 619
 620         TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);
 621         fOutput->Add(bkgDeltaInvMass);
 622
 623         TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);
 624         fOutput->Add(RecoPtDStar);
 625
 626         TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);
 627         fOutput->Add(RecoPtBkg);
 628
 629         TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);
 630         fOutput->Add(MCtagPtDStar);
 631
 632         TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);
 633         if(fselect==3) fOutput->Add(KZeroSpectra);
 634
 635         TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);
 636         if(fselect==3) fOutput->Add(KZeroSpectraifHF);
 637
 638         TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","NofTracksInPeak",500,0.5,500.5);
 639         fOutput->Add(NofTracksInPeak);
 640
 641         TH1D * NofTracksInSB = new TH1D("NofTracksInSB","NofTracksInSB",500,0.5,500.5);
 642         fOutput->Add(NofTracksInSB);
 643
 644         TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);
 645         if(fmixing) fOutput->Add(EventMixingCheck);
 646
 647
 648
 649
 650
 651         TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
 652         fOutput->Add(PhiEtaTrigger);
 653
 654         TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
 655         fOutput->Add(PhiEtaSideBand);
 656
 657         TH2F * PhiEtaPart = new TH2F("PhiEtaPart","#phi distribution of the associated particle",36,-0.5*Pi,1.5*Pi,18,-0.9,0.9);
 658         fOutput->Add(PhiEtaPart);
 659
 660
 661         //correlations histograms
 662         TString histoname1 = "DPhiDStar";
 663         if(fselect==1) histoname1 += "Hadron";
 664         if(fselect==2) histoname1 += "Kaon";
 665         if(fselect==3) histoname1 += "KZero";
 666
 667
 668         TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 669
 670         TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 671
 672         //side band background histograms
 673         TString histoname2 = "bkg";
 674         histoname2 += histoname1;
 675         TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 676         TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 677
 678
 679         fOutput->Add(DPhiDStar);
 680
 681         if(fselect==3){fOutput->Add(DPhiDStarKZero1);}
 682
 683         fOutput->Add(bkgDPhiDStar);
 684
 685         if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}
 686
 687
 688         // leading particle
 689         TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 690         TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 691
 692         fOutput->Add(leadingcand);
 693         fOutput->Add(leadingsidebands);
 694
 695         // ========================= histos for analysis on MC only
 696
 697         TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
 698         if(fmontecarlo) fOutputMC->Add(EventTypeMC);
 699
 700         TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);
 701         MCSources->GetXaxis()->SetBinLabel(1,"All ");
 702         MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");
 703         MCSources->GetXaxis()->SetBinLabel(3," from c->D");
 704         MCSources->GetXaxis()->SetBinLabel(4," from b->D");
 705         MCSources->GetXaxis()->SetBinLabel(5," from b->B");
 706         MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
 707         MCSources->GetXaxis()->SetBinLabel(7," from c");
 708         MCSources->GetXaxis()->SetBinLabel(8," from b");
 709
 710         if(fmontecarlo) fOutputMC->Add(MCSources);
 711
 712     // leading particle from mc source
 713     TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);
 714         LeadingMCSources->GetXaxis()->SetBinLabel(1,"All ");
 715         LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");
 716         LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");
 717         LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");
 718         LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");
 719         LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");
 720         LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");
 721         LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");
 722
 723         if(fmontecarlo) fOutputMC->Add(LeadingMCSources);
 724
 725     // all hadrons
 726         TString histoname3 = "MCTag";
 727         histoname3 += histoname1;
 728         TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 729
 730         TString histoname44 = "CharmDOrigin";
 731         histoname44 += histoname1;
 732         histoname44 += "MC";
 733
 734         TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 735
 736
 737         TString histoname54 = "BeautyDOrigin";
 738         histoname54 += histoname1;
 739         histoname54 += "MC";
 740         TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 741
 742         TString histoname55 = "BeautyBOrigin";
 743         histoname55 += histoname1;
 744         histoname55 += "MC";
 745         TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 746
 747         TString histoname4 = "CharmQuarkOrigin";
 748         histoname4 += histoname1;
 749         histoname4 += "MC";
 750         TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 751
 752         TString histoname5 = "BeautyQuarkOrigin";
 753         histoname5 += histoname1;
 754         histoname5 += "MC";
 755         TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 756
 757     if(fmontecarlo){
 758
 759         fOutputMC->Add(MCTagDPhiDStar);
 760         fOutputMC->Add(CharmDOriginDPhiDStar);
 761         fOutputMC->Add(BeautyDOriginDPhiDStar);
 762         fOutputMC->Add(BeautyBOriginDPhiDStar);
 763         fOutputMC->Add(CharmQuarkOriginDPhiDStar);
 764         fOutputMC->Add(BeautyQuarkOriginDPhiDStar);
 765
 766         }
 767
 768     // ========================= histos for analysis on MC
 769     // all leading hadron
 770         TString Leadinghistoname3 = "LeadingMCTag";
 771         Leadinghistoname3 += histoname1;
 772         TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 773
 774         TString Leadinghistoname44 = "LeadingCharmDOrigin";
 775         Leadinghistoname44 += histoname1;
 776         Leadinghistoname44 += "MC";
 777
 778         TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 779
 780
 781         TString Leadinghistoname54 = "LeadingBeautyDOrigin";
 782         Leadinghistoname54 += histoname1;
 783         Leadinghistoname54 += "MC";
 784         TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 785
 786         TString Leadinghistoname55 = "LeadingBeautyBOrigin";
 787         Leadinghistoname55 += histoname1;
 788         Leadinghistoname55 += "MC";
 789         TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 790
 791         TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";
 792         Leadinghistoname4 += histoname1;
 793         Leadinghistoname4 += "MC";
 794         TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 795
 796         TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";
 797         Leadinghistoname5 += histoname1;
 798         Leadinghistoname5 += "MC";
 799         TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-1.95,1.95);
 800
 801
 802
 803
 804         if(fmontecarlo){
 805
 806                 fOutputMC->Add(LeadingMCTagDPhiDStar);
 807                 fOutputMC->Add(LeadingCharmDOriginDPhiDStar);
 808                 fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);
 809                 fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);
 810                 fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);
 811                 fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);
 812
 813         }
 814
 815         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);
 816         MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");
 817         MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");
 818         MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");
 819         MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");
 820         MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");
 821         MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");
 822         if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);
 823
 824         // ============================= EVENT MIXING CHECKS ======================================
 825
 826         Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();
 827         Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();
 828         Int_t NofCentBins = fCuts2->GetNCentPoolBins();
 829         Double_t * CentBins = fCuts2->GetCentPoolBins();
 830         Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();
 831         Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();
 832
 833         Int_t k =0;
 834
 835         if(fSystem) k = 100; // PbPb centrality
 836         if(!fSystem) k = NofCentBins; // pp multiplicity
 837
 838
 839         Double_t minvalue = CentBins[0];
 840         Double_t maxvalue = CentBins[NofCentBins+1];
 841         Double_t Zminvalue = ZVrtxBins[0];
 842         Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];
 843
 844
 845
 846         Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};
 847         Double_t * events = Nevents;
 848
 849         TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,events);
 850         EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
 851         EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
 852         EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");
 853         if(fmixing) fOutput->Add(EventsPerPoolBin);
 854
 855         Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;
 856         Int_t Diff = MaxNofTracks-MinNofTracks;
 857
 858         Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};
 859         Double_t  * trackN = Ntracks;
 860
 861         TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,trackN);
 862         NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
 863         NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
 864         NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");
 865
 866         if(fmixing) fOutput->Add(NofTracksPerPoolBin);
 867
 868         TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);
 869         NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");
 870         NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");
 871         if(fmixing) fOutput->Add(NofPoolBinCalls);
 872
 873
 874
 875         TH2D * EventProps = new TH2D("EventProps","Number of tracks in bin pool",k,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);
 876         EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");
 877         EventProps->GetYaxis()->SetTitle("Z vertex [cm]");
 878         if(fmixing) fOutput->Add(EventProps);
 879
 880 }
 881
 882
 883
 884 //____________________________  Function for MC correlations ___________________________________________________
 885 void AliAnalysisTaskDStarCorrelations::FillMCTagCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Double_t ptTrack, Bool_t *mcSource){
 886
 887
 888
 889
 890
 891                 if(fselect==1) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
 892                 if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarKaon"))->Fill(DelPhi,ptTrig,DelEta);
 893                 if(fselect==3) ((TH3D*)fOutputMC->FindObject("MCTagDPhiDStarKZero"))->Fill(DelPhi,ptTrig,DelEta);
 894
 895
 896
 897                 ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(0);
 898
 899                 if(fDebugLevel){
 900                         std::cout << "MC source " << mcSource[0] << " "  << mcSource[1] << " " << mcSource[2] << " " << mcSource[3] << std::endl;
 901
 902                         if(mcSource[0]) std::cout << "mcSource 0 " << std::endl;
 903                         if(mcSource[1]) std::cout << "mcSource 1 " << std::endl;
 904                         if(mcSource[2]) std::cout << "mcSource 2 " << std::endl;
 905                         if(mcSource[3]) std::cout << "mcSource 3 " << std::endl;
 906
 907                 }
 908                 if(mcSource[0]){ // is from charm quark
 909                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(5); // all HF quarks
 910                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(6); //  charm quarks
 911                         if(fselect==1) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 912                         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
 913                         if(fselect==3) ((TH3D*)fOutputMC->FindObject("CharmQuarkOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
 914                 }
 915
 916                 if(mcSource[1]){ // is from b quark
 917                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(5); // all HF quarks
 918                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(7); // beauty quarks
 919                         if(fselect==1) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 920                         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
 921                         if(fselect==3) ((TH3D*)fOutputMC->FindObject("BeautyQuarkOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
 922
 923                 }
 924
 925                 if(mcSource[2]&&mcSource[0]){ // is from D meson and charm quark
 926                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(1); // all HF mesons
 927                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(2); //  charm + D
 928                         if(fselect==1) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 929                         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
 930                         if(fselect==3) ((TH3D*)fOutputMC->FindObject("CharmDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
 931                 }
 932
 933                 if(mcSource[2]&&mcSource[1]){ // is from D meson and b quark
 934                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(1); // all HF mesons
 935                         ((TH1F*)fOutputMC->FindObject("MCSources"))->Fill(3); //  beauty + D
 936                         if(fselect==1) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 937                         if(fselect==2 && ptTrack <1.5) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarKaonMC"))->Fill(DelPhi,ptTrig,DelEta);
 938                         if(fselect==3) ((TH3D*)fOutputMC->FindObject("BeautyDOriginDPhiDStarKZeroMC"))->Fill(DelPhi,ptTrig,DelEta);
 939                 }
 940
 941         return;
 942 }
 943
 944 //____________________________  Function for MC leading part correlations ___________________________________________________
 945 void AliAnalysisTaskDStarCorrelations::FillMCTagLeadingCorrelations(Double_t ptTrig, Double_t DelPhi,  Double_t DelEta, Bool_t *mcSource){
 946     // correlations with leading hadron on MC
 947
 948         if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingMCTagDPhiDStarHadron"))->Fill(DelPhi,ptTrig,DelEta);
 949
 950
 951
 952         ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(0);
 953
 954         if(fDebugLevel){ std::cout << "MC source " << mcSource[0] << " "  << mcSource[1] << " " << mcSource[2] << " " << mcSource[3] << std::endl;
 955
 956                 if(mcSource[0]) std::cout << "mcSource 0 " << std::endl;
 957                 if(mcSource[1]) std::cout << "mcSource 1 " << std::endl;
 958                 if(mcSource[2]) std::cout << "mcSource 2 " << std::endl;
 959                 if(mcSource[3]) std::cout << "mcSource 3 " << std::endl;
 960     }
 961
 962         if(mcSource[0]){ // is from charm quark
 963                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(5); // all HF quarks
 964                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(6); //  charm quarks
 965                 if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingCharmQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 966
 967         }
 968
 969         if(mcSource[1]){ // is from b quaLeadingrk
 970                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(5); // all HF quarks
 971                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(7); // beauty quarks
 972                 if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingBeautyQuarkOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 973
 974
 975         }
 976
 977         if(mcSource[2]&&mcSource[0]){ // is from D meson and charm quark
 978                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(1); // all HF mesons
 979                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(2); //  charm + D
 980                 if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingCharmDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 981
 982         }
 983
 984         if(mcSource[2]&&mcSource[1]){ // is from D meson and b quark
 985                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(1); // all HF mesons
 986                 ((TH1F*)fOutputMC->FindObject("LeadingMCSources"))->Fill(3); //  beauty + D
 987                 if(fselect==1) ((TH3D*)fOutputMC->FindObject("LeadingBeautyDOriginDPhiDStarHadronMC"))->Fill(DelPhi,ptTrig,DelEta);
 988
 989         }
 990
 991
 992         return;
 993 }
 994
 995
 996 //____________________________  Run checks on event mixing ___________________________________________________
 997 void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
 998
 999         AliCentrality *centralityObj = 0;
1000         Int_t multiplicity = -1;
1001         Double_t MultipOrCent = -1;
1002
1003         // get the pool for event mixing
1004         if(!fSystem){ // pp
1005                 multiplicity = AOD->GetNTracks();
1006                 MultipOrCent = multiplicity; // convert from Int_t to Double_t
1007         }
1008         if(fSystem){ // PbPb
1009
1010                 centralityObj = AOD->GetHeader()->GetCentralityP();
1011                 MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
1012                 AliInfo(Form("Centrality is %f", MultipOrCent));
1013         }
1014
1015         AliAODVertex *vtx = AOD->GetPrimaryVertex();
1016         Double_t zvertex = vtx->GetZ(); // zvertex
1017
1018
1019
1020
1021         AliEventPool * pool = fCorrelator->GetPool();
1022
1023
1024
1025
1026         ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
1027         ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
1028
1029         ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of events in the pool
1030         ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of calls of pool
1031 }
1032
1033
1034
1035
1036