PWGHF/correlationHF/AliAnalysisTaskSED0Correlations.cxx

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