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