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