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