PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2012, 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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for D0 candidates (2Prongs)
  21 // and hadrons correlations
  22 //
  23 // Authors: Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
  24 // Fabio Colamaria, fabio.colamaria@ba.infn.it (correlations)
  25 /////////////////////////////////////////////////////////////
  26
  27 #include <Riostream.h>
  28 #include <TClonesArray.h>
  29 #include <TCanvas.h>
  30 #include <TNtuple.h>
  31 #include <TTree.h>
  32 #include <TList.h>
  33 #include <TH1F.h>
  34 #include <TH2F.h>
  35 #include <THnSparse.h>
  36 #include <TDatabasePDG.h>
  37
  38 #include <AliAnalysisDataSlot.h>
  39 #include <AliAnalysisDataContainer.h>
  40 #include "AliAnalysisManager.h"
  41 #include "AliESDtrack.h"
  42 #include "AliVertexerTracks.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODRecoDecayHF2Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliAnalysisVertexingHF.h"
  52 #include "AliAnalysisTaskSE.h"
  53 #include "AliAnalysisTaskSED0Correlations.h"
  54 #include "AliNormalizationCounter.h"
  55
  56 using std::cout;
  57 using std::endl;
  58
  59 ClassImp(AliAnalysisTaskSED0Correlations)
  60
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations():
  64 AliAnalysisTaskSE(),
  65   fNPtBinsCorr(0),
  66   fBinLimsCorr(),
  67   fPtThreshLow(),
  68   fPtThreshUp(),
  69   fEvents(0),
  70   fAlreadyFilled(kFALSE),
  71   fOutputMass(0),
  72   fOutputCorr(0),
  73   fOutputStudy(0),
  74   fNentries(0),
  75   fCutsD0(0),
  76   fCutsTracks(0),
  77   fCorrelatorTr(0),
  78   fCorrelatorKc(0),
  79   fCorrelatorK0(0),
  80   fReadMC(0),
  81   fRecoTr(kTRUE),
  82   fRecoD0(kTRUE),
  83   fSelEvType(kFALSE),
  84   fMixing(kFALSE),
  85   fCounter(0),
  86   fNPtBins(1),
  87   fFillOnlyD0D0bar(0),
  88   fIsSelectedCandidate(0),
  89   fSys(0),
  90   fEtaForCorrel(0),
  91   fIsRejectSDDClusters(0),
  92   fFillGlobal(kTRUE)
  93 {
  94   // Default constructor
  95
  96 }
  97
  98 //________________________________________________________________________
  99 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const char *name,AliRDHFCutsD0toKpi* cutsD0, AliHFAssociatedTrackCuts* cutsTrk):
 100   AliAnalysisTaskSE(name),
 101   fNPtBinsCorr(0),
 102   fBinLimsCorr(),
 103   fPtThreshLow(),
 104   fPtThreshUp(),
 105   fEvents(0),
 106   fAlreadyFilled(kFALSE),
 107   fOutputMass(0),
 108   fOutputCorr(0),
 109   fOutputStudy(0),
 110   fNentries(0),
 111   fCutsD0(0),
 112   fCutsTracks(cutsTrk),
 113   fCorrelatorTr(0),
 114   fCorrelatorKc(0),
 115   fCorrelatorK0(0),
 116   fReadMC(0),
 117   fRecoTr(kTRUE),
 118   fRecoD0(kTRUE),
 119   fSelEvType(kFALSE),
 120   fMixing(kFALSE),
 121   fCounter(0),
 122   fNPtBins(1),
 123   fFillOnlyD0D0bar(0),
 124   fIsSelectedCandidate(0),
 125   fSys(0),
 126   fEtaForCorrel(0),
 127   fIsRejectSDDClusters(0),
 128   fFillGlobal(kTRUE)
 129 {
 130   // Default constructor
 131
 132   fNPtBins=cutsD0->GetNPtBins();
 133
 134   fCutsD0=cutsD0;
 135
 136   // Output slot #1 writes into a TList container (mass with cuts)
 137   DefineOutput(1,TList::Class());  //My private output
 138   // Output slot #2 writes into a TH1F container (number of events)
 139   DefineOutput(2,TH1F::Class());  //My private output
 140   // Output slot #3 writes into a AliRDHFD0toKpi container (cuts)
 141   DefineOutput(3,AliRDHFCutsD0toKpi::Class());  //My private output
 142   // Output slot #4 writes Normalization Counter
 143   DefineOutput(4,AliNormalizationCounter::Class());
 144   // Output slot #5 writes into a TList container (correl output)
 145   DefineOutput(5,TList::Class());  //My private output
 146   // Output slot #6 writes into a TList container (correl advanced)
 147   DefineOutput(6,TList::Class());  //My private output
 148   // Output slot #7 writes into a AliHFAssociatedTrackCuts container (cuts)
 149   DefineOutput(7,AliHFAssociatedTrackCuts::Class());  //My private output
 150 }
 151
 152 //________________________________________________________________________
 153 AliAnalysisTaskSED0Correlations::AliAnalysisTaskSED0Correlations(const AliAnalysisTaskSED0Correlations &source):
 154   AliAnalysisTaskSE(source),
 155   fNPtBinsCorr(source.fNPtBinsCorr),
 156   fBinLimsCorr(source.fBinLimsCorr),
 157   fPtThreshLow(source.fPtThreshLow),
 158   fPtThreshUp(source.fPtThreshUp),
 159   fEvents(source.fEvents),
 160   fAlreadyFilled(source.fAlreadyFilled),
 161   fOutputMass(source.fOutputMass),
 162   fOutputCorr(source.fOutputCorr),
 163   fOutputStudy(source.fOutputStudy),
 164   fNentries(source.fNentries),
 165   fCutsD0(source.fCutsD0),
 166   fCutsTracks(source.fCutsTracks),
 167   fCorrelatorTr(source.fCorrelatorTr),
 168   fCorrelatorKc(source.fCorrelatorKc),
 169   fCorrelatorK0(source.fCorrelatorK0),
 170   fReadMC(source.fReadMC),
 171   fRecoTr(source.fRecoTr),
 172   fRecoD0(source.fRecoD0),
 173   fSelEvType(source.fSelEvType),
 174   fMixing(source.fMixing),
 175   fCounter(source.fCounter),
 176   fNPtBins(source.fNPtBins),
 177   fFillOnlyD0D0bar(source.fFillOnlyD0D0bar),
 178   fIsSelectedCandidate(source.fIsSelectedCandidate),
 179   fSys(source.fSys),
 180   fEtaForCorrel(source.fEtaForCorrel),
 181   fIsRejectSDDClusters(source.fIsRejectSDDClusters),
 182   fFillGlobal(source.fFillGlobal)
 183 {
 184   // Copy constructor
 185 }
 186
 187 //________________________________________________________________________
 188 AliAnalysisTaskSED0Correlations::~AliAnalysisTaskSED0Correlations()
 189 {
 190   if (fOutputMass) {
 191     delete fOutputMass;
 192     fOutputMass = 0;
 193   }
 194   if (fOutputCorr) {
 195     delete fOutputCorr;
 196     fOutputCorr = 0;
 197   }
 198   if (fOutputStudy) {
 199     delete fOutputStudy;
 200     fOutputStudy = 0;
 201   }
 202   if (fCutsD0) {
 203     delete fCutsD0;
 204     fCutsD0 = 0;
 205   }
 206   if (fNentries){
 207     delete fNentries;
 208     fNentries = 0;
 209   }
 210   if (fCorrelatorTr) {
 211     delete fCorrelatorTr;
 212     fCorrelatorTr = 0;
 213   }
 214   if (fCorrelatorKc) {
 215     delete fCorrelatorKc;
 216     fCorrelatorKc = 0;
 217   }
 218   if (fCorrelatorK0) {
 219     delete fCorrelatorK0;
 220     fCorrelatorK0 = 0;
 221   }
 222   if (fCounter){
 223     delete fCounter;
 224     fCounter=0;
 225   }
 226 }
 227
 228 //______________________________________________________________________________
 229 AliAnalysisTaskSED0Correlations& AliAnalysisTaskSED0Correlations::operator=(const AliAnalysisTaskSED0Correlations& orig)
 230 {
 231 // Assignment
 232   if (&orig == this) return *this; //if address is the same (same object), returns itself
 233
 234   AliAnalysisTaskSE::operator=(orig); //Uses the AliAnalysisTaskSE operator to assign the inherited part of the class
 235   fNPtBinsCorr = orig.fNPtBinsCorr;
 236   fBinLimsCorr = orig.fBinLimsCorr;
 237   fPtThreshLow = orig.fPtThreshLow;
 238   fPtThreshUp = orig.fPtThreshUp;
 239   fEvents = orig.fEvents;
 240   fAlreadyFilled = orig.fAlreadyFilled;
 241   fOutputMass = orig.fOutputMass;
 242   fOutputCorr = orig.fOutputCorr;
 243   fOutputStudy = orig.fOutputStudy;
 244   fNentries = orig.fNentries;
 245   fCutsD0 = orig.fCutsD0;
 246   fCutsTracks = orig.fCutsTracks;
 247   fCorrelatorTr = orig.fCorrelatorTr;
 248   fCorrelatorKc = orig.fCorrelatorKc;
 249   fCorrelatorK0 = orig.fCorrelatorK0;
 250   fReadMC = orig.fReadMC;
 251   fRecoTr = orig.fRecoTr;
 252   fRecoD0 = orig.fRecoD0;
 253   fSelEvType = orig.fSelEvType;
 254   fMixing = orig.fMixing;
 255   fCounter = orig.fCounter;
 256   fNPtBins = orig.fNPtBins;
 257   fFillOnlyD0D0bar = orig.fFillOnlyD0D0bar;
 258   fIsSelectedCandidate = orig.fIsSelectedCandidate;
 259   fSys = orig.fSys;
 260   fEtaForCorrel = orig.fEtaForCorrel;
 261   fIsRejectSDDClusters = orig.fIsRejectSDDClusters;
 262   fFillGlobal = orig.fFillGlobal;
 263
 264   return *this; //returns pointer of the class
 265 }
 266
 267 //________________________________________________________________________
 268 void AliAnalysisTaskSED0Correlations::Init()
 269 {
 270   // Initialization
 271
 272   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::Init() \n");
 273
 274   //Copy of cuts objects
 275   AliRDHFCutsD0toKpi* copyfCutsD0 = new AliRDHFCutsD0toKpi(*fCutsD0);
 276   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 277   copyfCutsD0->SetName(nameoutput);
 278
 279   //needer to clear completely the objects inside with Clear() method
 280   // Post the data
 281   PostData(3,copyfCutsD0);
 282   PostData(7,fCutsTracks);
 283
 284   return;
 285 }
 286
 287 //________________________________________________________________________
 288 void AliAnalysisTaskSED0Correlations::UserCreateOutputObjects()
 289 {
 290
 291   // Create the output container
 292   //
 293   if(fDebug > 1) printf("AnalysisTaskSED0Correlations::UserCreateOutputObjects() \n");
 294
 295   //HFCorrelator creation and definition
 296   fCorrelatorTr = new AliHFCorrelator("CorrelatorTr",fCutsTracks,fSys);
 297   fCorrelatorKc = new AliHFCorrelator("CorrelatorKc",fCutsTracks,fSys);
 298   fCorrelatorK0 = new AliHFCorrelator("CorrelatorK0",fCutsTracks,fSys);
 299   fCorrelatorTr->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);// set the Delta Phi Interval you want (in this case -0.5Pi to 1.5 Pi)
 300   fCorrelatorKc->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
 301   fCorrelatorK0->SetDeltaPhiInterval(-TMath::Pi()/2,3*TMath::Pi()/2);
 302   fCorrelatorTr->SetEventMixing(fMixing);// sets the analysis on a single event (kFALSE) or mixed events (kTRUE)
 303   fCorrelatorKc->SetEventMixing(fMixing);
 304   fCorrelatorK0->SetEventMixing(fMixing);
 305   fCorrelatorTr->SetAssociatedParticleType(1);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 306   fCorrelatorKc->SetAssociatedParticleType(2);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 307   fCorrelatorK0->SetAssociatedParticleType(3);// set 1 for correlations with hadrons, 2 with kaons, 3 with KZeros
 308   fCorrelatorTr->SetApplyDisplacementCut(2); //0: don't calculate d0; 1: return d0; 2: return d0/d0err
 309   fCorrelatorKc->SetApplyDisplacementCut(2);
 310   fCorrelatorK0->SetApplyDisplacementCut(0);
 311   fCorrelatorTr->SetUseMC(fReadMC);// sets Montecarlo flag
 312   fCorrelatorKc->SetUseMC(fReadMC);
 313   fCorrelatorK0->SetUseMC(fReadMC);
 314   fCorrelatorTr->SetUseReco(fRecoTr);// sets (if MC analysis) wheter to analyze Reco or Kinem tracks
 315   fCorrelatorKc->SetUseReco(fRecoTr);
 316   fCorrelatorK0->SetUseReco(fRecoTr);
 317   fCorrelatorKc->SetPIDmode(2); //switch for K+/- PID option
 318   Bool_t pooldefTr = fCorrelatorTr->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 319   Bool_t pooldefKc = fCorrelatorKc->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 320   Bool_t pooldefK0 = fCorrelatorK0->DefineEventPool();// method that defines the properties ot the event mixing (zVtx and Multipl. bins)
 321   if(!pooldefTr) AliInfo("Warning:: Event pool not defined properly");
 322   if(!pooldefKc) AliInfo("Warning:: Event pool not defined properly");
 323   if(!pooldefK0) AliInfo("Warning:: Event pool not defined properly");
 324
 325   // Several histograms are more conveniently managed in a TList
 326   fOutputMass = new TList();
 327   fOutputMass->SetOwner();
 328   fOutputMass->SetName("listMass");
 329
 330   fOutputCorr = new TList();
 331   fOutputCorr->SetOwner();
 332   fOutputCorr->SetName("correlationslist");
 333
 334   fOutputStudy = new TList();
 335   fOutputStudy->SetOwner();
 336   fOutputStudy->SetName("MCstudyplots");
 337
 338   TString nameMass=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ";
 339
 340   for(Int_t i=0;i<fCutsD0->GetNPtBins();i++){
 341
 342     nameMass="histMass_";
 343     nameMass+=i;
 344     nameSgn="histSgn_";
 345     nameSgn+=i;
 346     nameBkg="histBkg_";
 347     nameBkg+=i;
 348     nameRfl="histRfl_";
 349     nameRfl+=i;
 350
 351     //histograms of invariant mass distributions
 352
 353     //MC signal
 354     if(fReadMC){
 355       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
 356       tmpSt->Sumw2();
 357
 358       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 359       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
 360       TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",120,1.5648,2.1648);
 361       tmpBt->Sumw2();
 362       tmpRt->Sumw2();
 363       fOutputMass->Add(tmpSt);
 364       fOutputMass->Add(tmpRt);
 365       fOutputMass->Add(tmpBt);
 366     }
 367
 368     //mass
 369     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",120,1.5648,2.1648);
 370     tmpMt->Sumw2();
 371     fOutputMass->Add(tmpMt);
 372   }
 373
 374   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 375
 376   fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 20,-0.5,19.5);
 377
 378   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 379   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 380   fReadMC ? fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected") : fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
 381   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 382   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 383   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 384   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 385   if(fReadMC && fSys==0){
 386     fNentries->GetXaxis()->SetBinLabel(12,"K");
 387     fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
 388   }
 389   fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
 390   fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
 391   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
 392   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
 393   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
 394   fNentries->GetXaxis()->SetBinLabel(19,"nEventsSelected");
 395   if(fReadMC) fNentries->GetXaxis()->SetBinLabel(20,"nEvsWithProdMech");
 396   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 397
 398   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(4)->GetContainer()->GetName()));
 399   fCounter->Init();
 400
 401   CreateCorrelationsObjs(); //creates histos for correlations analysis
 402
 403   // Post the data
 404   PostData(1,fOutputMass);
 405   PostData(2,fNentries);
 406   PostData(4,fCounter);
 407   PostData(5,fOutputCorr);
 408   PostData(6,fOutputStudy);
 409
 410   return;
 411 }
 412
 413 //________________________________________________________________________
 414 void AliAnalysisTaskSED0Correlations::UserExec(Option_t */*option*/)
 415 {
 416   // Execute analysis for current event:
 417   // heavy flavor candidates association to MC truth
 418   //cout<<"I'm in UserExec"<<endl;
 419
 420
 421   //cuts order
 422   //     printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 423   //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 424   //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 425   //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 426   //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 427   //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 428   //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 429   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 430   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 431
 432   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 433   fEvents++;
 434
 435   TString bname="D0toKpi";
 436
 437   TClonesArray *inputArray=0;
 438
 439   if(!aod && AODEvent() && IsStandardAOD()) {
 440     // In case there is an AOD handler writing a standard AOD, use the AOD
 441     // event in memory rather than the input (ESD) event.
 442     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 443     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 444     // have to taken from the AOD event hold by the AliAODExtension
 445     AliAODHandler* aodHandler = (AliAODHandler*)
 446       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 447
 448     if(aodHandler->GetExtensions()) {
 449       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 450       AliAODEvent* aodFromExt = ext->GetAOD();
 451       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 452     }
 453   } else if(aod) {
 454     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 455   }
 456
 457   if(!inputArray || !aod) {
 458     printf("AliAnalysisTaskSED0Correlations::UserExec: input branch not found!\n");
 459     return;
 460   }
 461
 462   // fix for temporary bug in ESDfilter
 463   // the AODs with null vertex pointer didn't pass the PhysSel
 464   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 465
 466   TClonesArray *mcArray = 0;
 467   AliAODMCHeader *mcHeader = 0;
 468
 469   if(fReadMC) {
 470     // load MC particles
 471     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 472     if(!mcArray) {
 473       printf("AliAnalysisTaskSED0Correlations::UserExec: MC particles branch not found!\n");
 474       return;
 475     }
 476
 477     // load MC header
 478     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 479     if(!mcHeader) {
 480       printf("AliAnalysisTaskSED0Correlations::UserExec: MC header branch not found!\n");
 481       return;
 482     }
 483   }
 484
 485   //histogram filled with 1 for every AOD
 486   fNentries->Fill(0);
 487   fCounter->StoreEvent(aod,fCutsD0,fReadMC);
 488
 489   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 490   TString trigclass=aod->GetFiredTriggerClasses();
 491   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
 492
 493   if(!fCutsD0->IsEventSelected(aod)) {
 494     if(fCutsD0->GetWhyRejection()==1) // rejected for pileup
 495       fNentries->Fill(13);
 496     if(fSys==1 && (fCutsD0->GetWhyRejection()==2 || fCutsD0->GetWhyRejection()==3)) fNentries->Fill(15);
 497     if(fCutsD0->GetWhyRejection()==7) fNentries->Fill(17);
 498     return;
 499   }
 500
 501   fNentries->Fill(18); //event selected after selection
 502
 503   //Setting PIDResponse for associated tracks
 504   fCorrelatorTr->SetPidAssociated();
 505   fCorrelatorKc->SetPidAssociated();
 506   fCorrelatorK0->SetPidAssociated();
 507
 508   //Selection on production type (MC)
 509   if(fReadMC && fSelEvType){
 510
 511     Bool_t isMCeventgood = kFALSE;
 512
 513     Int_t eventType = mcHeader->GetEventType();
 514     Int_t NMCevents = fCutsTracks->GetNofMCEventType();
 515
 516     for(Int_t k=0; k<NMCevents; k++){
 517       Int_t * MCEventType = fCutsTracks->GetMCEventType();
 518
 519       if(eventType == MCEventType[k]) isMCeventgood= kTRUE;
 520       ((TH1D*)fOutputStudy->FindObject("EventTypeMC"))->Fill(eventType);
 521     }
 522
 523     if(NMCevents && !isMCeventgood){
 524       std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;
 525       return;
 526     }
 527     fNentries->Fill(19); //event with particular production type
 528
 529   } //end of selection
 530
 531
 532   // Check the Nb of SDD clusters
 533   if (fIsRejectSDDClusters) {
 534     Bool_t skipEvent = kFALSE;
 535     Int_t ntracks = 0;
 536     if (aod) ntracks = aod->GetNTracks();
 537     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
 538       //    ... get the track
 539       AliAODTrack * track = aod->GetTrack(itrack);
 540       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 541         skipEvent=kTRUE;
 542         fNentries->Fill(16);
 543         break;
 544       }
 545     }
 546     if (skipEvent) return;
 547   }
 548
 549   //HFCorrelators initialization (for this event)
 550   fCorrelatorTr->SetAODEvent(aod); // set the AOD event from which you are processing
 551   fCorrelatorKc->SetAODEvent(aod);
 552   fCorrelatorK0->SetAODEvent(aod);
 553   Bool_t correlatorONTr = fCorrelatorTr->Initialize(); // initialize the pool for event mixing
 554   Bool_t correlatorONKc = fCorrelatorKc->Initialize();
 555   Bool_t correlatorONK0 = fCorrelatorK0->Initialize();
 556   if(!correlatorONTr) {AliInfo("AliHFCorrelator (tracks) didn't initialize the pool correctly or processed a bad event"); return;}
 557   if(!correlatorONKc) {AliInfo("AliHFCorrelator (charged K) didn't initialize the pool correctly or processed a bad event"); return;}
 558   if(!correlatorONK0) {AliInfo("AliHFCorrelator (K0) didn't initialize the pool correctly or processed a bad event"); return;}
 559
 560   if(fReadMC) {
 561     fCorrelatorTr->SetMCArray(mcArray); // set the TClonesArray *fmcArray for analysis on monte carlo
 562     fCorrelatorKc->SetMCArray(mcArray);
 563     fCorrelatorK0->SetMCArray(mcArray);
 564   }
 565
 566   // AOD primary vertex
 567   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 568
 569   //vtx1->Print();
 570   TString primTitle = vtx1->GetTitle();
 571   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 572     fNentries->Fill(3);
 573   }
 574
 575   //Reset flag for tracks distributions fill
 576   fAlreadyFilled=kFALSE;
 577
 578   //***** Loop over D0 candidates *****
 579   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 580   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 581
 582   if(fFillGlobal) { //loop on V0 and tracks for each event, to fill Pt distr. and InvMass distr.
 583
 584     TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
 585     Int_t pdgCodes[2] = {211,211};
 586     Int_t idArrayV0[v0array->GetEntriesFast()][2];
 587     for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
 588       for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
 589       AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
 590       if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
 591         if(fReadMC && fRecoTr && (v0->MatchToMC(310,mcArray,2,pdgCodes)<0)) continue; //310 = K0s, 311 = K0 generico!!
 592         ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
 593         ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
 594       }
 595     }
 596
 597     for(Int_t itrack=0; itrack<aod->GetNTracks(); itrack++) { // loop on tacks
 598       AliAODTrack * track = aod->GetTrack(itrack);
 599       //rejection of tracks
 600       if(track->GetID() < 0) continue; //discard negative ID tracks
 601       if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
 602       if(!fCutsTracks->IsHadronSelected(track) || !fCutsTracks->CheckHadronKinematic(track->Pt(),0.1)) continue; //0.1 = dummy (d0 value, no cut on it for me)
 603       //pT distribution (in all events), charg and hadr cases
 604       ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt());
 605       if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
 606     }
 607
 608   } //end of loops for global plot fill
 609
 610   Int_t nSelectedloose=0,nSelectedtight=0;
 611
 612   //RecoD0 case ************************************************
 613   if(fRecoD0) {
 614
 615     for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 616       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 617
 618       if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 619         fNentries->Fill(2);
 620         continue; //skip the D0 from Dstar
 621       }
 622
 623       if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
 624         nSelectedloose++;
 625         nSelectedtight++;
 626         if(fSys==0){
 627           if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
 628         }
 629         Int_t ptbin=fCutsD0->PtBin(d->Pt());
 630         if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 631
 632         fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
 633         if(!fIsSelectedCandidate) continue;
 634
 635         //D0 infos
 636         Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(d->Phi());
 637                  phiD0 = fCorrelatorKc->SetCorrectPhiRange(d->Phi());  //bad usage, but returns a Double_t...
 638                  phiD0 = fCorrelatorK0->SetCorrectPhiRange(d->Phi());
 639         fCorrelatorTr->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta()); // sets the parameters of the trigger particles that are needed
 640         fCorrelatorKc->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 641         fCorrelatorK0->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 642         fCorrelatorTr->SetD0Properties(d,fIsSelectedCandidate); //sets special properties for D0
 643         fCorrelatorKc->SetD0Properties(d,fIsSelectedCandidate);
 644         fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
 645
 646         if(!fReadMC) {
 647           if (TMath::Abs(d->Eta())<fEtaForCorrel) CalculateCorrelations(d); //correlations on real data
 648         } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
 649           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
 650             Int_t pdgDgD0toKpi[2]={321,211};
 651             Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
 652             if (labD0>-1) CalculateCorrelations(d,labD0,mcArray);
 653           }
 654         }
 655
 656         FillMassHists(d,mcArray,fCutsD0,fOutputMass);
 657       }
 658     }
 659   }
 660   //End RecoD0 case ************************************************
 661
 662   //MCKineD0 case ************************************************
 663   if(fReadMC && !fRecoD0) {
 664
 665     for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { //Loop over all the tracks of MCArray
 666       AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 667       if (!mcPart) {
 668         AliWarning("Particle not found in tree, skipping");
 669         continue;
 670       }
 671
 672       if(TMath::Abs(mcPart->GetPdgCode()) == 421){  // THIS IS A D0
 673         if (fCutsD0->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y()) ) {
 674           nSelectedloose++;
 675           nSelectedtight++;
 676           if(fSys==0) fNentries->Fill(6);
 677           Int_t ptbin=fCutsD0->PtBin(mcPart->Pt());
 678           if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 679
 680           //D0 infos
 681           Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(mcPart->Phi());
 682                    phiD0 = fCorrelatorKc->SetCorrectPhiRange(mcPart->Phi());  //bad usage, but returns a Double_t...
 683                    phiD0 = fCorrelatorK0->SetCorrectPhiRange(mcPart->Phi());
 684           fCorrelatorTr->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta()); // sets the parameters of the trigger particles that are needed
 685           fCorrelatorKc->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 686           fCorrelatorK0->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 687           //fCorrelatorTr->SetD0Properties(mcPart,fIsSelectedCandidate); //needed for D* soft pions rejection, useless in MCKine
 688           //fCorrelatorKc->SetD0Properties(mcPart,fIsSelectedCandidate);
 689           //fCorrelatorK0->SetD0Properties(mcPart,fIsSelectedCandidate);
 690
 691           if (TMath::Abs(mcPart->Eta())<fEtaForCorrel) {
 692
 693             //Removal of D0 from D* feeddown! This solves also the problem of soft pions, now excluded
 694          /*   Int_t mother = mcPart->GetMother();
 695             AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
 696             if(!mcMoth) continue;
 697             if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;
 698         */
 699             if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
 700             else fIsSelectedCandidate = 2;
 701
 702             TString fillthis="histSgn_"; fillthis+=ptbin;
 703             ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
 704
 705             CalculateCorrelationsMCKine(mcPart,mcArray);
 706           }
 707         }
 708       }
 709     }
 710   }
 711   //End MCKineD0 case ************************************************
 712
 713   if(fMixing /* && fAlreadyFilled*/) { // update the pool for Event Mixing, if: enabled,  event is ok, at least a SelD0 found! (fAlreadyFilled's role!)
 714     Bool_t updatedTr = fCorrelatorTr->PoolUpdate();
 715     Bool_t updatedKc = fCorrelatorKc->PoolUpdate();
 716     Bool_t updatedK0 = fCorrelatorK0->PoolUpdate();
 717     if(!updatedTr || !updatedKc || !updatedK0) AliInfo("Pool was not updated");
 718   }
 719
 720   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 721   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 722
 723   // Post the data
 724   PostData(1,fOutputMass);
 725   PostData(2,fNentries);
 726   PostData(4,fCounter);
 727   PostData(5,fOutputCorr);
 728   PostData(6,fOutputStudy);
 729
 730   return;
 731 }
 732
 733 //____________________________________________________________________________
 734 void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
 735   //
 736   // function used in UserExec to fill mass histograms:
 737   //
 738   if (!fIsSelectedCandidate) return;
 739
 740   if(fDebug>2)  cout<<"Candidate selected"<<endl;
 741
 742   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 743   Int_t ptbin = cuts->PtBin(part->Pt());
 744
 745   TString fillthis="";
 746   Int_t pdgDgD0toKpi[2]={321,211};
 747   Int_t labD0=-1;
 748   if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 749
 750   //count candidates selected by cuts
 751   fNentries->Fill(1);
 752   //count true D0 selected by cuts
 753   if (fReadMC && labD0>=0) fNentries->Fill(2);
 754
 755   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
 756
 757     if(fReadMC){
 758       if(labD0>=0) {
 759         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 760         Int_t pdgD0 = partD0->GetPdgCode();
 761         if (pdgD0==421){ //D0
 762           fillthis="histSgn_";
 763           fillthis+=ptbin;
 764           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 765         } else{ //it was a D0bar
 766           fillthis="histRfl_";
 767           fillthis+=ptbin;
 768           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 769         }
 770       } else {//background
 771         fillthis="histBkg_";
 772         fillthis+=ptbin;
 773         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 774       }
 775     }else{
 776       fillthis="histMass_";
 777       fillthis+=ptbin;
 778       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 779     }
 780
 781   }
 782   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
 783
 784     if(fReadMC){
 785       if(labD0>=0) {
 786         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 787         Int_t pdgD0 = partD0->GetPdgCode();
 788
 789         if (pdgD0==-421){ //D0bar
 790           fillthis="histSgn_";
 791           fillthis+=ptbin;
 792           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 793         } else{
 794           fillthis="histRfl_";
 795           fillthis+=ptbin;
 796           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 797         }
 798       } else {//background or LS
 799         fillthis="histBkg_";
 800         fillthis+=ptbin;
 801         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 802       }
 803     }else{
 804       fillthis="histMass_";
 805       fillthis+=ptbin;
 806       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
 807     }
 808   }
 809
 810   return;
 811 }
 812
 813 //________________________________________________________________________
 814 void AliAnalysisTaskSED0Correlations::Terminate(Option_t */*option*/)
 815 {
 816   // Terminate analysis
 817   //
 818   if(fDebug > 1) printf("AnalysisTaskSED0Correlations: Terminate() \n");
 819
 820   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
 821   if (!fOutputMass) {
 822     printf("ERROR: fOutputMass not available\n");
 823     return;
 824   }
 825
 826   fNentries = dynamic_cast<TH1F*>(GetOutputData(2));
 827
 828   if(!fNentries){
 829     printf("ERROR: fNEntries not available\n");
 830     return;
 831   }
 832
 833   fCutsD0 = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(3));
 834   if(!fCutsD0){
 835     printf("ERROR: fCuts not available\n");
 836     return;
 837   }
 838
 839   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(4));
 840   if (!fCounter) {
 841     printf("ERROR: fCounter not available\n");
 842     return;
 843   }
 844   fOutputCorr = dynamic_cast<TList*> (GetOutputData(5));
 845   if (!fOutputCorr) {
 846     printf("ERROR: fOutputCorr not available\n");
 847     return;
 848   }
 849   fOutputStudy = dynamic_cast<TList*> (GetOutputData(6));
 850   if (!fOutputStudy) {
 851     printf("ERROR: fOutputStudy not available\n");
 852     return;
 853   }
 854   fCutsTracks = dynamic_cast<AliHFAssociatedTrackCuts*>(GetOutputData(7));
 855   if(!fCutsTracks){
 856     printf("ERROR: fCutsTracks not available\n");
 857     return;
 858   }
 859
 860   return;
 861 }
 862
 863 //_________________________________________________________________________________________________
 864 Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 865   //
 866   // checking whether the mother of the particles come from a charm or a bottom quark
 867   //
 868   printf("AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
 869
 870   Int_t pdgGranma = 0;
 871   Int_t mother = 0;
 872   mother = mcPartCandidate->GetMother();
 873   Int_t abspdgGranma =0;
 874   Bool_t isFromB=kFALSE;
 875   Bool_t isQuarkFound=kFALSE;
 876
 877   while (mother > 0){
 878     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 879     if (mcMoth){
 880       pdgGranma = mcMoth->GetPdgCode();
 881       abspdgGranma = TMath::Abs(pdgGranma);
 882       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 883         isFromB=kTRUE;
 884       }
 885       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 886       mother = mcMoth->GetMother();
 887     }else{
 888       AliError("Failed casting the mother particle!");
 889       break;
 890     }
 891   }
 892
 893   if(isQuarkFound) {
 894     if(isFromB) return 5;
 895     else return 4;
 896   }
 897   else return 1;
 898 }
 899
 900 //________________________________________________________________________
 901 void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
 902 //
 903
 904   TString namePlot = "";
 905
 906   //These for limits in THnSparse (one per bin, same limits).
 907   //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
 908   Int_t nBinsPhi[5] = {32,150,6,3,16};
 909   Double_t binMinPhi[5] = {-TMath::Pi()/2,1.6,0.,0.,-1.6};  //is the minimum for all the bins
 910   Double_t binMaxPhi[5] = {3*TMath::Pi()/2,2.2,3.0,3.,1.6};  //is the maximum for all the bins
 911
 912   //Vars: DeltaPhi, InvMass, DeltaEta
 913   Int_t nBinsMix[3] = {32,150,16};
 914   Double_t binMinMix[3] = {-TMath::Pi()/2,1.6,-1.6};  //is the minimum for all the bins
 915   Double_t binMaxMix[3] = {3*TMath::Pi()/2,2.2,1.6};  //is the maximum for all the bins
 916
 917   for(Int_t i=0;i<fNPtBinsCorr;i++){
 918
 919     if(!fMixing) {
 920       //THnSparse plots: correlations for various invariant mass (MC and data)
 921       namePlot="hPhi_K0_Bin";
 922       namePlot+=i;
 923
 924       THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 925       hPhiK->Sumw2();
 926       fOutputCorr->Add(hPhiK);
 927
 928       namePlot="hPhi_Kcharg_Bin";
 929       namePlot+=i;
 930
 931       THnSparseF *hPhiH = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 932       hPhiH->Sumw2();
 933       fOutputCorr->Add(hPhiH);
 934
 935       namePlot="hPhi_Charg_Bin";
 936       namePlot+=i;
 937
 938       THnSparseF *hPhiC = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 939       hPhiC->Sumw2();
 940       fOutputCorr->Add(hPhiC);
 941
 942       //histos for c/b origin for D0 (MC only)
 943       if (fReadMC) {
 944
 945         //generic origin for tracks
 946         namePlot="hPhi_K0_From_c_Bin";
 947         namePlot+=i;
 948
 949         THnSparseF *hPhiK_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 950         hPhiK_c->Sumw2();
 951         fOutputCorr->Add(hPhiK_c);
 952
 953         namePlot="hPhi_Kcharg_From_c_Bin";
 954         namePlot+=i;
 955
 956         THnSparseF *hPhiH_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 957         hPhiH_c->Sumw2();
 958         fOutputCorr->Add(hPhiH_c);
 959
 960         namePlot="hPhi_Charg_From_c_Bin";
 961         namePlot+=i;
 962
 963         THnSparseF *hPhiC_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 964         hPhiC_c->Sumw2();
 965         fOutputCorr->Add(hPhiC_c);
 966
 967         namePlot="hPhi_K0_From_b_Bin";
 968         namePlot+=i;
 969
 970         THnSparseF *hPhiK_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 971         hPhiK_b->Sumw2();
 972         fOutputCorr->Add(hPhiK_b);
 973
 974         namePlot="hPhi_Kcharg_From_b_Bin";
 975         namePlot+=i;
 976
 977         THnSparseF *hPhiH_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 978         hPhiH_b->Sumw2();
 979         fOutputCorr->Add(hPhiH_b);
 980
 981         namePlot="hPhi_Charg_From_b_Bin";
 982         namePlot+=i;
 983
 984         THnSparseF *hPhiC_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 985         hPhiC_b->Sumw2();
 986         fOutputCorr->Add(hPhiC_b);
 987
 988         //HF-only tracks (c for c->D0, b for b->D0)
 989         namePlot="hPhi_K0_HF_From_c_Bin";
 990         namePlot+=i;
 991
 992         THnSparseF *hPhiK_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
 993         hPhiK_HF_c->Sumw2();
 994         fOutputCorr->Add(hPhiK_HF_c);
 995
 996         namePlot="hPhi_Kcharg_HF_From_c_Bin";
 997         namePlot+=i;
 998
 999         THnSparseF *hPhiH_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1000         hPhiH_HF_c->Sumw2();
1001         fOutputCorr->Add(hPhiH_HF_c);
1002
1003         namePlot="hPhi_Charg_HF_From_c_Bin";
1004         namePlot+=i;
1005
1006         THnSparseF *hPhiC_HF_c = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1007         hPhiC_HF_c->Sumw2();
1008         fOutputCorr->Add(hPhiC_HF_c);
1009
1010         namePlot="hPhi_K0_HF_From_b_Bin";
1011         namePlot+=i;
1012
1013         THnSparseF *hPhiK_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1014         hPhiK_HF_b->Sumw2();
1015         fOutputCorr->Add(hPhiK_HF_b);
1016
1017         namePlot="hPhi_Kcharg_HF_From_b_Bin";
1018         namePlot+=i;
1019
1020         THnSparseF *hPhiH_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1021         hPhiH_HF_b->Sumw2();
1022         fOutputCorr->Add(hPhiH_HF_b);
1023
1024         namePlot="hPhi_Charg_HF_From_b_Bin";
1025         namePlot+=i;
1026
1027         THnSparseF *hPhiC_HF_b = new THnSparseF(namePlot.Data(), "Azimuthal correlation HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1028         hPhiC_HF_b->Sumw2();
1029         fOutputCorr->Add(hPhiC_HF_b);
1030
1031         namePlot="hPhi_K0_NonHF_Bin";
1032         namePlot+=i;
1033
1034         THnSparseF *hPhiK_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1035         hPhiK_Non->Sumw2();
1036         fOutputCorr->Add(hPhiK_Non);
1037
1038         namePlot="hPhi_Kcharg_NonHF_Bin";
1039         namePlot+=i;
1040
1041         THnSparseF *hPhiH_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1042         hPhiH_Non->Sumw2();
1043         fOutputCorr->Add(hPhiH_Non);
1044
1045         namePlot="hPhi_Charg_NonHF_Bin";
1046         namePlot+=i;
1047
1048         THnSparseF *hPhiC_Non = new THnSparseF(namePlot.Data(), "Azimuthal correlation - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1049         hPhiC_Non->Sumw2();
1050         fOutputCorr->Add(hPhiC_Non);
1051       }
1052
1053       //leading hadron correlations
1054       namePlot="hPhi_Lead_Bin";
1055       namePlot+=i;
1056
1057       THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1058       hCorrLead->Sumw2();
1059       fOutputCorr->Add(hCorrLead);
1060
1061       if (fReadMC) {
1062         namePlot="hPhi_Lead_From_c_Bin";
1063         namePlot+=i;
1064
1065         THnSparseF *hCorrLead_c = new THnSparseF(namePlot.Data(), "Leading particle correlations - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1066         hCorrLead_c->Sumw2();
1067         fOutputCorr->Add(hCorrLead_c);
1068
1069         namePlot="hPhi_Lead_From_b_Bin";
1070         namePlot+=i;
1071
1072         THnSparseF *hCorrLead_b = new THnSparseF(namePlot.Data(), "Leading particle correlations - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1073         hCorrLead_b->Sumw2();
1074         fOutputCorr->Add(hCorrLead_b);
1075
1076         namePlot="hPhi_Lead_HF_From_c_Bin";
1077         namePlot+=i;
1078
1079         THnSparseF *hCorrLead_HF_c = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1080         hCorrLead_HF_c->Sumw2();
1081         fOutputCorr->Add(hCorrLead_HF_c);
1082
1083         namePlot="hPhi_Lead_HF_From_b_Bin";
1084         namePlot+=i;
1085
1086         THnSparseF *hCorrLead_HF_b = new THnSparseF(namePlot.Data(), "Leading particle correlations HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1087         hCorrLead_HF_b->Sumw2();
1088         fOutputCorr->Add(hCorrLead_HF_b);
1089
1090         namePlot="hPhi_Lead_NonHF_Bin";
1091         namePlot+=i;
1092
1093         THnSparseF *hCorrLead_Non = new THnSparseF(namePlot.Data(), "Leading particle correlations - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1094         hCorrLead_Non->Sumw2();
1095         fOutputCorr->Add(hCorrLead_Non);
1096       }
1097
1098       //pT weighted correlations
1099       namePlot="hPhi_Weig_Bin";
1100       namePlot+=i;
1101
1102       THnSparseF *hCorrWeig = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted); #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1103       fOutputCorr->Add(hCorrWeig);
1104
1105       if (fReadMC) {
1106         namePlot="hPhi_Weig_From_c_Bin";
1107         namePlot+=i;
1108
1109         THnSparseF *hCorrWeig_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1110         fOutputCorr->Add(hCorrWeig_c);
1111
1112         namePlot="hPhi_Weig_From_b_Bin";
1113         namePlot+=i;
1114
1115         THnSparseF *hCorrWeig_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1116         fOutputCorr->Add(hCorrWeig_b);
1117
1118         namePlot="hPhi_Weig_HF_From_c_Bin";
1119         namePlot+=i;
1120
1121         THnSparseF *hCorrWeig_HF_c = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - c origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1122         fOutputCorr->Add(hCorrWeig_HF_c);
1123
1124         namePlot="hPhi_Weig_HF_From_b_Bin";
1125         namePlot+=i;
1126
1127         THnSparseF *hCorrWeig_HF_b = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) HF - b origin; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1128         fOutputCorr->Add(hCorrWeig_HF_b);
1129
1130         namePlot="hPhi_Weig_NonHF_Bin";
1131         namePlot+=i;
1132
1133         THnSparseF *hCorrWeig_Non = new THnSparseF(namePlot.Data(), "Charged particle correlations (pT weighted) - Non HF; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1134         fOutputCorr->Add(hCorrWeig_Non);
1135       }
1136
1137     //pT distribution histos
1138     namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
1139     TH1F *hPtC = new TH1F(namePlot.Data(), "Charged track pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1140     hPtC->SetMinimum(0);
1141     fOutputStudy->Add(hPtC);
1142
1143     namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
1144     TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1145     hPtH->SetMinimum(0);
1146     fOutputStudy->Add(hPtH);
1147
1148     namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
1149     TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1150     hPtK->SetMinimum(0);
1151     fOutputStudy->Add(hPtK);
1152
1153     //D* feeddown pions rejection histos
1154     namePlot = "hDstarPions_Bin"; namePlot+=i;
1155     TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,150,1.6,2.2);
1156     hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
1157     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
1158     hDstarPions->SetMinimum(0);
1159     fOutputStudy->Add(hDstarPions);
1160
1161     }
1162
1163     if(fMixing) {
1164       //THnSparse plots for event mixing!
1165       namePlot="hPhi_K0_Bin";
1166       namePlot+=i;namePlot+="_EvMix";
1167
1168       THnSparseF *hPhiK_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1169       hPhiK_EvMix->Sumw2();
1170       fOutputCorr->Add(hPhiK_EvMix);
1171
1172       namePlot="hPhi_Kcharg_Bin";
1173       namePlot+=i;namePlot+="_EvMix";
1174
1175       THnSparseF *hPhiH_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1176       hPhiH_EvMix->Sumw2();
1177       fOutputCorr->Add(hPhiH_EvMix);
1178
1179       namePlot="hPhi_Charg_Bin";
1180       namePlot+=i;namePlot+="_EvMix";
1181
1182       THnSparseF *hPhiC_EvMix = new THnSparseF(namePlot.Data(), "Az. corr. EvMix; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",3,nBinsMix,binMinMix,binMaxMix);
1183       hPhiC_EvMix->Sumw2();
1184       fOutputCorr->Add(hPhiC_EvMix);
1185     }
1186   }
1187
1188   //out of bin loop
1189   if(!fMixing) {
1190     TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1191     hCountC->SetMinimum(0);
1192     fOutputStudy->Add(hCountC);
1193
1194     TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1195     hCountH->SetMinimum(0);
1196     fOutputStudy->Add(hCountH);
1197
1198     TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1199     hCountK->SetMinimum(0);
1200     fOutputStudy->Add(hCountK);
1201   }
1202
1203   if (fReadMC) {
1204     TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1205     fOutputStudy->Add(hEventTypeMC);
1206   }
1207
1208   if (fFillGlobal) { //all-events plots
1209     //pt distributions
1210     TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1211     hPtCAll->SetMinimum(0);
1212     fOutputStudy->Add(hPtCAll);
1213
1214     TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1215     hPtHAll->SetMinimum(0);
1216     fOutputStudy->Add(hPtHAll);
1217
1218     TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1219     hPtKAll->SetMinimum(0);
1220     fOutputStudy->Add(hPtKAll);
1221
1222     if(!fMixing) {
1223       //phi distributions
1224       TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1225       hPhiDistCAll->SetMinimum(0);
1226       fOutputStudy->Add(hPhiDistCAll);
1227
1228       TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1229       hPhiDistHAll->SetMinimum(0);
1230       fOutputStudy->Add(hPhiDistHAll);
1231
1232       TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1233       hPhiDistKAll->SetMinimum(0);
1234       fOutputStudy->Add(hPhiDistKAll);
1235
1236       TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1237       hPhiDistDAll->SetMinimum(0);
1238       fOutputStudy->Add(hPhiDistDAll);
1239     }
1240
1241     //K0 Invariant Mass plots
1242     TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1243     hK0MassInv->SetMinimum(0);
1244     fOutputStudy->Add(hK0MassInv);
1245   }
1246
1247   //for MC analysis only
1248   for(Int_t i=0;i<fNPtBinsCorr;i++) {
1249
1250     if (fReadMC && !fMixing) {
1251
1252       //displacement histos
1253       namePlot="histDispl_K0_Bin"; namePlot+=i;
1254       TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1255       hDisplK->SetMinimum(0);
1256       fOutputStudy->Add(hDisplK);
1257
1258       namePlot="histDispl_K0_HF_Bin";  namePlot+=i;
1259       TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1260       hDisplK_HF->SetMinimum(0);
1261       fOutputStudy->Add(hDisplK_HF);
1262
1263       namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1264       TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1265       hDisplHadr->SetMinimum(0);
1266       fOutputStudy->Add(hDisplHadr);
1267
1268       namePlot="histDispl_Kcharg_HF_Bin";  namePlot+=i;
1269       TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1270       hDisplHadr_HF->SetMinimum(0);
1271       fOutputStudy->Add(hDisplHadr_HF);
1272
1273       namePlot="histDispl_Charg_Bin"; namePlot+=i;
1274       TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1275       hDisplCharg->SetMinimum(0);
1276       fOutputStudy->Add(hDisplCharg);
1277
1278       namePlot="histDispl_Charg_HF_Bin";  namePlot+=i;
1279       TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1280       hDisplCharg_HF->SetMinimum(0);
1281       fOutputStudy->Add(hDisplCharg_HF);
1282
1283       namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
1284       TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1285       hDisplK_c->SetMinimum(0);
1286       fOutputStudy->Add(hDisplK_c);
1287
1288       namePlot="histDispl_K0_HF_From_c_Bin";  namePlot+=i;
1289       TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1290       hDisplK_HF_c->SetMinimum(0);
1291       fOutputStudy->Add(hDisplK_HF_c);
1292
1293       namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1294       TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1295       hDisplHadr_c->SetMinimum(0);
1296       fOutputStudy->Add(hDisplHadr_c);
1297
1298       namePlot="histDispl_Kcharg_HF_From_c_Bin";  namePlot+=i;
1299       TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1300       hDisplHadr_HF_c->SetMinimum(0);
1301       fOutputStudy->Add(hDisplHadr_HF_c);
1302
1303       namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1304       TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1305       hDisplCharg_c->Sumw2();
1306       hDisplCharg_c->SetMinimum(0);
1307       fOutputStudy->Add(hDisplCharg_c);
1308
1309       namePlot="histDispl_Charg_HF_From_c_Bin";  namePlot+=i;
1310       TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1311       hDisplCharg_HF_c->SetMinimum(0);
1312       fOutputStudy->Add(hDisplCharg_HF_c);
1313
1314       namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
1315       TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1316       hDisplK_b->SetMinimum(0);
1317       fOutputStudy->Add(hDisplK_b);
1318
1319       namePlot="histDispl_K0_HF_From_b_Bin";  namePlot+=i;
1320       TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1321       hDisplK_HF_b->SetMinimum(0);
1322       fOutputStudy->Add(hDisplK_HF_b);
1323
1324       namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1325       TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1326       hDisplHadr_b->SetMinimum(0);
1327       fOutputStudy->Add(hDisplHadr_b);
1328
1329       namePlot="histDispl_Kcharg_HF_From_b_Bin";  namePlot+=i;
1330       TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1331       hDisplHadr_HF_b->SetMinimum(0);
1332       fOutputStudy->Add(hDisplHadr_HF_b);
1333
1334       namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1335       TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1336       hDisplCharg_b->SetMinimum(0);
1337       fOutputStudy->Add(hDisplCharg_b);
1338
1339       namePlot="histDispl_Charg_HF_From_b_Bin";  namePlot+=i;
1340       TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1341       hDisplCharg_HF_b->SetMinimum(0);
1342       fOutputStudy->Add(hDisplCharg_HF_b);
1343
1344       //origin of tracks histos
1345       namePlot="histOrig_Charg_Bin";  namePlot+=i;
1346       TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1347       hOrigin_Charm->SetMinimum(0);
1348       hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1349       hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1350       hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1351       hOrigin_Charm->GetXaxis()->SetBinLabel(4,"c hadr.");
1352       hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->#");
1353       hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1354       hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->#");
1355       hOrigin_Charm->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1356       hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1357       fOutputStudy->Add(hOrigin_Charm);
1358
1359       namePlot="histOrig_Kcharg_Bin";  namePlot+=i;
1360       TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1361       hOrigin_Kcharg->SetMinimum(0);
1362       hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1363       hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1364       hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1365       hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"c hadr.");
1366       hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->#");
1367       hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1368       hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->#");
1369       hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1370       hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1371       fOutputStudy->Add(hOrigin_Kcharg);
1372
1373       namePlot="histOrig_K0_Bin";  namePlot+=i;
1374       TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1375       hOrigin_K->SetMinimum(0);
1376       hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1377       hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1378       hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1379       hOrigin_K->GetXaxis()->SetBinLabel(4,"c hadr.");
1380       hOrigin_K->GetXaxis()->SetBinLabel(5,"B->#");
1381       hOrigin_K->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1382       hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->#");
1383       hOrigin_K->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1384       hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1385       fOutputStudy->Add(hOrigin_K);
1386     }
1387
1388     if (fReadMC) {
1389       //origin of D0 histos
1390       namePlot="histOrig_D0_Bin";  namePlot+=i;
1391       TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1392       hOrigin_D0->SetMinimum(0);
1393       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1394       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1395       fOutputStudy->Add(hOrigin_D0);
1396     }
1397   }
1398 }
1399
1400 //________________________________________________________________________
1401 void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1402 //
1403 // Method for correlations D0-hadrons study
1404 //
1405   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1406   Double_t mD0, mD0bar;
1407   Int_t origD0 = 0, PDGD0 = 0, ptbin = 0;
1408   d->InvMassD0(mD0,mD0bar);
1409   Double_t mInv[2] = {mD0, mD0bar};
1410   ptbin = PtBinCorr(d->Pt());
1411
1412   if(ptbin < 0) return;
1413
1414   //Fill of D0 phi distribution
1415   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1416
1417   //Origin of D0
1418   TString orig="";
1419   if(fReadMC) {
1420     origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1421     PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
1422     switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
1423       case 4:
1424         orig = "_From_c";
1425         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1426         break;
1427       case 5:
1428         orig = "_From_b";
1429         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1430         break;
1431       default:
1432         return;
1433     }
1434   }
1435
1436   Double_t highPt = 0; Double_t lead[4] = {0,0,0,1};  //infos for leading particle (pt,deltaphi)
1437
1438   //loop over the tracks in the pool
1439   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1440   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1441   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1442
1443   Int_t NofEventsinPool = 1;
1444   if(fMixing) {
1445     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1446     if(!execPoolTr) {
1447       AliInfo("Mixed event analysis: track pool is not ready");
1448       NofEventsinPool = 0;
1449     }
1450   }
1451
1452   //Charged tracks
1453   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1454     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1455     if(!analyzetracksTr) {
1456       AliInfo("AliHFCorrelator::Cannot process the track array");
1457       continue;
1458     }
1459
1460     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1461
1462       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1463       if(!runcorrelation) continue;
1464
1465       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1466
1467       if(!fMixing) {
1468         /*if(!track->CheckSoftPi()) { //removal of soft pions
1469           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1470           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1471           continue;
1472         } else { //not a soft pion
1473           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1474           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1475         }*/
1476         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1477         if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1478       }
1479       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1480
1481       Double_t eff = track->GetWeight(); //extract track efficiency
1482
1483       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks, weight = 1./eff
1484
1485       if(!fMixing) N_Charg++;
1486
1487       //retrieving leading info...
1488       if(track->Pt() > highPt) {
1489         if(fReadMC && track->GetLabel()<1) continue;
1490         if(fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) continue;
1491         lead[0] = fCorrelatorTr->GetDeltaPhi();
1492         lead[1] = fCorrelatorTr->GetDeltaEta();
1493         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1494         lead[3] = 1./track->GetWeight(); //weight is 1./efficiency
1495         highPt = track->Pt();
1496       }
1497
1498     } // end of tracks loop
1499   } //end of event loop for fCorrelatorTr
1500
1501   if(fMixing) {
1502     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1503     if(!execPoolKc) {
1504       AliInfo("Mixed event analysis: K+/- pool is not ready");
1505       NofEventsinPool = 0;
1506     }
1507   }
1508
1509   //Charged Kaons loop
1510   for (Int_t jMix = 0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1511     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1512     if(!analyzetracksKc) {
1513       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1514       continue;
1515     }
1516
1517     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1518
1519       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1520       if(!runcorrelation) continue;
1521
1522       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1523
1524       if(!fMixing) {
1525         if(!kCharg->CheckSoftPi()) { //removal of soft pions
1526           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1527           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1528           continue;
1529         } else {
1530           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1531           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1532         }
1533         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1534         if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1535       }
1536       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1537
1538       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1539
1540       if(!fMixing) N_KCharg++;
1541
1542     } // end of charged kaons loop
1543   } //end of event loop for fCorrelatorKc
1544
1545   if(fMixing) {
1546     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1547     if(!execPoolK0) {
1548       AliInfo("Mixed event analysis: K0 pool is not ready");
1549       NofEventsinPool = 0;
1550     }
1551   }
1552
1553   //K0 loop
1554   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1555     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1556     if(!analyzetracksK0) {
1557       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1558       continue;
1559     }
1560
1561     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1562
1563       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1564       if(!runcorrelation) continue;
1565
1566       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1567
1568       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1569
1570       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1571
1572       if(!fMixing) N_Kaons++;
1573
1574     } // end of charged kaons loop
1575   } //end of event loop for fCorrelatorK0
1576
1577   Double_t fillSpLeadD0[3] = {lead[0],mD0,lead[1]};
1578   Double_t fillSpLeadD0bar[3] = {lead[0],mD0bar,lead[1]};
1579
1580   //leading track correlations fill
1581   if(!fMixing) {
1582     if(fReadMC) {
1583       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1584         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
1585         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
1586         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1587         if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]);
1588         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);  //non HF
1589       }
1590       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1591         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1592         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
1593         if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1594         if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1595         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  //non HF
1596       }
1597     } else {
1598         if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);
1599         if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1600     }
1601
1602     //Fill of count histograms
1603     if (!fAlreadyFilled) {
1604       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1605       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1606       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1607     }
1608   }
1609
1610   fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1611
1612 }
1613
1614 //________________________________________________________________________
1615 void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCParticle* d, TClonesArray* mcArray) {
1616 //
1617 // Method for correlations D0-hadrons study
1618 //
1619   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1620   Double_t mD0 = 1.864, mD0bar = 1.864;
1621   Double_t mInv[2] = {mD0, mD0bar};
1622   Int_t origD0 = 0, PDGD0 = 0;
1623   Int_t ptbin = PtBinCorr(d->Pt());
1624
1625   if(ptbin < 0) return;
1626
1627   //Fill of D0 phi distribution
1628   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1629
1630   //Origin of D0
1631   TString orig="";
1632   origD0=CheckD0Origin(mcArray,d);
1633   PDGD0 = d->GetPdgCode();
1634   switch (CheckD0Origin(mcArray,d)) {
1635     case 4:
1636       orig = "_From_c";
1637       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1638       break;
1639     case 5:
1640       orig = "_From_b";
1641       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1642       break;
1643     default:
1644       return;
1645   }
1646
1647   Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
1648
1649   //loop over the tracks in the pool
1650   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1651   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1652   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1653
1654   Int_t NofEventsinPool = 1;
1655   if(fMixing) {
1656     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1657     if(!execPoolTr) {
1658       AliInfo("Mixed event analysis: track pool is not ready");
1659       NofEventsinPool = 0;
1660     }
1661   }
1662
1663   //Charged tracks
1664   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1665
1666     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1667     if(!analyzetracksTr) {
1668       AliInfo("AliHFCorrelator::Cannot process the track array");
1669       continue;
1670     }
1671
1672     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1673
1674       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1675       if(!runcorrelation) continue;
1676
1677       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1678       if(track->GetLabel()<0) continue;
1679       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1680       if(track->Pt() < 0.3 || TMath::Abs(track->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1681       if(!fMixing) N_Charg++;
1682
1683       AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1684       if(!trkMC) continue;
1685
1686       if (!trkMC->IsPhysicalPrimary()) continue; //reject material budget, or other fake tracks
1687       if (IsDDaughter(d,trkMC,mcArray)) continue;
1688       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
1689
1690       //retrieving leading info...
1691       if(track->Pt() > highPt) {
1692         lead[0] = fCorrelatorTr->GetDeltaPhi();
1693         lead[1] = fCorrelatorTr->GetDeltaEta();
1694         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,trkMC) : 0;
1695         highPt = track->Pt();
1696       }
1697
1698     } // end of tracks loop
1699   } //end of event loop for fCorrelatorTr
1700
1701   if(fMixing) {
1702     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1703     if(!execPoolKc) {
1704       AliInfo("Mixed event analysis: K+/- pool is not ready");
1705       NofEventsinPool = 0;
1706     }
1707   }
1708
1709   //Charged Kaons loop
1710   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1711     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1712     if(!analyzetracksKc) {
1713       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1714       continue;
1715     }
1716
1717     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1718
1719       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1720       if(!runcorrelation) continue;
1721
1722       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1723       if(kCharg->GetLabel()<1) continue;
1724       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1725       if(TMath::Abs(kCharg->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1726       if(!fMixing) N_KCharg++;
1727
1728       AliAODMCParticle *kChargMC = (AliAODMCParticle*)mcArray->At(kCharg->GetLabel());
1729       if(!kChargMC) continue;
1730
1731       if (IsDDaughter(d,kChargMC,mcArray)) continue;
1732       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1733
1734     } // end of charged kaons loop
1735   } //end of event loop for fCorrelatorKc
1736
1737   if(fMixing) {
1738     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1739     if(!execPoolK0) {
1740       AliInfo("Mixed event analysis: K0 pool is not ready");
1741       NofEventsinPool = 0;
1742     }
1743   }
1744
1745   //K0 loop
1746   for (Int_t jMix =0; jMix < NofEventsinPool; jMix++) {// loop on events in the pool; if it is SE analysis, stops at one (index not needed there)
1747     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1748     if(!analyzetracksK0) {
1749       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1750       continue;
1751     }
1752
1753     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1754
1755       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1756       if(!runcorrelation) continue;
1757
1758       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1759       if(k0->GetLabel()<1) continue;
1760       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1761       if(TMath::Abs(k0->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1762
1763       AliAODMCParticle *k0MC = (AliAODMCParticle*)mcArray->At(k0->GetLabel());
1764       if(!k0MC) continue;
1765
1766       if (IsDDaughter(d,k0MC,mcArray)) continue;
1767       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1768
1769       if(!fMixing) N_Kaons++;
1770
1771     } // end of charged kaons loop
1772   } //end of event loop for fCorrelatorK0
1773
1774   Double_t fillSpLeadMC[3] = {lead[0],mD0,lead[1]}; //mD0 = mD0bar = 1.864
1775
1776   //leading track correlations fill
1777   if(!fMixing) {
1778     if(d->GetPdgCode()==421) { //D0
1779       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
1780       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
1781       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1782       if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1783       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
1784     }
1785     if(d->GetPdgCode()==-421) { //D0bar
1786       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
1787       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
1788       if(origD0==4&&(int)lead[2]>=1&&(int)lead[2]<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1789       if(origD0==5&&(int)lead[2]>=4&&(int)lead[2]<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC);
1790       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
1791     }
1792
1793     //Fill of count histograms
1794     if (!fAlreadyFilled) {
1795       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1796       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1797       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1798     }
1799
1800     fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1801   }
1802
1803 }
1804
1805 //________________________________________________________________________
1806 void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, Double_t mInv[], Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t ptbin, Int_t type, Double_t wg) {
1807   //
1808   //fills the THnSparse for correlations, calculating the variables
1809   //
1810
1811   //Initialization of variables
1812   Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
1813   mD0 = mInv[0];
1814   mD0bar = mInv[1];
1815
1816   if (fReadMC && track->GetLabel()<1) return;
1817   if (fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) return;
1818   Double_t ptTrack = track->Pt();
1819   Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
1820   Double_t phiTr = track->Phi();
1821   Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1822
1823   TString part = "", orig = "";
1824
1825   switch (type) {
1826     case(kTrack): {
1827       part = "Charg";
1828       deltaphi = fCorrelatorTr->GetDeltaPhi();
1829       deltaeta = fCorrelatorTr->GetDeltaEta();
1830       break;
1831     }
1832     case(kKCharg): {
1833       part = "Kcharg";
1834       deltaphi = fCorrelatorKc->GetDeltaPhi();
1835       deltaeta = fCorrelatorKc->GetDeltaEta();
1836       break;
1837     }
1838     case(kK0): {
1839       part = "K0";
1840       deltaphi = fCorrelatorK0->GetDeltaPhi();
1841       deltaeta = fCorrelatorK0->GetDeltaEta();
1842       break;
1843     }
1844   }
1845
1846   if(fMixing == kSE) {
1847
1848     //Fixes limits; needed to include overflow into THnSparse projections!
1849     Double_t pTorig = track->Pt();
1850     Double_t d0orig = track->GetImpPar();
1851     Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1852     Double_t displLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
1853     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
1854     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
1855     if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
1856     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1857     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1858
1859     //variables for filling histos
1860     Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
1861     Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
1862     Double_t fillSpWeigD0[3] = {deltaphi,mD0,deltaeta};
1863     Double_t fillSpWeigD0bar[3] = {deltaphi,mD0bar,deltaeta};
1864
1865     if(fReadMC == 0) {
1866       //sparse fill for data (tracks, K+-, K0) + weighted
1867       if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1868         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1869         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1870       }
1871       if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1872         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1873         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1874       }
1875       if(!fAlreadyFilled) {
1876         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1877         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1878       }
1879     }
1880
1881     if(fReadMC) {
1882
1883       if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
1884
1885       //sparse fill for data (tracks, K+-, K0) + weighted
1886       if(PdgD0==421) { //D0 (from MCTruth)
1887          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1888          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1889          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1890          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1891          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1892          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1893          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1894          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1895          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1896          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
1897       }
1898       if(PdgD0==-421) { //D0bar
1899          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1900          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1901          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1902          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1903          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1904          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1905          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1906          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1907          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1908          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
1909       }
1910       if(!fAlreadyFilled) {
1911         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1912         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
1913         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
1914         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1915         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1916         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
1917         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1918       }
1919     }//end MC case
1920
1921   } //end of SE fill
1922
1923   if(fMixing == kME) {
1924
1925     //Fixes limits; needed to include overflow into THnSparse projections!
1926     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
1927     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1928     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1929
1930     //variables for filling histos
1931     Double_t fillSpPhiD0[3] = {deltaphi,mD0,deltaeta};
1932     Double_t fillSpPhiD0bar[3] = {deltaphi,mD0bar,deltaeta};
1933
1934     if(fReadMC == 0) {
1935       //sparse fill for data (tracks, K+-, K0)
1936       if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1937       if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1938     }
1939     if(fReadMC == 1) {
1940       //sparse fill for data (tracks, K+-, K0)
1941       if(PdgD0==421)  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
1942       if(PdgD0==-421) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
1943     }//end MC case
1944
1945   } //end of ME fill
1946
1947   return;
1948 }
1949
1950 //_________________________________________________________________________________________________
1951 Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1952   //
1953   // checks on particle (#) origin:
1954   // 0) Not HF
1955   // 1) D->#
1956   // 2) D->X->#
1957   // 3) c hadronization
1958   // 4) B->#
1959   // 5) B->X-># (X!=D)
1960   // 6) B->D->#
1961   // 7) B->D->X->#
1962   // 8) b hadronization
1963   //
1964   if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
1965
1966   Int_t pdgGranma = 0;
1967   Int_t mother = 0;
1968   mother = mcPartCandidate->GetMother();
1969   Int_t istep = 0;
1970   Int_t abspdgGranma =0;
1971   Bool_t isFromB=kFALSE;
1972   Bool_t isDdaugh=kFALSE;
1973   Bool_t isDchaindaugh=kFALSE;
1974   Bool_t isBdaugh=kFALSE;
1975   Bool_t isBchaindaugh=kFALSE;
1976   Bool_t isQuarkFound=kFALSE;
1977
1978   if (mother<0) return -1;
1979   while (mother >= 0){
1980     istep++;
1981     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1982     if (mcMoth){
1983       pdgGranma = mcMoth->GetPdgCode();
1984       abspdgGranma = TMath::Abs(pdgGranma);
1985       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1986         isBchaindaugh=kTRUE;
1987         if(istep==1) isBdaugh=kTRUE;
1988       }
1989       if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
1990         isDchaindaugh=kTRUE;
1991         if(istep==1) isDdaugh=kTRUE;
1992       }
1993       if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
1994       mother = mcMoth->GetMother();
1995     }else{
1996       AliError("Failed casting the mother particle!");
1997       return -1;
1998     }
1999   }
2000
2001   //decides what to return based on the flag status
2002   if(isQuarkFound) {
2003     if(!isFromB) {  //charm
2004       if(isDdaugh) return 1; //charm immediate
2005       else if(isDchaindaugh) return 2; //charm chain
2006       else return 3; //charm hadronization
2007     }
2008     else { //beauty
2009       if(isBdaugh) return 4; //b immediate
2010       else if(isBchaindaugh) { //b chain
2011         if(isDchaindaugh) {
2012           if(isDdaugh) return 6; //d immediate
2013           return 7; //d chain
2014           }
2015         else return 5; //b, not d
2016       }
2017       else return 8; //b hadronization
2018     }
2019   }
2020   else if(!isDdaugh && !isDchaindaugh && !isBdaugh && !isBchaindaugh) return 0; //no HF decay
2021      //in this case, it's !isQuarkFound, but not in 100% cases it's a non HF particle!
2022      //rarely you can find a D/B meson which comes from a -1! It isn't a Non-HF, in that case! And I'll return -1...
2023
2024   return -1; //some problem spotted
2025 }
2026 //________________________________________________________________________
2027 Bool_t AliAnalysisTaskSED0Correlations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
2028
2029   //Daughter removal in MCKine case
2030   Bool_t isDaughter = kFALSE;
2031   Int_t labelD0 = d->GetLabel();
2032
2033   Int_t mother = track->GetMother();
2034
2035   //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
2036   while (mother > 0){
2037     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother)); //it's the mother of the track!
2038     if (mcMoth){
2039       if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
2040       mother = mcMoth->GetMother(); //goes back by one
2041     } else{
2042       AliError("Failed casting the mother particle!");
2043       break;
2044     }
2045   }
2046
2047   return isDaughter;
2048 }
2049
2050 //________________________________________________________________________
2051 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
2052   //
2053   //give the pt bin where the pt lies.
2054   //
2055   Int_t ptbin=-1;
2056   if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
2057
2058   Int_t i = 0;
2059   while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
2060
2061   return ptbin;
2062 }
2063
2064 //---------------------------------------------------------------------------
2065 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
2066 {
2067   //
2068   // Selection for K0 hypotheses
2069   // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
2070   //          2 = no previous selections
2071
2072   if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
2073
2074   AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
2075   AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
2076
2077   if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
2078     if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
2079   }
2080
2081   //This part removes double counting for swapped tracks!
2082   Int_t i = 0;  //while loop (until the last-written entry pair of ID!
2083   while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
2084     if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
2085        (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
2086     i++;
2087   }
2088   idArrayV0[i][0]=v0Daug1->GetID();
2089   idArrayV0[i][1]=v0Daug2->GetID();
2090
2091   return kTRUE;
2092 }
2093
2094 //________________________________________________________________________
2095 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
2096
2097   cout << "--------------------------\n";
2098   cout << "PtBins = " << fNPtBinsCorr << "\n";
2099   cout << "PtBin limits--------------\n";
2100   for (int i=0; i<fNPtBinsCorr; i++) {
2101     cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
2102   }
2103   cout << "\n--------------------------\n";
2104   cout << "PtBin tresh. tracks low---\n";
2105   for (int i=0; i<fNPtBinsCorr; i++) {
2106     cout << fPtThreshLow.at(i) << "; ";
2107   }
2108   cout << "PtBin tresh. tracks up----\n";
2109   for (int i=0; i<fNPtBinsCorr; i++) {
2110     cout << fPtThreshUp.at(i) << "; ";
2111   }
2112   cout << "\n--------------------------\n";
2113   cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
2114   cout << "--------------------------\n";
2115   cout << "MC Truth = "<<fReadMC<<" - MC Reco: Trk = "<<fRecoTr<<", D0 = "<<fRecoD0<<"\n";
2116   cout << "--------------------------\n";
2117   cout << "Sel of Event tpye (PP,GS,FE,...)= "<<fSelEvType<<"\n";
2118   cout << "--------------------------\n";
2119   cout << "Ev Mixing = "<<fMixing<<"\n";
2120   cout << "--------------------------\n";
2121 }
2122