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->GetNumberOfTracks();
 611     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
 612       //    ... get the track
 613       AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
 614       if(!track) AliFatal("Not a standard AOD");
 615       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 616         skipEvent=kTRUE;
 617         fNentries->Fill(16);
 618         break;
 619       }
 620     }
 621     if (skipEvent) return;
 622   }
 623
 624   //HFCorrelators initialization (for this event)
 625   fCorrelatorTr->SetAODEvent(aod); // set the AOD event from which you are processing
 626   fCorrelatorKc->SetAODEvent(aod);
 627   fCorrelatorK0->SetAODEvent(aod);
 628   Bool_t correlatorONTr = fCorrelatorTr->Initialize(); // initialize the pool for event mixing
 629   Bool_t correlatorONKc = fCorrelatorKc->Initialize();
 630   Bool_t correlatorONK0 = fCorrelatorK0->Initialize();
 631   if(!correlatorONTr) {AliInfo("AliHFCorrelator (tracks) didn't initialize the pool correctly or processed a bad event"); return;}
 632   if(!correlatorONKc) {AliInfo("AliHFCorrelator (charged K) didn't initialize the pool correctly or processed a bad event"); return;}
 633   if(!correlatorONK0) {AliInfo("AliHFCorrelator (K0) didn't initialize the pool correctly or processed a bad event"); return;}
 634
 635   if(fReadMC) {
 636     fCorrelatorTr->SetMCArray(mcArray); // set the TClonesArray *fmcArray for analysis on monte carlo
 637     fCorrelatorKc->SetMCArray(mcArray);
 638     fCorrelatorK0->SetMCArray(mcArray);
 639   }
 640
 641   // AOD primary vertex
 642   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 643
 644   //vtx1->Print();
 645   TString primTitle = vtx1->GetTitle();
 646   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 647     fNentries->Fill(3);
 648   }
 649
 650   //Reset flag for tracks distributions fill
 651   fAlreadyFilled=kFALSE;
 652
 653   //***** Loop over D0 candidates *****
 654   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 655   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 656
 657   if(fFillGlobal) { //loop on V0 and tracks for each event, to fill Pt distr. and InvMass distr.
 658
 659     TClonesArray *v0array = (TClonesArray*)aod->GetList()->FindObject("v0s");
 660     Int_t pdgCodes[2] = {211,211};
 661     Int_t idArrayV0[v0array->GetEntriesFast()][2];
 662     for(int iV0=0; iV0<v0array->GetEntriesFast(); iV0++) {
 663       for(int j=0; j<2; j++) {idArrayV0[iV0][j]=-2;}
 664       AliAODv0 *v0 = (AliAODv0*)v0array->UncheckedAt(iV0);
 665       if(SelectV0(v0,vtx1,2,idArrayV0)) { //option 2 = for mass inv plots only
 666         if(fReadMC && fRecoTr && (v0->MatchToMC(310,mcArray,2,pdgCodes)<0)) continue; //310 = K0s, 311 = K0 generico!!
 667         ((TH2F*)fOutputStudy->FindObject("hK0MassInv"))->Fill(v0->MassK0Short(),v0->Pt()); //invariant mass plot
 668         ((TH1F*)fOutputStudy->FindObject("hist_Pt_K0_AllEv"))->Fill(v0->Pt()); //pT distribution (in all events), K0 case
 669       }
 670     }
 671
 672     for(Int_t itrack=0; itrack<aod->GetNumberOfTracks(); itrack++) { // loop on tacks
 673       AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
 674       if(!track) AliFatal("Not a standard AOD");
 675       //rejection of tracks
 676       if(track->GetID() < 0) continue; //discard negative ID tracks
 677       if(track->Pt() < fPtThreshLow.at(0) || track->Pt() > fPtThreshUp.at(0)) continue; //discard tracks outside pt range for hadrons/K
 678       if(!fCutsTracks->IsHadronSelected(track) || !fCutsTracks->CheckHadronKinematic(track->Pt(),0.1)) continue; //0.1 = dummy (d0 value, no cut on it for me)
 679       //pT distribution (in all events), charg and hadr cases
 680       ((TH1F*)fOutputStudy->FindObject("hist_Pt_Charg_AllEv"))->Fill(track->Pt());
 681       if(fCutsTracks->CheckKaonCompatibility(track,kFALSE,0,2)) ((TH1F*)fOutputStudy->FindObject("hist_Pt_Kcharg_AllEv"))->Fill(track->Pt());
 682     }
 683
 684   } //end of loops for global plot fill
 685
 686   Int_t nSelectedloose=0,nSelectedtight=0;
 687
 688   //Fill Event Multiplicity (needed only in Reco)
 689   fMultEv = (Double_t)(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.));
 690
 691   //RecoD0 case ************************************************
 692   if(fRecoD0) {
 693
 694     for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 695       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 696
 697       if(d->Pt()<2.) continue; //to save time and merging memory...
 698
 699       if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 700         fNentries->Fill(2);
 701         continue; //skip the D0 from Dstar
 702       }
 703
 704       if (fCutsD0->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
 705         nSelectedloose++;
 706         nSelectedtight++;
 707         if(fSys==0){
 708           if(fCutsD0->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
 709         }
 710         Int_t ptbin=fCutsD0->PtBin(d->Pt());
 711         if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 712
 713         fIsSelectedCandidate=fCutsD0->IsSelected(d,AliRDHFCuts::kAll,aod); //D0 selected
 714         if(!fIsSelectedCandidate) continue;
 715
 716         //D0 infos
 717         Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(d->Phi());
 718                  phiD0 = fCorrelatorKc->SetCorrectPhiRange(d->Phi());  //bad usage, but returns a Double_t...
 719                  phiD0 = fCorrelatorK0->SetCorrectPhiRange(d->Phi());
 720         fCorrelatorTr->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta()); // sets the parameters of the trigger particles that are needed
 721         fCorrelatorKc->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 722         fCorrelatorK0->SetTriggerParticleProperties(d->Pt(),phiD0,d->Eta());
 723         fCorrelatorTr->SetD0Properties(d,fIsSelectedCandidate); //sets special properties for D0
 724         fCorrelatorKc->SetD0Properties(d,fIsSelectedCandidate);
 725         fCorrelatorK0->SetD0Properties(d,fIsSelectedCandidate);
 726
 727         if(!fReadMC) {
 728           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
 729             CalculateCorrelations(d); //correlations on real data
 730             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv);
 731           }
 732         } else { //correlations on MC -> association of selected D0 to MCinfo with MCtruth
 733           if (TMath::Abs(d->Eta())<fEtaForCorrel) {
 734             Int_t pdgDgD0toKpi[2]={321,211};
 735             Int_t labD0 = d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not
 736             if (labD0>-1) {
 737               CalculateCorrelations(d,labD0,mcArray);
 738               if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
 739             }
 740           }
 741         }
 742
 743         FillMassHists(d,mcArray,fCutsD0,fOutputMass);
 744       }
 745     }
 746   }
 747   //End RecoD0 case ************************************************
 748
 749   //MCKineD0 case ************************************************
 750   if(fReadMC && !fRecoD0) {
 751
 752     for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { //Loop over all the tracks of MCArray
 753       AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 754       if (!mcPart) {
 755         AliWarning("Particle not found in tree, skipping");
 756         continue;
 757       }
 758
 759       if(TMath::Abs(mcPart->GetPdgCode()) == 421){  // THIS IS A D0
 760         if (fCutsD0->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y()) ) {
 761           nSelectedloose++;
 762           nSelectedtight++;
 763
 764           //Removal of cases in which D0 decay is not in Kpi!
 765           if(mcPart->GetNDaughters()!=2) continue;
 766           AliAODMCParticle* mcDau1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
 767           AliAODMCParticle* mcDau2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)));
 768           if(!mcDau1 || !mcDau2) continue;
 769           Int_t pdg1 = TMath::Abs(mcDau1->GetPdgCode());
 770           Int_t pdg2 = TMath::Abs(mcDau2->GetPdgCode());
 771           if(!((pdg1 == 211 && pdg2 == 321) || (pdg2 == 211 && pdg1 == 321))) continue;
 772           if(TMath::Abs(mcDau1->Eta())>0.8||TMath::Abs(mcDau2->Eta())>0.8) continue;
 773             //Check momentum conservation (to exclude 4-prong decays with tracks outside y=1.5)
 774             Double_t p1[3]  = {mcDau1->Px(),mcDau1->Py(),mcDau1->Pz()};
 775             Double_t p2[3]  = {mcDau2->Px(),mcDau2->Py(),mcDau2->Pz()};
 776             Double_t pD0[3] = {mcPart->Px(),mcPart->Py(),mcPart->Pz()};
 777             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;
 778
 779           if(fSys==0) fNentries->Fill(6);
 780           Int_t ptbin=fCutsD0->PtBin(mcPart->Pt());
 781           if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 782
 783           //D0 infos
 784           Double_t phiD0 = fCorrelatorTr->SetCorrectPhiRange(mcPart->Phi());
 785                    phiD0 = fCorrelatorKc->SetCorrectPhiRange(mcPart->Phi());  //bad usage, but returns a Double_t...
 786                    phiD0 = fCorrelatorK0->SetCorrectPhiRange(mcPart->Phi());
 787           fCorrelatorTr->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta()); // sets the parameters of the trigger particles that are needed
 788           fCorrelatorKc->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 789           fCorrelatorK0->SetTriggerParticleProperties(mcPart->Pt(),phiD0,mcPart->Eta());
 790           //fCorrelatorTr->SetD0Properties(mcPart,fIsSelectedCandidate); //needed for D* soft pions rejection, useless in MCKine
 791           //fCorrelatorKc->SetD0Properties(mcPart,fIsSelectedCandidate);
 792           //fCorrelatorK0->SetD0Properties(mcPart,fIsSelectedCandidate);
 793
 794           if (TMath::Abs(mcPart->Eta())<fEtaForCorrel) {
 795
 796             //Removal of D0 from D* feeddown! This solves also the problem of soft pions, now excluded
 797          /*   Int_t mother = mcPart->GetMother();
 798             AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
 799             if(!mcMoth) continue;
 800             if(TMath::Abs(mcMoth->GetPdgCode())==413) continue;*/
 801
 802             if (mcPart->GetPdgCode()==421) fIsSelectedCandidate = 1;
 803             else fIsSelectedCandidate = 2;
 804
 805             TString fillthis="histSgn_";
 806             if(CheckD0Origin(mcArray,mcPart)==4) fillthis+="c_";
 807             else if(CheckD0Origin(mcArray,mcPart)==5) fillthis+="b_";
 808             else continue;
 809             fillthis+=ptbin;
 810             ((TH1F*)(fOutputMass->FindObject(fillthis)))->Fill(1.864);
 811
 812             CalculateCorrelationsMCKine(mcPart,mcArray);
 813             if(!fMixing) ((TH1F*)fOutputStudy->FindObject(Form("hMultiplEvt_Bin%d",ptbin)))->Fill(fMultEv); //Fill multiplicity histo
 814           }
 815         }
 816       }
 817     }
 818   }
 819   //End MCKineD0 case ************************************************
 820
 821   if(fMixing /* && fAlreadyFilled*/) { // update the pool for Event Mixing, if: enabled,  event is ok, at least a SelD0 found! (fAlreadyFilled's role!)
 822     Bool_t updatedTr = fCorrelatorTr->PoolUpdate();
 823     Bool_t updatedKc = fCorrelatorKc->PoolUpdate();
 824     Bool_t updatedK0 = fCorrelatorK0->PoolUpdate();
 825     if(!updatedTr || !updatedKc || !updatedK0) AliInfo("Pool was not updated");
 826   }
 827
 828   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
 829   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
 830
 831   // Post the data
 832   PostData(1,fOutputMass);
 833   PostData(2,fNentries);
 834   PostData(4,fCounter);
 835   PostData(5,fOutputCorr);
 836   PostData(6,fOutputStudy);
 837
 838   return;
 839 }
 840
 841 //____________________________________________________________________________
 842 void AliAnalysisTaskSED0Correlations::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
 843   //
 844   // function used in UserExec to fill mass histograms:
 845   //
 846   if (!fIsSelectedCandidate) return;
 847
 848   if(fDebug>2)  cout<<"Candidate selected"<<endl;
 849
 850   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 851   Int_t ptbin = cuts->PtBin(part->Pt());
 852
 853   TString fillthis="";
 854   Int_t pdgDgD0toKpi[2]={321,211};
 855   Int_t labD0=-1;
 856   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)
 857
 858   //count candidates selected by cuts
 859   fNentries->Fill(1);
 860   //count true D0 selected by cuts
 861   if (fReadMC && labD0>=0) fNentries->Fill(2);
 862
 863   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
 864
 865     if(fReadMC){
 866       if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
 867         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 868         Int_t pdgD0 = partD0->GetPdgCode();
 869         if (pdgD0==421){ //D0
 870           fillthis="histSgn_c_";
 871           fillthis+=ptbin;
 872           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 873           fillthis="histSgn_WeigD0Eff_c_";
 874           fillthis+=ptbin;
 875           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 876           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 877         } else{ //it was a D0bar
 878           fillthis="histRfl_c_";
 879           fillthis+=ptbin;
 880           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 881           fillthis="histRfl_WeigD0Eff_c_";
 882           fillthis+=ptbin;
 883           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 884           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 885         }
 886       } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
 887           AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 888           Int_t pdgD0 = partD0->GetPdgCode();
 889           if (pdgD0==421){ //D0
 890             fillthis="histSgn_b_";
 891             fillthis+=ptbin;
 892             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 893             fillthis="histSgn_WeigD0Eff_b_";
 894             fillthis+=ptbin;
 895             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 896             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 897           } else{ //it was a D0bar
 898             fillthis="histRfl_b_";
 899             fillthis+=ptbin;
 900             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 901             fillthis="histRfl_WeigD0Eff_b_";
 902             fillthis+=ptbin;
 903             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 904             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 905           }
 906       } else {//background
 907         fillthis="histBkg_c_";
 908         fillthis+=ptbin;
 909         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 910         fillthis="histBkg_WeigD0Eff_c_";
 911         fillthis+=ptbin;
 912         Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 913         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 914       }
 915     }else{
 916       fillthis="histMass_";
 917       fillthis+=ptbin;
 918       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 919       fillthis="histMass_WeigD0Eff_";
 920       fillthis+=ptbin;
 921       Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 922       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,1./effD0);
 923     }
 924
 925   }
 926   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
 927
 928     if(fReadMC){
 929       if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==4) {
 930         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 931         Int_t pdgD0 = partD0->GetPdgCode();
 932         if (pdgD0==-421){ //D0
 933           fillthis="histSgn_c_";
 934           fillthis+=ptbin;
 935           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 936           fillthis="histSgn_WeigD0Eff_c_";
 937           fillthis+=ptbin;
 938           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 939           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 940         } else{ //it was a D0bar
 941           fillthis="histRfl_c_";
 942           fillthis+=ptbin;
 943           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 944           fillthis="histRfl_WeigD0Eff_c_";
 945           fillthis+=ptbin;
 946           Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 947           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 948         }
 949       } else if(labD0>=0 && CheckD0Origin(arrMC,(AliAODMCParticle*)arrMC->At(labD0))==5) {
 950           AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 951           Int_t pdgD0 = partD0->GetPdgCode();
 952           if (pdgD0==-421){ //D0
 953             fillthis="histSgn_b_";
 954             fillthis+=ptbin;
 955             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 956             fillthis="histSgn_WeigD0Eff_b_";
 957             fillthis+=ptbin;
 958             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 959             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 960           } else{ //it was a D0bar
 961             fillthis="histRfl_b_";
 962             fillthis+=ptbin;
 963             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 964             fillthis="histRfl_WeigD0Eff_b_";
 965             fillthis+=ptbin;
 966             Double_t effD0 = fCutsTracks->GetTrigWeightB(part->Pt(),fMultEv);
 967             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 968           }
 969       } else {//background
 970         fillthis="histBkg_c_";
 971         fillthis+=ptbin;
 972         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 973         fillthis="histBkg_WeigD0Eff_c_";
 974         fillthis+=ptbin;
 975         Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 976         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 977       }
 978     }else{
 979       fillthis="histMass_";
 980       fillthis+=ptbin;
 981       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 982       fillthis="histMass_WeigD0Eff_";
 983       fillthis+=ptbin;
 984       Double_t effD0 = fCutsTracks->GetTrigWeight(part->Pt(),fMultEv);
 985       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,1./effD0);
 986     }
 987
 988   }
 989
 990   return;
 991 }
 992
 993 //________________________________________________________________________
 994 void AliAnalysisTaskSED0Correlations::Terminate(Option_t */*option*/)
 995 {
 996   // Terminate analysis
 997   //
 998   if(fDebug > 1) printf("AnalysisTaskSED0Correlations: Terminate() \n");
 999
1000   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1001   if (!fOutputMass) {
1002     printf("ERROR: fOutputMass not available\n");
1003     return;
1004   }
1005
1006   fNentries = dynamic_cast<TH1F*>(GetOutputData(2));
1007
1008   if(!fNentries){
1009     printf("ERROR: fNEntries not available\n");
1010     return;
1011   }
1012
1013   fCutsD0 = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(3));
1014   if(!fCutsD0){
1015     printf("ERROR: fCuts not available\n");
1016     return;
1017   }
1018
1019   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(4));
1020   if (!fCounter) {
1021     printf("ERROR: fCounter not available\n");
1022     return;
1023   }
1024   fOutputCorr = dynamic_cast<TList*> (GetOutputData(5));
1025   if (!fOutputCorr) {
1026     printf("ERROR: fOutputCorr not available\n");
1027     return;
1028   }
1029   fOutputStudy = dynamic_cast<TList*> (GetOutputData(6));
1030   if (!fOutputStudy) {
1031     printf("ERROR: fOutputStudy not available\n");
1032     return;
1033   }
1034   fCutsTracks = dynamic_cast<AliHFAssociatedTrackCuts*>(GetOutputData(7));
1035   if(!fCutsTracks){
1036     printf("ERROR: fCutsTracks not available\n");
1037     return;
1038   }
1039
1040   return;
1041 }
1042
1043 //_________________________________________________________________________________________________
1044 Int_t AliAnalysisTaskSED0Correlations::CheckD0Origin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1045   //
1046   // checking whether the mother of the particles come from a charm or a bottom quark
1047   //
1048   printf("AliAnalysisTaskSED0Correlations::CheckD0Origin() \n");
1049
1050   Int_t pdgGranma = 0;
1051   Int_t mother = 0;
1052   mother = mcPartCandidate->GetMother();
1053   Int_t abspdgGranma =0;
1054   Bool_t isFromB=kFALSE;
1055   Bool_t isQuarkFound=kFALSE;
1056
1057   while (mother > 0){
1058     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1059     if (mcMoth){
1060       pdgGranma = mcMoth->GetPdgCode();
1061       abspdgGranma = TMath::Abs(pdgGranma);
1062       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1063         isFromB=kTRUE;
1064       }
1065       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1066       mother = mcMoth->GetMother();
1067     }else{
1068       AliError("Failed casting the mother particle!");
1069       break;
1070     }
1071   }
1072
1073   if(isQuarkFound) {
1074     if(isFromB) return 5;
1075     else return 4;
1076   }
1077   else return 1;
1078 }
1079
1080 //________________________________________________________________________
1081 void AliAnalysisTaskSED0Correlations::CreateCorrelationsObjs() {
1082 //
1083
1084   TString namePlot = "";
1085
1086   //These for limits in THnSparse (one per bin, same limits).
1087   //Vars: DeltaPhi, InvMass, PtTrack, Displacement, DeltaEta
1088   Int_t nBinsPhi[5] = {32,150,6,3,16};
1089   Double_t binMinPhi[5] = {-TMath::Pi()/2.,1.6,0.,0.,-1.6};  //is the minimum for all the bins
1090   Double_t binMaxPhi[5] = {3.*TMath::Pi()/2.,2.2,3.0,3.,1.6};  //is the maximum for all the bins
1091
1092   //Vars: DeltaPhi, InvMass, DeltaEta
1093   Int_t nBinsMix[4] = {32,150,16,6};
1094   Double_t binMinMix[4] = {-TMath::Pi()/2.,1.6,-1.6,0.};  //is the minimum for all the bins
1095   Double_t binMaxMix[4] = {3.*TMath::Pi()/2.,2.2,1.6,3.};  //is the maximum for all the bins
1096
1097   for(Int_t i=0;i<fNPtBinsCorr;i++){
1098
1099     if(!fMixing) {
1100       //THnSparse plots: correlations for various invariant mass (MC and data)
1101       namePlot="hPhi_K0_Bin";
1102       namePlot+=i;
1103
1104       THnSparseF *hPhiK = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1105       hPhiK->Sumw2();
1106       fOutputCorr->Add(hPhiK);
1107
1108       namePlot="hPhi_Kcharg_Bin";
1109       namePlot+=i;
1110
1111       THnSparseF *hPhiH = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1112       hPhiH->Sumw2();
1113       fOutputCorr->Add(hPhiH);
1114
1115       namePlot="hPhi_Charg_Bin";
1116       namePlot+=i;
1117
1118       THnSparseF *hPhiC = new THnSparseF(namePlot.Data(), "Azimuthal correlation; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",5,nBinsPhi,binMinPhi,binMaxPhi);
1119       hPhiC->Sumw2();
1120       fOutputCorr->Add(hPhiC);
1121
1122       //histos for c/b origin for D0 (MC only)
1123       if (fReadMC) {
1124
1125         //generic origin for tracks
1126         namePlot="hPhi_K0_From_c_Bin";
1127         namePlot+=i;
1128
1129         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);
1130         hPhiK_c->Sumw2();
1131         fOutputCorr->Add(hPhiK_c);
1132
1133         namePlot="hPhi_Kcharg_From_c_Bin";
1134         namePlot+=i;
1135
1136         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);
1137         hPhiH_c->Sumw2();
1138         fOutputCorr->Add(hPhiH_c);
1139
1140         namePlot="hPhi_Charg_From_c_Bin";
1141         namePlot+=i;
1142
1143         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);
1144         hPhiC_c->Sumw2();
1145         fOutputCorr->Add(hPhiC_c);
1146
1147         namePlot="hPhi_K0_From_b_Bin";
1148         namePlot+=i;
1149
1150         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);
1151         hPhiK_b->Sumw2();
1152         fOutputCorr->Add(hPhiK_b);
1153
1154         namePlot="hPhi_Kcharg_From_b_Bin";
1155         namePlot+=i;
1156
1157         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);
1158         hPhiH_b->Sumw2();
1159         fOutputCorr->Add(hPhiH_b);
1160
1161         namePlot="hPhi_Charg_From_b_Bin";
1162         namePlot+=i;
1163
1164         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);
1165         hPhiC_b->Sumw2();
1166         fOutputCorr->Add(hPhiC_b);
1167
1168         //HF-only tracks (c for c->D0, b for b->D0)
1169         namePlot="hPhi_K0_HF_From_c_Bin";
1170         namePlot+=i;
1171
1172         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);
1173         hPhiK_HF_c->Sumw2();
1174         fOutputCorr->Add(hPhiK_HF_c);
1175
1176         namePlot="hPhi_Kcharg_HF_From_c_Bin";
1177         namePlot+=i;
1178
1179         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);
1180         hPhiH_HF_c->Sumw2();
1181         fOutputCorr->Add(hPhiH_HF_c);
1182
1183         namePlot="hPhi_Charg_HF_From_c_Bin";
1184         namePlot+=i;
1185
1186         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);
1187         hPhiC_HF_c->Sumw2();
1188         fOutputCorr->Add(hPhiC_HF_c);
1189
1190         namePlot="hPhi_K0_HF_From_b_Bin";
1191         namePlot+=i;
1192
1193         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);
1194         hPhiK_HF_b->Sumw2();
1195         fOutputCorr->Add(hPhiK_HF_b);
1196
1197         namePlot="hPhi_Kcharg_HF_From_b_Bin";
1198         namePlot+=i;
1199
1200         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);
1201         hPhiH_HF_b->Sumw2();
1202         fOutputCorr->Add(hPhiH_HF_b);
1203
1204         namePlot="hPhi_Charg_HF_From_b_Bin";
1205         namePlot+=i;
1206
1207         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);
1208         hPhiC_HF_b->Sumw2();
1209         fOutputCorr->Add(hPhiC_HF_b);
1210
1211         namePlot="hPhi_K0_NonHF_Bin";
1212         namePlot+=i;
1213
1214         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);
1215         hPhiK_Non->Sumw2();
1216         fOutputCorr->Add(hPhiK_Non);
1217
1218         namePlot="hPhi_Kcharg_NonHF_Bin";
1219         namePlot+=i;
1220
1221         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);
1222         hPhiH_Non->Sumw2();
1223         fOutputCorr->Add(hPhiH_Non);
1224
1225         namePlot="hPhi_Charg_NonHF_Bin";
1226         namePlot+=i;
1227
1228         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);
1229         hPhiC_Non->Sumw2();
1230         fOutputCorr->Add(hPhiC_Non);
1231       }
1232
1233       //leading hadron correlations
1234       namePlot="hPhi_Lead_Bin";
1235       namePlot+=i;
1236
1237       THnSparseF *hCorrLead = new THnSparseF(namePlot.Data(), "Leading particle correlations; #Delta#phi; Inv. Mass (GeV/c^{2}); p_{t} (GeV/c)",4,nBinsMix,binMinMix,binMaxMix);
1238       hCorrLead->Sumw2();
1239       fOutputCorr->Add(hCorrLead);
1240
1241       if (fReadMC) {
1242         namePlot="hPhi_Lead_From_c_Bin";
1243         namePlot+=i;
1244
1245         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);
1246         hCorrLead_c->Sumw2();
1247         fOutputCorr->Add(hCorrLead_c);
1248
1249         namePlot="hPhi_Lead_From_b_Bin";
1250         namePlot+=i;
1251
1252         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);
1253         hCorrLead_b->Sumw2();
1254         fOutputCorr->Add(hCorrLead_b);
1255
1256         namePlot="hPhi_Lead_HF_From_c_Bin";
1257         namePlot+=i;
1258
1259         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);
1260         hCorrLead_HF_c->Sumw2();
1261         fOutputCorr->Add(hCorrLead_HF_c);
1262
1263         namePlot="hPhi_Lead_HF_From_b_Bin";
1264         namePlot+=i;
1265
1266         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);
1267         hCorrLead_HF_b->Sumw2();
1268         fOutputCorr->Add(hCorrLead_HF_b);
1269
1270         namePlot="hPhi_Lead_NonHF_Bin";
1271         namePlot+=i;
1272
1273         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);
1274         hCorrLead_Non->Sumw2();
1275         fOutputCorr->Add(hCorrLead_Non);
1276       }
1277
1278       //pT weighted correlations
1279       namePlot="hPhi_Weig_Bin";
1280       namePlot+=i;
1281
1282       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);
1283       fOutputCorr->Add(hCorrWeig);
1284
1285       if (fReadMC) {
1286         namePlot="hPhi_Weig_From_c_Bin";
1287         namePlot+=i;
1288
1289         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);
1290         fOutputCorr->Add(hCorrWeig_c);
1291
1292         namePlot="hPhi_Weig_From_b_Bin";
1293         namePlot+=i;
1294
1295         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);
1296         fOutputCorr->Add(hCorrWeig_b);
1297
1298         namePlot="hPhi_Weig_HF_From_c_Bin";
1299         namePlot+=i;
1300
1301         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);
1302         fOutputCorr->Add(hCorrWeig_HF_c);
1303
1304         namePlot="hPhi_Weig_HF_From_b_Bin";
1305         namePlot+=i;
1306
1307         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);
1308         fOutputCorr->Add(hCorrWeig_HF_b);
1309
1310         namePlot="hPhi_Weig_NonHF_Bin";
1311         namePlot+=i;
1312
1313         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);
1314         fOutputCorr->Add(hCorrWeig_Non);
1315       }
1316
1317     //pT distribution histos
1318     namePlot = "hist_Pt_Charg_Bin"; namePlot+=i;
1319     TH1F *hPtC = new TH1F(namePlot.Data(), "Charged track pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1320     hPtC->SetMinimum(0);
1321     fOutputStudy->Add(hPtC);
1322
1323     namePlot = "hist_Pt_Kcharg_Bin"; namePlot+=i;
1324     TH1F *hPtH = new TH1F(namePlot.Data(), "Hadrons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1325     hPtH->SetMinimum(0);
1326     fOutputStudy->Add(hPtH);
1327
1328     namePlot = "hist_Pt_K0_Bin"; namePlot+=i;
1329     TH1F *hPtK = new TH1F(namePlot.Data(), "Kaons pT (in D0 evs); p_{T} (GeV/c)",240,0.,12.);
1330     hPtK->SetMinimum(0);
1331     fOutputStudy->Add(hPtK);
1332
1333     //D* feeddown pions rejection histos
1334     namePlot = "hDstarPions_Bin"; namePlot+=i;
1335     TH2F *hDstarPions = new TH2F(namePlot.Data(), "Tracks rejected for D* inv.mass cut; # Tracks",2,0.,2.,150,1.6,2.2);
1336     hDstarPions->GetXaxis()->SetBinLabel(1,"Not rejected");
1337     hDstarPions->GetXaxis()->SetBinLabel(2,"Rejected");
1338     hDstarPions->SetMinimum(0);
1339     fOutputStudy->Add(hDstarPions);
1340
1341     //Events multiplicity
1342     Double_t yAxisMult[13] = {0, 4, 8, 12, 16, 20, 28, 36, 44, 100};
1343     namePlot = "hMultiplEvt_Bin"; namePlot+=i;
1344     TH1F *hMultEv = new TH1F(namePlot.Data(), "Event multiplicity",9,yAxisMult);
1345     hMultEv->SetMinimum(0);
1346     fOutputStudy->Add(hMultEv);
1347
1348     }
1349
1350     if(fMixing) {
1351       //THnSparse plots for event mixing!
1352       namePlot="hPhi_K0_Bin";
1353       namePlot+=i;namePlot+="_EvMix";
1354
1355       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);
1356       hPhiK_EvMix->Sumw2();
1357       fOutputCorr->Add(hPhiK_EvMix);
1358
1359       namePlot="hPhi_Kcharg_Bin";
1360       namePlot+=i;namePlot+="_EvMix";
1361
1362       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);
1363       hPhiH_EvMix->Sumw2();
1364       fOutputCorr->Add(hPhiH_EvMix);
1365
1366       namePlot="hPhi_Charg_Bin";
1367       namePlot+=i;namePlot+="_EvMix";
1368
1369       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);
1370       hPhiC_EvMix->Sumw2();
1371       fOutputCorr->Add(hPhiC_EvMix);
1372     }
1373   }
1374
1375   //out of bin loop
1376   if(!fMixing) {
1377     TH1F *hCountC = new TH1F("hist_Count_Charg", "Charged track counter; # Tracks",100,0.,100.);
1378     hCountC->SetMinimum(0);
1379     fOutputStudy->Add(hCountC);
1380
1381     TH1F *hCountH = new TH1F("hist_Count_Kcharg", "Hadrons counter; # Tracks",20,0.,20.);
1382     hCountH->SetMinimum(0);
1383     fOutputStudy->Add(hCountH);
1384
1385     TH1F *hCountK = new TH1F("hist_Count_K0", "Kaons counter; # Tracks",20,0.,20.);
1386     hCountK->SetMinimum(0);
1387     fOutputStudy->Add(hCountK);
1388   }
1389
1390   if (fReadMC) {
1391     TH1D *hEventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);
1392     fOutputStudy->Add(hEventTypeMC);
1393   }
1394
1395   if (fFillGlobal) { //all-events plots
1396     //pt distributions
1397     TH1F *hPtCAll = new TH1F("hist_Pt_Charg_AllEv", "Charged track pT (All); p_{T} (GeV/c)",240,0.,12.);
1398     hPtCAll->SetMinimum(0);
1399     fOutputStudy->Add(hPtCAll);
1400
1401     TH1F *hPtHAll = new TH1F("hist_Pt_Kcharg_AllEv", "Hadrons pT (All); p_{T} (GeV/c)",240,0.,12.);
1402     hPtHAll->SetMinimum(0);
1403     fOutputStudy->Add(hPtHAll);
1404
1405     TH1F *hPtKAll = new TH1F("hist_Pt_K0_AllEv", "Kaons pT (All); p_{T} (GeV/c)",240,0.,12.);
1406     hPtKAll->SetMinimum(0);
1407     fOutputStudy->Add(hPtKAll);
1408
1409     //K0 Invariant Mass plots
1410     TH2F *hK0MassInv = new TH2F("hK0MassInv", "K0 invariant mass; Invariant mass (MeV/c^{2}); pT (GeV/c)",200,0.4,0.6,100,0.,10.);
1411     hK0MassInv->SetMinimum(0);
1412     fOutputStudy->Add(hK0MassInv);
1413   }
1414
1415   if(!fMixing) {
1416     //phi distributions
1417     TH1F *hPhiDistCAll = new TH1F("hist_PhiDistr_Charg", "Charged track phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1418     hPhiDistCAll->SetMinimum(0);
1419     fOutputStudy->Add(hPhiDistCAll);
1420
1421     TH1F *hPhiDistHAll = new TH1F("hist_PhiDistr_Kcharg", "Hadrons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1422     hPhiDistHAll->SetMinimum(0);
1423     fOutputStudy->Add(hPhiDistHAll);
1424
1425     TH1F *hPhiDistKAll = new TH1F("hist_PhiDistr_K0", "Kaons phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1426     hPhiDistKAll->SetMinimum(0);
1427     fOutputStudy->Add(hPhiDistKAll);
1428
1429     TH1F *hPhiDistDAll = new TH1F("hist_PhiDistr_D0", "D^{0} phi distr. (All); p_{T} (GeV/c)",64,0,6.283);
1430     hPhiDistDAll->SetMinimum(0);
1431     fOutputStudy->Add(hPhiDistDAll);
1432   }
1433
1434   //for MC analysis only
1435   for(Int_t i=0;i<fNPtBinsCorr;i++) {
1436
1437     if (fReadMC && !fMixing) {
1438
1439       //displacement histos
1440       namePlot="histDispl_K0_Bin"; namePlot+=i;
1441       TH1F *hDisplK = new TH1F(namePlot.Data(), "Kaons Displacement; DCA",150,0.,0.15);
1442       hDisplK->SetMinimum(0);
1443       fOutputStudy->Add(hDisplK);
1444
1445       namePlot="histDispl_K0_HF_Bin";  namePlot+=i;
1446       TH1F *hDisplK_HF = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only); DCA",150,0.,0.15);
1447       hDisplK_HF->SetMinimum(0);
1448       fOutputStudy->Add(hDisplK_HF);
1449
1450       namePlot="histDispl_Kcharg_Bin"; namePlot+=i;
1451       TH1F *hDisplHadr = new TH1F(namePlot.Data(), "Hadrons Displacement; DCA",150,0.,0.15);
1452       hDisplHadr->SetMinimum(0);
1453       fOutputStudy->Add(hDisplHadr);
1454
1455       namePlot="histDispl_Kcharg_HF_Bin";  namePlot+=i;
1456       TH1F *hDisplHadr_HF = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only); DCA",150,0.,0.15);
1457       hDisplHadr_HF->SetMinimum(0);
1458       fOutputStudy->Add(hDisplHadr_HF);
1459
1460       namePlot="histDispl_Charg_Bin"; namePlot+=i;
1461       TH1F *hDisplCharg = new TH1F(namePlot.Data(), "Charged tracks Displacement; DCA",150,0.,0.15);
1462       hDisplCharg->SetMinimum(0);
1463       fOutputStudy->Add(hDisplCharg);
1464
1465       namePlot="histDispl_Charg_HF_Bin";  namePlot+=i;
1466       TH1F *hDisplCharg_HF = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only); DCA",150,0.,0.15);
1467       hDisplCharg_HF->SetMinimum(0);
1468       fOutputStudy->Add(hDisplCharg_HF);
1469
1470       namePlot="histDispl_K0_From_c_Bin"; namePlot+=i;
1471       TH1F *hDisplK_c = new TH1F(namePlot.Data(), "Kaons Displacement - c origin; DCA",150,0.,0.15);
1472       hDisplK_c->SetMinimum(0);
1473       fOutputStudy->Add(hDisplK_c);
1474
1475       namePlot="histDispl_K0_HF_From_c_Bin";  namePlot+=i;
1476       TH1F *hDisplK_HF_c = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1477       hDisplK_HF_c->SetMinimum(0);
1478       fOutputStudy->Add(hDisplK_HF_c);
1479
1480       namePlot="histDispl_Kcharg_From_c_Bin"; namePlot+=i;
1481       TH1F *hDisplHadr_c = new TH1F(namePlot.Data(), "Hadrons Displacement - c origin; DCA",150,0.,0.15);
1482       hDisplHadr_c->SetMinimum(0);
1483       fOutputStudy->Add(hDisplHadr_c);
1484
1485       namePlot="histDispl_Kcharg_HF_From_c_Bin";  namePlot+=i;
1486       TH1F *hDisplHadr_HF_c = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1487       hDisplHadr_HF_c->SetMinimum(0);
1488       fOutputStudy->Add(hDisplHadr_HF_c);
1489
1490       namePlot="histDispl_Charg_From_c_Bin"; namePlot+=i;
1491       TH1F *hDisplCharg_c = new TH1F(namePlot.Data(), "Charged tracks Displacement - c origin; DCA",150,0.,0.15);
1492       hDisplCharg_c->Sumw2();
1493       hDisplCharg_c->SetMinimum(0);
1494       fOutputStudy->Add(hDisplCharg_c);
1495
1496       namePlot="histDispl_Charg_HF_From_c_Bin";  namePlot+=i;
1497       TH1F *hDisplCharg_HF_c = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - c origin; DCA",150,0.,0.15);
1498       hDisplCharg_HF_c->SetMinimum(0);
1499       fOutputStudy->Add(hDisplCharg_HF_c);
1500
1501       namePlot="histDispl_K0_From_b_Bin"; namePlot+=i;
1502       TH1F *hDisplK_b = new TH1F(namePlot.Data(), "Kaons Displacement - b origin; DCA",150,0.,0.15);
1503       hDisplK_b->SetMinimum(0);
1504       fOutputStudy->Add(hDisplK_b);
1505
1506       namePlot="histDispl_K0_HF_From_b_Bin";  namePlot+=i;
1507       TH1F *hDisplK_HF_b = new TH1F(namePlot.Data(), "Kaons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1508       hDisplK_HF_b->SetMinimum(0);
1509       fOutputStudy->Add(hDisplK_HF_b);
1510
1511       namePlot="histDispl_Kcharg_From_b_Bin"; namePlot+=i;
1512       TH1F *hDisplHadr_b = new TH1F(namePlot.Data(), "Hadrons Displacement - b origin; DCA",150,0.,0.15);
1513       hDisplHadr_b->SetMinimum(0);
1514       fOutputStudy->Add(hDisplHadr_b);
1515
1516       namePlot="histDispl_Kcharg_HF_From_b_Bin";  namePlot+=i;
1517       TH1F *hDisplHadr_HF_b = new TH1F(namePlot.Data(), "Hadrons Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1518       hDisplHadr_HF_b->SetMinimum(0);
1519       fOutputStudy->Add(hDisplHadr_HF_b);
1520
1521       namePlot="histDispl_Charg_From_b_Bin"; namePlot+=i;
1522       TH1F *hDisplCharg_b = new TH1F(namePlot.Data(), "Charged tracks Displacement - b origin; DCA",150,0.,0.15);
1523       hDisplCharg_b->SetMinimum(0);
1524       fOutputStudy->Add(hDisplCharg_b);
1525
1526       namePlot="histDispl_Charg_HF_From_b_Bin";  namePlot+=i;
1527       TH1F *hDisplCharg_HF_b = new TH1F(namePlot.Data(), "Charged tracks Displacement (from HF decay only) - b origin; DCA",150,0.,0.15);
1528       hDisplCharg_HF_b->SetMinimum(0);
1529       fOutputStudy->Add(hDisplCharg_HF_b);
1530
1531       //origin of tracks histos
1532       namePlot="histOrig_Charg_Bin";  namePlot+=i;
1533       TH1F *hOrigin_Charm = new TH1F(namePlot.Data(), "Origin of charged tracks",9,0.,9.);
1534       hOrigin_Charm->SetMinimum(0);
1535       hOrigin_Charm->GetXaxis()->SetBinLabel(1,"Not HF");
1536       hOrigin_Charm->GetXaxis()->SetBinLabel(2,"D->#");
1537       hOrigin_Charm->GetXaxis()->SetBinLabel(3,"D->X->#");
1538       hOrigin_Charm->GetXaxis()->SetBinLabel(4,"c hadr.");
1539       hOrigin_Charm->GetXaxis()->SetBinLabel(5,"B->#");
1540       hOrigin_Charm->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1541       hOrigin_Charm->GetXaxis()->SetBinLabel(7,"B->D->#");
1542       hOrigin_Charm->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1543       hOrigin_Charm->GetXaxis()->SetBinLabel(9,"b hadr.");
1544       fOutputStudy->Add(hOrigin_Charm);
1545
1546       namePlot="histOrig_Kcharg_Bin";  namePlot+=i;
1547       TH1F *hOrigin_Kcharg = new TH1F(namePlot.Data(), "Origin of hadrons",9,0.,9.);
1548       hOrigin_Kcharg->SetMinimum(0);
1549       hOrigin_Kcharg->GetXaxis()->SetBinLabel(1,"Not HF");
1550       hOrigin_Kcharg->GetXaxis()->SetBinLabel(2,"D->#");
1551       hOrigin_Kcharg->GetXaxis()->SetBinLabel(3,"D->X->#");
1552       hOrigin_Kcharg->GetXaxis()->SetBinLabel(4,"c hadr.");
1553       hOrigin_Kcharg->GetXaxis()->SetBinLabel(5,"B->#");
1554       hOrigin_Kcharg->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1555       hOrigin_Kcharg->GetXaxis()->SetBinLabel(7,"B->D->#");
1556       hOrigin_Kcharg->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1557       hOrigin_Kcharg->GetXaxis()->SetBinLabel(9,"b hadr.");
1558       fOutputStudy->Add(hOrigin_Kcharg);
1559
1560       namePlot="histOrig_K0_Bin";  namePlot+=i;
1561       TH1F *hOrigin_K = new TH1F(namePlot.Data(), "Origin of kaons",9,0.,9.);
1562       hOrigin_K->SetMinimum(0);
1563       hOrigin_K->GetXaxis()->SetBinLabel(1,"Not HF");
1564       hOrigin_K->GetXaxis()->SetBinLabel(2,"D->#");
1565       hOrigin_K->GetXaxis()->SetBinLabel(3,"D->X->#");
1566       hOrigin_K->GetXaxis()->SetBinLabel(4,"c hadr.");
1567       hOrigin_K->GetXaxis()->SetBinLabel(5,"B->#");
1568       hOrigin_K->GetXaxis()->SetBinLabel(6,"B->X-># (X!=D)");
1569       hOrigin_K->GetXaxis()->SetBinLabel(7,"B->D->#");
1570       hOrigin_K->GetXaxis()->SetBinLabel(8,"B->D->X->#");
1571       hOrigin_K->GetXaxis()->SetBinLabel(9,"b hadr.");
1572       fOutputStudy->Add(hOrigin_K);
1573     }
1574
1575     if (fReadMC) {
1576       //origin of D0 histos
1577       namePlot="histOrig_D0_Bin";  namePlot+=i;
1578       TH1F *hOrigin_D0 = new TH1F(namePlot.Data(), "Origin of D0",2,0.,2.);
1579       hOrigin_D0->SetMinimum(0);
1580       hOrigin_D0->GetXaxis()->SetBinLabel(1,"From c");
1581       hOrigin_D0->GetXaxis()->SetBinLabel(2,"From b");
1582       fOutputStudy->Add(hOrigin_D0);
1583
1584       //primary tracks (Kine & Reco)
1585       namePlot="hPhysPrim_Bin";  namePlot+=i;
1586       TH1F *hPhysPrim = new TH1F(namePlot.Data(), "Origin of hadrons",2,0.,2.);
1587       hPhysPrim->SetMinimum(0);
1588       hPhysPrim->GetXaxis()->SetBinLabel(1,"OK");
1589       hPhysPrim->GetXaxis()->SetBinLabel(2,"NO");
1590       fOutputStudy->Add(hPhysPrim);
1591     }
1592   }
1593 }
1594
1595 //________________________________________________________________________
1596 void AliAnalysisTaskSED0Correlations::CalculateCorrelations(AliAODRecoDecayHF2Prong* d, Int_t labD0, TClonesArray* mcArray) {
1597 //
1598 // Method for correlations D0-hadrons study
1599 //
1600   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1601   Double_t mD0, mD0bar;
1602   Int_t origD0 = 0, PDGD0 = 0, ptbin = 0;
1603   d->InvMassD0(mD0,mD0bar);
1604   Double_t mInv[2] = {mD0, mD0bar};
1605   ptbin = PtBinCorr(d->Pt());
1606
1607   if(ptbin < 0) return;
1608
1609   //Fill of D0 phi distribution
1610   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1611
1612   //Origin of D0
1613   TString orig="";
1614   if(fReadMC) {
1615     origD0=CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0));
1616     PDGD0 = ((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode();
1617     switch (CheckD0Origin(mcArray,(AliAODMCParticle*)mcArray->At(labD0))) {
1618       case 4:
1619         orig = "_From_c";
1620         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1621         break;
1622       case 5:
1623         orig = "_From_b";
1624         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1625         break;
1626       default:
1627         return;
1628     }
1629   }
1630
1631   Double_t highPt = 0; Double_t lead[4] = {0,0,0,1};  //infos for leading particle (pt,deltaphi)
1632
1633   //loop over the tracks in the pool
1634   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1635   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1636   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1637
1638   Int_t NofEventsinPool = 1;
1639   if(fMixing) {
1640     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1641     if(!execPoolTr) {
1642       AliInfo("Mixed event analysis: track pool is not ready");
1643       NofEventsinPool = 0;
1644     }
1645   }
1646
1647   //Charged tracks
1648   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)
1649     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1650     if(!analyzetracksTr) {
1651       AliInfo("AliHFCorrelator::Cannot process the track array");
1652       continue;
1653     }
1654
1655     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1656
1657       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1658       if(!runcorrelation) continue;
1659
1660       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1661
1662       if(!fMixing) {
1663         if(fSoftPiCut && !track->CheckSoftPi()) { //removal of soft pions
1664           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1665           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1666           continue;
1667         } else { //not a soft pion
1668           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1669           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1670         }
1671         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1672         if(track->GetID() == idDaughs[0] || track->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1673       }
1674       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1675
1676       if(fReadMC) {
1677         AliAODMCParticle* trkKine = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1678         if (!trkKine) continue;
1679         if (!trkKine->IsPhysicalPrimary()) {
1680           ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);
1681           continue; //reject the Reco track if correspondent Kine track is not primary
1682         } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
1683       }
1684
1685       Double_t effTr = track->GetWeight(); //extract track efficiency
1686       Double_t effD0 = 1.;
1687       if(fReadMC) {
1688         if(origD0==4) effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
1689         if(origD0==5) effD0 = fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv);
1690       } else effD0 = fCutsTracks->GetTrigWeight(d->Pt(),fMultEv);
1691       Double_t eff = effTr*effD0;
1692
1693       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack,1./eff); //fills for charged tracks
1694
1695       if(!fMixing) N_Charg++;
1696
1697       //retrieving leading info...
1698       if(track->Pt() > highPt) {
1699         if(fReadMC && track->GetLabel()<1) continue;
1700         if(fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) continue;
1701         lead[0] = fCorrelatorTr->GetDeltaPhi();
1702         lead[1] = fCorrelatorTr->GetDeltaEta();
1703         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
1704         if(fReadMC) {
1705           if(origD0==4) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv)); //weight is 1./efficiency
1706           if(origD0==5) lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeightB(d->Pt(),fMultEv)); //weight is 1./efficiency
1707         } else lead[3] = 1./(track->GetWeight()*fCutsTracks->GetTrigWeight(d->Pt(),fMultEv));
1708         highPt = track->Pt();
1709       }
1710
1711     } // end of tracks loop
1712   } //end of event loop for fCorrelatorTr
1713
1714  if(fKaonCorr) { //loops for Kcharg and K0
1715
1716   if(fMixing) {
1717     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1718     if(!execPoolKc) {
1719       AliInfo("Mixed event analysis: K+/- pool is not ready");
1720       NofEventsinPool = 0;
1721     }
1722   }
1723
1724   //Charged Kaons loop
1725   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)
1726     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1727     if(!analyzetracksKc) {
1728       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1729       continue;
1730     }
1731
1732     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1733
1734       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1735       if(!runcorrelation) continue;
1736
1737       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1738
1739       if(!fMixing) {
1740         if(fSoftPiCut && !kCharg->CheckSoftPi()) { //removal of soft pions
1741           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0);
1742           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(1.,mD0bar);
1743           continue;
1744         } else {
1745           if (fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0);
1746           if (fIsSelectedCandidate >= 2) ((TH2F*)fOutputStudy->FindObject(Form("hDstarPions_Bin%d",ptbin)))->Fill(0.,mD0bar);
1747         }
1748         Int_t idDaughs[2] = {((AliVTrack*)d->GetDaughter(0))->GetID(),((AliVTrack*)d->GetDaughter(1))->GetID()}; //IDs of daughters to be skipped
1749         if(kCharg->GetID() == idDaughs[0] || kCharg->GetID() == idDaughs[1]) continue; //discards daughters of candidate
1750       }
1751       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1752
1753       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1754
1755       if(!fMixing) N_KCharg++;
1756
1757     } // end of charged kaons loop
1758   } //end of event loop for fCorrelatorKc
1759
1760   if(fMixing) {
1761     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1762     if(!execPoolK0) {
1763       AliInfo("Mixed event analysis: K0 pool is not ready");
1764       NofEventsinPool = 0;
1765     }
1766   }
1767
1768   //K0 loop
1769   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)
1770     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1771     if(!analyzetracksK0) {
1772       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1773       continue;
1774     }
1775
1776     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1777
1778       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1779       if(!runcorrelation) continue;
1780
1781       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1782
1783       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1784
1785       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1786
1787       if(!fMixing) N_Kaons++;
1788
1789     } // end of charged kaons loop
1790   } //end of event loop for fCorrelatorK0
1791
1792  } //end of 'if(fKaonCorr)'
1793
1794   Double_t fillSpLeadD0[4] = {lead[0],mD0,lead[1],0.4}; //dummy value for threshold of leading!
1795   Double_t fillSpLeadD0bar[4] = {lead[0],mD0bar,lead[1],0.4};
1796
1797   //leading track correlations fill
1798   if(!fMixing) {
1799     if(fReadMC) {
1800       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
1801         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]); //c and b D0
1802         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0,lead[3]); //c or b D0
1803         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]);
1804         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]);
1805         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);  //non HF
1806       }
1807       if(((AliAODMCParticle*)mcArray->At(labD0))->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
1808         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1809         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadD0bar,lead[3]); //c or b D0
1810         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]);
1811         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]);
1812         if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);  //non HF
1813       }
1814     } else {
1815         if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0,lead[3]);
1816         if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadD0bar,lead[3]);
1817     }
1818
1819     //Fill of count histograms
1820     if (!fAlreadyFilled) {
1821       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
1822       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
1823       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
1824     }
1825   }
1826
1827   fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
1828
1829 }
1830
1831 //________________________________________________________________________
1832 void AliAnalysisTaskSED0Correlations::CalculateCorrelationsMCKine(AliAODMCParticle* d, TClonesArray* mcArray) {
1833 //
1834 // Method for correlations D0-hadrons study
1835 //
1836   Int_t N_Charg = 0, N_KCharg = 0, N_Kaons = 0;
1837   Double_t mD0 = 1.864, mD0bar = 1.864;
1838   Double_t mInv[2] = {mD0, mD0bar};
1839   Int_t origD0 = 0, PDGD0 = 0;
1840   Int_t ptbin = PtBinCorr(d->Pt());
1841
1842   if(ptbin < 0) return;
1843
1844   //Fill of D0 phi distribution
1845   if (!fMixing) ((TH1F*)fOutputStudy->FindObject("hist_PhiDistr_D0"))->Fill(d->Phi());
1846
1847   //Origin of D0
1848   TString orig="";
1849   origD0=CheckD0Origin(mcArray,d);
1850   PDGD0 = d->GetPdgCode();
1851   switch (CheckD0Origin(mcArray,d)) {
1852     case 4:
1853       orig = "_From_c";
1854       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(0.);
1855       break;
1856     case 5:
1857       orig = "_From_b";
1858       ((TH1F*)fOutputStudy->FindObject(Form("histOrig_D0_Bin%d",ptbin)))->Fill(1.);
1859       break;
1860     default:
1861       return;
1862   }
1863
1864   Double_t highPt = 0; Double_t lead[3] = {0,0,0};  //infos for leading particle (pt,deltaphi)
1865
1866   //loop over the tracks in the pool
1867   Bool_t execPoolTr = fCorrelatorTr->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1868   Bool_t execPoolKc = fCorrelatorKc->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1869   Bool_t execPoolK0 = fCorrelatorK0->ProcessEventPool(); //pool is ready? (only in ME, in SE returns kFALSE)
1870
1871   Int_t NofEventsinPool = 1;
1872   if(fMixing) {
1873     NofEventsinPool = fCorrelatorTr->GetNofEventsInPool();
1874     if(!execPoolTr) {
1875       AliInfo("Mixed event analysis: track pool is not ready");
1876       NofEventsinPool = 0;
1877     }
1878   }
1879
1880   //Charged tracks
1881   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)
1882
1883     Bool_t analyzetracksTr = fCorrelatorTr->ProcessAssociatedTracks(jMix);// process all the tracks in the aodEvent, by applying the selection cuts
1884     if(!analyzetracksTr) {
1885       AliInfo("AliHFCorrelator::Cannot process the track array");
1886       continue;
1887     }
1888
1889     for(Int_t iTrack = 0; iTrack<fCorrelatorTr->GetNofTracks(); iTrack++){ // looping on track candidates
1890
1891       Bool_t runcorrelation = fCorrelatorTr->Correlate(iTrack);
1892       if(!runcorrelation) continue;
1893
1894       AliReducedParticle* track = fCorrelatorTr->GetAssociatedParticle();
1895       if(track->GetLabel()<0) continue;
1896       if(track->Pt() < fPtThreshLow.at(ptbin) || track->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1897       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
1898       if(!fMixing) N_Charg++;
1899
1900       AliAODMCParticle *trkMC = (AliAODMCParticle*)mcArray->At(track->GetLabel());
1901       if(!trkMC) continue;
1902
1903       if (!trkMC->IsPhysicalPrimary()) {  //reject material budget, or other fake tracks
1904         ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(1.);
1905         continue;
1906       } else ((TH1F*)fOutputStudy->FindObject(Form("hPhysPrim_Bin%d",ptbin)))->Fill(0.);
1907
1908       if (IsDDaughter(d,trkMC,mcArray)) continue;
1909       if (fSoftPiCut && IsSoftPion_MCKine(d,trkMC,mcArray)) continue; //remove soft pions (if requestes, e.g. for templates)
1910
1911       FillSparsePlots(mcArray,mInv,origD0,PDGD0,track,ptbin,kTrack); //fills for charged tracks
1912
1913       //retrieving leading info...
1914       if(track->Pt() > highPt) {
1915         lead[0] = fCorrelatorTr->GetDeltaPhi();
1916         lead[1] = fCorrelatorTr->GetDeltaEta();
1917         lead[2] = fReadMC ? CheckTrackOrigin(mcArray,trkMC) : 0;
1918         highPt = track->Pt();
1919       }
1920
1921     } // end of tracks loop
1922   } //end of event loop for fCorrelatorTr
1923
1924  if(fKaonCorr) { //loops for Kcharg and K0
1925
1926   if(fMixing) {
1927     NofEventsinPool = fCorrelatorKc->GetNofEventsInPool();
1928     if(!execPoolKc) {
1929       AliInfo("Mixed event analysis: K+/- pool is not ready");
1930       NofEventsinPool = 0;
1931     }
1932   }
1933
1934   //Charged Kaons loop
1935   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)
1936     Bool_t analyzetracksKc = fCorrelatorKc->ProcessAssociatedTracks(jMix);
1937     if(!analyzetracksKc) {
1938       AliInfo("AliHFCorrelator::Cannot process the K+/- array");
1939       continue;
1940     }
1941
1942     for(Int_t iTrack = 0; iTrack<fCorrelatorKc->GetNofTracks(); iTrack++){ // looping on charged kaons candidates
1943
1944       Bool_t runcorrelation = fCorrelatorKc->Correlate(iTrack);
1945       if(!runcorrelation) continue;
1946
1947       AliReducedParticle* kCharg = fCorrelatorKc->GetAssociatedParticle();
1948       if(kCharg->GetLabel()<1) continue;
1949       if(kCharg->Pt() < fPtThreshLow.at(ptbin) || kCharg->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1950       if(TMath::Abs(kCharg->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1951       if(!fMixing) N_KCharg++;
1952
1953       AliAODMCParticle *kChargMC = (AliAODMCParticle*)mcArray->At(kCharg->GetLabel());
1954       if(!kChargMC) continue;
1955
1956       if (IsDDaughter(d,kChargMC,mcArray)) continue;
1957       FillSparsePlots(mcArray,mInv,origD0,PDGD0,kCharg,ptbin,kKCharg); //fills for charged tracks
1958
1959     } // end of charged kaons loop
1960   } //end of event loop for fCorrelatorKc
1961
1962   if(fMixing) {
1963     NofEventsinPool = fCorrelatorK0->GetNofEventsInPool();
1964     if(!execPoolK0) {
1965       AliInfo("Mixed event analysis: K0 pool is not ready");
1966       NofEventsinPool = 0;
1967     }
1968   }
1969
1970   //K0 loop
1971   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)
1972     Bool_t analyzetracksK0 = fCorrelatorK0->ProcessAssociatedTracks(jMix);
1973     if(!analyzetracksK0) {
1974       AliInfo("AliHFCorrelator::Cannot process the K0 array");
1975       continue;
1976     }
1977
1978     for(Int_t iTrack = 0; iTrack<fCorrelatorK0->GetNofTracks(); iTrack++){ // looping on k0 candidates
1979
1980       Bool_t runcorrelation = fCorrelatorK0->Correlate(iTrack);
1981       if(!runcorrelation) continue;
1982
1983       AliReducedParticle* k0 = fCorrelatorK0->GetAssociatedParticle();
1984       if(k0->GetLabel()<1) continue;
1985       if(k0->Pt() < fPtThreshLow.at(ptbin) || k0->Pt() > fPtThreshUp.at(ptbin)) continue; //discard tracks outside pt range for hadrons/K
1986       if(TMath::Abs(k0->Eta())>0.8) continue; //discard tracks outside barrel (since it's kinematic MC and produces tracks all over rapidity region
1987
1988       AliAODMCParticle *k0MC = (AliAODMCParticle*)mcArray->At(k0->GetLabel());
1989       if(!k0MC) continue;
1990
1991       if (IsDDaughter(d,k0MC,mcArray)) continue;
1992       FillSparsePlots(mcArray,mInv,origD0,PDGD0,k0,ptbin,kK0); //fills for charged tracks
1993
1994       if(!fMixing) N_Kaons++;
1995
1996     } // end of charged kaons loop
1997   } //end of event loop for fCorrelatorK0
1998
1999  } //end of 'if(fKaonCorr)'
2000
2001   Double_t fillSpLeadMC[4] = {lead[0],mD0,lead[1],0.4}; //mD0 = mD0bar = 1.864
2002
2003   //leading track correlations fill
2004   if(!fMixing) {
2005     if(d->GetPdgCode()==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0
2006       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC); //c and b D0
2007       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
2008       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);
2009       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);
2010       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
2011     }
2012     if(d->GetPdgCode()==-421 && fIsSelectedCandidate>1) { //D0bar
2013       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_Bin%d",ptbin)))->Fill(fillSpLeadMC);
2014       ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpLeadMC); //c or b D0
2015       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);
2016       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);
2017       if((int)lead[2]==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Lead_NonHF_Bin%d",ptbin)))->Fill(fillSpLeadMC);  //non HF
2018     }
2019
2020     //Fill of count histograms
2021     if (!fAlreadyFilled) {
2022       ((TH1F*)fOutputStudy->FindObject("hist_Count_Charg"))->Fill(N_Charg);
2023       ((TH1F*)fOutputStudy->FindObject("hist_Count_Kcharg"))->Fill(N_KCharg);
2024       ((TH1F*)fOutputStudy->FindObject("hist_Count_K0"))->Fill(N_Kaons);
2025     }
2026
2027     fAlreadyFilled=kTRUE; //at least a D0 analyzed in the event; distribution plots already filled
2028   }
2029
2030 }
2031
2032 //________________________________________________________________________
2033 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) {
2034   //
2035   //fills the THnSparse for correlations, calculating the variables
2036   //
2037
2038   //Initialization of variables
2039   Double_t mD0, mD0bar, deltaphi = 0., deltaeta = 0.;
2040   mD0 = mInv[0];
2041   mD0bar = mInv[1];
2042
2043   if (fReadMC && track->GetLabel()<1) return;
2044   if (fReadMC && !(AliAODMCParticle*)mcArray->At(track->GetLabel())) return;
2045   Double_t ptTrack = track->Pt();
2046   Double_t d0Track = type!=kK0 ? track->GetImpPar() : 0.;
2047   Double_t phiTr = track->Phi();
2048   Double_t origTr = fReadMC ? CheckTrackOrigin(mcArray,(AliAODMCParticle*)mcArray->At(track->GetLabel())) : 0;
2049
2050   TString part = "", orig = "";
2051
2052   switch (type) {
2053     case(kTrack): {
2054       part = "Charg";
2055       deltaphi = fCorrelatorTr->GetDeltaPhi();
2056       deltaeta = fCorrelatorTr->GetDeltaEta();
2057       break;
2058     }
2059     case(kKCharg): {
2060       part = "Kcharg";
2061       deltaphi = fCorrelatorKc->GetDeltaPhi();
2062       deltaeta = fCorrelatorKc->GetDeltaEta();
2063       break;
2064     }
2065     case(kK0): {
2066       part = "K0";
2067       deltaphi = fCorrelatorK0->GetDeltaPhi();
2068       deltaeta = fCorrelatorK0->GetDeltaEta();
2069       break;
2070     }
2071   }
2072
2073   if(fMixing == kSE) {
2074
2075     //Fixes limits; needed to include overflow into THnSparse projections!
2076     Double_t pTorig = track->Pt();
2077     Double_t d0orig = track->GetImpPar();
2078     Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(2)->GetXmax(); //all plots have same axes...
2079     Double_t displLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(3)->GetXmax();
2080     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d",ptbin)))->GetAxis(4)->GetXmax();
2081     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
2082     if(d0Track > displLim_Sparse) d0Track = (displLim_Sparse-0.001);
2083     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
2084     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
2085
2086     //variables for filling histos
2087     Double_t fillSpPhiD0[5] = {deltaphi,mD0,ptTrack,d0Track,deltaeta};
2088     Double_t fillSpPhiD0bar[5] = {deltaphi,mD0bar,ptTrack,d0Track,deltaeta};
2089     Double_t fillSpWeigD0[4] = {deltaphi,mD0,deltaeta,ptTrack};
2090     Double_t fillSpWeigD0bar[4] = {deltaphi,mD0bar,deltaeta,ptTrack};
2091
2092     if(fReadMC == 0) {
2093       //sparse fill for data (tracks, K+-, K0) + weighted
2094       if(fIsSelectedCandidate == 1 || fIsSelectedCandidate == 3) { //D0
2095         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2096         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2097       }
2098       if(fIsSelectedCandidate == 2 || fIsSelectedCandidate == 3) { //D0bar
2099         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2100         ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2101       }
2102       if(!fAlreadyFilled) {
2103         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
2104         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
2105       }
2106     }
2107
2108     if(fReadMC) {
2109
2110       if(origD0==4) {orig = "_From_c";} else {orig = "_From_b";}
2111
2112       //sparse fill for data (tracks, K+-, K0) + weighted
2113       if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3)) { //D0 (from MCTruth)
2114          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2115          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2116          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);
2117          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);
2118          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2119          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2120          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2121          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2122          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2123          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0,pTorig*wg);
2124       }
2125       if(PdgD0==-421 && fIsSelectedCandidate>1) { //D0bar
2126          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2127          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2128          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);
2129          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);
2130          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_NonHF_Bin%d",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2131          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2132          ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2133          if(origD0==4&&origTr>=1&&origTr<=3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2134          if(origD0==5&&origTr>=4&&origTr<=8) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_HF%s_Bin%d",orig.Data(),ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2135          if(origTr==0) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Weig_NonHF_Bin%d",ptbin)))->Fill(fillSpWeigD0bar,pTorig*wg);
2136       }
2137       if(!fAlreadyFilled) {
2138         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_Bin%d",part.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
2139         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF_Bin%d",part.Data(),ptbin)))->Fill(d0orig);
2140         if (origTr>=1&&origTr<=8) ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s_HF%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig);
2141         ((TH1F*)fOutputStudy->FindObject(Form("histDispl_%s%s_Bin%d",part.Data(),orig.Data(),ptbin)))->Fill(d0orig); //Fills displacement histos
2142         ((TH1F*)fOutputStudy->FindObject(Form("hist_Pt_%s_Bin%d",part.Data(),ptbin)))->Fill(pTorig);
2143         ((TH1F*)fOutputStudy->FindObject(Form("histOrig_%s_Bin%d",part.Data(),ptbin)))->Fill(origTr);
2144         ((TH1F*)fOutputStudy->FindObject(Form("hist_PhiDistr_%s",part.Data())))->Fill(phiTr);
2145       }
2146     }//end MC case
2147
2148   } //end of SE fill
2149
2150   if(fMixing == kME) {
2151
2152     //Fixes limits; needed to include overflow into THnSparse projections!
2153     Double_t EtaLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(2)->GetXmax();
2154     if(deltaeta > EtaLim_Sparse) deltaeta = EtaLim_Sparse-0.01;
2155     if(deltaeta < -EtaLim_Sparse) deltaeta = -EtaLim_Sparse+0.01;
2156     Double_t ptLim_Sparse = ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_Charg_Bin%d_EvMix",ptbin)))->GetAxis(3)->GetXmax(); //all plots have same axes...
2157     if(ptTrack > ptLim_Sparse) ptTrack = ptLim_Sparse-0.01;
2158
2159     //variables for filling histos
2160     Double_t fillSpPhiD0[4] = {deltaphi,mD0,deltaeta,0.4}; //dummy for ME threshold! unless explicitly set by flag...
2161     Double_t fillSpPhiD0bar[4] = {deltaphi,mD0bar,deltaeta,0.4};
2162     if(fMEAxisThresh) {
2163       fillSpPhiD0[3] = ptTrack;
2164       fillSpPhiD0bar[3] = ptTrack;
2165     }
2166
2167     if(fReadMC == 0) {
2168       //sparse fill for data (tracks, K+-, K0)
2169       if(fIsSelectedCandidate == 1||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2170       if(fIsSelectedCandidate == 2||fIsSelectedCandidate == 3) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2171     }
2172     if(fReadMC == 1) {
2173       //sparse fill for data (tracks, K+-, K0)
2174       if(PdgD0==421 && (fIsSelectedCandidate==1||fIsSelectedCandidate==3))  ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0,wg);
2175       if(PdgD0==-421 && fIsSelectedCandidate>1) ((THnSparseF*)fOutputCorr->FindObject(Form("hPhi_%s_Bin%d_EvMix",part.Data(),ptbin)))->Fill(fillSpPhiD0bar,wg);
2176     }//end MC case
2177
2178   } //end of ME fill
2179
2180   return;
2181 }
2182
2183 //_________________________________________________________________________________________________
2184 Int_t AliAnalysisTaskSED0Correlations::CheckTrackOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2185   //
2186   // checks on particle (#) origin:
2187   // 0) Not HF
2188   // 1) D->#
2189   // 2) D->X->#
2190   // 3) c hadronization
2191   // 4) B->#
2192   // 5) B->X-># (X!=D)
2193   // 6) B->D->#
2194   // 7) B->D->X->#
2195   // 8) b hadronization
2196   //
2197   if(fDebug>2) printf("AliAnalysisTaskSED0Correlations::CheckTrkOrigin() \n");
2198
2199   Int_t pdgGranma = 0;
2200   Int_t mother = 0;
2201   mother = mcPartCandidate->GetMother();
2202   Int_t istep = 0;
2203   Int_t abspdgGranma =0;
2204   Bool_t isFromB=kFALSE;
2205   Bool_t isDdaugh=kFALSE;
2206   Bool_t isDchaindaugh=kFALSE;
2207   Bool_t isBdaugh=kFALSE;
2208   Bool_t isBchaindaugh=kFALSE;
2209   Bool_t isQuarkFound=kFALSE;
2210
2211   if (mother<0) return -1;
2212   while (mother >= 0){
2213     istep++;
2214     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2215     if (mcMoth){
2216       pdgGranma = mcMoth->GetPdgCode();
2217       abspdgGranma = TMath::Abs(pdgGranma);
2218       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2219         isBchaindaugh=kTRUE;
2220         if(istep==1) isBdaugh=kTRUE;
2221       }
2222       if ((abspdgGranma > 400 && abspdgGranma < 500) || (abspdgGranma > 4000 && abspdgGranma < 5000)){
2223         isDchaindaugh=kTRUE;
2224         if(istep==1) isDdaugh=kTRUE;
2225       }
2226       if(abspdgGranma==4 || abspdgGranma==5) {isQuarkFound=kTRUE; if(abspdgGranma==5) isFromB = kTRUE;}
2227       mother = mcMoth->GetMother();
2228     }else{
2229       AliError("Failed casting the mother particle!");
2230       return -1;
2231     }
2232   }
2233
2234   //decides what to return based on the flag status
2235   if(isQuarkFound) {
2236     if(!isFromB) {  //charm
2237       if(isDdaugh) return 1; //charm immediate
2238       else if(isDchaindaugh) return 2; //charm chain
2239       else return 3; //charm hadronization
2240     }
2241     else { //beauty
2242       if(isBdaugh) return 4; //b immediate
2243       else if(isBchaindaugh) { //b chain
2244         if(isDchaindaugh) {
2245           if(isDdaugh) return 6; //d immediate
2246           return 7; //d chain
2247           }
2248         else return 5; //b, not d
2249       }
2250       else return 8; //b hadronization
2251     }
2252   }
2253   else if(!isDdaugh && !isDchaindaugh && !isBdaugh && !isBchaindaugh) return 0; //no HF decay
2254      //in this case, it's !isQuarkFound, but not in 100% cases it's a non HF particle!
2255      //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...
2256
2257   return -1; //some problem spotted
2258 }
2259 //________________________________________________________________________
2260 Bool_t AliAnalysisTaskSED0Correlations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* mcArray) const {
2261
2262   //Daughter removal in MCKine case
2263   Bool_t isDaughter = kFALSE;
2264   Int_t labelD0 = d->GetLabel();
2265
2266   Int_t mother = track->GetMother();
2267
2268   //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)
2269   while (mother > 0){
2270     AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother)); //it's the mother of the track!
2271     if (mcMoth){
2272       if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;
2273       mother = mcMoth->GetMother(); //goes back by one
2274     } else{
2275       AliError("Failed casting the mother particle!");
2276       break;
2277     }
2278   }
2279
2280   return isDaughter;
2281 }
2282
2283 //________________________________________________________________________
2284 Int_t AliAnalysisTaskSED0Correlations::PtBinCorr(Double_t pt) const {
2285   //
2286   //give the pt bin where the pt lies.
2287   //
2288   Int_t ptbin=-1;
2289   if(pt<fBinLimsCorr.at(0)) return ptbin; //out of bounds
2290
2291   Int_t i = 0;
2292   while(pt>fBinLimsCorr.at(i)) {ptbin=i; i++;}
2293
2294   return ptbin;
2295 }
2296
2297 //---------------------------------------------------------------------------
2298 Bool_t AliAnalysisTaskSED0Correlations::SelectV0(AliAODv0* v0, AliAODVertex *vtx, Int_t opt, Int_t idArrayV0[][2]) const
2299 {
2300   //
2301   // Selection for K0 hypotheses
2302   // options: 1 = selects mass invariant about 3 sigma inside the peak + threshold of 0.3 GeV
2303   //          2 = no previous selections
2304
2305   if(!fCutsTracks->IsKZeroSelected(v0,vtx)) return kFALSE;
2306
2307   AliAODTrack *v0Daug1 = (AliAODTrack*)v0->GetDaughter(0);
2308   AliAODTrack *v0Daug2 = (AliAODTrack*)v0->GetDaughter(1);
2309
2310   if(opt==1) { //additional cuts for correlations (V0 has to be closer than 3 sigma from K0 mass)
2311     if(TMath::Abs(v0->MassK0Short()-0.4976) > 3*0.004) return kFALSE;
2312   }
2313
2314   //This part removes double counting for swapped tracks!
2315   Int_t i = 0;  //while loop (until the last-written entry pair of ID!
2316   while(idArrayV0[i][0]!=-2 && idArrayV0[i][1]!=-2) {
2317     if((v0Daug1->GetID()==idArrayV0[i][0] && v0Daug2->GetID()==idArrayV0[i][1])||
2318        (v0Daug1->GetID()==idArrayV0[i][1] && v0Daug2->GetID()==idArrayV0[i][0])) return kFALSE;
2319     i++;
2320   }
2321   idArrayV0[i][0]=v0Daug1->GetID();
2322   idArrayV0[i][1]=v0Daug2->GetID();
2323
2324   return kTRUE;
2325 }
2326
2327 //---------------------------------------------------------------------------
2328 Bool_t AliAnalysisTaskSED0Correlations::IsSoftPion_MCKine(AliAODMCParticle* d, AliAODMCParticle* track, TClonesArray* arrayMC) const
2329 {
2330   //
2331   // Removes soft pions in Kine
2332
2333   //Daughter removal in MCKine case
2334   Bool_t isSoftPi = kFALSE;
2335   Int_t labelD0 = d->GetLabel();
2336
2337   Int_t mother = track->GetMother();
2338   if(mother<0) return isSoftPi; //safety check
2339
2340   AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother)); //it's the mother of the track!
2341   if(!mcMoth){
2342     return isSoftPi;
2343   }
2344   if(TMath::Abs(mcMoth->GetPdgCode())==413 && mcMoth->GetNDaughters()==2) { //mother is D* with 2 daughs
2345     Int_t labdau1 = mcMoth->GetDaughter(0);
2346     Int_t labdau2 = mcMoth->GetDaughter(1);
2347     AliAODMCParticle* dau1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau1));
2348     AliAODMCParticle* dau2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(labdau2));
2349     if(!dau1 || !dau2) return isSoftPi; //safety check
2350     if(dau1->GetLabel()==labelD0 || dau2->GetLabel()==labelD0) { //one of the daughs is the D0 trigger
2351       if((TMath::Abs(dau1->GetPdgCode())==421 && TMath::Abs(dau2->GetPdgCode())==211)||(TMath::Abs(dau1->GetPdgCode())==211 && TMath::Abs(dau2->GetPdgCode())==421)) {
2352         isSoftPi = kTRUE; //ok, soft pion was found
2353         return isSoftPi;
2354       }
2355     }
2356   }
2357
2358   return isSoftPi;
2359 }
2360
2361 //________________________________________________________________________
2362 void AliAnalysisTaskSED0Correlations::PrintBinsAndLimits() {
2363
2364   cout << "--------------------------\n";
2365   cout << "PtBins = " << fNPtBinsCorr << "\n";
2366   cout << "PtBin limits--------------\n";
2367   for (int i=0; i<fNPtBinsCorr; i++) {
2368     cout << "Bin "<<i+1<<" = "<<fBinLimsCorr.at(i)<<" to "<<fBinLimsCorr.at(i+1)<<"\n";
2369   }
2370   cout << "\n--------------------------\n";
2371   cout << "PtBin tresh. tracks low---\n";
2372   for (int i=0; i<fNPtBinsCorr; i++) {
2373     cout << fPtThreshLow.at(i) << "; ";
2374   }
2375   cout << "PtBin tresh. tracks up----\n";
2376   for (int i=0; i<fNPtBinsCorr; i++) {
2377     cout << fPtThreshUp.at(i) << "; ";
2378   }
2379   cout << "\n--------------------------\n";
2380   cout << "D0 Eta cut for Correl = "<<fEtaForCorrel<<"\n";
2381   cout << "--------------------------\n";
2382   cout << "MC Truth = "<<fReadMC<<" - MC Reco: Trk = "<<fRecoTr<<", D0 = "<<fRecoD0<<"\n";
2383   cout << "--------------------------\n";
2384   cout << "Sel of Event tpye (PP,GS,FE,...)= "<<fSelEvType<<"\n";
2385   cout << "--------------------------\n";
2386   cout << "Ev Mixing = "<<fMixing<<"\n";
2387   cout << "--------------------------\n";
2388   cout << "ME thresh axis = "<<fMEAxisThresh<<"\n";
2389   cout << "--------------------------\n";
2390   cout << "Soft Pi Cut = "<<fSoftPiCut<<"\n";
2391   cout << "--------------------------\n";
2392 }
2393