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   TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1051   hCountC->SetMinimum(0);
1052   fOutputStudy->Add(hCountC);
1053
1054   TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1055   hCountH->SetMinimum(0);
1056   fOutputStudy->Add(hCountH);
1057
1058   TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1059   hCountK->SetMinimum(0);
1060   fOutputStudy->Add(hCountK);
1061
1062   if (fFillGlobal) { //all-events plots
1063     //pt distributions
1064     TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1065     hPtCAll->SetMinimum(0);
1066     fOutputStudy->Add(hPtCAll);
1067
1068     TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1069     hPtHAll->SetMinimum(0);
1070     fOutputStudy->Add(hPtHAll);
1071
1072     TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1073     hPtKAll->SetMinimum(0);
1074     fOutputStudy->Add(hPtKAll);
1075
1076     //phi distributions
1077     TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
1078     hPhiDistCAll->SetMinimum(0);
1079     fOutputStudy->Add(hPhiDistCAll);
1080
1081     TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
1082     hPhiDistHAll->SetMinimum(0);
1083     fOutputStudy->Add(hPhiDistHAll);
1084
1085     TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
1086     hPhiDistKAll->SetMinimum(0);
1087     fOutputStudy->Add(hPhiDistKAll);
1088
1089     TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",32,-1.6,4.8);
1090     hPhiDistDAll->SetMinimum(0);
1091     fOutputStudy->Add(hPhiDistDAll);
1092
1093     //K0 Invariant Mass plots
1094     TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1095     hK0MassInv->SetMinimum(0);
1096     fOutputStudy->Add(hK0MassInv);
1097   }
1098
1099   //for MC analysis only
1100   if (fReadMC) {
1101
1102     for(Int_t i=0;i<fNPtBinsCorr;i++){
1103
1104       //displacement histos
1105       namePlot="histDispl_K0_Bin"; namePlot+=i;
1106       TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1107       hDisplK->SetMinimum(0);
1108       fOutputStudy->Add(hDisplK);
1109
1110       namePlot="histDispl_K0_HF_Bin";  namePlot+=i;
1111       TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1112       hDisplK_HF->SetMinimum(0);
1113       fOutputStudy->Add(hDisplK_HF);
1114
1115       namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1116       TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1117       hDisplHadr->SetMinimum(0);
1118       fOutputStudy->Add(hDisplHadr);
1119
1120       namePlot="histDispl_Kcharg_HF_Bin";  namePlot+=i;
1121       TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1122       hDisplHadr_HF->SetMinimum(0);
1123       fOutputStudy->Add(hDisplHadr_HF);
1124
1125       namePlot="histDispl_Charg_Bin"; namePlot+=i;
1126       TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1127       hDisplCharg->SetMinimum(0);
1128       fOutputStudy->Add(hDisplCharg);
1129
1130       namePlot="histDispl_Charg_HF_Bin";  namePlot+=i;
1131       TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1132       hDisplCharg_HF->SetMinimum(0);
1133       fOutputStudy->Add(hDisplCharg_HF);
1134
1135       namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
1136       TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1137       hDisplK_c->SetMinimum(0);
1138       fOutputStudy->Add(hDisplK_c);
1139
1140       namePlot="histDispl_K0_HF_From_c_Bin";  namePlot+=i;
1141       TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1142       hDisplK_HF_c->SetMinimum(0);
1143       fOutputStudy->Add(hDisplK_HF_c);
1144
1145       namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1146       TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1147       hDisplHadr_c->SetMinimum(0);
1148       fOutputStudy->Add(hDisplHadr_c);
1149
1150       namePlot="histDispl_Kcharg_HF_From_c_Bin";  namePlot+=i;
1151       TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1152       hDisplHadr_HF_c->SetMinimum(0);
1153       fOutputStudy->Add(hDisplHadr_HF_c);
1154
1155       namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1156       TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1157       hDisplCharg_c->Sumw2();
1158       hDisplCharg_c->SetMinimum(0);
1159       fOutputStudy->Add(hDisplCharg_c);
1160
1161       namePlot="histDispl_Charg_HF_From_c_Bin";  namePlot+=i;
1162       TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1163       hDisplCharg_HF_c->SetMinimum(0);
1164       fOutputStudy->Add(hDisplCharg_HF_c);
1165
1166       namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
1167       TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1168       hDisplK_b->SetMinimum(0);
1169       fOutputStudy->Add(hDisplK_b);
1170
1171       namePlot="histDispl_K0_HF_From_b_Bin";  namePlot+=i;
1172       TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1173       hDisplK_HF_b->SetMinimum(0);
1174       fOutputStudy->Add(hDisplK_HF_b);
1175
1176       namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1177       TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1178       hDisplHadr_b->SetMinimum(0);
1179       fOutputStudy->Add(hDisplHadr_b);
1180
1181       namePlot="histDispl_Kcharg_HF_From_b_Bin";  namePlot+=i;
1182       TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1183       hDisplHadr_HF_b->SetMinimum(0);
1184       fOutputStudy->Add(hDisplHadr_HF_b);
1185
1186       namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1187       TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1188       hDisplCharg_b->SetMinimum(0);
1189       fOutputStudy->Add(hDisplCharg_b);
1190
1191       namePlot="histDispl_Charg_HF_From_b_Bin";  namePlot+=i;
1192       TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1193       hDisplCharg_HF_b->SetMinimum(0);
1194       fOutputStudy->Add(hDisplCharg_HF_b);
1195
1196       //origin of tracks histos
1197       namePlot="histOrig_Charg_Bin";  namePlot+=i;
1198       TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1199       hOrigin_Charm->SetMinimum(0);
1200       hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1201       hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1202       hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1203       hOrigin_Charm->GetXaxis()->SetBinLabel(4,"B->#");
1204       hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1205       hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->D->#");
1206       hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1207       hOrigin_Charm->GetXaxis()->SetBinLabel(8,"c hadr.");
1208       hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1209       fOutputStudy->Add(hOrigin_Charm);
1210
1211       namePlot="histOrig_Kcharg_Bin";  namePlot+=i;
1212       TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1213       hOrigin_Kcharg->SetMinimum(0);
1214       hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1215       hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1216       hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1217       hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"B->#");
1218       hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1219       hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->D->#");
1220       hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1221       hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"c hadr.");
1222       hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1223       fOutputStudy->Add(hOrigin_Kcharg);
1224
1225       namePlot="histOrig_K0_Bin";  namePlot+=i;
1226       TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1227       hOrigin_K->SetMinimum(0);
1228       hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1229       hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1230       hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1231       hOrigin_K->GetXaxis()->SetBinLabel(4,"B->#");
1232       hOrigin_K->GetXaxis()->SetBinLabel(5,"B->X-># (X!=D)");
1233       hOrigin_K->GetXaxis()->SetBinLabel(6,"B->D->#");
1234       hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->X->#");
1235       hOrigin_K->GetXaxis()->SetBinLabel(8,"c hadr.");
1236       hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1237       fOutputStudy->Add(hOrigin_K);
1238
1239       //origin of D0 histos
1240       namePlot="histOrig_D0_Bin";  namePlot+=i;
1241       TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1242       hOrigin_D0->SetMinimum(0);
1243       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1244       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1245       fOutputStudy->Add(hOrigin_D0);
1246     }
1247   }
1248
1249 }
1250
1251 //________________________________________________________________________
1252 void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1253 //
1254 // Method for correlations D0-hadrons study
1255 //
1256   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1257   Double_t mD0, mD0bar;
1258   Int_t origD0 = 0, PDGD0 = 0;
1259   d->InvMassD0(mD0,mD0bar);
1260   Int_t ptbin = PtBinCorr(d->Pt());
1261   if(ptbin < 0) return;
1262
1263   //Origin of D0
1264   TString orig="";
1265   if(fReadMC) {
1266     origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1267     PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
1268     switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
1269       case 4:
1270         orig = "_From_c";
1271         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1272         break;
1273       case 5:
1274         orig = "_From_b";
1275         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1276         break;
1277       default:
1278         return;
1279     }
1280   }
1281
1282   Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
1283
1284   //loop over the tracks in the pool
1285   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1286   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1287   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1288
1289   Int_t NofEventsinPool = 1;
1290   if(fMixing) {
1291     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1292     if(!execPoolTr) {
1293       AliInfo("Mixed event analysis: track pool is not ready");
1294       NofEventsinPool = 0;
1295     }
1296   }
1297
1298   //Charged tracks
1299   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)
1300     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1301     if(!analyzetracksTr) {
1302       AliInfo("AliHFCorrelator::Cannot process the track array");
1303       continue;
1304     }
1305
1306     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1307
1308       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1309       if(!runcorrelation) continue;
1310
1311       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1312
1313       if(!fMixing) {
1314         if(!track->CheckSoftPi()) { //removal of soft pions
1315           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1316           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1317           continue;
1318         } else { //not a soft pion
1319           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1320           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1321         }
1322       Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1323       if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1324       }
1325       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1326
1327       FillSparsePlots(mcArray,d,origD0,PDGD0,track,kTrack); //fills for charged tracks
1328
1329       if(!fMixing) N_Charg++;
1330
1331       //retrieving leading info...
1332       if(track->Pt() > highPt) {
1333         lead[0] = fCorrelatorTr->GetDeltaPhi();
1334         lead[1] = fCorrelatorTr->GetDeltaEta();
1335         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1336         highPt = track->Pt();
1337       }
1338
1339     } // end of tracks loop
1340   } //end of event loop for fCorrelatorTr
1341
1342   if(fMixing) {
1343     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1344     if(!execPoolKc) {
1345       AliInfo("Mixed event analysis: K+/- pool is not ready");
1346       NofEventsinPool = 0;
1347     }
1348   }
1349
1350   //Charged Kaons loop
1351   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)
1352     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1353     if(!analyzetracksKc) {
1354       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1355       continue;
1356     }
1357
1358     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1359
1360       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1361       if(!runcorrelation) continue;
1362
1363       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1364
1365       if(!fMixing) {
1366         if(!kCharg->CheckSoftPi()) { //removal of soft pions
1367           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1368           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1369           continue;
1370         } else {
1371           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1372           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1373         }
1374       Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1375       if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1376       }
1377       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1378
1379       FillSparsePlots(mcArray,d,origD0,PDGD0,kCharg,kKCharg); //fills for charged tracks
1380
1381       if(!fMixing) N_KCharg++;
1382
1383     } // end of charged kaons loop
1384   } //end of event loop for fCorrelatorKc
1385
1386   if(fMixing) {
1387     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1388     if(!execPoolK0) {
1389       AliInfo("Mixed event analysis: K0 pool is not ready");
1390       NofEventsinPool = 0;
1391     }
1392   }
1393
1394   //K0 loop
1395   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)
1396     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1397     if(!analyzetracksK0) {
1398       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1399       continue;
1400     }
1401
1402     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1403
1404       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1405       if(!runcorrelation) continue;
1406
1407       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1408
1409       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1410
1411       FillSparsePlots(mcArray,d,origD0,PDGD0,k0,kK0); //fills for charged tracks
1412
1413       if(!fMixing) N_Kaons++;
1414
1415     } // end of charged kaons loop
1416   } //end of event loop for fCorrelatorK0
1417
1418   //leading track correlations fill
1419   if(!fMixing) {
1420     if(fReadMC) {
1421       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421) { //D0
1422         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0); //c and b D0
1423         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[0],mD0); //c or b D0
1424         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);
1425         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);
1426       }
1427       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421) { //D0bar
1428         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[1],mD0bar);
1429         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(lead[0],mD0bar); //c or b D0
1430         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);
1431         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);
1432       }
1433     } else {
1434         if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[0],mD0); //c and b D0
1435         if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(lead[0],mD0bar);
1436     }
1437
1438     //Fill of count histograms
1439     if (!fAlreadyFilled) {
1440       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1441       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1442       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1443     }
1444   }
1445
1446   fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1447
1448 }
1449
1450 //________________________________________________________________________
1451 void AliAnalysisTaskSED0Correlations::FillSparsePlots(TClonesArray* mcArray, AliAODRecoDecayHF2Prong *d, Int_t origD0, Int_t PdgD0, AliReducedParticle* track, Int_t type) {
1452   //
1453   //fills the THnSparse for correlations, calculating the variables
1454   //
1455
1456   //Initialization of variables
1457   Int_t ptbin = PtBinCorr(d->Pt());
1458   if(ptbin < 0) return;
1459   Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
1460   d->InvMassD0(mD0,mD0bar);
1461
1462   Double_t ptTrack = track->Pt();
1463   Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
1464   Double_t phiTr = track->Phi();
1465   Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1466
1467   TString part = "", orig = "";
1468
1469   switch (type) {
1470     case(kTrack): {
1471       part = "Charg";
1472       deltaphi = fCorrelatorTr->GetDeltaPhi();
1473       deltaeta = fCorrelatorTr->GetDeltaEta();
1474       break;
1475     }
1476     case(kKCharg): {
1477       part = "Kcharg";
1478       deltaphi = fCorrelatorKc->GetDeltaPhi();
1479       deltaeta = fCorrelatorKc->GetDeltaEta();
1480       break;
1481     }
1482     case(kK0): {
1483       part = "K0";
1484       deltaphi = fCorrelatorK0->GetDeltaPhi();
1485       deltaeta = fCorrelatorK0->GetDeltaEta();
1486       break;
1487     }
1488   }
1489
1490   if(fMixing == kSE) {
1491
1492     //Fixes limits; needed to include overflow into THnSparse projections!
1493     Double_t pTorig = track->Pt();
1494     Double_t d0orig = track->GetImpPar();
1495     Double_t ptLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
1496     Double_t displLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
1497     Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
1498     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
1499     if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
1500     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1501     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1502
1503     //variables for filling histos
1504     Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
1505     Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
1506
1507     if(fReadMC == 0) {
1508       //sparse fill for data (tracks, K+-, K0) + weighted
1509       if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
1510         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0);
1511         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pTorig);
1512       }
1513       if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
1514         ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1515         ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pTorig);
1516       }
1517       if(!fAlreadyFilled) {
1518         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1519         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1520       }
1521     }
1522
1523     if(fReadMC) {
1524
1525       if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
1526
1527       //sparse fill for data (tracks, K+-, K0) + weighted
1528       if(PdgD0==421) { //D0 (from MCTruth)
1529          ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0);
1530          ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1531          if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1532          if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0);
1533          ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0,pTorig);
1534          ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
1535          if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
1536          if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0,pTorig);
1537       }
1538       if(PdgD0==-421) { //D0bar
1539          ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1540          ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1541          if(origD0==4&&origTr>=1&&origTr<=2) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1542          if(origD0==5&&origTr>=3&&origTr<=6) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1543          ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(deltaphi,mD0bar,pTorig);
1544          ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
1545          if(origD0==4&&origTr>=1&&origTr<=2) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
1546          if(origD0==5&&origTr>=3&&origTr<=6) ((TH2F*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(deltaphi,mD0bar,pTorig);
1547       }
1548       if(!fAlreadyFilled) {
1549         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1550         if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
1551         if (origTr>=1&&origTr<=6) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
1552         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
1553         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
1554         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
1555         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
1556       }
1557     }//end MC case
1558
1559   } //end of SE fill
1560
1561   if(fMixing == kME) {
1562
1563     //Fixes limits; needed to include overflow into THnSparse projections!
1564     Double_t EtaLim_Sparse = ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
1565     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
1566     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
1567
1568     //variables for filling histos
1569     Double_t fillSpPhiD0[3] = {deltaphi,mD0,deltaeta};
1570     Double_t fillSpPhiD0bar[3] = {deltaphi,mD0bar,deltaeta};
1571
1572     if(fReadMC == 0) {
1573       //sparse fill for data (tracks, K+-, K0)
1574       if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0);
1575       if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1576     }
1577     if(fReadMC == 1) {
1578       //sparse fill for data (tracks, K+-, K0)
1579       if(PdgD0==421)  ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0);
1580       if(PdgD0==-421) ((THnSparseI*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar);
1581     }//end MC case
1582
1583   } //end of ME fill
1584
1585   return;
1586 }
1587
1588 //_________________________________________________________________________________________________
1589 Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1590   //
1591   // checks on particle (#) origin:
1592   // 0) Not HF
1593   // 1) D->#
1594   // 2) D->X->#
1595   // 3) B->#
1596   // 4) B->X-># (X!=D)
1597   // 5) B->D->#
1598   // 6) B->D->X->#
1599   // 7) c hadronization
1600   // 8) b hadronization
1601   //
1602   if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
1603
1604   Int_t pdgGranma = 0;
1605   Int_t mother = 0;
1606   mother = mcPartCandidate->GetMother();
1607   Int_t istep = 0;
1608   Int_t abspdgGranma =0;
1609   Bool_t isFromB=kFALSE;
1610   Bool_t isDdaugh=kFALSE;
1611   Bool_t isDchaindaugh=kFALSE;
1612   Bool_t isBdaugh=kFALSE;
1613   Bool_t isBchaindaugh=kFALSE;
1614   Bool_t isQuarkFound=kFALSE;
1615
1616   while (mother > 0){
1617     istep++;
1618     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1619     if (mcGranma){
1620       pdgGranma = mcGranma->GetPdgCode();
1621       abspdgGranma = TMath::Abs(pdgGranma);
1622       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1623         isBchaindaugh=kTRUE;
1624         if(istep==1) isBdaugh=kTRUE;
1625       }
1626       if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
1627         isDchaindaugh=kTRUE;
1628         if(istep==1) isDdaugh=kTRUE;
1629       }
1630       if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
1631       mother = mcGranma->GetMother();
1632     }else{
1633       AliError("Failed casting the mother particle!");
1634       break;
1635     }
1636   }
1637
1638   //decides what to return based on the flag status
1639   if(isQuarkFound) {
1640     if(!isFromB) {  //charm
1641       if(isDdaugh) return 1; //charm immediate
1642       else if(isDchaindaugh) return 2; //charm chain
1643       else return 7; //charm hadronization
1644     }
1645     else { //beauty
1646       if(isBdaugh) return 3; //b immediate
1647       else if(isBchaindaugh) { //b chain
1648         if(isDchaindaugh) {
1649           if(isDdaugh) return 5; //d immediate
1650           return 6; //d chain
1651           }
1652         else return 4; //b, not d
1653       }
1654       else return 8; //b hadronization
1655     }
1656   }
1657   else return 0; //no HF quark
1658 }
1659
1660 //________________________________________________________________________
1661 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
1662   //
1663   //give the pt bin where the pt lies.
1664   //
1665   Int_t ptbin=-1;
1666   if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
1667
1668   Int_t i = 0;
1669   while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
1670
1671   return ptbin;
1672 }
1673
1674 //---------------------------------------------------------------------------
1675 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
1676 {
1677   //
1678   // Selection for K0 hypotheses
1679   // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
1680   //          2 = no previous selections
1681
1682   if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
1683
1684   AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
1685   AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
1686
1687   if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
1688     if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
1689   }
1690
1691   //This part removes double counting for swapped tracks!
1692   Int_t i = 0;  //while loop (until the last-written entry pair of ID!
1693   while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
1694     if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
1695        (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
1696     i++;
1697   }
1698   idArrayV0[i][0]=v0Daug1->GetID();
1699   idArrayV0[i][1]=v0Daug2->GetID();
1700
1701   return kTRUE;
1702 }
1703
1704 //________________________________________________________________________
1705 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
1706
1707   cout << "--------------------------\n";
1708   cout << "PtBins = " << fNPtBinsCorr << "\n";
1709   cout << "PtBin limits--------------\n";
1710   for (int i=0; i<fNPtBinsCorr; i++) {
1711     cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
1712   }
1713   cout << "\n--------------------------\n";
1714   cout << "PtBin tresh. tracks low---\n";
1715   for (int i=0; i<fNPtBinsCorr; i++) {
1716     cout << fPtThreshLow.at(i) << "; ";
1717   }
1718   cout << "PtBin tresh. tracks up----\n";
1719   for (int i=0; i<fNPtBinsCorr; i++) {
1720     cout << fPtThreshUp.at(i) << "; ";
1721   }
1722   cout << "\n--------------------------\n";
1723   cout << "MC Truth = "<<fReadMC<<"\n";
1724   cout << "--------------------------\n";
1725   cout << "Ev Mixing = "<<fMixing<<"\n";
1726   cout << "--------------------------\n";
1727 }
1728