PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskV0sInJetsEmcal.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // task for analysis of V0s (K0S, (anti-)Lambda) in charged jets
  17 // fork of AliAnalysisTaskV0sInJets for the EMCal framework
  18 // Author: Vit Kucera (vit.kucera@cern.ch)
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH1D.h"
  23 #include "TH2D.h"
  24 #include "THnSparse.h"
  25 #include "TCanvas.h"
  26
  27 #include "AliAnalysisTask.h"
  28 #include "AliAnalysisManager.h"
  29
  30 #include "AliESDEvent.h"
  31 #include "AliAODEvent.h"
  32 #include "AliAODTrack.h"
  33 #include <TDatabasePDG.h>
  34 #include <TPDGCode.h>
  35 #include "AliPIDResponse.h"
  36 #include "AliInputEventHandler.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODMCParticle.h"
  39 #include "TClonesArray.h"
  40 #include "TRandom3.h"
  41
  42 #include "AliVCluster.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliESDCaloCluster.h"
  45 #include "AliVTrack.h"
  46 #include "AliEmcalJet.h"
  47 #include "AliRhoParameter.h"
  48 #include "AliLog.h"
  49 #include "AliJetContainer.h"
  50 #include "AliParticleContainer.h"
  51 #include "AliClusterContainer.h"
  52 #include "AliPicoTrack.h"
  53
  54 #include "AliAnalysisTaskV0sInJetsEmcal.h"
  55
  56 ClassImp(AliAnalysisTaskV0sInJetsEmcal)
  57
  58 // upper edges of centrality bins
  59 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsCent] = {10, 30, 50, 80}; // Alice Zimmermann
  60 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsCent] = {10, 20, 40, 60, 80}; // Vit Kucera, initial binning
  61 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsCent] = {5, 10, 20, 40, 60, 80}; // Iouri Belikov, LF analysis
  62 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsCent] = {10}; // only central
  63
  64 // axis: pT of V0
  65 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtV0[2] = {0, 12};
  66 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtV0 = sizeof(AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtV0) / sizeof((AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtV0)[0]) - 1;
  67 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtV0Init = int(((AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtV0)[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtV0] - (AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtV0)[0]) / 0.1); // bin width 0.1 GeV/c
  68 // axis: pT of jets
  69 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet[2] = {0, 100};
  70 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJet = sizeof(AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet) / sizeof(AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet[0]) - 1;
  71 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJet] - (AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[0]) / 5.); // bin width 5 GeV/c
  72 // axis: K0S invariant mass
  73 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsMassK0s = 300;
  74 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdMassK0sMin = 0.35; // [GeV/c^2]
  75 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdMassK0sMax = 0.65; // [GeV/c^2]
  76 // axis: Lambda invariant mass
  77 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsMassLambda = 200;
  78 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdMassLambdaMin = 1.05; // [GeV/c^2]
  79 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdMassLambdaMax = 1.25; // [GeV/c^2]
  80
  81
  82 // Default constructor
  83 AliAnalysisTaskV0sInJetsEmcal::AliAnalysisTaskV0sInJetsEmcal():
  84   AliAnalysisTaskEmcalJet(),
  85   fAODIn(0),
  86   fAODOut(0),
  87   fOutputListStd(0),
  88   fOutputListQA(0),
  89   fOutputListCuts(0),
  90   fOutputListMC(0),
  91   fbIsPbPb(1),
  92
  93   fdCutDCAToPrimVtxMin(0.1),
  94   fdCutDCADaughtersMax(1.),
  95   fdCutNSigmadEdxMax(3),
  96   fdCutCPAMin(0.998),
  97   fdCutNTauMax(5),
  98
  99   fdCutPtJetMin(0),
 100   fdCutPtTrackMin(5),
 101   fdRadiusJet(0.4),
 102   fbJetSelection(0),
 103   fbMCAnalysis(0),
 104   fRandom(0),
 105
 106   fJetsCont(0),
 107   fJetsBgCont(0),
 108   fTracksCont(0),
 109   fCaloClustersCont(0),
 110
 111   fdCutVertexZ(10),
 112   fdCutVertexR2(1),
 113   fdCutCentLow(0),
 114   fdCutCentHigh(80),
 115
 116   fdCentrality(0),
 117   fh1EventCounterCut(0),
 118   fh1EventCent(0),
 119   fh1EventCent2(0),
 120   fh1EventCent2Jets(0),
 121   fh1EventCent2NoJets(0),
 122   fh2EventCentTracks(0),
 123   fh1V0CandPerEvent(0),
 124   fh1NRndConeCent(0),
 125   fh1NMedConeCent(0),
 126   fh1AreaExcluded(0),
 127
 128   fh2CCK0s(0),
 129   fh2CCLambda(0),
 130   fh3CCMassCorrelBoth(0),
 131   fh3CCMassCorrelKNotL(0),
 132   fh3CCMassCorrelLNotK(0)
 133 {
 134   for(Int_t i = 0; i < fgkiNQAIndeces; i++)
 135   {
 136     fh1QAV0Status[i] = 0;
 137     fh1QAV0TPCRefit[i] = 0;
 138     fh1QAV0TPCRows[i] = 0;
 139     fh1QAV0TPCFindable[i] = 0;
 140     fh1QAV0TPCRowsFind[i] = 0;
 141     fh1QAV0Eta[i] = 0;
 142     fh2QAV0EtaRows[i] = 0;
 143     fh2QAV0PtRows[i] = 0;
 144     fh2QAV0PhiRows[i] = 0;
 145     fh2QAV0NClRows[i] = 0;
 146     fh2QAV0EtaNCl[i] = 0;
 147
 148     fh2QAV0EtaPtK0sPeak[i] = 0;
 149     fh2QAV0EtaEtaK0s[i] = 0;
 150     fh2QAV0PhiPhiK0s[i] = 0;
 151     fh1QAV0RapK0s[i] = 0;
 152     fh2QAV0PtPtK0sPeak[i] = 0;
 153     fh2ArmPodK0s[i] = 0;
 154
 155     fh2QAV0EtaPtLambdaPeak[i] = 0;
 156     fh2QAV0EtaEtaLambda[i] = 0;
 157     fh2QAV0PhiPhiLambda[i] = 0;
 158     fh1QAV0RapLambda[i] = 0;
 159     fh2QAV0PtPtLambdaPeak[i] = 0;
 160     fh2ArmPodLambda[i] = 0;
 161
 162     fh2QAV0EtaPtALambdaPeak[i] = 0;
 163     fh2QAV0EtaEtaALambda[i] = 0;
 164     fh2QAV0PhiPhiALambda[i] = 0;
 165     fh1QAV0RapALambda[i] = 0;
 166     fh2QAV0PtPtALambdaPeak[i] = 0;
 167     fh2ArmPodALambda[i] = 0;
 168
 169     fh1QAV0Pt[i] = 0;
 170     fh1QAV0Charge[i] = 0;
 171     fh1QAV0DCAVtx[i] = 0;
 172     fh1QAV0DCAV0[i] = 0;
 173     fh1QAV0Cos[i] = 0;
 174     fh1QAV0R[i] = 0;
 175     fh1QACTau2D[i] = 0;
 176     fh1QACTau3D[i] = 0;
 177
 178     fh2ArmPod[i] = 0;
 179
 180     /*
 181     fh2CutTPCRowsK0s[i] = 0;
 182     fh2CutTPCRowsLambda[i] = 0;
 183     fh2CutPtPosK0s[i] = 0;
 184     fh2CutPtNegK0s[i] = 0;
 185     fh2CutPtPosLambda[i] = 0;
 186     fh2CutPtNegLambda[i] = 0;
 187     fh2CutDCAVtx[i] = 0;
 188     fh2CutDCAV0[i] = 0;
 189     fh2CutCos[i] = 0;
 190     fh2CutR[i] = 0;
 191     fh2CutEtaK0s[i] = 0;
 192     fh2CutEtaLambda[i] = 0;
 193     fh2CutRapK0s[i] = 0;
 194     fh2CutRapLambda[i] = 0;
 195     fh2CutCTauK0s[i] = 0;
 196     fh2CutCTauLambda[i] = 0;
 197     fh2CutPIDPosK0s[i] = 0;
 198     fh2CutPIDNegK0s[i] = 0;
 199     fh2CutPIDPosLambda[i] = 0;
 200     fh2CutPIDNegLambda[i] = 0;
 201
 202     fh2Tau3DVs2D[i] = 0;
 203     */
 204   }
 205   for(Int_t i = 0; i < fgkiNCategV0; i++)
 206   {
 207     fh1V0InvMassK0sAll[i] = 0;
 208     fh1V0InvMassLambdaAll[i] = 0;
 209     fh1V0InvMassALambdaAll[i] = 0;
 210   }
 211   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 212   {
 213     fh1EventCounterCutCent[i] = 0;
 214     fh1V0CounterCentK0s[i] = 0;
 215     fh1V0CounterCentLambda[i] = 0;
 216     fh1V0CounterCentALambda[i] = 0;
 217     fh1V0CandPerEventCentK0s[i] = 0;
 218     fh1V0CandPerEventCentLambda[i] = 0;
 219     fh1V0CandPerEventCentALambda[i] = 0;
 220     fh1V0InvMassK0sCent[i] = 0;
 221     fh1V0InvMassLambdaCent[i] = 0;
 222     fh1V0InvMassALambdaCent[i] = 0;
 223     fh1V0K0sPtMCGen[i] = 0;
 224     fh2V0K0sPtMassMCRec[i] = 0;
 225     fh1V0K0sPtMCRecFalse[i] = 0;
 226     fh2V0K0sEtaPtMCGen[i] = 0;
 227     fh3V0K0sEtaPtMassMCRec[i] = 0;
 228     fh2V0K0sInJetPtMCGen[i] = 0;
 229     fh3V0K0sInJetPtMassMCRec[i] = 0;
 230     fh3V0K0sInJetEtaPtMCGen[i] = 0;
 231     fh4V0K0sInJetEtaPtMassMCRec[i] = 0;
 232     fh2V0K0sMCResolMPt[i] = 0;
 233     fh2V0K0sMCPtGenPtRec[i] = 0;
 234     fh1V0LambdaPtMCGen[i] = 0;
 235     fh2V0LambdaPtMassMCRec[i] = 0;
 236     fh1V0LambdaPtMCRecFalse[i] = 0;
 237     fh2V0LambdaEtaPtMCGen[i] = 0;
 238     fh3V0LambdaEtaPtMassMCRec[i] = 0;
 239     fh2V0LambdaInJetPtMCGen[i] = 0;
 240     fh3V0LambdaInJetPtMassMCRec[i] = 0;
 241     fh3V0LambdaInJetEtaPtMCGen[i] = 0;
 242     fh4V0LambdaInJetEtaPtMassMCRec[i] = 0;
 243     fh2V0LambdaMCResolMPt[i] = 0;
 244     fh2V0LambdaMCPtGenPtRec[i] = 0;
 245     fhnV0LambdaInclMCFD[i] = 0;
 246     fhnV0LambdaInJetsMCFD[i] = 0;
 247     fhnV0LambdaBulkMCFD[i] = 0;
 248     fh1V0XiPtMCGen[i] = 0;
 249     fh1V0ALambdaPt[i] = 0;
 250     fh1V0ALambdaPtMCGen[i] = 0;
 251     fh1V0ALambdaPtMCRec[i] = 0;
 252     fh2V0ALambdaPtMassMCRec[i] = 0;
 253     fh1V0ALambdaPtMCRecFalse[i] = 0;
 254     fh2V0ALambdaEtaPtMCGen[i] = 0;
 255     fh3V0ALambdaEtaPtMassMCRec[i] = 0;
 256     fh2V0ALambdaInJetPtMCGen[i] = 0;
 257     fh2V0ALambdaInJetPtMCRec[i] = 0;
 258     fh3V0ALambdaInJetPtMassMCRec[i] = 0;
 259     fh3V0ALambdaInJetEtaPtMCGen[i] = 0;
 260     fh4V0ALambdaInJetEtaPtMassMCRec[i] = 0;
 261     fh2V0ALambdaMCResolMPt[i] = 0;
 262     fh2V0ALambdaMCPtGenPtRec[i] = 0;
 263     fhnV0ALambdaInclMCFD[i] = 0;
 264     fhnV0ALambdaInJetsMCFD[i] = 0;
 265     fhnV0ALambdaBulkMCFD[i] = 0;
 266     fh1V0AXiPtMCGen[i] = 0;
 267
 268     // eta daughters
 269 //      fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
 270     fhnV0K0sInclDaughterEtaPtPtMCRec[i] = 0;
 271 //      fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
 272     fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = 0;
 273 //      fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
 274     fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = 0;
 275 //      fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
 276     fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
 277 //      fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
 278     fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = 0;
 279 //      fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
 280     fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
 281
 282     // Inclusive
 283     fhnV0InclusiveK0s[i] = 0;
 284     fhnV0InclusiveLambda[i] = 0;
 285     fhnV0InclusiveALambda[i] = 0;
 286     // Cones
 287     fhnV0InJetK0s[i] = 0;
 288     fhnV0InPerpK0s[i] = 0;
 289     fhnV0InRndK0s[i] = 0;
 290     fhnV0InMedK0s[i] = 0;
 291     fhnV0OutJetK0s[i] = 0;
 292     fhnV0NoJetK0s[i] = 0;
 293     fhnV0InJetLambda[i] = 0;
 294     fhnV0InPerpLambda[i] = 0;
 295     fhnV0InRndLambda[i] = 0;
 296     fhnV0InMedLambda[i] = 0;
 297     fhnV0OutJetLambda[i] = 0;
 298     fhnV0NoJetLambda[i] = 0;
 299     fhnV0InJetALambda[i] = 0;
 300     fhnV0InPerpALambda[i] = 0;
 301     fhnV0InRndALambda[i] = 0;
 302     fhnV0InMedALambda[i] = 0;
 303     fhnV0OutJetALambda[i] = 0;
 304     fhnV0NoJetALambda[i] = 0;
 305
 306     fh2V0PtJetAngleK0s[i] = 0;
 307     fh2V0PtJetAngleLambda[i] = 0;
 308     fh2V0PtJetAngleALambda[i] = 0;
 309     fh1DCAInK0s[i] = 0;
 310     fh1DCAInLambda[i] = 0;
 311     fh1DCAInALambda[i] = 0;
 312     fh1DCAOutK0s[i] = 0;
 313     fh1DCAOutLambda[i] = 0;
 314     fh1DCAOutALambda[i] = 0;
 315
 316     fh1PtJet[i] = 0;
 317     fh1EtaJet[i] = 0;
 318     fh2EtaPtJet[i] = 0;
 319     fh1PhiJet[i] = 0;
 320     fh1PtJetTrackLeading[i] = 0;
 321     fh1NJetPerEvent[i] = 0;
 322     fh2EtaPhiRndCone[i] = 0;
 323     fh2EtaPhiMedCone[i] = 0;
 324
 325     fh1VtxZ[i] = 0;
 326     fh2VtxXY[i] = 0;
 327   }
 328 }
 329
 330 // Constructor
 331 AliAnalysisTaskV0sInJetsEmcal::AliAnalysisTaskV0sInJetsEmcal(const char* name):
 332   AliAnalysisTaskEmcalJet(name, kTRUE),
 333   fAODIn(0),
 334   fAODOut(0),
 335   fOutputListStd(0),
 336   fOutputListQA(0),
 337   fOutputListCuts(0),
 338   fOutputListMC(0),
 339   fbIsPbPb(1),
 340
 341   fdCutDCAToPrimVtxMin(0.1),
 342   fdCutDCADaughtersMax(1.),
 343   fdCutNSigmadEdxMax(3),
 344   fdCutCPAMin(0.998),
 345   fdCutNTauMax(5),
 346
 347   fdCutPtJetMin(0),
 348   fdCutPtTrackMin(5),
 349   fdRadiusJet(0.4),
 350   fbJetSelection(0),
 351   fbMCAnalysis(0),
 352   fRandom(0),
 353
 354   fJetsCont(0),
 355   fJetsBgCont(0),
 356   fTracksCont(0),
 357   fCaloClustersCont(0),
 358
 359   fdCutVertexZ(10),
 360   fdCutVertexR2(1),
 361   fdCutCentLow(0),
 362   fdCutCentHigh(80),
 363
 364   fdCentrality(0),
 365   fh1EventCounterCut(0),
 366   fh1EventCent(0),
 367   fh1EventCent2(0),
 368   fh1EventCent2Jets(0),
 369   fh1EventCent2NoJets(0),
 370   fh2EventCentTracks(0),
 371   fh1V0CandPerEvent(0),
 372   fh1NRndConeCent(0),
 373   fh1NMedConeCent(0),
 374   fh1AreaExcluded(0),
 375
 376   fh2CCK0s(0),
 377   fh2CCLambda(0),
 378   fh3CCMassCorrelBoth(0),
 379   fh3CCMassCorrelKNotL(0),
 380   fh3CCMassCorrelLNotK(0)
 381 {
 382   for(Int_t i = 0; i < fgkiNQAIndeces; i++)
 383   {
 384     fh1QAV0Status[i] = 0;
 385     fh1QAV0TPCRefit[i] = 0;
 386     fh1QAV0TPCRows[i] = 0;
 387     fh1QAV0TPCFindable[i] = 0;
 388     fh1QAV0TPCRowsFind[i] = 0;
 389     fh1QAV0Eta[i] = 0;
 390     fh2QAV0EtaRows[i] = 0;
 391     fh2QAV0PtRows[i] = 0;
 392     fh2QAV0PhiRows[i] = 0;
 393     fh2QAV0NClRows[i] = 0;
 394     fh2QAV0EtaNCl[i] = 0;
 395
 396     fh2QAV0EtaPtK0sPeak[i] = 0;
 397     fh2QAV0EtaEtaK0s[i] = 0;
 398     fh2QAV0PhiPhiK0s[i] = 0;
 399     fh1QAV0RapK0s[i] = 0;
 400     fh2QAV0PtPtK0sPeak[i] = 0;
 401     fh2ArmPodK0s[i] = 0;
 402
 403     fh2QAV0EtaPtLambdaPeak[i] = 0;
 404     fh2QAV0EtaEtaLambda[i] = 0;
 405     fh2QAV0PhiPhiLambda[i] = 0;
 406     fh1QAV0RapLambda[i] = 0;
 407     fh2QAV0PtPtLambdaPeak[i] = 0;
 408     fh2ArmPodLambda[i] = 0;
 409
 410     fh2QAV0EtaPtALambdaPeak[i] = 0;
 411     fh2QAV0EtaEtaALambda[i] = 0;
 412     fh2QAV0PhiPhiALambda[i] = 0;
 413     fh1QAV0RapALambda[i] = 0;
 414     fh2QAV0PtPtALambdaPeak[i] = 0;
 415     fh2ArmPodALambda[i] = 0;
 416
 417     fh1QAV0Pt[i] = 0;
 418     fh1QAV0Charge[i] = 0;
 419     fh1QAV0DCAVtx[i] = 0;
 420     fh1QAV0DCAV0[i] = 0;
 421     fh1QAV0Cos[i] = 0;
 422     fh1QAV0R[i] = 0;
 423     fh1QACTau2D[i] = 0;
 424     fh1QACTau3D[i] = 0;
 425
 426     fh2ArmPod[i] = 0;
 427
 428     /*
 429     fh2CutTPCRowsK0s[i] = 0;
 430     fh2CutTPCRowsLambda[i] = 0;
 431     fh2CutPtPosK0s[i] = 0;
 432     fh2CutPtNegK0s[i] = 0;
 433     fh2CutPtPosLambda[i] = 0;
 434     fh2CutPtNegLambda[i] = 0;
 435     fh2CutDCAVtx[i] = 0;
 436     fh2CutDCAV0[i] = 0;
 437     fh2CutCos[i] = 0;
 438     fh2CutR[i] = 0;
 439     fh2CutEtaK0s[i] = 0;
 440     fh2CutEtaLambda[i] = 0;
 441     fh2CutRapK0s[i] = 0;
 442     fh2CutRapLambda[i] = 0;
 443     fh2CutCTauK0s[i] = 0;
 444     fh2CutCTauLambda[i] = 0;
 445     fh2CutPIDPosK0s[i] = 0;
 446     fh2CutPIDNegK0s[i] = 0;
 447     fh2CutPIDPosLambda[i] = 0;
 448     fh2CutPIDNegLambda[i] = 0;
 449
 450     fh2Tau3DVs2D[i] = 0;
 451     */
 452   }
 453   for(Int_t i = 0; i < fgkiNCategV0; i++)
 454   {
 455     fh1V0InvMassK0sAll[i] = 0;
 456     fh1V0InvMassLambdaAll[i] = 0;
 457     fh1V0InvMassALambdaAll[i] = 0;
 458   }
 459   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 460   {
 461     fh1EventCounterCutCent[i] = 0;
 462     fh1V0CounterCentK0s[i] = 0;
 463     fh1V0CounterCentLambda[i] = 0;
 464     fh1V0CounterCentALambda[i] = 0;
 465     fh1V0CandPerEventCentK0s[i] = 0;
 466     fh1V0CandPerEventCentLambda[i] = 0;
 467     fh1V0CandPerEventCentALambda[i] = 0;
 468     fh1V0InvMassK0sCent[i] = 0;
 469     fh1V0InvMassLambdaCent[i] = 0;
 470     fh1V0InvMassALambdaCent[i] = 0;
 471     fh1V0K0sPtMCGen[i] = 0;
 472     fh2V0K0sPtMassMCRec[i] = 0;
 473     fh1V0K0sPtMCRecFalse[i] = 0;
 474     fh2V0K0sEtaPtMCGen[i] = 0;
 475     fh3V0K0sEtaPtMassMCRec[i] = 0;
 476     fh2V0K0sInJetPtMCGen[i] = 0;
 477     fh3V0K0sInJetPtMassMCRec[i] = 0;
 478     fh3V0K0sInJetEtaPtMCGen[i] = 0;
 479     fh4V0K0sInJetEtaPtMassMCRec[i] = 0;
 480     fh2V0K0sMCResolMPt[i] = 0;
 481     fh2V0K0sMCPtGenPtRec[i] = 0;
 482     fh1V0LambdaPtMCGen[i] = 0;
 483     fh2V0LambdaPtMassMCRec[i] = 0;
 484     fh1V0LambdaPtMCRecFalse[i] = 0;
 485     fh2V0LambdaEtaPtMCGen[i] = 0;
 486     fh3V0LambdaEtaPtMassMCRec[i] = 0;
 487     fh2V0LambdaInJetPtMCGen[i] = 0;
 488     fh3V0LambdaInJetPtMassMCRec[i] = 0;
 489     fh3V0LambdaInJetEtaPtMCGen[i] = 0;
 490     fh4V0LambdaInJetEtaPtMassMCRec[i] = 0;
 491     fh2V0LambdaMCResolMPt[i] = 0;
 492     fh2V0LambdaMCPtGenPtRec[i] = 0;
 493     fhnV0LambdaInclMCFD[i] = 0;
 494     fhnV0LambdaInJetsMCFD[i] = 0;
 495     fhnV0LambdaBulkMCFD[i] = 0;
 496     fh1V0XiPtMCGen[i] = 0;
 497     fh1V0ALambdaPt[i] = 0;
 498     fh1V0ALambdaPtMCGen[i] = 0;
 499     fh1V0ALambdaPtMCRec[i] = 0;
 500     fh2V0ALambdaPtMassMCRec[i] = 0;
 501     fh1V0ALambdaPtMCRecFalse[i] = 0;
 502     fh2V0ALambdaEtaPtMCGen[i] = 0;
 503     fh3V0ALambdaEtaPtMassMCRec[i] = 0;
 504     fh2V0ALambdaInJetPtMCGen[i] = 0;
 505     fh2V0ALambdaInJetPtMCRec[i] = 0;
 506     fh3V0ALambdaInJetPtMassMCRec[i] = 0;
 507     fh3V0ALambdaInJetEtaPtMCGen[i] = 0;
 508     fh4V0ALambdaInJetEtaPtMassMCRec[i] = 0;
 509     fh2V0ALambdaMCResolMPt[i] = 0;
 510     fh2V0ALambdaMCPtGenPtRec[i] = 0;
 511     fhnV0ALambdaInclMCFD[i] = 0;
 512     fhnV0ALambdaInJetsMCFD[i] = 0;
 513     fhnV0ALambdaBulkMCFD[i] = 0;
 514     fh1V0AXiPtMCGen[i] = 0;
 515
 516     // eta daughters
 517 //      fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
 518     fhnV0K0sInclDaughterEtaPtPtMCRec[i] = 0;
 519 //      fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
 520     fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = 0;
 521 //      fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
 522     fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = 0;
 523 //      fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
 524     fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
 525 //      fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
 526     fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = 0;
 527 //      fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
 528     fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
 529
 530     // Inclusive
 531     fhnV0InclusiveK0s[i] = 0;
 532     fhnV0InclusiveLambda[i] = 0;
 533     fhnV0InclusiveALambda[i] = 0;
 534     // Cones
 535     fhnV0InJetK0s[i] = 0;
 536     fhnV0InPerpK0s[i] = 0;
 537     fhnV0InRndK0s[i] = 0;
 538     fhnV0InMedK0s[i] = 0;
 539     fhnV0OutJetK0s[i] = 0;
 540     fhnV0NoJetK0s[i] = 0;
 541     fhnV0InJetLambda[i] = 0;
 542     fhnV0InPerpLambda[i] = 0;
 543     fhnV0InRndLambda[i] = 0;
 544     fhnV0InMedLambda[i] = 0;
 545     fhnV0OutJetLambda[i] = 0;
 546     fhnV0NoJetLambda[i] = 0;
 547     fhnV0InJetALambda[i] = 0;
 548     fhnV0InPerpALambda[i] = 0;
 549     fhnV0InRndALambda[i] = 0;
 550     fhnV0InMedALambda[i] = 0;
 551     fhnV0OutJetALambda[i] = 0;
 552     fhnV0NoJetALambda[i] = 0;
 553
 554     fh2V0PtJetAngleK0s[i] = 0;
 555     fh2V0PtJetAngleLambda[i] = 0;
 556     fh2V0PtJetAngleALambda[i] = 0;
 557     fh1DCAInK0s[i] = 0;
 558     fh1DCAInLambda[i] = 0;
 559     fh1DCAInALambda[i] = 0;
 560     fh1DCAOutK0s[i] = 0;
 561     fh1DCAOutLambda[i] = 0;
 562     fh1DCAOutALambda[i] = 0;
 563
 564     fh1PtJet[i] = 0;
 565     fh1EtaJet[i] = 0;
 566     fh2EtaPtJet[i] = 0;
 567     fh1PhiJet[i] = 0;
 568     fh1PtJetTrackLeading[i] = 0;
 569     fh1NJetPerEvent[i] = 0;
 570     fh2EtaPhiRndCone[i] = 0;
 571     fh2EtaPhiMedCone[i] = 0;
 572
 573     fh1VtxZ[i] = 0;
 574     fh2VtxXY[i] = 0;
 575   }
 576   // Define input and output slots here
 577   // Input slot #0 works with a TChain
 578   DefineInput(0, TChain::Class());
 579   // Output slot #0 id reserved by the base class for AOD
 580   // Output slot #1 writes into a TList container
 581   DefineOutput(1, TList::Class());
 582   DefineOutput(2, TList::Class());
 583   DefineOutput(3, TList::Class());
 584   DefineOutput(4, TList::Class());
 585   DefineOutput(5, TTree::Class());
 586 }
 587
 588 AliAnalysisTaskV0sInJetsEmcal::~AliAnalysisTaskV0sInJetsEmcal()
 589 {
 590   delete fRandom;
 591   fRandom = 0;
 592 }
 593
 594 void AliAnalysisTaskV0sInJetsEmcal::ExecOnce()
 595 {
 596   AliAnalysisTaskEmcalJet::ExecOnce();
 597 //  printf("AliAnalysisTaskV0sInJetsEmcal: ExecOnce\n");
 598
 599   if(fJetsCont && fJetsCont->GetArray() == 0)
 600     fJetsCont = 0;
 601   if(fJetsBgCont && fJetsBgCont->GetArray() == 0)
 602     fJetsBgCont = 0;
 603   if(fTracksCont && fTracksCont->GetArray() == 0)
 604     fTracksCont = 0;
 605   if(fCaloClustersCont && fCaloClustersCont->GetArray() == 0)
 606     fCaloClustersCont = 0;
 607 }
 608
 609 Bool_t AliAnalysisTaskV0sInJetsEmcal::Run()
 610 {
 611   // Run analysis code here, if needed. It will be executed before FillHistograms().
 612 //  printf("AliAnalysisTaskV0sInJetsEmcal: Run\n");
 613   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 614 }
 615
 616 void AliAnalysisTaskV0sInJetsEmcal::UserCreateOutputObjects()
 617 {
 618   // Called once
 619
 620   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 621 //  printf("AliAnalysisTaskV0sInJetsEmcal: UserCreateOutputObjects\n");
 622
 623   fJetsCont = GetJetContainer(0);
 624   fJetsBgCont = GetJetContainer(1);
 625   if(fJetsCont) //get particles and clusters connected to jets
 626   {
 627     fTracksCont = fJetsCont->GetParticleContainer();
 628     fCaloClustersCont = fJetsCont->GetClusterContainer();
 629   }
 630   else //no jets, just analysis tracks and clusters
 631   {
 632     fTracksCont = GetParticleContainer(0);
 633     fCaloClustersCont = GetClusterContainer(0);
 634   }
 635   fTracksCont->SetClassName("AliVTrack");
 636   fCaloClustersCont->SetClassName("AliAODCaloCluster");
 637
 638   // Initialise random-number generator
 639   fRandom = new TRandom3(0);
 640
 641   // Create histograms
 642
 643   fOutputListStd = new TList();
 644   fOutputListStd->SetOwner();
 645   fOutputListQA = new TList();
 646   fOutputListQA->SetOwner();
 647   fOutputListCuts = new TList();
 648   fOutputListCuts->SetOwner();
 649   fOutputListMC = new TList();
 650   fOutputListMC->SetOwner();
 651
 652   // event categories
 653   const Int_t iNCategEvent = 6;
 654   TString categEvent[iNCategEvent] = {"coll. candid.", "AOD OK", "vtx & cent", "with V0", "with jets", "jet selection"};
 655   // labels for stages of V0 selection
 656   TString categV0[fgkiNCategV0] = {"all"/*0*/, "mass range"/*1*/, "rec. method"/*2*/, "tracks TPC"/*3*/, "track pt"/*4*/, "DCA prim v"/*5*/, "DCA daughters"/*6*/, "CPA"/*7*/, "volume"/*8*/, "track #it{#eta}"/*9*/, "V0 #it{y} & #it{#eta}"/*10*/, "lifetime"/*11*/, "PID"/*12*/, "Arm.-Pod."/*13*/, "inclusive"/*14*/, "in jet event"/*15*/, "in jet"/*16*/};
 657
 658   fh1EventCounterCut = new TH1D("fh1EventCounterCut", "Number of events after filtering;selection filter;counts", iNCategEvent, 0, iNCategEvent);
 659   for(Int_t i = 0; i < iNCategEvent; i++)
 660     fh1EventCounterCut->GetXaxis()->SetBinLabel(i + 1, categEvent[i].Data());
 661   fh1EventCent2 = new TH1D("fh1EventCent2", "Number of events vs centrality;centrality;counts", 100, 0, 100);
 662   fh1EventCent2Jets = new TH1D("fh1EventCent2Jets", "Number of sel.-jet events vs centrality;centrality;counts", 100, 0, 100);
 663   fh1EventCent2NoJets = new TH1D("fh1EventCent2NoJets", "Number of no-jet events vs centrality;centrality;counts", 100, 0, 100);
 664   fh2EventCentTracks = new TH2D("fh2EventCentTracks", "Number of tracks vs centrality;centrality;tracks;counts", 100, 0, 100, 150, 0, 15e3);
 665   fh1EventCent = new TH1D("fh1EventCent", "Number of events in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
 666   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 667     fh1EventCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
 668   fh1NRndConeCent = new TH1D("fh1NRndConeCent", "Number of rnd. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
 669   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 670     fh1NRndConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
 671   fh1NMedConeCent = new TH1D("fh1NMedConeCent", "Number of med.-cl. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
 672   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 673     fh1NMedConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
 674   fh1AreaExcluded = new TH1D("fh1AreaExcluded", "Area of excluded cones in centrality bins;centrality;area", fgkiNBinsCent, 0, fgkiNBinsCent);
 675   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 676     fh1AreaExcluded->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
 677   fOutputListStd->Add(fh1EventCounterCut);
 678   fOutputListStd->Add(fh1EventCent);
 679   fOutputListStd->Add(fh1EventCent2);
 680   fOutputListStd->Add(fh1EventCent2Jets);
 681   fOutputListStd->Add(fh1EventCent2NoJets);
 682   fOutputListStd->Add(fh1NRndConeCent);
 683   fOutputListStd->Add(fh1NMedConeCent);
 684   fOutputListStd->Add(fh1AreaExcluded);
 685   fOutputListStd->Add(fh2EventCentTracks);
 686
 687   fh1V0CandPerEvent = new TH1D("fh1V0CandPerEvent", "Number of all V0 candidates per event;candidates;events", 1000, 0, 1000);
 688   fOutputListStd->Add(fh1V0CandPerEvent);
 689
 690   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 691   {
 692     fh1EventCounterCutCent[i] = new TH1D(Form("fh1EventCounterCutCent_%d", i), Form("Number of events after filtering, cent %s;selection filter;counts", GetCentBinLabel(i).Data()), iNCategEvent, 0, iNCategEvent);
 693     for(Int_t j = 0; j < iNCategEvent; j++)
 694       fh1EventCounterCutCent[i]->GetXaxis()->SetBinLabel(j + 1, categEvent[j].Data());
 695     fh1V0CandPerEventCentK0s[i] = new TH1D(Form("fh1V0CandPerEventCentK0s_%d", i), Form("Number of selected K0s candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
 696     fh1V0CandPerEventCentLambda[i] = new TH1D(Form("fh1V0CandPerEventCentLambda_%d", i), Form("Number of selected Lambda candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
 697     fh1V0CandPerEventCentALambda[i] = new TH1D(Form("fh1V0CandPerEventCentALambda_%d", i), Form("Number of selected ALambda candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
 698     fh1V0CounterCentK0s[i] = new TH1D(Form("fh1V0CounterCentK0s_%d", i), Form("Number of K0s candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
 699     fh1V0CounterCentLambda[i] = new TH1D(Form("fh1V0CounterCentLambda_%d", i), Form("Number of Lambda candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
 700     fh1V0CounterCentALambda[i] = new TH1D(Form("fh1V0CounterCentALambda_%d", i), Form("Number of ALambda candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
 701     for(Int_t j = 0; j < fgkiNCategV0; j++)
 702     {
 703       fh1V0CounterCentK0s[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
 704       fh1V0CounterCentLambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
 705       fh1V0CounterCentALambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
 706     }
 707     fOutputListStd->Add(fh1EventCounterCutCent[i]);
 708     fOutputListStd->Add(fh1V0CandPerEventCentK0s[i]);
 709     fOutputListStd->Add(fh1V0CandPerEventCentLambda[i]);
 710     fOutputListStd->Add(fh1V0CandPerEventCentALambda[i]);
 711     fOutputListStd->Add(fh1V0CounterCentK0s[i]);
 712     fOutputListStd->Add(fh1V0CounterCentLambda[i]);
 713     fOutputListStd->Add(fh1V0CounterCentALambda[i]);
 714   }
 715   // pt binning for V0 and jets
 716   Int_t iNBinsPtV0 = fgkiNBinsPtV0Init;
 717   Double_t dPtV0Min = fgkdBinsPtV0[0];
 718   Double_t dPtV0Max = fgkdBinsPtV0[fgkiNBinsPtV0];
 719   Int_t iNJetPtBins = fgkiNBinsPtJetInit;
 720   Double_t dJetPtMin = fgkdBinsPtJet[0];
 721   Double_t dJetPtMax = fgkdBinsPtJet[fgkiNBinsPtJet];
 722
 723   fh2CCK0s = new TH2D("fh2CCK0s", "K0s candidates in Lambda peak", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
 724   fh2CCLambda = new TH2D("fh2CCLambda", "Lambda candidates in K0s peak", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
 725   Int_t binsCorrel[3] = {fgkiNBinsMassK0s, fgkiNBinsMassLambda, iNBinsPtV0};
 726   Double_t xminCorrel[3] = {fgkdMassK0sMin, fgkdMassLambdaMin, dPtV0Min};
 727   Double_t xmaxCorrel[3] = {fgkdMassK0sMax, fgkdMassLambdaMax, dPtV0Max};
 728 //  Int_t binsCorrel[3] = {200, 200, iNBinsPtV0};
 729 //  Double_t xminCorrel[3] = {0, 0, dPtV0Min};
 730 //  Double_t xmaxCorrel[3] = {2, 2, dPtV0Max};
 731   fh3CCMassCorrelBoth = new THnSparseD("fh3CCMassCorrelBoth", "Mass correlation: K0S && Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
 732   fh3CCMassCorrelKNotL = new THnSparseD("fh3CCMassCorrelKNotL", "Mass correlation: K0S, not Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
 733   fh3CCMassCorrelLNotK = new THnSparseD("fh3CCMassCorrelLNotK", "Mass correlation: Lambda, not K0S;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
 734   fOutputListQA->Add(fh2CCK0s);
 735   fOutputListQA->Add(fh2CCLambda);
 736   fOutputListQA->Add(fh3CCMassCorrelBoth);
 737   fOutputListQA->Add(fh3CCMassCorrelKNotL);
 738   fOutputListQA->Add(fh3CCMassCorrelLNotK);
 739
 740   Double_t dStepEtaV0 = 0.025;
 741   Double_t dRangeEtaV0Max = 0.8;
 742   const Int_t iNBinsEtaV0 = 2 * Int_t(dRangeEtaV0Max / dStepEtaV0);
 743   // inclusive
 744   const Int_t iNDimIncl = 3;
 745   Int_t binsKIncl[iNDimIncl] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
 746   Double_t xminKIncl[iNDimIncl] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
 747   Double_t xmaxKIncl[iNDimIncl] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
 748   Int_t binsLIncl[iNDimIncl] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
 749   Double_t xminLIncl[iNDimIncl] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
 750   Double_t xmaxLIncl[iNDimIncl] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
 751   // binning in jets
 752   const Int_t iNDimInJC = 4;
 753   Int_t binsKInJC[iNDimInJC] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
 754   Double_t xminKInJC[iNDimInJC] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
 755   Double_t xmaxKInJC[iNDimInJC] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
 756   Int_t binsLInJC[iNDimInJC] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
 757   Double_t xminLInJC[iNDimInJC] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
 758   Double_t xmaxLInJC[iNDimInJC] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
 759
 760   // binning eff inclusive vs eta-pT
 761   Double_t dStepDeltaEta = 0.1;
 762   Double_t dRangeDeltaEtaMax = 0.5;
 763   const Int_t iNBinsDeltaEta = 2 * Int_t(dRangeDeltaEtaMax / dStepDeltaEta);
 764   Int_t binsEtaK[3] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
 765   Double_t xminEtaK[3] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
 766   Double_t xmaxEtaK[3] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
 767   Int_t binsEtaL[3] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
 768   Double_t xminEtaL[3] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
 769   Double_t xmaxEtaL[3] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
 770   // binning eff in jets vs eta-pT
 771   // associated
 772   Int_t binsEtaKInRec[5] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
 773   Double_t xminEtaKInRec[5] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
 774   Double_t xmaxEtaKInRec[5] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
 775   Int_t binsEtaLInRec[5] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
 776   Double_t xminEtaLInRec[5] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
 777   Double_t xmaxEtaLInRec[5] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
 778   // generated
 779   Int_t binsEtaInGen[4] = {iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
 780   Double_t xminEtaInGen[4] = {dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
 781   Double_t xmaxEtaInGen[4] = {dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
 782   // daughter eta: charge-etaD-ptD-etaV0-ptV0-ptJet
 783   const Int_t iNDimEtaD = 6;
 784   Int_t binsEtaDaughter[iNDimEtaD] = {2, 20, iNBinsPtV0, iNBinsEtaV0, iNBinsPtV0, iNJetPtBins};
 785   Double_t xminEtaDaughter[iNDimEtaD] = {0, -1, dPtV0Min, -dRangeEtaV0Max, dPtV0Min, dJetPtMin};
 786   Double_t xmaxEtaDaughter[iNDimEtaD] = {2, 1, dPtV0Max, dRangeEtaV0Max, dPtV0Max, dJetPtMax};
 787
 788   for(Int_t i = 0; i < fgkiNBinsCent; i++)
 789   {
 790     fh1V0InvMassK0sCent[i] = new TH1D(Form("fh1V0InvMassK0sCent_%d", i), Form("K0s: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
 791     fh1V0InvMassLambdaCent[i] = new TH1D(Form("fh1V0InvMassLambdaCent_%d", i), Form("Lambda: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
 792     fh1V0InvMassALambdaCent[i] = new TH1D(Form("fh1V0InvMassALambdaCent_%d", i), Form("ALambda: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
 793     fOutputListStd->Add(fh1V0InvMassK0sCent[i]);
 794     fOutputListStd->Add(fh1V0InvMassLambdaCent[i]);
 795     fOutputListStd->Add(fh1V0InvMassALambdaCent[i]);
 796     // Inclusive
 797     fhnV0InclusiveK0s[i] = new THnSparseD(Form("fhnV0InclusiveK0s_C%d", i), "K0s: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
 798     fhnV0InclusiveLambda[i] = new THnSparseD(Form("fhnV0InclusiveLambda_C%d", i), "Lambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 799     fhnV0InclusiveALambda[i] = new THnSparseD(Form("fhnV0InclusiveALambda_C%d", i), "ALambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 800     fOutputListStd->Add(fhnV0InclusiveK0s[i]);
 801     fOutputListStd->Add(fhnV0InclusiveLambda[i]);
 802     fOutputListStd->Add(fhnV0InclusiveALambda[i]);
 803     // In cones
 804     fhnV0InJetK0s[i] = new THnSparseD(Form("fhnV0InJetK0s_%d", i), Form("K0s: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
 805     fOutputListStd->Add(fhnV0InJetK0s[i]);
 806     fhnV0InPerpK0s[i] = new THnSparseD(Form("fhnV0InPerpK0s_%d", i), Form("K0s: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
 807     fOutputListStd->Add(fhnV0InPerpK0s[i]);
 808     fhnV0InRndK0s[i] = new THnSparseD(Form("fhnV0InRndK0s_%d", i), Form("K0s: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
 809     fOutputListStd->Add(fhnV0InRndK0s[i]);
 810     fhnV0InMedK0s[i] = new THnSparseD(Form("fhnV0InMedK0s_%d", i), Form("K0s: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
 811     fOutputListStd->Add(fhnV0InMedK0s[i]);
 812     fhnV0OutJetK0s[i] = new THnSparseD(Form("fhnV0OutJetK0s_%d", i), Form("K0s: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
 813     fOutputListStd->Add(fhnV0OutJetK0s[i]);
 814     fhnV0NoJetK0s[i] = new THnSparseD(Form("fhnV0NoJetK0s_%d", i), Form("K0s: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
 815     fOutputListStd->Add(fhnV0NoJetK0s[i]);
 816     fhnV0InJetLambda[i] = new THnSparseD(Form("fhnV0InJetLambda_%d", i), Form("Lambda: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 817     fOutputListStd->Add(fhnV0InJetLambda[i]);
 818     fhnV0InPerpLambda[i] = new THnSparseD(Form("fhnV0InPerpLambda_%d", i), Form("Lambda: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 819     fOutputListStd->Add(fhnV0InPerpLambda[i]);
 820     fhnV0InRndLambda[i] = new THnSparseD(Form("fhnV0InRndLambda_%d", i), Form("Lambda: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 821     fOutputListStd->Add(fhnV0InRndLambda[i]);
 822     fhnV0InMedLambda[i] = new THnSparseD(Form("fhnV0InMedLambda_%d", i), Form("Lambda: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 823     fOutputListStd->Add(fhnV0InMedLambda[i]);
 824     fhnV0OutJetLambda[i] = new THnSparseD(Form("fhnV0OutJetLambda_%d", i), Form("Lambda: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 825     fOutputListStd->Add(fhnV0OutJetLambda[i]);
 826     fhnV0NoJetLambda[i] = new THnSparseD(Form("fhnV0NoJetLambda_%d", i), Form("Lambda: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 827     fOutputListStd->Add(fhnV0NoJetLambda[i]);
 828     fhnV0InJetALambda[i] = new THnSparseD(Form("fhnV0InJetALambda_%d", i), Form("ALambda: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 829     fOutputListStd->Add(fhnV0InJetALambda[i]);
 830     fhnV0InPerpALambda[i] = new THnSparseD(Form("fhnV0InPerpALambda_%d", i), Form("ALambda: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 831     fOutputListStd->Add(fhnV0InPerpALambda[i]);
 832     fhnV0InRndALambda[i] = new THnSparseD(Form("fhnV0InRndALambda_%d", i), Form("ALambda: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 833     fOutputListStd->Add(fhnV0InRndALambda[i]);
 834     fhnV0InMedALambda[i] = new THnSparseD(Form("fhnV0InMedALambda_%d", i), Form("ALambda: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 835     fOutputListStd->Add(fhnV0InMedALambda[i]);
 836     fhnV0OutJetALambda[i] = new THnSparseD(Form("fhnV0OutJetALambda_%d", i), Form("ALambda: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 837     fOutputListStd->Add(fhnV0OutJetALambda[i]);
 838     fhnV0NoJetALambda[i] = new THnSparseD(Form("fhnV0NoJetALambda_%d", i), Form("ALambda: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
 839     fOutputListStd->Add(fhnV0NoJetALambda[i]);
 840
 841     fh2V0PtJetAngleK0s[i] = new TH2D(Form("fh2V0PtJetAngleK0s_%d", i), Form("K0s: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
 842     fOutputListStd->Add(fh2V0PtJetAngleK0s[i]);
 843     fh2V0PtJetAngleLambda[i] = new TH2D(Form("fh2V0PtJetAngleLambda_%d", i), Form("Lambda: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
 844     fOutputListStd->Add(fh2V0PtJetAngleLambda[i]);
 845     fh2V0PtJetAngleALambda[i] = new TH2D(Form("fh2V0PtJetAngleALambda_%d", i), Form("ALambda: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
 846     fOutputListStd->Add(fh2V0PtJetAngleALambda[i]);
 847
 848     fh1DCAInK0s[i] = new TH1D(Form("fh1DCAInK0s_%d", i), Form("K0s in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 849     fOutputListQA->Add(fh1DCAInK0s[i]);
 850     fh1DCAInLambda[i] = new TH1D(Form("fh1DCAInLambda_%d", i), Form("Lambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 851     fOutputListQA->Add(fh1DCAInLambda[i]);
 852     fh1DCAInALambda[i] = new TH1D(Form("fh1DCAInALambda_%d", i), Form("ALambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 853     fOutputListQA->Add(fh1DCAInALambda[i]);
 854
 855     fh1DCAOutK0s[i] = new TH1D(Form("fh1DCAOutK0s_%d", i), Form("K0s outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 856     fOutputListQA->Add(fh1DCAOutK0s[i]);
 857     fh1DCAOutLambda[i] = new TH1D(Form("fh1DCAOutLambda_%d", i), Form("Lambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 858     fOutputListQA->Add(fh1DCAOutLambda[i]);
 859     fh1DCAOutALambda[i] = new TH1D(Form("fh1DCAOutALambda_%d", i), Form("ALambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
 860     fOutputListQA->Add(fh1DCAOutALambda[i]);
 861
 862     // jet histograms
 863     fh1PtJet[i] = new TH1D(Form("fh1PtJet_%d", i), Form("Jet pt spectrum, cent: %s;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax);
 864     fOutputListStd->Add(fh1PtJet[i]);
 865     fh1EtaJet[i] = new TH1D(Form("fh1EtaJet_%d", i), Form("Jet eta spectrum, cent: %s;#it{#eta} jet", GetCentBinLabel(i).Data()), 80, -1., 1.);
 866     fOutputListStd->Add(fh1EtaJet[i]);
 867     fh2EtaPtJet[i] = new TH2D(Form("fh2EtaPtJet_%d", i), Form("Jet eta vs pT spectrum, cent: %s;#it{#eta} jet;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), 80, -1., 1., iNJetPtBins, dJetPtMin, dJetPtMax);
 868     fOutputListStd->Add(fh2EtaPtJet[i]);
 869     fh2EtaPhiRndCone[i] = new TH2D(Form("fh2EtaPhiRndCone_%d", i), Form("Rnd. cones: eta vs phi, cent: %s;#it{#eta} cone;#it{#phi} cone", GetCentBinLabel(i).Data()), 80, -1., 1., 100, 0., TMath::TwoPi());
 870     fOutputListStd->Add(fh2EtaPhiRndCone[i]);
 871     fh2EtaPhiMedCone[i] = new TH2D(Form("fh2EtaPhiMedCone_%d", i), Form("Med.-cl. cones: eta vs phi, cent: %s;#it{#eta} cone;#it{#phi} cone", GetCentBinLabel(i).Data()), 80, -1., 1., 100, 0., TMath::TwoPi());
 872     fOutputListStd->Add(fh2EtaPhiMedCone[i]);
 873     fh1PhiJet[i] = new TH1D(Form("fh1PhiJet_%d", i), Form("Jet phi spectrum, cent: %s;#it{#phi} jet", GetCentBinLabel(i).Data()), 100, 0., TMath::TwoPi());
 874     fOutputListStd->Add(fh1PhiJet[i]);
 875     fh1PtJetTrackLeading[i] = new TH1D(Form("fh1PtJetTrackLeading_%d", i), Form("Leading track pt spectrum, cent: %s;#it{p}_{T} leading track (GeV/#it{c})", GetCentBinLabel(i).Data()), 200, 0., 20);
 876     fOutputListStd->Add(fh1PtJetTrackLeading[i]);
 877     fh1NJetPerEvent[i] = new TH1D(Form("fh1NJetPerEvent_%d", i), Form("Number of selected jets per event, cent: %s;# jets;# events", GetCentBinLabel(i).Data()), 100, 0., 100.);
 878     fOutputListStd->Add(fh1NJetPerEvent[i]);
 879     // event histograms
 880     fh1VtxZ[i] = new TH1D(Form("fh1VtxZ_%d", i), Form("#it{z} coordinate of the primary vertex, cent: %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 150, -1.5 * fdCutVertexZ, 1.5 * fdCutVertexZ);
 881     fOutputListQA->Add(fh1VtxZ[i]);
 882     fh2VtxXY[i] = new TH2D(Form("fh2VtxXY_%d", i), Form("#it{xy} coordinate of the primary vertex, cent: %s;#it{x} (cm);#it{y} (cm)", GetCentBinLabel(i).Data()), 200, -0.2, 0.2, 500, -0.5, 0.5);
 883     fOutputListQA->Add(fh2VtxXY[i]);
 884     // fOutputListStd->Add([i]);
 885     if(fbMCAnalysis)
 886     {
 887       // inclusive pt
 888       fh1V0K0sPtMCGen[i] = new TH1D(Form("fh1V0K0sPtMCGen_%d", i), Form("MC K0s generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 889       fOutputListMC->Add(fh1V0K0sPtMCGen[i]);
 890       fh2V0K0sPtMassMCRec[i] = new TH2D(Form("fh2V0K0sPtMassMCRec_%d", i), Form("MC K0s associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
 891       fOutputListMC->Add(fh2V0K0sPtMassMCRec[i]);
 892       fh1V0K0sPtMCRecFalse[i] = new TH1D(Form("fh1V0K0sPtMCRecFalse_%d", i), Form("MC K0s false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 893       fOutputListMC->Add(fh1V0K0sPtMCRecFalse[i]);
 894       // inclusive pt-eta
 895       fh2V0K0sEtaPtMCGen[i] = new TH2D(Form("fh2V0K0sEtaPtMCGen_%d", i), Form("MC K0s generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
 896       fOutputListMC->Add(fh2V0K0sEtaPtMCGen[i]);
 897       fh3V0K0sEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0K0sEtaPtMassMCRec_%d", i), Form("MC K0s associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaK, xminEtaK, xmaxEtaK);
 898       fOutputListMC->Add(fh3V0K0sEtaPtMassMCRec[i]);
 899       // in jet pt
 900       fh2V0K0sInJetPtMCGen[i] = new TH2D(Form("fh2V0K0sInJetPtMCGen_%d", i), Form("MC K0s in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
 901       fOutputListMC->Add(fh2V0K0sInJetPtMCGen[i]);
 902       fh3V0K0sInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0K0sInJetPtMassMCRec_%d", i), Form("MC K0s in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
 903       fOutputListMC->Add(fh3V0K0sInJetPtMassMCRec[i]);
 904       // in jet pt-eta
 905       fh3V0K0sInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0K0sInJetEtaPtMCGen_%d", i), Form("MC K0s generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
 906       fOutputListMC->Add(fh3V0K0sInJetEtaPtMCGen[i]);
 907       fh4V0K0sInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0K0sInJetEtaPtMassMCRec_%d", i), Form("MC K0s associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaKInRec, xminEtaKInRec, xmaxEtaKInRec);
 908       fOutputListMC->Add(fh4V0K0sInJetEtaPtMassMCRec[i]);
 909
 910       fh2V0K0sMCResolMPt[i] = new TH2D(Form("fh2V0K0sMCResolMPt_%d", i), Form("MC K0s associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
 911       fOutputListMC->Add(fh2V0K0sMCResolMPt[i]);
 912       fh2V0K0sMCPtGenPtRec[i] = new TH2D(Form("fh2V0K0sMCPtGenPtRec_%d", i), Form("MC K0s associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
 913       fOutputListMC->Add(fh2V0K0sMCPtGenPtRec[i]);
 914
 915       // inclusive pt
 916       fh1V0LambdaPtMCGen[i] = new TH1D(Form("fh1V0LambdaPtMCGen_%d", i), Form("MC Lambda generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 917       fOutputListMC->Add(fh1V0LambdaPtMCGen[i]);
 918       fh2V0LambdaPtMassMCRec[i] = new TH2D(Form("fh2V0LambdaPtMassMCRec_%d", i), Form("MC Lambda associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
 919       fOutputListMC->Add(fh2V0LambdaPtMassMCRec[i]);
 920       fh1V0LambdaPtMCRecFalse[i] = new TH1D(Form("fh1V0LambdaPtMCRecFalse_%d", i), Form("MC Lambda false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 921       fOutputListMC->Add(fh1V0LambdaPtMCRecFalse[i]);
 922       // inclusive pt-eta
 923       fh2V0LambdaEtaPtMCGen[i] = new TH2D(Form("fh2V0LambdaEtaPtMCGen_%d", i), Form("MC Lambda generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
 924       fOutputListMC->Add(fh2V0LambdaEtaPtMCGen[i]);
 925       fh3V0LambdaEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0LambdaEtaPtMassMCRec_%d", i), Form("MC Lambda associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaL, xminEtaL, xmaxEtaL);
 926       fOutputListMC->Add(fh3V0LambdaEtaPtMassMCRec[i]);
 927       // in jet pt
 928       fh2V0LambdaInJetPtMCGen[i] = new TH2D(Form("fh2V0LambdaInJetPtMCGen_%d", i), Form("MC Lambda in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
 929       fOutputListMC->Add(fh2V0LambdaInJetPtMCGen[i]);
 930       fh3V0LambdaInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0LambdaInJetPtMassMCRec_%d", i), Form("MC Lambda in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 931       fOutputListMC->Add(fh3V0LambdaInJetPtMassMCRec[i]);
 932       // in jet pt-eta
 933       fh3V0LambdaInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0LambdaInJetEtaPtMCGen_%d", i), Form("MC Lambda generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
 934       fOutputListMC->Add(fh3V0LambdaInJetEtaPtMCGen[i]);
 935       fh4V0LambdaInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0LambdaInJetEtaPtMassMCRec_%d", i), Form("MC Lambda associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaLInRec, xminEtaLInRec, xmaxEtaLInRec);
 936       fOutputListMC->Add(fh4V0LambdaInJetEtaPtMassMCRec[i]);
 937
 938       fh2V0LambdaMCResolMPt[i] = new TH2D(Form("fh2V0LambdaMCResolMPt_%d", i), Form("MC Lambda associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
 939       fOutputListMC->Add(fh2V0LambdaMCResolMPt[i]);
 940       fh2V0LambdaMCPtGenPtRec[i] = new TH2D(Form("fh2V0LambdaMCPtGenPtRec_%d", i), Form("MC Lambda associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
 941       fOutputListMC->Add(fh2V0LambdaMCPtGenPtRec[i]);
 942
 943       // inclusive pt
 944       fh1V0ALambdaPtMCGen[i] = new TH1D(Form("fh1V0ALambdaPtMCGen_%d", i), Form("MC ALambda generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 945       fOutputListMC->Add(fh1V0ALambdaPtMCGen[i]);
 946       fh2V0ALambdaPtMassMCRec[i] = new TH2D(Form("fh2V0ALambdaPtMassMCRec_%d", i), Form("MC ALambda associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
 947       fOutputListMC->Add(fh2V0ALambdaPtMassMCRec[i]);
 948       fh1V0ALambdaPtMCRecFalse[i] = new TH1D(Form("fh1V0ALambdaPtMCRecFalse_%d", i), Form("MC ALambda false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
 949       fOutputListMC->Add(fh1V0ALambdaPtMCRecFalse[i]);
 950       // inclusive pt-eta
 951       fh2V0ALambdaEtaPtMCGen[i] = new TH2D(Form("fh2V0ALambdaEtaPtMCGen_%d", i), Form("MC ALambda generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
 952       fOutputListMC->Add(fh2V0ALambdaEtaPtMCGen[i]);
 953       fh3V0ALambdaEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0ALambdaEtaPtMassMCRec_%d", i), Form("MC ALambda associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaL, xminEtaL, xmaxEtaL);
 954       fOutputListMC->Add(fh3V0ALambdaEtaPtMassMCRec[i]);
 955       // in jet pt
 956       fh2V0ALambdaInJetPtMCGen[i] = new TH2D(Form("fh2V0ALambdaInJetPtMCGen_%d", i), Form("MC ALambda in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
 957       fOutputListMC->Add(fh2V0ALambdaInJetPtMCGen[i]);
 958       fh3V0ALambdaInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0ALambdaInJetPtMassMCRec_%d", i), Form("MC ALambda in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
 959       fOutputListMC->Add(fh3V0ALambdaInJetPtMassMCRec[i]);
 960       // in jet pt-eta
 961       fh3V0ALambdaInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0ALambdaInJetEtaPtMCGen_%d", i), Form("MC ALambda generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
 962       fOutputListMC->Add(fh3V0ALambdaInJetEtaPtMCGen[i]);
 963       fh4V0ALambdaInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0ALambdaInJetEtaPtMassMCRec_%d", i), Form("MC ALambda associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaLInRec, xminEtaLInRec, xmaxEtaLInRec);
 964       fOutputListMC->Add(fh4V0ALambdaInJetEtaPtMassMCRec[i]);
 965
 966       fh2V0ALambdaMCResolMPt[i] = new TH2D(Form("fh2V0ALambdaMCResolMPt_%d", i), Form("MC ALambda associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
 967       fOutputListMC->Add(fh2V0ALambdaMCResolMPt[i]);
 968       fh2V0ALambdaMCPtGenPtRec[i] = new TH2D(Form("fh2V0ALambdaMCPtGenPtRec_%d", i), Form("MC ALambda associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
 969       fOutputListMC->Add(fh2V0ALambdaMCPtGenPtRec[i]);
 970
 971       Int_t iNBinsPtXi = 80;
 972       Double_t dPtXiMin = 0;
 973       Double_t dPtXiMax = 8;
 974       const Int_t iNDimFD = 3;
 975       Int_t binsFD[iNDimFD] = {iNBinsPtV0, iNBinsPtXi, iNJetPtBins};
 976       Double_t xminFD[iNDimFD] = {dPtV0Min, dPtXiMin, dJetPtMin};
 977       Double_t xmaxFD[iNDimFD] = {dPtV0Max, dPtXiMax, dJetPtMax};
 978       fhnV0LambdaInclMCFD[i] = new THnSparseD(Form("fhnV0LambdaInclMCFD_%d", i), Form("MC Lambda associated, inclusive, from Xi: pt-pt, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 979       fOutputListMC->Add(fhnV0LambdaInclMCFD[i]);
 980       fhnV0LambdaInJetsMCFD[i] = new THnSparseD(Form("fhnV0LambdaInJetsMCFD_%d", i), Form("MC Lambda associated, in JC, from Xi: pt-pt-ptJet, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 981       fOutputListMC->Add(fhnV0LambdaInJetsMCFD[i]);
 982       fhnV0LambdaBulkMCFD[i] = new THnSparseD(Form("fhnV0LambdaBulkMCFD_%d", i), Form("MC Lambda associated, in no jet events, from Xi: pt-pt, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 983       fOutputListMC->Add(fhnV0LambdaBulkMCFD[i]);
 984       fh1V0XiPtMCGen[i] = new TH1D(Form("fh1V0XiPtMCGen_%d", i), Form("MC Xi^{-} generated: Pt spectrum, cent %s;#it{p}_{T}^{#Xi^{-},gen.} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtXi, dPtXiMin, dPtXiMax);
 985       fOutputListMC->Add(fh1V0XiPtMCGen[i]);
 986       fhnV0ALambdaInclMCFD[i] = new THnSparseD(Form("fhnV0ALambdaInclMCFD_%d", i), Form("MC ALambda associated, from AXi: pt-pt, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 987       fOutputListMC->Add(fhnV0ALambdaInclMCFD[i]);
 988       fhnV0ALambdaInJetsMCFD[i] = new THnSparseD(Form("fhnV0ALambdaInJetsMCFD_%d", i), Form("MC ALambda associated, in JC, from AXi: pt-pt-ptJet, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 989       fOutputListMC->Add(fhnV0ALambdaInJetsMCFD[i]);
 990       fhnV0ALambdaBulkMCFD[i] = new THnSparseD(Form("fhnV0ALambdaBulkMCFD_%d", i), Form("MC ALambda associated, in no jet events, from AXi: pt-pt-ptJet, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
 991       fOutputListMC->Add(fhnV0ALambdaBulkMCFD[i]);
 992       fh1V0AXiPtMCGen[i] = new TH1D(Form("fh1V0AXiPtMCGen_%d", i), Form("MC AXi^{-} generated: Pt spectrum, cent %s;#it{p}_{T}^{A#Xi^{-},gen.} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtXi, dPtXiMin, dPtXiMax);
 993       fOutputListMC->Add(fh1V0AXiPtMCGen[i]);
 994
 995       // daughter eta
 996 //          fhnV0K0sInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
 997       fhnV0K0sInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
 998 //          fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
 999       fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1000 //          fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1001       fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1002 //          fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1003       fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1004 //          fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1005       fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1006 //          fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1007       fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1008
1009 //          fOutputListMC->Add(fhnV0K0sInclDaughterEtaPtPtMCGen[i]);
1010       fOutputListMC->Add(fhnV0K0sInclDaughterEtaPtPtMCRec[i]);
1011 //          fOutputListMC->Add(fhnV0K0sInJetsDaughterEtaPtPtMCGen[i]);
1012       fOutputListMC->Add(fhnV0K0sInJetsDaughterEtaPtPtMCRec[i]);
1013 //          fOutputListMC->Add(fhnV0LambdaInclDaughterEtaPtPtMCGen[i]);
1014       fOutputListMC->Add(fhnV0LambdaInclDaughterEtaPtPtMCRec[i]);
1015 //          fOutputListMC->Add(fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i]);
1016       fOutputListMC->Add(fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i]);
1017 //          fOutputListMC->Add(fhnV0ALambdaInclDaughterEtaPtPtMCGen[i]);
1018       fOutputListMC->Add(fhnV0ALambdaInclDaughterEtaPtPtMCRec[i]);
1019 //          fOutputListMC->Add(fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i]);
1020       fOutputListMC->Add(fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i]);
1021     }
1022   }
1023
1024   // QA Histograms
1025   for(Int_t i = 0; i < fgkiNQAIndeces; i++)
1026   {
1027 //    [i] = new TH1D(Form("%d",i),";;Counts",,,);
1028     fh1QAV0Status[i] = new TH1D(Form("fh1QAV0Status_%d", i), "QA: V0 status", 2, 0, 2);
1029     fh1QAV0TPCRefit[i] = new TH1D(Form("fh1QAV0TPCRefit_%d", i), "QA: TPC refit", 2, 0, 2);
1030     fh1QAV0TPCRows[i] = new TH1D(Form("fh1QAV0TPCRows_%d", i), "QA: TPC Rows", 160, 0, 160);
1031     fh1QAV0TPCFindable[i] = new TH1D(Form("fh1QAV0TPCFindable_%d", i), "QA: TPC Findable", 160, 0, 160);
1032     fh1QAV0TPCRowsFind[i] = new TH1D(Form("fh1QAV0TPCRowsFind_%d", i), "QA: TPC Rows/Findable", 100, 0, 2);
1033     fh1QAV0Eta[i] = new TH1D(Form("fh1QAV0Eta_%d", i), "QA: Daughter Eta", 200, -2, 2);
1034     fh2QAV0EtaRows[i] = new TH2D(Form("fh2QAV0EtaRows_%d", i), "QA: Daughter Eta vs TPC rows;#eta;TPC rows", 200, -2, 2, 160, 0, 160);
1035     fh2QAV0PtRows[i] = new TH2D(Form("fh2QAV0PtRows_%d", i), "QA: Daughter Pt vs TPC rows;pt;TPC rows", 100, 0, 10, 160, 0, 160);
1036     fh2QAV0PhiRows[i] = new TH2D(Form("fh2QAV0PhiRows_%d", i), "QA: Daughter Phi vs TPC rows;#phi;TPC rows", 100, 0, TMath::TwoPi(), 160, 0, 160);
1037     fh2QAV0NClRows[i] = new TH2D(Form("fh2QAV0NClRows_%d", i), "QA: Daughter NCl vs TPC rows;findable clusters;TPC rows", 100, 0, 160, 160, 0, 160);
1038     fh2QAV0EtaNCl[i] = new TH2D(Form("fh2QAV0EtaNCl_%d", i), "QA: Daughter Eta vs NCl;#eta;findable clusters", 200, -2, 2, 160, 0, 160);
1039
1040     fh2QAV0EtaPtK0sPeak[i] = new TH2D(Form("fh2QAV0EtaPtK0sPeak_%d", i), "QA: K0s: Daughter Eta vs V0 pt, peak;track eta;V0 pt", 200, -2, 2, iNBinsPtV0, dPtV0Min, dPtV0Max);
1041     fh2QAV0EtaEtaK0s[i] = new TH2D(Form("fh2QAV0EtaEtaK0s_%d", i), "QA: K0s: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1042     fh2QAV0PhiPhiK0s[i] = new TH2D(Form("fh2QAV0PhiPhiK0s_%d", i), "QA: K0s: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1043     fh1QAV0RapK0s[i] = new TH1D(Form("fh1QAV0RapK0s_%d", i), "QA: K0s: V0 Rapidity", 200, -2, 2);
1044     fh2QAV0PtPtK0sPeak[i] = new TH2D(Form("fh2QAV0PtPtK0sPeak_%d", i), "QA: K0s: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1045
1046     fh2QAV0EtaPtLambdaPeak[i] = new TH2D(Form("fh2QAV0EtaPtLambdaPeak_%d", i), "QA: Lambda: Daughter Eta vs V0 pt, peak;track eta;V0 pt", 200, -2, 2, iNBinsPtV0, dPtV0Min, dPtV0Max);
1047     fh2QAV0EtaEtaLambda[i] = new TH2D(Form("fh2QAV0EtaEtaLambda_%d", i), "QA: Lambda: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1048     fh2QAV0PhiPhiLambda[i] = new TH2D(Form("fh2QAV0PhiPhiLambda_%d", i), "QA: Lambda: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1049     fh1QAV0RapLambda[i] = new TH1D(Form("fh1QAV0RapLambda_%d", i), "QA: Lambda: V0 Rapidity", 200, -2, 2);
1050     fh2QAV0PtPtLambdaPeak[i] = new TH2D(Form("fh2QAV0PtPtLambdaPeak_%d", i), "QA: Lambda: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1051
1052     fh1QAV0Pt[i] = new TH1D(Form("fh1QAV0Pt_%d", i), "QA: Daughter Pt", 100, 0, 5);
1053     fh1QAV0Charge[i] = new TH1D(Form("fh1QAV0Charge_%d", i), "QA: V0 Charge", 3, -1, 2);
1054     fh1QAV0DCAVtx[i] = new TH1D(Form("fh1QAV0DCAVtx_%d", i), "QA: DCA daughters to primary vertex", 100, 0, 10);
1055     fh1QAV0DCAV0[i] = new TH1D(Form("fh1QAV0DCAV0_%d", i), "QA: DCA daughters", 100, 0, 2);
1056     fh1QAV0Cos[i] = new TH1D(Form("fh1QAV0Cos_%d", i), "QA: CPA", 10000, 0.9, 1);
1057     fh1QAV0R[i] = new TH1D(Form("fh1QAV0R_%d", i), "QA: R", 1500, 0, 150);
1058     fh1QACTau2D[i] = new TH1D(Form("fh1QACTau2D_%d", i), "QA: K0s: c#tau 2D;mR/pt#tau", 100, 0, 10);
1059     fh1QACTau3D[i] = new TH1D(Form("fh1QACTau3D_%d", i), "QA: K0s: c#tau 3D;mL/p#tau", 100, 0, 10);
1060
1061     fh2ArmPod[i] = new TH2D(Form("fh2ArmPod_%d", i), "Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1062     fh2ArmPodK0s[i] = new TH2D(Form("fh2ArmPodK0s_%d", i), "K0s: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1063     fh2ArmPodLambda[i] = new TH2D(Form("fh2ArmPodLambda_%d", i), "Lambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1064     fh2ArmPodALambda[i] = new TH2D(Form("fh2ArmPodALambda_%d", i), "ALambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1065
1066     fOutputListQA->Add(fh1QAV0Status[i]);
1067     fOutputListQA->Add(fh1QAV0TPCRefit[i]);
1068     fOutputListQA->Add(fh1QAV0TPCRows[i]);
1069     fOutputListQA->Add(fh1QAV0TPCFindable[i]);
1070     fOutputListQA->Add(fh1QAV0TPCRowsFind[i]);
1071     fOutputListQA->Add(fh1QAV0Eta[i]);
1072     fOutputListQA->Add(fh2QAV0EtaRows[i]);
1073     fOutputListQA->Add(fh2QAV0PtRows[i]);
1074     fOutputListQA->Add(fh2QAV0PhiRows[i]);
1075     fOutputListQA->Add(fh2QAV0NClRows[i]);
1076     fOutputListQA->Add(fh2QAV0EtaNCl[i]);
1077
1078     fOutputListQA->Add(fh2QAV0EtaPtK0sPeak[i]);
1079     fOutputListQA->Add(fh2QAV0EtaEtaK0s[i]);
1080     fOutputListQA->Add(fh2QAV0PhiPhiK0s[i]);
1081     fOutputListQA->Add(fh1QAV0RapK0s[i]);
1082     fOutputListQA->Add(fh2QAV0PtPtK0sPeak[i]);
1083
1084     fOutputListQA->Add(fh2QAV0EtaPtLambdaPeak[i]);
1085     fOutputListQA->Add(fh2QAV0EtaEtaLambda[i]);
1086     fOutputListQA->Add(fh2QAV0PhiPhiLambda[i]);
1087     fOutputListQA->Add(fh1QAV0RapLambda[i]);
1088     fOutputListQA->Add(fh2QAV0PtPtLambdaPeak[i]);
1089
1090     fOutputListQA->Add(fh1QAV0Pt[i]);
1091     fOutputListQA->Add(fh1QAV0Charge[i]);
1092     fOutputListQA->Add(fh1QAV0DCAVtx[i]);
1093     fOutputListQA->Add(fh1QAV0DCAV0[i]);
1094     fOutputListQA->Add(fh1QAV0Cos[i]);
1095     fOutputListQA->Add(fh1QAV0R[i]);
1096     fOutputListQA->Add(fh1QACTau2D[i]);
1097     fOutputListQA->Add(fh1QACTau3D[i]);
1098
1099     fOutputListQA->Add(fh2ArmPod[i]);
1100     fOutputListQA->Add(fh2ArmPodK0s[i]);
1101     fOutputListQA->Add(fh2ArmPodLambda[i]);
1102     fOutputListQA->Add(fh2ArmPodALambda[i]);
1103
1104     /*
1105     fh2CutTPCRowsK0s[i] = new TH2D(Form("fh2CutTPCRowsK0s_%d", i), "Cuts: K0s: TPC Rows vs mass;#it{m}_{inv} (GeV/#it{c}^{2});TPC rows", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 160, 0, 160);
1106     fh2CutTPCRowsLambda[i] = new TH2D(Form("fh2CutTPCRowsLambda_%d", i), "Cuts: Lambda: TPC Rows vs mass;#it{m}_{inv} (GeV/#it{c}^{2});TPC rows", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 160, 0, 160);
1107     fh2CutPtPosK0s[i] = new TH2D(Form("fh2CutPtPosK0s_%d", i), "Cuts: K0s: Pt pos;#it{m}_{inv} (GeV/#it{c}^{2});pt pos", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 5);
1108     fh2CutPtNegK0s[i] = new TH2D(Form("fh2CutPtNegK0s_%d", i), "Cuts: K0s: Pt neg;#it{m}_{inv} (GeV/#it{c}^{2});pt neg", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 5);
1109     fh2CutPtPosLambda[i] = new TH2D(Form("fh2CutPtPosLambda_%d", i), "Cuts: Lambda: Pt pos;#it{m}_{inv} (GeV/#it{c}^{2});pt pos", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 5);
1110     fh2CutPtNegLambda[i] = new TH2D(Form("fh2CutPtNegLambda_%d", i), "Cuts: Lambda: Pt neg;#it{m}_{inv} (GeV/#it{c}^{2});pt neg", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 5);
1111     fh2CutDCAVtx[i] = new TH2D(Form("fh2CutDCAVtx_%d", i), "Cuts: DCA daughters to prim. vtx.;#it{m}_{inv} (GeV/#it{c}^{2});DCA daughter to prim. vtx. (cm)", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1112     fh2CutDCAV0[i] = new TH2D(Form("fh2CutDCAV0_%d", i), "Cuts: DCA daughters;#it{m}_{inv} (GeV/#it{c}^{2});DCA daughters / #sigma_{TPC}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 2);
1113     fh2CutCos[i] = new TH2D(Form("fh2CutCos_%d", i), "Cuts: CPA;#it{m}_{inv} (GeV/#it{c}^{2});CPA", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 10000, 0.9, 1);
1114     fh2CutR[i] = new TH2D(Form("fh2CutR_%d", i), "Cuts: R;#it{m}_{inv} (GeV/#it{c}^{2});R (cm)", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 1500, 0, 150);
1115     fh2CutEtaK0s[i] = new TH2D(Form("fh2CutEtaK0s_%d", i), "Cuts: K0s: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 200, -2, 2);
1116     fh2CutEtaLambda[i] = new TH2D(Form("fh2CutEtaLambda_%d", i), "Cuts: Lambda: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 200, -2, 2);
1117     fh2CutRapK0s[i] = new TH2D(Form("fh2CutRapK0s_%d", i), "Cuts: K0s: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 200, -2, 2);
1118     fh2CutRapLambda[i] = new TH2D(Form("fh2CutRapLambda_%d", i), "Cuts: Lambda: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 200, -2, 2);
1119     fh2CutCTauK0s[i] = new TH2D(Form("fh2CutCTauK0s_%d", i), "Cuts: K0s: #it{c#tau};#it{m}_{inv} (GeV/#it{c}^{2});#it{mL/p#tau}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1120     fh2CutCTauLambda[i] = new TH2D(Form("fh2CutCTauLambda_%d", i), "Cuts: Lambda: #it{c#tau};#it{m}_{inv} (GeV/#it{c}^{2});#it{mL/p#tau}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1121     fh2CutPIDPosK0s[i] = new TH2D(Form("fh2CutPIDPosK0s_%d", i), "Cuts: K0s: PID pos;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1122     fh2CutPIDNegK0s[i] = new TH2D(Form("fh2CutPIDNegK0s_%d", i), "Cuts: K0s: PID neg;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1123     fh2CutPIDPosLambda[i] = new TH2D(Form("fh2CutPIDPosLambda_%d", i), "Cuts: Lambda: PID pos;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1124     fh2CutPIDNegLambda[i] = new TH2D(Form("fh2CutPIDNegLambda_%d", i), "Cuts: Lambda: PID neg;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1125     fh2Tau3DVs2D[i] = new TH2D(Form("fh2Tau3DVs2D_%d", i), "Decay 3D vs 2D;pt;3D/2D", 100, 0, 10, 200, 0.5, 1.5);
1126
1127     fOutputListCuts->Add(fh2CutTPCRowsK0s[i]);
1128     fOutputListCuts->Add(fh2CutTPCRowsLambda[i]);
1129     fOutputListCuts->Add(fh2CutPtPosK0s[i]);
1130     fOutputListCuts->Add(fh2CutPtNegK0s[i]);
1131     fOutputListCuts->Add(fh2CutPtPosLambda[i]);
1132     fOutputListCuts->Add(fh2CutPtNegLambda[i]);
1133     fOutputListCuts->Add(fh2CutDCAVtx[i]);
1134     fOutputListCuts->Add(fh2CutDCAV0[i]);
1135     fOutputListCuts->Add(fh2CutCos[i]);
1136     fOutputListCuts->Add(fh2CutR[i]);
1137     fOutputListCuts->Add(fh2CutEtaK0s[i]);
1138     fOutputListCuts->Add(fh2CutEtaLambda[i]);
1139     fOutputListCuts->Add(fh2CutRapK0s[i]);
1140     fOutputListCuts->Add(fh2CutRapLambda[i]);
1141     fOutputListCuts->Add(fh2CutCTauK0s[i]);
1142     fOutputListCuts->Add(fh2CutCTauLambda[i]);
1143     fOutputListCuts->Add(fh2CutPIDPosK0s[i]);
1144     fOutputListCuts->Add(fh2CutPIDNegK0s[i]);
1145     fOutputListCuts->Add(fh2CutPIDPosLambda[i]);
1146     fOutputListCuts->Add(fh2CutPIDNegLambda[i]);
1147     fOutputListCuts->Add(fh2Tau3DVs2D[i]);
1148     */
1149   }
1150
1151   for(Int_t i = 0; i < fgkiNCategV0; i++)
1152   {
1153     fh1V0InvMassK0sAll[i] = new TH1D(Form("fh1V0InvMassK0sAll_%d", i), Form("K0s: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
1154     fh1V0InvMassLambdaAll[i] = new TH1D(Form("fh1V0InvMassLambdaAll_%d", i), Form("Lambda: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
1155     fh1V0InvMassALambdaAll[i] = new TH1D(Form("fh1V0InvMassALambdaAll_%d", i), Form("ALambda: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
1156     fOutputListStd->Add(fh1V0InvMassK0sAll[i]);
1157     fOutputListStd->Add(fh1V0InvMassLambdaAll[i]);
1158     fOutputListStd->Add(fh1V0InvMassALambdaAll[i]);
1159   }
1160
1161   for(Int_t i = 0; i < fOutputListStd->GetEntries(); ++i)
1162   {
1163     TH1* h1 = dynamic_cast<TH1*>(fOutputListStd->At(i));
1164     if(h1)
1165     {
1166       h1->Sumw2();
1167       continue;
1168     }
1169     THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListStd->At(i));
1170     if(hn) hn->Sumw2();
1171   }
1172   for(Int_t i = 0; i < fOutputListQA->GetEntries(); ++i)
1173   {
1174     TH1* h1 = dynamic_cast<TH1*>(fOutputListQA->At(i));
1175     if(h1)
1176     {
1177       h1->Sumw2();
1178       continue;
1179     }
1180     THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListQA->At(i));
1181     if(hn) hn->Sumw2();
1182   }
1183   for(Int_t i = 0; i < fOutputListCuts->GetEntries(); ++i)
1184   {
1185     TH1* h1 = dynamic_cast<TH1*>(fOutputListCuts->At(i));
1186     if(h1)
1187     {
1188       h1->Sumw2();
1189       continue;
1190     }
1191     THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListCuts->At(i));
1192     if(hn) hn->Sumw2();
1193   }
1194   for(Int_t i = 0; i < fOutputListMC->GetEntries(); ++i)
1195   {
1196     TH1* h1 = dynamic_cast<TH1*>(fOutputListMC->At(i));
1197     if(h1)
1198     {
1199       h1->Sumw2();
1200       continue;
1201     }
1202     THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListMC->At(i));
1203     if(hn) hn->Sumw2();
1204   }
1205
1206   PostData(1, fOutputListStd);
1207   PostData(2, fOutputListQA);
1208   PostData(3, fOutputListCuts);
1209   PostData(4, fOutputListMC);
1210 }
1211
1212 Bool_t AliAnalysisTaskV0sInJetsEmcal::FillHistograms()
1213 {
1214   // Main loop, called for each event
1215   if(fDebug > 5) printf("TaskV0sInJetsEmcal: FillHistograms: Start\n");
1216
1217   if(fDebug > 2) printf("TaskV0sInJetsEmcal: AOD analysis\n");
1218   fh1EventCounterCut->Fill(0); // all available selected events (collision candidates)
1219
1220   if(fDebug > 5) printf("TaskV0sInJetsEmcal: FillHistograms: Loading AOD\n");
1221   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); // input AOD
1222   if(!fAODIn)
1223   {
1224     if(fDebug > 0) printf("TaskV0sInJetsEmcal: No input AOD found\n");
1225     return kFALSE;
1226   }
1227   if(fDebug > 5) printf("TaskV0sInJetsEmcal: FillHistograms: Loading AOD OK\n");
1228
1229   TClonesArray* arrayMC = 0; // array particles in the MC event
1230   AliAODMCHeader* headerMC = 0; // MC header
1231   Int_t iNTracksMC = 0; // number of MC tracks
1232   Double_t dPrimVtxMCX = 0., dPrimVtxMCY = 0., dPrimVtxMCZ = 0.; // position of the MC primary vertex
1233
1234   // Simulation info
1235   if(fbMCAnalysis)
1236   {
1237     arrayMC = (TClonesArray*)fAODIn->FindListObject(AliAODMCParticle::StdBranchName());
1238     if(!arrayMC)
1239     {
1240       if(fDebug > 0) printf("TaskV0sInJetsEmcal: No MC array found\n");
1241       return kFALSE;
1242     }
1243     if(fDebug > 5) printf("TaskV0sInJetsEmcal: MC array found\n");
1244     iNTracksMC = arrayMC->GetEntriesFast();
1245     if(fDebug > 5) printf("TaskV0sInJetsEmcal: There are %d MC tracks in this event\n", iNTracksMC);
1246     headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
1247     if(!headerMC)
1248     {
1249       if(fDebug > 0) printf("TaskV0sInJetsEmcal: No MC header found\n");
1250       return kFALSE;
1251     }
1252     // get position of the MC primary vertex
1253     dPrimVtxMCX = headerMC->GetVtxX();
1254     dPrimVtxMCY = headerMC->GetVtxY();
1255     dPrimVtxMCZ = headerMC->GetVtxZ();
1256   }
1257
1258   // PID Response Task object
1259   AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1260   AliInputEventHandler* inputHandler = (AliInputEventHandler*)mgr->GetInputEventHandler();
1261   AliPIDResponse* fPIDResponse = inputHandler->GetPIDResponse();
1262   if(!fPIDResponse)
1263   {
1264     if(fDebug > 0) printf("TaskV0sInJetsEmcal: No PID response object found\n");
1265     return kFALSE;
1266   }
1267
1268   // AOD files are OK
1269   fh1EventCounterCut->Fill(1);
1270
1271   // Event selection
1272   if(!IsSelectedForJets(fAODIn, fdCutVertexZ, fdCutVertexR2, fdCutCentLow, fdCutCentHigh, 1, 0.1)) // cut on |delta z| in 2011 data between SPD vertex and nominal primary vertex
1273 //  if (!IsSelectedForJets(fAODIn,fdCutVertexZ,fdCutVertexR2,fdCutCentLow,fdCutCentHigh)) // no need for cutting in 2010 data
1274   {
1275     if(fDebug > 5) printf("TaskV0sInJetsEmcal: Event rejected\n");
1276     return kFALSE;
1277   }
1278
1279 //  fdCentrality = fAODIn->GetHeader()->GetCentrality(); // event centrality
1280   fdCentrality = fAODIn->GetHeader()->GetCentralityP()->GetCentralityPercentile("V0M"); // event centrality
1281   if(!fbIsPbPb)
1282     fdCentrality = 0.;
1283   Int_t iCentIndex = GetCentralityBinIndex(fdCentrality); // get index of centrality bin
1284   if(iCentIndex < 0)
1285   {
1286     if(fDebug > 5) printf("TaskV0sInJetsEmcal: Event is out of histogram range\n");
1287     return kFALSE;
1288   }
1289   fh1EventCounterCut->Fill(2); // selected events (vertex, centrality)
1290   fh1EventCounterCutCent[iCentIndex]->Fill(2);
1291
1292   UInt_t iNTracks = fAODIn->GetNumberOfTracks(); // get number of tracks in event
1293   if(fDebug > 5) printf("TaskV0sInJetsEmcal: There are %d tracks in this event\n", iNTracks);
1294
1295   Int_t iNV0s = fAODIn->GetNumberOfV0s(); // get the number of V0 candidates
1296   if(!iNV0s)
1297   {
1298     if(fDebug > 2) printf("TaskV0sInJetsEmcal: No V0s found in event\n");
1299   }
1300
1301   //===== Event is OK for the analysis =====
1302   fh1EventCent->Fill(iCentIndex);
1303   fh1EventCent2->Fill(fdCentrality);
1304   fh2EventCentTracks->Fill(fdCentrality, iNTracks);
1305
1306   if(iNV0s)
1307   {
1308     fh1EventCounterCut->Fill(3); // events with V0s
1309     fh1EventCounterCutCent[iCentIndex]->Fill(3);
1310   }
1311
1312   AliAODv0* v0 = 0; // pointer to V0 candidates
1313   TVector3 vecV0Momentum; // 3D vector of V0 momentum
1314   Double_t dMassV0K0s = 0; // invariant mass of the K0s candidate
1315   Double_t dMassV0Lambda = 0; // invariant mass of the Lambda candidate
1316   Double_t dMassV0ALambda = 0; // invariant mass of the Lambda candidate
1317   Int_t iNV0CandTot = 0; // counter of all V0 candidates at the beginning
1318   Int_t iNV0CandK0s = 0; // counter of K0s candidates at the end
1319   Int_t iNV0CandLambda = 0; // counter of Lambda candidates at the end
1320   Int_t iNV0CandALambda = 0; // counter of Lambda candidates at the end
1321
1322   Bool_t bUseOldCuts = 0; // old reconstruction cuts
1323   Bool_t bUseAliceCuts = 0; // cuts used by Alice Zimmermann
1324   Bool_t bUseIouriCuts = 0; // cuts used by Iouri
1325   Bool_t bPrintCuts = 0; // print out which cuts are applied
1326   Bool_t bPrintJetSelection = 0; // print out which jets are selected
1327
1328   // Values of V0 reconstruction cuts:
1329   // Daughter tracks
1330   Int_t iRefit = AliAODTrack::kTPCrefit; // TPC refit for daughter tracks
1331   Double_t dDCAToPrimVtxMin = fdCutDCAToPrimVtxMin; // 0.1; // [cm] min DCA of daughters to the prim vtx
1332   Double_t dDCADaughtersMax = fdCutDCADaughtersMax; // 1.; // [sigma of TPC tracking] max DCA between daughters
1333   Double_t dEtaDaughterMax = 0.8; // max |pseudorapidity| of daughter tracks
1334   Double_t dNSigmadEdxMax = fdCutNSigmadEdxMax;// 3.; // [sigma dE/dx] max difference between measured and expected signal of dE/dx in the TPC
1335   Double_t dPtProtonPIDMax = 1.; // [GeV/c] maxium pT of proton for applying PID cut
1336   // V0 candidate
1337   Bool_t bOnFly = 0; // on-the-fly (yes) or offline (no) reconstructed
1338   Double_t dCPAMin = fdCutCPAMin;// 0.998; // min cosine of the pointing angle
1339   Double_t dRadiusDecayMin = 5.; // [cm] min radial distance of the decay vertex
1340   Double_t dRadiusDecayMax = 100.; // [cm] max radial distance of the decay vertex
1341   Double_t dEtaMax = 0.7; // max |pseudorapidity| of V0
1342   Double_t dNTauMax = fdCutNTauMax; // 5.0; // [tau] max proper lifetime in multiples of the mean lifetime
1343
1344   // Old cuts Start
1345   Double_t dNCrossedRowsTPCMin = 70.; // min number of crossed TPC rows (turned off)
1346 //  Double_t dCrossedRowsOverFindMin = 0.8; // min ratio crossed rows / findable clusters (turned off)
1347 //  Double_t dCrossedRowsOverFindMax = 1e3; // max ratio crossed rows / findable clusters (turned off)
1348   Double_t dPtDaughterMin = 0.150; // [GeV/c] min transverse momentum of daughter tracks (turned off)
1349   Double_t dRapMax = 0.75; // max |rapidity| of V0 (turned off)
1350   // Old cuts End
1351
1352   // Other cuts
1353   Double_t dNSigmaMassMax = 3.; // [sigma m] max difference between candidate mass and real particle mass (used only for mass peak method of signal extraction)
1354   Double_t dDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
1355
1356   // Selection of active cuts
1357   Bool_t bCutEtaDaughter = 1; // daughter pseudorapidity
1358   Bool_t bCutRapV0 = 0; // V0 rapidity
1359   Bool_t bCutEtaV0 = 1; // V0 pseudorapidity
1360   Bool_t bCutTau = 1; // V0 lifetime
1361   Bool_t bCutPid = 1; // PID (TPC dE/dx)
1362   Bool_t bCutArmPod = 1; // Armenteros-Podolanski for K0S
1363 //  Bool_t bCutCross = 0; // cross contamination
1364
1365   if(bUseOldCuts)
1366   {
1367     bCutRapV0 = 1;
1368     dEtaMax = 0.75;
1369     dNTauMax = 3.0;
1370   }
1371   else if(bUseAliceCuts)
1372   {
1373 //      bOnFly = 1;
1374     dEtaMax = 0.75;
1375     dNTauMax = 5.0;
1376   }
1377   else if(bUseIouriCuts)
1378   {
1379     bCutRapV0 = 1;
1380     bCutEtaV0 = 0;
1381     dNTauMax = 3.0;
1382     dRapMax = 0.5;
1383   }
1384
1385   Double_t dCTauK0s = 2.6844; // [cm] c tau of K0S
1386   Double_t dCTauLambda = 7.89; // [cm] c tau of Lambda
1387
1388   // Load PDG values of particle masses
1389   Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1390   Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1391
1392   // PDG codes of used particles
1393   Int_t iPdgCodePion = 211;
1394   Int_t iPdgCodeProton = 2212;
1395   Int_t iPdgCodeK0s = 310;
1396   Int_t iPdgCodeLambda = 3122;
1397
1398   // Jet selection: fdCutPtJetMin, fdCutPtTrackMin
1399   Double_t dJetEtaWindow = dEtaMax - fdRadiusJet; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1400   Double_t dCutJetAreaMin = 0.6 * TMath::Pi() * fdRadiusJet * fdRadiusJet; // minimum jet area
1401   Double_t dRadiusExcludeCone = 2 * fdRadiusJet; // radius of cones around jets excluded for V0 outside jets
1402   Bool_t bLeadingJetOnly = 0;
1403
1404   if(bUseAliceCuts)
1405   {
1406     fdCutPtJetMin = 5;
1407     fdCutPtTrackMin = 5;
1408     dCutJetAreaMin = 0;
1409     bLeadingJetOnly = 0;
1410   }
1411
1412   if(fJetsCont)
1413   {
1414 //    fJetsCont->SetJetPtCut(fdCutPtJetMin); // needs to be applied on the pt after bg subtraction
1415     fJetsCont->SetPtBiasJetTrack(fdCutPtTrackMin);
1416     fJetsCont->SetPercAreaCut(0.6);
1417     fJetsCont->SetJetEtaLimits(-dJetEtaWindow, dJetEtaWindow);
1418   }
1419
1420   Int_t iNJet = 0; // number of reconstructed jets in the input
1421   TClonesArray* jetArraySel = new TClonesArray("AliAODJet", 0); // object where the selected jets are copied
1422   Int_t iNJetSel = 0; // number of selected reconstructed jets
1423   TClonesArray* jetArrayPerp = new TClonesArray("AliAODJet", 0); // object where the perp. cones are stored
1424   Int_t iNJetPerp = 0; // number of perpendicular cones
1425
1426   AliAODJet* jet = 0; // pointer to a jet
1427   AliAODJet* jetPerp = 0; // pointer to a perp. cone
1428   AliAODJet* jetRnd = 0; // pointer to a rand. cone
1429   AliEmcalJet* jetMed = 0; // pointer to a median cluster
1430   TVector3 vecJetMomentum; // 3D vector of jet momentum
1431   Bool_t bJetEventGood = kTRUE; // indicator of good jet events
1432   Double_t dRho = 0; // average bg pt density
1433   TLorentzVector vecJetSel; // 4-momentum of selected jet
1434   TLorentzVector vecPerpPlus; // 4-momentum of perpendicular cone plus
1435   TLorentzVector vecPerpMinus; // 4-momentum of perpendicular cone minus
1436
1437   if(fbJetSelection)  // analysis of V0s in jets is switched on
1438   {
1439     if(!fJetsCont)
1440     {
1441       if(fDebug > 0) printf("TaskV0sInJetsEmcal: No jet container\n");
1442       bJetEventGood = kFALSE;
1443     }
1444     if(bJetEventGood)
1445       iNJet = fJetsCont->GetNJets();
1446     if(bJetEventGood && !iNJet)  // check whether there are some jets
1447     {
1448       if(fDebug > 2) printf("TaskV0sInJetsEmcal: No jets in array\n");
1449       bJetEventGood = kFALSE;
1450     }
1451     if(bJetEventGood && !fJetsBgCont)
1452     {
1453       if(fDebug > 0) printf("TaskV0sInJetsEmcal: No bg jet container\n");
1454 //        bJetEventGood = kFALSE;
1455     }
1456   }
1457   else // no in-jet analysis
1458     bJetEventGood = kFALSE;
1459
1460   // select good jets and copy them to another array
1461   if(bJetEventGood)
1462   {
1463     dRho = fJetsCont->GetRhoVal();
1464 //    printf("TaskV0sInJetsEmcal: Loaded rho value: %g\n",dRho);
1465     if(bLeadingJetOnly)
1466       iNJet = 1; // only leading jets
1467     if(fDebug > 5) printf("TaskV0sInJetsEmcal: Jet selection for %d jets\n", iNJet);
1468     for(Int_t iJet = 0; iJet < iNJet; iJet++)
1469     {
1470       AliEmcalJet* jetSel = (AliEmcalJet*)(fJetsCont->GetAcceptJet(iJet)); // load a jet in the list
1471       if(bLeadingJetOnly)
1472         jetSel = fJetsCont->GetLeadingJet();
1473       if(!jetSel)
1474       {
1475         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Jet %d not accepted in container\n", iJet);
1476         continue;
1477       }
1478       Double_t dPtJetCorr = jetSel->Pt() - dRho * jetSel->Area();
1479       if(bPrintJetSelection)
1480         if(fDebug > 7) printf("jet: i = %d, pT = %g, eta = %g, phi = %g, pt lead tr = %g, pt corr = %g ", iJet, jetSel->Pt(), jetSel->Eta(), jetSel->Phi(), fJetsCont->GetLeadingHadronPt(jetSel), dPtJetCorr);
1481 //          printf("TaskV0sInJetsEmcal: Checking pt > %.2f for jet %d with pt %.2f\n",fdCutPtJetMin,iJet,jetSel->Pt());
1482       if(dPtJetCorr < fdCutPtJetMin)  // selection of high-pt jets, needs to be applied on the pt after bg subtraction
1483       {
1484         if(bPrintJetSelection)
1485           if(fDebug > 7) printf("rejected (pt)\n");
1486         continue;
1487       }
1488 //          printf("TaskV0sInJetsEmcal: Checking |eta| < %.2f for jet %d with |eta| %.2f\n",dJetEtaWindow,iJet,TMath::Abs(jetSel->Eta()));
1489       if(TMath::Abs(jetSel->Eta()) > dJetEtaWindow)  // selection of jets in the chosen pseudorapidity range
1490       {
1491         if(bPrintJetSelection)
1492           if(fDebug > 7) printf("rejected (eta)\n");
1493         continue;
1494       }
1495       if(!bUseOldCuts)
1496       {
1497         if(jetSel->Area() < dCutJetAreaMin)
1498         {
1499           if(bPrintJetSelection)
1500             if(fDebug > 7) printf("rejected (area)\n");
1501           continue;
1502         }
1503       }
1504       Double_t dPtTrack = fJetsCont->GetLeadingHadronPt(jetSel);
1505       if(fdCutPtTrackMin > 0)  // a positive min leading track pt is set
1506       {
1507         if(dPtTrack < fdCutPtTrackMin)  // selection of high-pt jet-track events
1508         {
1509           if(bPrintJetSelection)
1510             if(fDebug > 7) printf("rejected (track pt)\n");
1511           continue;
1512         }
1513       }
1514       if(bPrintJetSelection)
1515         if(fDebug > 7) printf("accepted\n");
1516       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Jet %d with pt %.2f passed selection\n", iJet, dPtJetCorr);
1517
1518       vecJetSel.SetPtEtaPhiM(dPtJetCorr, jetSel->Eta(), jetSel->Phi(), 0.);
1519       vecPerpPlus.SetPtEtaPhiM(dPtJetCorr, jetSel->Eta(), jetSel->Phi(), 0.);
1520       vecPerpMinus.SetPtEtaPhiM(dPtJetCorr, jetSel->Eta(), jetSel->Phi(), 0.);
1521       vecPerpPlus.RotateZ(TMath::Pi() / 2.); // rotate vector by +90 deg around z
1522       vecPerpMinus.RotateZ(-TMath::Pi() / 2.); // rotate vector by -90 deg around z
1523       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Adding perp. cones number %d, %d\n", iNJetPerp, iNJetPerp + 1);
1524       new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpPlus); // write perp. cone to the array
1525       new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpMinus); // write perp. cone to the array
1526       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Adding jet number %d\n", iNJetSel);
1527       new((*jetArraySel)[iNJetSel++]) AliAODJet(vecJetSel); // copy selected jet to the array
1528       fh1PtJetTrackLeading[iCentIndex]->Fill(dPtTrack); // pt of leading jet track
1529     }
1530     if(fDebug > 5) printf("TaskV0sInJetsEmcal: Added jets: %d\n", iNJetSel);
1531     iNJetSel = jetArraySel->GetEntriesFast();
1532     if(fDebug > 2) printf("TaskV0sInJetsEmcal: Selected jets in array: %d\n", iNJetSel);
1533     fh1NJetPerEvent[iCentIndex]->Fill(iNJetSel);
1534     // fill jet spectra
1535     for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1536     {
1537       jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
1538       fh1PtJet[iCentIndex]->Fill(jet->Pt()); // pt spectrum of selected jets
1539       fh1EtaJet[iCentIndex]->Fill(jet->Eta()); // eta spectrum of selected jets
1540       fh2EtaPtJet[iCentIndex]->Fill(jet->Eta(), jet->Pt()); // eta-pT spectrum of selected jets
1541       fh1PhiJet[iCentIndex]->Fill(jet->Phi()); // phi spectrum of selected jets
1542       Double_t dAreaExcluded = TMath::Pi() * dRadiusExcludeCone * dRadiusExcludeCone; // area of the cone
1543       dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax - jet->Eta()); // positive eta overhang
1544       dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax + jet->Eta()); // negative eta overhang
1545       fh1AreaExcluded->Fill(iCentIndex, dAreaExcluded);
1546     }
1547     jet = 0;
1548   }
1549
1550   if(bJetEventGood)  // there should be some reconstructed jets
1551   {
1552     fh1EventCounterCut->Fill(4); // events with jet(s)
1553     fh1EventCounterCutCent[iCentIndex]->Fill(4); // events with jet(s)
1554     if(iNJetSel)
1555     {
1556       fh1EventCounterCut->Fill(5); // events with selected jets
1557       fh1EventCounterCutCent[iCentIndex]->Fill(5);
1558     }
1559   }
1560   if(iNJetSel)
1561     fh1EventCent2Jets->Fill(fdCentrality);
1562   else
1563     fh1EventCent2NoJets->Fill(fdCentrality);
1564
1565   if(iNJetSel)
1566   {
1567     jetRnd = GetRandomCone(jetArraySel, dJetEtaWindow, 2 * fdRadiusJet);
1568     if(jetRnd)
1569     {
1570       fh1NRndConeCent->Fill(iCentIndex);
1571       fh2EtaPhiRndCone[iCentIndex]->Fill(jetRnd->Eta(), jetRnd->Phi());
1572     }
1573     jetMed = GetMedianCluster(fJetsBgCont, dJetEtaWindow);
1574     if(jetMed)
1575     {
1576       fh1NMedConeCent->Fill(iCentIndex);
1577       fh2EtaPhiMedCone[iCentIndex]->Fill(jetMed->Eta(), jetMed->Phi());
1578     }
1579   }
1580
1581   // Loading primary vertex info
1582   AliAODVertex* primVtx = fAODIn->GetPrimaryVertex(); // get the primary vertex
1583   Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
1584   primVtx->GetXYZ(dPrimVtxPos);
1585   fh1VtxZ[iCentIndex]->Fill(dPrimVtxPos[2]);
1586   fh2VtxXY[iCentIndex]->Fill(dPrimVtxPos[0], dPrimVtxPos[1]);
1587
1588   //===== Start of loop over V0 candidates =====
1589   if(fDebug > 2) printf("TaskV0sInJetsEmcal: Start of V0 loop\n");
1590   for(Int_t iV0 = 0; iV0 < iNV0s; iV0++)
1591   {
1592     v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1593     if(!v0)
1594       continue;
1595
1596     iNV0CandTot++;
1597
1598     // Initialization of status indicators
1599     Bool_t bIsCandidateK0s = kTRUE; // candidate for K0s
1600     Bool_t bIsCandidateLambda = kTRUE; // candidate for Lambda
1601     Bool_t bIsCandidateALambda = kTRUE; // candidate for anti-Lambda
1602     Bool_t bIsInPeakK0s = kFALSE; // candidate within the K0s mass peak
1603     Bool_t bIsInPeakLambda = kFALSE; // candidate within the Lambda mass peak
1604     Bool_t bIsInConeJet = kFALSE; // candidate within the jet cones
1605     Bool_t bIsInConePerp = kFALSE; // candidate within the perpendicular cone
1606     Bool_t bIsInConeRnd = kFALSE; // candidate within the random cone
1607     Bool_t bIsInConeMed = kFALSE; // candidate within the median-cluster cone
1608     Bool_t bIsOutsideCones = kFALSE; // candidate outside excluded cones
1609
1610     // Invariant mass calculation
1611     dMassV0K0s = v0->MassK0Short();
1612     dMassV0Lambda = v0->MassLambda();
1613     dMassV0ALambda = v0->MassAntiLambda();
1614
1615     Int_t iCutIndex = 0; // indicator of current selection step
1616     // 0
1617     // All V0 candidates
1618     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1619     iCutIndex++;
1620
1621     // Skip candidates outside the histogram range
1622     if((dMassV0K0s < fgkdMassK0sMin) || (dMassV0K0s >= fgkdMassK0sMax))
1623       bIsCandidateK0s = kFALSE;
1624     if((dMassV0Lambda < fgkdMassLambdaMin) || (dMassV0Lambda >= fgkdMassLambdaMax))
1625       bIsCandidateLambda = kFALSE;
1626     if((dMassV0ALambda < fgkdMassLambdaMin) || (dMassV0ALambda >= fgkdMassLambdaMax))
1627       bIsCandidateALambda = kFALSE;
1628     if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1629       continue;
1630
1631     Double_t dPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
1632     vecV0Momentum = TVector3(v0->Px(), v0->Py(), v0->Pz()); // set the vector of V0 momentum
1633
1634     // Sigma of the mass peak window
1635     Double_t dMassPeakWindowK0s = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 0);
1636     Double_t dMassPeakWindowLambda = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 1);
1637
1638     // Invariant mass peak selection
1639     if(TMath::Abs(dMassV0K0s - dMassPDGK0s) < dMassPeakWindowK0s)
1640       bIsInPeakK0s = kTRUE;
1641     if(TMath::Abs(dMassV0Lambda - dMassPDGLambda) < dMassPeakWindowLambda)
1642       bIsInPeakLambda = kTRUE;
1643
1644     // Retrieving all relevant properties of the V0 candidate
1645     Bool_t bOnFlyStatus = v0->GetOnFlyStatus(); // online (on fly) reconstructed vs offline reconstructed
1646     const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1647     const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
1648     Double_t dPtDaughterPos = trackPos->Pt(); // transverse momentum of a daughter track
1649     Double_t dPtDaughterNeg = trackNeg->Pt();
1650     Double_t dNRowsPos = trackPos->GetTPCClusterInfo(2, 1); // crossed TPC pad rows of a daughter track
1651     Double_t dNRowsNeg = trackNeg->GetTPCClusterInfo(2, 1);
1652     Double_t dDCAToPrimVtxPos = TMath::Abs(v0->DcaPosToPrimVertex()); // dca of a daughter to the primary vertex
1653     Double_t dDCAToPrimVtxNeg = TMath::Abs(v0->DcaNegToPrimVertex());
1654     Double_t dDCADaughters = v0->DcaV0Daughters(); // dca between daughters
1655     Double_t dCPA = v0->CosPointingAngle(primVtx); // cosine of the pointing angle
1656     Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
1657 //      Double_t dSecVtxPos[3] = {v0->DecayVertexV0X(),v0->DecayVertexV0Y(),v0->DecayVertexV0Z()}; // V0 vertex position
1658     v0->GetSecondaryVtx(dSecVtxPos);
1659     Double_t dRadiusDecay = TMath::Sqrt(dSecVtxPos[0] * dSecVtxPos[0] + dSecVtxPos[1] * dSecVtxPos[1]); // distance of the V0 vertex from the z-axis
1660     Double_t dEtaDaughterNeg = trackNeg->Eta(); // = v0->EtaProng(1), pseudorapidity of a daughter track
1661     Double_t dEtaDaughterPos = trackPos->Eta(); // = v0->EtaProng(0)
1662     Double_t dRapK0s = v0->RapK0Short(); // rapidity calculated for K0s assumption
1663     Double_t dRapLambda = v0->RapLambda(); // rapidity calculated for Lambda assumption
1664     Double_t dEtaV0 = v0->Eta(); // V0 pseudorapidity
1665 //      Double_t dPhiV0 = v0->Phi(); // V0 pseudorapidity
1666     Double_t dDecayPath[3];
1667     for(Int_t iPos = 0; iPos < 3; iPos++)
1668       dDecayPath[iPos] = dSecVtxPos[iPos] - dPrimVtxPos[iPos]; // vector of the V0 path
1669     Double_t dDecLen = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]); // path length L
1670     Double_t dDecLen2D = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]); // transverse path length R
1671     Double_t dLOverP = dDecLen / v0->P(); // L/p
1672     Double_t dROverPt = dDecLen2D / dPtV0; // R/pT
1673     Double_t dMLOverPK0s = dMassPDGK0s * dLOverP; // m*L/p = c*(proper lifetime)
1674 //      Double_t dMLOverPLambda = dMassPDGLambda*dLOverP; // m*L/p
1675     Double_t dMROverPtK0s = dMassPDGK0s * dROverPt; // m*R/pT
1676     Double_t dMROverPtLambda = dMassPDGLambda * dROverPt; // m*R/pT
1677     Double_t dNSigmaPosPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kPion)); // difference between measured and expected signal of the dE/dx in the TPC
1678     Double_t dNSigmaPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kProton));
1679     Double_t dNSigmaNegPion   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kPion));
1680     Double_t dNSigmaNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kProton));
1681     Double_t dAlpha = v0->AlphaV0(); // Armenteros-Podolanski alpha
1682     Double_t dPtArm = v0->PtArmV0(); // Armenteros-Podolanski pT
1683     AliAODVertex* prodVtxDaughterPos = (AliAODVertex*)(trackPos->GetProdVertex()); // production vertex of the positive daughter track
1684     Char_t cTypeVtxProdPos = prodVtxDaughterPos->GetType(); // type of the production vertex
1685     AliAODVertex* prodVtxDaughterNeg = (AliAODVertex*)(trackNeg->GetProdVertex()); // production vertex of the negative daughter track
1686     Char_t cTypeVtxProdNeg = prodVtxDaughterNeg->GetType(); // type of the production vertex
1687
1688 //    fh2Tau3DVs2D[0]->Fill(dPtV0, dLOverP / dROverPt);
1689
1690     // QA histograms before cuts
1691     FillQAHistogramV0(primVtx, v0, 0, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
1692     // Cut vs mass histograms before cuts
1693     /*
1694     if(bIsCandidateK0s)
1695     {
1696       fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsPos);
1697       fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsNeg);
1698       fh2CutPtPosK0s[0]->Fill(dMassV0K0s, dPtDaughterPos);
1699       fh2CutPtNegK0s[0]->Fill(dMassV0K0s, dPtDaughterNeg);
1700       fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
1701       fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
1702       fh2CutDCAV0[0]->Fill(dMassV0K0s, dDCADaughters);
1703       fh2CutCos[0]->Fill(dMassV0K0s, dCPA);
1704       fh2CutR[0]->Fill(dMassV0K0s, dRadiusDecay);
1705       fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterPos);
1706       fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterNeg);
1707       fh2CutRapK0s[0]->Fill(dMassV0K0s, dRapK0s);
1708       fh2CutCTauK0s[0]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
1709       fh2CutPIDPosK0s[0]->Fill(dMassV0K0s, dNSigmaPosPion);
1710       fh2CutPIDNegK0s[0]->Fill(dMassV0K0s, dNSigmaNegPion);
1711     }
1712     if(bIsCandidateLambda)
1713     {
1714       fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsPos);
1715       fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsNeg);
1716       fh2CutPtPosLambda[0]->Fill(dMassV0Lambda, dPtDaughterPos);
1717       fh2CutPtNegLambda[0]->Fill(dMassV0Lambda, dPtDaughterNeg);
1718       fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterPos);
1719       fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterNeg);
1720       fh2CutRapLambda[0]->Fill(dMassV0Lambda, dRapLambda);
1721       fh2CutCTauLambda[0]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
1722       fh2CutPIDPosLambda[0]->Fill(dMassV0Lambda, dNSigmaPosProton);
1723       fh2CutPIDNegLambda[0]->Fill(dMassV0Lambda, dNSigmaNegPion);
1724     }
1725     */
1726
1727     //===== Start of reconstruction cutting =====
1728
1729     // 1
1730     // All V0 candidates
1731     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1732     iCutIndex++;
1733
1734     // Start of global cuts
1735     // 2
1736     // Reconstruction method
1737     if(bPrintCuts) printf("Rec: Applying cut: Reconstruction method: on-the-fly? %s\n", (bOnFly ? "yes" : "no"));
1738     if(bOnFlyStatus != bOnFly)
1739       continue;
1740     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1741     iCutIndex++;
1742
1743     // 3
1744     // Tracks TPC OK
1745     if(bPrintCuts) printf("Rec: Applying cut: Correct charge of daughters\n");
1746     if(!trackNeg || !trackPos)
1747       continue;
1748     if(trackNeg->Charge() == trackPos->Charge())  // daughters have different charge?
1749       continue;
1750     if(trackNeg->Charge() != -1)  // daughters have expected charge?
1751       continue;
1752     if(trackPos->Charge() != 1)  // daughters have expected charge?
1753       continue;
1754
1755     if(bPrintCuts) printf("Rec: Applying cut: TPC refit: %d\n", iRefit);
1756     if(!trackNeg->IsOn(iRefit))  // TPC refit is ON?
1757       continue;
1758     if(bPrintCuts) printf("Rec: Applying cut: Type of production vertex of daughter: Not %d\n", AliAODVertex::kKink);
1759     if(cTypeVtxProdNeg == AliAODVertex::kKink) // kink daughter rejection
1760       continue;
1761     // Old cuts Start
1762     if(bUseOldCuts)
1763     {
1764       if(bPrintCuts) printf("Rec: Applying cut: Number of TPC rows: > %f\n", dNCrossedRowsTPCMin);
1765       if(dNRowsNeg < dNCrossedRowsTPCMin)  // Crossed TPC padrows
1766         continue;
1767 //      Int_t findable = trackNeg->GetTPCNclsF(); // Findable clusters
1768 //      if (findable <= 0)
1769 //        continue;
1770 //      if (dNRowsNeg/findable < dCrossedRowsOverFindMin)
1771 //        continue;
1772 //      if (dNRowsNeg/findable > dCrossedRowsOverFindMax)
1773 //        continue;
1774     }
1775     // Old cuts End
1776
1777     if(!trackPos->IsOn(iRefit))
1778       continue;
1779     if(cTypeVtxProdPos == AliAODVertex::kKink) // kink daughter rejection
1780       continue;
1781     // Old cuts Start
1782     if(bUseOldCuts)
1783     {
1784       if(dNRowsPos < dNCrossedRowsTPCMin)
1785         continue;
1786 //      findable = trackPos->GetTPCNclsF();
1787 //      if (findable <= 0)
1788 //        continue;
1789 //      if (dNRowsPos/findable < dCrossedRowsOverFindMin)
1790 //        continue;
1791 //      if (dNRowsPos/findable > dCrossedRowsOverFindMax)
1792 //        continue;
1793     }
1794     // Old cuts End
1795
1796     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1797     iCutIndex++;
1798
1799     // 4
1800     // Daughters: transverse momentum cut
1801     if(bUseOldCuts)
1802     {
1803       if(bPrintCuts) printf("Rec: Applying cut: Daughter pT: > %f\n", dPtDaughterMin);
1804       if((dPtDaughterNeg < dPtDaughterMin) || (dPtDaughterPos < dPtDaughterMin))
1805         continue;
1806       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1807     }
1808     iCutIndex++;
1809
1810     // 5
1811     // Daughters: Impact parameter of daughters to prim vtx
1812     if(bPrintCuts) printf("Rec: Applying cut: Daughter DCA to prim vtx: > %f\n", dDCAToPrimVtxMin);
1813     if((dDCAToPrimVtxNeg < dDCAToPrimVtxMin) || (dDCAToPrimVtxPos < dDCAToPrimVtxMin))
1814       continue;
1815     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1816     iCutIndex++;
1817
1818     // 6
1819     // Daughters: DCA
1820     if(bPrintCuts) printf("Rec: Applying cut: DCA between daughters: < %f\n", dDCADaughtersMax);
1821     if(dDCADaughters > dDCADaughtersMax)
1822       continue;
1823     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1824     iCutIndex++;
1825
1826     // 7
1827     // V0: Cosine of the pointing angle
1828     if(bPrintCuts) printf("Rec: Applying cut: CPA: > %f\n", dCPAMin);
1829     if(dCPA < dCPAMin)
1830       continue;
1831     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1832     iCutIndex++;
1833
1834     // 8
1835     // V0: Fiducial volume
1836     if(bPrintCuts) printf("Rec: Applying cut: Decay radius: > %f, < %f\n", dRadiusDecayMin, dRadiusDecayMax);
1837     if((dRadiusDecay < dRadiusDecayMin) || (dRadiusDecay > dRadiusDecayMax))
1838       continue;
1839     FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1840     iCutIndex++;
1841
1842     // 9
1843     // Daughters: pseudorapidity cut
1844     if(bCutEtaDaughter)
1845     {
1846       if(bPrintCuts) printf("Rec: Applying cut: Daughter |eta|: < %f\n", dEtaDaughterMax);
1847       if((TMath::Abs(dEtaDaughterNeg) > dEtaDaughterMax) || (TMath::Abs(dEtaDaughterPos) > dEtaDaughterMax))
1848         continue;
1849       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1850     }
1851     iCutIndex++;
1852     // End of global cuts
1853
1854     // Start of particle-dependent cuts
1855     // 10
1856     // V0: rapidity cut & pseudorapidity cut
1857     if(bCutRapV0)
1858     {
1859       if(bPrintCuts) printf("Rec: Applying cut: V0 |y|: < %f\n", dRapMax);
1860       if(TMath::Abs(dRapK0s) > dRapMax)
1861         bIsCandidateK0s = kFALSE;
1862       if(TMath::Abs(dRapLambda) > dRapMax)
1863       {
1864         bIsCandidateLambda = kFALSE;
1865         bIsCandidateALambda = kFALSE;
1866       }
1867     }
1868     if(bCutEtaV0)
1869     {
1870       if(bPrintCuts) printf("Rec: Applying cut: V0 |eta|: < %f\n", dEtaMax);
1871       if(TMath::Abs(dEtaV0) > dEtaMax)
1872       {
1873         bIsCandidateK0s = kFALSE;
1874         bIsCandidateLambda = kFALSE;
1875         bIsCandidateALambda = kFALSE;
1876       }
1877       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1878     }
1879     iCutIndex++;
1880
1881     // 11
1882     // Lifetime cut
1883     if(bCutTau)
1884     {
1885       if(bPrintCuts) printf("Rec: Applying cut: Proper lifetime: < %f\n", dNTauMax);
1886       if(dMROverPtK0s > dNTauMax * dCTauK0s)
1887         bIsCandidateK0s = kFALSE;
1888       if(dMROverPtLambda > dNTauMax * dCTauLambda)
1889       {
1890         bIsCandidateLambda = kFALSE;
1891         bIsCandidateALambda = kFALSE;
1892       }
1893       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1894     }
1895     iCutIndex++;
1896
1897     // 12
1898     // Daughter PID
1899     if(bCutPid)
1900     {
1901       if(bUseOldCuts)
1902       {
1903         if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (both daughters): < %f\n", dNSigmadEdxMax);
1904         if(dNSigmaPosPion > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax)  // pi+, pi-
1905           bIsCandidateK0s = kFALSE;
1906         if(dNSigmaPosProton > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax)  // p+, pi-
1907           bIsCandidateLambda = kFALSE;
1908         if(dNSigmaNegProton > dNSigmadEdxMax || dNSigmaPosPion > dNSigmadEdxMax)  // p-, pi+
1909           bIsCandidateALambda = kFALSE;
1910       }
1911       else
1912       {
1913         if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (proton below %f GeV/c): < %f\n", dPtProtonPIDMax, dNSigmadEdxMax);
1914         if((dPtDaughterPos < dPtProtonPIDMax) && (dNSigmaPosProton > dNSigmadEdxMax))    // p+
1915           bIsCandidateLambda = kFALSE;
1916         if((dPtDaughterNeg < dPtProtonPIDMax) && (dNSigmaNegProton > dNSigmadEdxMax))    // p-
1917           bIsCandidateALambda = kFALSE;
1918       }
1919       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1920     }
1921     iCutIndex++;
1922
1923     Double_t valueCorrel[3] = {dMassV0K0s, dMassV0Lambda, dPtV0};
1924     if(bIsCandidateK0s && bIsCandidateLambda)
1925       fh3CCMassCorrelBoth->Fill(valueCorrel); // correlation of mass distribution of candidates selected as both K0s and Lambda
1926     if(bIsCandidateK0s && !bIsCandidateLambda)
1927       fh3CCMassCorrelKNotL->Fill(valueCorrel); // correlation of mass distribution of candidates selected as K0s and not Lambda
1928     if(!bIsCandidateK0s && bIsCandidateLambda)
1929       fh3CCMassCorrelLNotK->Fill(valueCorrel); // correlation of mass distribution of candidates selected as not K0s and Lambda
1930
1931     // 13
1932     // Armenteros-Podolanski cut
1933     if(bCutArmPod)
1934     {
1935       if(bPrintCuts) printf("Rec: Applying cut: Armenteros-Podolanski (K0S): pT > %f * |alpha|\n", 0.2);
1936       if(dPtArm < TMath::Abs(0.2 * dAlpha))
1937         bIsCandidateK0s = kFALSE;
1938       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1939     }
1940     iCutIndex++;
1941
1942     // Cross contamination
1943     if(bIsInPeakK0s)
1944     {
1945       if(bIsCandidateLambda)  // Lambda candidates in K0s peak, excluded from Lambda candidates by CC cut
1946         fh2CCLambda->Fill(dMassV0Lambda, dPtV0);
1947     }
1948     if(bIsInPeakLambda)
1949     {
1950       if(bIsCandidateK0s)  // K0s candidates in Lambda peak, excluded from K0s candidates by CC cut
1951         fh2CCK0s->Fill(dMassV0K0s, dPtV0);
1952     }
1953 //      if (bCutCross)
1954 //        {
1955 //          if (bIsInPeakK0s)
1956 //            bIsCandidateLambda = kFALSE;
1957 //          if (bIsInPeakLambda)
1958 //            bIsCandidateK0s = kFALSE;
1959 //          FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1960 //        }
1961 //      iCutIndex++;
1962
1963     // End of particle-dependent cuts
1964
1965     //===== End of reconstruction cutting =====
1966
1967     if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1968       continue;
1969
1970     // Selection of V0s in jet cones, perpendicular cones, random cones, outside cones
1971     if(bJetEventGood && iNJetSel && (bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda))
1972     {
1973       // Selection of V0s in jet cones
1974       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for V0 %d %d in %d jet cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetSel);
1975       for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1976       {
1977         jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
1978         vecJetMomentum.SetXYZ(jet->Px(), jet->Py(), jet->Pz()); // set the vector of jet momentum
1979         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Checking if V0 %d %d in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
1980         if(IsParticleInCone(v0, jet, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
1981         {
1982           if(fDebug > 5) printf("TaskV0sInJetsEmcal: V0 %d %d found in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
1983           bIsInConeJet = kTRUE;
1984           break;
1985         }
1986       }
1987       // Selection of V0s in perp. cones
1988       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for V0 %d %d in %d perp. cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetPerp);
1989       for(Int_t iJet = 0; iJet < iNJetPerp; iJet++)
1990       {
1991         jetPerp = (AliAODJet*)jetArrayPerp->At(iJet); // load a jet in the list
1992         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Checking if V0 %d %d in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
1993         if(IsParticleInCone(v0, jetPerp, fdRadiusJet)) // V0 in perp. cone
1994         {
1995           if(fDebug > 5) printf("TaskV0sInJetsEmcal: V0 %d %d found in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
1996           bIsInConePerp = kTRUE;
1997           break;
1998         }
1999       }
2000       // Selection of V0s in random cones
2001       if(jetRnd)
2002       {
2003         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for V0 %d %d in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2004         if(IsParticleInCone(v0, jetRnd, fdRadiusJet)) // V0 in rnd. cone?
2005         {
2006           if(fDebug > 5) printf("TaskV0sInJetsEmcal: V0 %d %d found in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2007           bIsInConeRnd = kTRUE;
2008         }
2009       }
2010       // Selection of V0s in median-cluster cones
2011       if(jetMed)
2012       {
2013         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for V0 %d %d in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2014         if(IsParticleInCone(v0, jetMed, fdRadiusJet)) // V0 in med. cone?
2015         {
2016           if(fDebug > 5) printf("TaskV0sInJetsEmcal: V0 %d %d found in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2017           bIsInConeMed = kTRUE;
2018         }
2019       }
2020       // Selection of V0s outside jet cones
2021       if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for V0 %d %d outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2022       if(!OverlapWithJets(jetArraySel, v0, dRadiusExcludeCone)) // V0 oustide jet cones
2023       {
2024         if(fDebug > 5) printf("TaskV0sInJetsEmcal: V0 %d %d found outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2025         bIsOutsideCones = kTRUE;
2026       }
2027     }
2028
2029     // QA histograms after cuts
2030     FillQAHistogramV0(primVtx, v0, 1, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
2031     // Cut vs mass histograms after cuts
2032     /*
2033     if(bIsCandidateK0s)
2034     {
2035       fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsPos);
2036       fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsNeg);
2037       fh2CutPtPosK0s[1]->Fill(dMassV0K0s, dPtDaughterPos);
2038       fh2CutPtNegK0s[1]->Fill(dMassV0K0s, dPtDaughterNeg);
2039       fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
2040       fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
2041       fh2CutDCAV0[1]->Fill(dMassV0K0s, dDCADaughters);
2042       fh2CutCos[1]->Fill(dMassV0K0s, dCPA);
2043       fh2CutR[1]->Fill(dMassV0K0s, dRadiusDecay);
2044       fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterPos);
2045       fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterNeg);
2046       fh2CutRapK0s[1]->Fill(dMassV0K0s, dRapK0s);
2047       fh2CutCTauK0s[1]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
2048       fh2CutPIDPosK0s[1]->Fill(dMassV0K0s, dNSigmaPosPion);
2049       fh2CutPIDNegK0s[1]->Fill(dMassV0K0s, dNSigmaNegPion);
2050     }
2051     if(bIsCandidateLambda)
2052     {
2053       fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsPos);
2054       fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsNeg);
2055       fh2CutPtPosLambda[1]->Fill(dMassV0Lambda, dPtDaughterPos);
2056       fh2CutPtNegLambda[1]->Fill(dMassV0Lambda, dPtDaughterNeg);
2057       fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterPos);
2058       fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterNeg);
2059       fh2CutRapLambda[1]->Fill(dMassV0Lambda, dRapLambda);
2060       fh2CutCTauLambda[1]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
2061       fh2CutPIDPosLambda[1]->Fill(dMassV0Lambda, dNSigmaPosProton);
2062       fh2CutPIDNegLambda[1]->Fill(dMassV0Lambda, dNSigmaNegPion);
2063     }
2064     */
2065
2066     //===== Start of filling V0 spectra =====
2067
2068     Double_t dAngle = TMath::Pi(); // angle between V0 momentum and jet momentum
2069     if(bIsInConeJet)
2070     {
2071       dAngle = vecV0Momentum.Angle(vecJetMomentum);
2072     }
2073
2074     // iCutIndex = 14
2075     if(bIsCandidateK0s)
2076     {
2077       // 14 K0s candidates after cuts
2078 //          printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0K0s,dPtV0,dEtaV0,dPhiV0);
2079       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex, iCentIndex);
2080       Double_t valueKIncl[3] = {dMassV0K0s, dPtV0, dEtaV0};
2081       fhnV0InclusiveK0s[iCentIndex]->Fill(valueKIncl);
2082       fh1V0InvMassK0sCent[iCentIndex]->Fill(dMassV0K0s);
2083
2084       fh1QACTau2D[1]->Fill(dMROverPtK0s / dCTauK0s);
2085       fh1QACTau3D[1]->Fill(dMLOverPK0s / dCTauK0s);
2086 //      fh2Tau3DVs2D[1]->Fill(dPtV0, dLOverP / dROverPt);
2087
2088       if(iNJetSel)
2089       {
2090         // 15 K0s in jet events
2091         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 1, iCentIndex);
2092       }
2093       if(bIsInConeJet)
2094       {
2095         // 16 K0s in jets
2096         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 2, iCentIndex);
2097         Double_t valueKInJC[4] = {dMassV0K0s, dPtV0, dEtaV0, jet->Pt()};
2098         fhnV0InJetK0s[iCentIndex]->Fill(valueKInJC);
2099         fh2V0PtJetAngleK0s[iCentIndex]->Fill(jet->Pt(), dAngle);
2100       }
2101       if(bIsOutsideCones)
2102       {
2103         Double_t valueKOutJC[3] = {dMassV0K0s, dPtV0, dEtaV0};
2104         fhnV0OutJetK0s[iCentIndex]->Fill(valueKOutJC);
2105       }
2106       if(bIsInConePerp)
2107       {
2108         Double_t valueKInPC[4] = {dMassV0K0s, dPtV0, dEtaV0, jetPerp->Pt()};
2109         fhnV0InPerpK0s[iCentIndex]->Fill(valueKInPC);
2110       }
2111       if(bIsInConeRnd)
2112       {
2113         Double_t valueKInRnd[3] = {dMassV0K0s, dPtV0, dEtaV0};
2114         fhnV0InRndK0s[iCentIndex]->Fill(valueKInRnd);
2115       }
2116       if(bIsInConeMed)
2117       {
2118         Double_t valueKInMed[3] = {dMassV0K0s, dPtV0, dEtaV0};
2119         fhnV0InMedK0s[iCentIndex]->Fill(valueKInMed);
2120       }
2121       if(!iNJetSel)
2122       {
2123         Double_t valueKNoJet[3] = {dMassV0K0s, dPtV0, dEtaV0};
2124         fhnV0NoJetK0s[iCentIndex]->Fill(valueKNoJet);
2125       }
2126       iNV0CandK0s++;
2127     }
2128     if(bIsCandidateLambda)
2129     {
2130       // 14 Lambda candidates after cuts
2131 //          printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0Lambda,dPtV0,dEtaV0,dPhiV0);
2132       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex, iCentIndex);
2133       Double_t valueLIncl[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2134       fhnV0InclusiveLambda[iCentIndex]->Fill(valueLIncl);
2135       fh1V0InvMassLambdaCent[iCentIndex]->Fill(dMassV0Lambda);
2136       if(iNJetSel)
2137       {
2138         // 15 Lambda in jet events
2139         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 1, iCentIndex);
2140       }
2141       if(bIsInConeJet)
2142       {
2143         // 16 Lambda in jets
2144         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 2, iCentIndex);
2145         Double_t valueLInJC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jet->Pt()};
2146         fhnV0InJetLambda[iCentIndex]->Fill(valueLInJC);
2147         fh2V0PtJetAngleLambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2148       }
2149       if(bIsOutsideCones)
2150       {
2151         Double_t valueLOutJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2152         fhnV0OutJetLambda[iCentIndex]->Fill(valueLOutJet);
2153       }
2154       if(bIsInConePerp)
2155       {
2156         Double_t valueLInPC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jetPerp->Pt()};
2157         fhnV0InPerpLambda[iCentIndex]->Fill(valueLInPC);
2158       }
2159       if(bIsInConeRnd)
2160       {
2161         Double_t valueLInRnd[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2162         fhnV0InRndLambda[iCentIndex]->Fill(valueLInRnd);
2163       }
2164       if(bIsInConeMed)
2165       {
2166         Double_t valueLInMed[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2167         fhnV0InMedLambda[iCentIndex]->Fill(valueLInMed);
2168       }
2169       if(!iNJetSel)
2170       {
2171         Double_t valueLNoJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2172         fhnV0NoJetLambda[iCentIndex]->Fill(valueLNoJet);
2173       }
2174       iNV0CandLambda++;
2175     }
2176     if(bIsCandidateALambda)
2177     {
2178       // 14 ALambda candidates after cuts
2179 //          printf("AL: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0ALambda,dPtV0,dEtaV0,dPhiV0);
2180       FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex, iCentIndex);
2181       Double_t valueALIncl[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2182       fhnV0InclusiveALambda[iCentIndex]->Fill(valueALIncl);
2183       fh1V0InvMassALambdaCent[iCentIndex]->Fill(dMassV0ALambda);
2184       if(iNJetSel)
2185       {
2186         // 15 ALambda in jet events
2187         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 1, iCentIndex);
2188       }
2189       if(bIsInConeJet)
2190       {
2191         // 16 ALambda in jets
2192         FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 2, iCentIndex);
2193         Double_t valueLInJC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jet->Pt()};
2194         fhnV0InJetALambda[iCentIndex]->Fill(valueLInJC);
2195         fh2V0PtJetAngleALambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2196       }
2197       if(bIsOutsideCones)
2198       {
2199         Double_t valueALOutJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2200         fhnV0OutJetALambda[iCentIndex]->Fill(valueALOutJet);
2201       }
2202       if(bIsInConePerp)
2203       {
2204         Double_t valueLInPC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jetPerp->Pt()};
2205         fhnV0InPerpALambda[iCentIndex]->Fill(valueLInPC);
2206       }
2207       if(bIsInConeRnd)
2208       {
2209         Double_t valueALInRnd[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2210         fhnV0InRndALambda[iCentIndex]->Fill(valueALInRnd);
2211       }
2212       if(bIsInConeMed)
2213       {
2214         Double_t valueALInMed[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2215         fhnV0InMedALambda[iCentIndex]->Fill(valueALInMed);
2216       }
2217       if(!iNJetSel)
2218       {
2219         Double_t valueALNoJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2220         fhnV0NoJetALambda[iCentIndex]->Fill(valueALNoJet);
2221       }
2222       iNV0CandALambda++;
2223     }
2224     //===== End of filling V0 spectra =====
2225
2226
2227     //===== Association of reconstructed V0 candidates with MC particles =====
2228     if(fbMCAnalysis)
2229     {
2230       // Associate selected candidates only
2231 //          if ( !(bIsCandidateK0s && bIsInPeakK0s) && !(bIsCandidateLambda && bIsInPeakLambda) ) // signal candidates
2232       if(!(bIsCandidateK0s) && !(bIsCandidateLambda)  && !(bIsCandidateALambda))    // chosen candidates with any mass
2233         continue;
2234
2235       // Get MC labels of reconstructed daughter tracks
2236       Int_t iLabelPos = TMath::Abs(trackPos->GetLabel());
2237       Int_t iLabelNeg = TMath::Abs(trackNeg->GetLabel());
2238
2239       // Make sure MC daughters are in the array range
2240       if((iLabelNeg < 0) || (iLabelNeg >= iNTracksMC) || (iLabelPos < 0) || (iLabelPos >= iNTracksMC))
2241         continue;
2242
2243       // Get MC particles corresponding to reconstructed daughter tracks
2244       AliAODMCParticle* particleMCDaughterNeg = (AliAODMCParticle*)arrayMC->At(iLabelNeg);
2245       AliAODMCParticle* particleMCDaughterPos = (AliAODMCParticle*)arrayMC->At(iLabelPos);
2246       if(!particleMCDaughterNeg || !particleMCDaughterPos)
2247         continue;
2248
2249       // Make sure MC daughter particles are not physical primary
2250       if((particleMCDaughterNeg->IsPhysicalPrimary()) || (particleMCDaughterPos->IsPhysicalPrimary()))
2251         continue;
2252
2253       // Get identities of MC daughter particles
2254       Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
2255       Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
2256
2257       // Get index of the mother particle for each MC daughter particle
2258       Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
2259       Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
2260
2261       if((iIndexMotherNeg < 0) || (iIndexMotherNeg >= iNTracksMC) || (iIndexMotherPos < 0) || (iIndexMotherPos >= iNTracksMC))
2262         continue;
2263
2264       // Check whether MC daughter particles have the same mother
2265       if(iIndexMotherNeg != iIndexMotherPos)
2266         continue;
2267
2268       // Get the MC mother particle of both MC daughter particles
2269       AliAODMCParticle* particleMCMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherPos);
2270       if(!particleMCMother)
2271         continue;
2272
2273       // Get identity of the MC mother particle
2274       Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
2275
2276       // Skip not interesting particles
2277       if((iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda))
2278         continue;
2279
2280       // Check identity of the MC mother particle and the decay channel
2281       // Is MC mother particle K0S?
2282       Bool_t bV0MCIsK0s = ((iPdgCodeMother == iPdgCodeK0s) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2283       // Is MC mother particle Lambda?
2284       Bool_t bV0MCIsLambda = ((iPdgCodeMother == +iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodeProton) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2285       // Is MC mother particle anti-Lambda?
2286       Bool_t bV0MCIsALambda = ((iPdgCodeMother == -iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodeProton));
2287
2288       Double_t dPtV0Gen = particleMCMother->Pt();
2289 //          Double_t dRapV0MC = particleMCMother->Y();
2290       Double_t dEtaV0Gen = particleMCMother->Eta();
2291 //          Double_t dPhiV0Gen = particleMCMother->Phi();
2292
2293       // Is MC mother particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2294 //          Bool_t bV0MCIsPrimary = particleMCMother->IsPhysicalPrimary();
2295       // Get the MC mother particle of the MC mother particle
2296       Int_t iIndexMotherOfMother = particleMCMother->GetMother();
2297       AliAODMCParticle* particleMCMotherOfMother = 0;
2298       if(iIndexMotherOfMother >= 0)
2299         particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2300       // Get identity of the MC mother particle of the MC mother particle if it exists
2301       Int_t iPdgCodeMotherOfMother = 0;
2302       if(particleMCMotherOfMother)
2303         iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2304       // Check if the MC mother particle of the MC mother particle is a physical primary Sigma (3212 - Sigma0, 3224 - Sigma*+, 3214 - Sigma*0, 3114 - Sigma*-)
2305 //          Bool_t bV0MCComesFromSigma = kFALSE; // Is MC mother particle daughter of a Sigma?
2306 //          if ( (particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && ( (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) ) )
2307 //            bV0MCComesFromSigma = kTRUE;
2308       // Should MC mother particle be considered as primary when it is Lambda?
2309 //          Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2310       // Check if the MC mother particle of the MC mother particle is a Xi (3322 - Xi0, 3312 - Xi-)
2311       Bool_t bV0MCComesFromXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == 3322) || (iPdgCodeMotherOfMother == 3312))); // Is MC mother particle daughter of a Xi?
2312       Bool_t bV0MCComesFromAXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == -3322) || (iPdgCodeMotherOfMother == -3312))); // Is MC mother particle daughter of a anti-Xi?
2313
2314       // Get the distance between production point of the MC mother particle and the primary vertex
2315       Double_t dx = dPrimVtxMCX - particleMCMother->Xv();
2316       Double_t dy = dPrimVtxMCY - particleMCMother->Yv();
2317       Double_t dz = dPrimVtxMCZ - particleMCMother->Zv();
2318       Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2319       Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2320
2321       // phi, eta resolution for V0-reconstruction
2322 //          Double_t dResolutionV0Eta = particleMCMother->Eta()-v0->Eta();
2323 //          Double_t dResolutionV0Phi = particleMCMother->Phi()-v0->Phi();
2324
2325       // K0s
2326 //          if (bIsCandidateK0s && bIsInPeakK0s) // selected candidates in peak
2327       if(bIsCandidateK0s)  // selected candidates with any mass
2328       {
2329 //              if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2330         if(bV0MCIsK0s && bV0MCIsPrimaryDist)  // well reconstructed candidates
2331         {
2332           fh2V0K0sPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0K0s);
2333           Double_t valueEtaK[3] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen};
2334           fh3V0K0sEtaPtMassMCRec[iCentIndex]->Fill(valueEtaK);
2335
2336           Double_t valueEtaDKNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2337           fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKNeg);
2338           Double_t valueEtaDKPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2339           fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKPos);
2340
2341           fh2V0K0sMCResolMPt[iCentIndex]->Fill(dMassV0K0s - dMassPDGK0s, dPtV0);
2342           fh2V0K0sMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2343           if(bIsInConeJet)  // true V0 associated to a candidate in jet
2344           {
2345             Double_t valueKInJCMC[4] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2346             fh3V0K0sInJetPtMassMCRec[iCentIndex]->Fill(valueKInJCMC);
2347             Double_t valueEtaKIn[5] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2348             fh4V0K0sInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaKIn);
2349
2350             Double_t valueEtaDKJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2351             fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCNeg);
2352             Double_t valueEtaDKJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2353             fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCPos);
2354           }
2355         }
2356         if(bV0MCIsK0s && !bV0MCIsPrimaryDist)  // not primary K0s
2357         {
2358           fh1V0K0sPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2359         }
2360       }
2361       // Lambda
2362 //          if (bIsCandidateLambda && bIsInPeakLambda) // selected candidates in peak
2363       if(bIsCandidateLambda)  // selected candidates with any mass
2364       {
2365 //              if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2366         if(bV0MCIsLambda && bV0MCIsPrimaryDist)  // well reconstructed candidates
2367         {
2368           fh2V0LambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0Lambda);
2369           Double_t valueEtaL[3] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen};
2370           fh3V0LambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaL);
2371
2372           Double_t valueEtaDLNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2373           fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLNeg);
2374           Double_t valueEtaDLPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2375           fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLPos);
2376
2377           fh2V0LambdaMCResolMPt[iCentIndex]->Fill(dMassV0Lambda - dMassPDGLambda, dPtV0);
2378           fh2V0LambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2379           if(bIsInConeJet)  // true V0 associated to a reconstructed candidate in jet
2380           {
2381             Double_t valueLInJCMC[4] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2382             fh3V0LambdaInJetPtMassMCRec[iCentIndex]->Fill(valueLInJCMC);
2383             Double_t valueEtaLIn[5] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2384             fh4V0LambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaLIn);
2385
2386             Double_t valueEtaDLJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2387             fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCNeg);
2388             Double_t valueEtaDLJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2389             fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCPos);
2390           }
2391         }
2392         // Fill the feed-down histograms
2393         if(bV0MCIsLambda && bV0MCComesFromXi)
2394         {
2395           Double_t valueFDLIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2396           fhnV0LambdaInclMCFD[iCentIndex]->Fill(valueFDLIncl);
2397           if(bIsInConeRnd)
2398           {
2399             fhnV0LambdaBulkMCFD[iCentIndex]->Fill(valueFDLIncl);
2400           }
2401           if(bIsInConeJet)
2402           {
2403             Double_t valueFDLInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2404             fhnV0LambdaInJetsMCFD[iCentIndex]->Fill(valueFDLInJets);
2405           }
2406         }
2407         if(bV0MCIsLambda && !bV0MCIsPrimaryDist && !bV0MCComesFromXi)  // not primary Lambda
2408         {
2409           fh1V0LambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2410         }
2411       }
2412       // anti-Lambda
2413 //          if (bIsCandidateALambda && bIsInPeakALambda) // selected candidates in peak
2414       if(bIsCandidateALambda)  // selected candidates with any mass
2415       {
2416 //              if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2417         if(bV0MCIsALambda && bV0MCIsPrimaryDist)  // well reconstructed candidates
2418         {
2419           fh2V0ALambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0ALambda);
2420           Double_t valueEtaAL[3] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen};
2421           fh3V0ALambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaAL);
2422
2423           Double_t valueEtaDALNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2424           fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALNeg);
2425           Double_t valueEtaDALPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2426           fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALPos);
2427
2428           fh2V0ALambdaMCResolMPt[iCentIndex]->Fill(dMassV0ALambda - dMassPDGLambda, dPtV0);
2429           fh2V0ALambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2430           if(bIsInConeJet)  // true V0 associated to a reconstructed candidate in jet
2431           {
2432             Double_t valueALInJCMC[4] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2433             fh3V0ALambdaInJetPtMassMCRec[iCentIndex]->Fill(valueALInJCMC);
2434             Double_t valueEtaALIn[5] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2435             fh4V0ALambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaALIn);
2436
2437             Double_t valueEtaDALJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2438             fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCNeg);
2439             Double_t valueEtaDALJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2440             fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCPos);
2441           }
2442         }
2443         // Fill the feed-down histograms
2444         if(bV0MCIsALambda && bV0MCComesFromAXi)
2445         {
2446           Double_t valueFDALIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2447           fhnV0ALambdaInclMCFD[iCentIndex]->Fill(valueFDALIncl);
2448           if(bIsInConeRnd)
2449           {
2450             fhnV0ALambdaBulkMCFD[iCentIndex]->Fill(valueFDALIncl);
2451           }
2452           if(bIsInConeJet)
2453           {
2454             Double_t valueFDALInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2455             fhnV0ALambdaInJetsMCFD[iCentIndex]->Fill(valueFDALInJets);
2456           }
2457         }
2458         if(bV0MCIsALambda && !bV0MCIsPrimaryDist && !bV0MCComesFromAXi)  // not primary anti-Lambda
2459         {
2460           fh1V0ALambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2461         }
2462       }
2463     }
2464     //===== End Association of reconstructed V0 candidates with MC particles =====
2465   }
2466   //===== End of V0 loop =====
2467
2468   fh1V0CandPerEvent->Fill(iNV0CandTot);
2469   fh1V0CandPerEventCentK0s[iCentIndex]->Fill(iNV0CandK0s);
2470   fh1V0CandPerEventCentLambda[iCentIndex]->Fill(iNV0CandLambda);
2471   fh1V0CandPerEventCentALambda[iCentIndex]->Fill(iNV0CandALambda);
2472
2473   if(fDebug > 2) printf("TaskV0sInJetsEmcal: End of V0 loop\n");
2474
2475   // Spectra of generated particles
2476   if(fbMCAnalysis)
2477   {
2478     for(Int_t iPartMC = 0; iPartMC < iNTracksMC; iPartMC++)
2479     {
2480       // Get MC particle
2481       AliAODMCParticle* particleMC = (AliAODMCParticle*)arrayMC->At(iPartMC);
2482       if(!particleMC)
2483         continue;
2484
2485       // Get identity of MC particle
2486       Int_t iPdgCodeParticleMC = particleMC->GetPdgCode();
2487       // Fill Xi spectrum (3322 - Xi0, 3312 - Xi-)
2488 //          if ( (iPdgCodeParticleMC==3322) || (iPdgCodeParticleMC==3312) )
2489       if((iPdgCodeParticleMC == 3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2490       {
2491         fh1V0XiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2492       }
2493       if((iPdgCodeParticleMC == -3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2494       {
2495         fh1V0AXiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2496       }
2497       // Skip not interesting particles
2498       if((iPdgCodeParticleMC != iPdgCodeK0s) && (TMath::Abs(iPdgCodeParticleMC) != iPdgCodeLambda))
2499         continue;
2500
2501       // Check identity of the MC V0 particle
2502       // Is MC V0 particle K0S?
2503       Bool_t bV0MCIsK0s = (iPdgCodeParticleMC == iPdgCodeK0s);
2504       // Is MC V0 particle Lambda?
2505       Bool_t bV0MCIsLambda = (iPdgCodeParticleMC == +iPdgCodeLambda);
2506       // Is MC V0 particle anti-Lambda?
2507       Bool_t bV0MCIsALambda = (iPdgCodeParticleMC == -iPdgCodeLambda);
2508
2509       Double_t dPtV0Gen = particleMC->Pt();
2510       Double_t dRapV0Gen = particleMC->Y();
2511       Double_t dEtaV0Gen = particleMC->Eta();
2512
2513       // V0 rapidity cut
2514       if(bCutRapV0)
2515       {
2516         if(bPrintCuts) printf("Gen: Applying cut: V0 |y|: < %f\n", dRapMax);
2517         if((TMath::Abs(dRapV0Gen) > dRapMax))
2518           continue;
2519       }
2520       // V0 pseudorapidity cut
2521       if(bCutEtaV0)
2522       {
2523         if(bPrintCuts) printf("Gen: Applying cut: V0 |eta|: < %f\n", dEtaMax);
2524         if((TMath::Abs(dEtaV0Gen) > dEtaMax))
2525           continue;
2526       }
2527       /*
2528                 // Is MC V0 particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2529                 Bool_t bV0MCIsPrimary = particleMC->IsPhysicalPrimary();
2530
2531                 // Get the MC mother particle of the MC V0 particle
2532                 Int_t iIndexMotherOfMother = particleMC->GetMother();
2533                 AliAODMCParticle* particleMCMotherOfMother = 0;
2534                 if (iIndexMotherOfMother >= 0)
2535                   particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2536                 // Get identity of the MC mother particle of the MC V0 particle if it exists
2537                 Int_t iPdgCodeMotherOfMother = 0;
2538                 if (particleMCMotherOfMother)
2539                   iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2540                 // Check if the MC mother particle is a physical primary Sigma
2541                 Bool_t bV0MCComesFromSigma = kFALSE;
2542                 if ((particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) )
2543                   bV0MCComesFromSigma = kTRUE;
2544                 // Should the MC V0 particle be considered as primary when it is Lambda?
2545                 Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2546       */
2547       // Reject non primary particles
2548 //          if (!bV0MCIsPrimaryLambda)
2549 //            continue;
2550
2551       // Get the distance between the production point of the MC V0 particle and the primary vertex
2552       Double_t dx = dPrimVtxMCX - particleMC->Xv();
2553       Double_t dy = dPrimVtxMCY - particleMC->Yv();
2554       Double_t dz = dPrimVtxMCZ - particleMC->Zv();
2555       Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2556       Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2557
2558       // Check whether the MC V0 particle is in a MC jet
2559       AliAODJet* jetMC = 0;
2560       Bool_t bIsMCV0InJet = kFALSE;
2561       if(iNJetSel)
2562       {
2563         if(fDebug > 5) printf("TaskV0sInJetsEmcal: Searching for gen V0 in %d MC jets\n", iNJetSel);
2564         for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2565         {
2566           jetMC = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2567           if(fDebug > 5) printf("TaskV0sInJetsEmcal: Checking if gen V0 in MC jet %d\n", iJet);
2568           if(IsParticleInCone(particleMC, jetMC, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2569           {
2570             if(fDebug > 5) printf("TaskV0sInJetsEmcal: gen V0 found in MC jet %d\n", iJet);
2571             bIsMCV0InJet = kTRUE;
2572             break;
2573           }
2574         }
2575       }
2576
2577       // Select only primary-like MC V0 particles
2578       // K0s
2579 //          if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2580       if(bV0MCIsK0s && bV0MCIsPrimaryDist)  // well reconstructed candidates
2581       {
2582         fh1V0K0sPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2583         fh2V0K0sEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2584         if(bIsMCV0InJet)
2585         {
2586           fh2V0K0sInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2587           Double_t valueEtaKInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2588           fh3V0K0sInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaKInGen);
2589         }
2590       }
2591       // Lambda
2592 //          if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2593       if(bV0MCIsLambda && bV0MCIsPrimaryDist)  // well reconstructed candidates
2594       {
2595         fh1V0LambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2596         fh2V0LambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2597         if(bIsMCV0InJet)
2598         {
2599           fh2V0LambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2600           Double_t valueEtaLInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2601           fh3V0LambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaLInGen);
2602         }
2603       }
2604       // anti-Lambda
2605 //          if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2606       if(bV0MCIsALambda && bV0MCIsPrimaryDist)  // well reconstructed candidates
2607       {
2608         fh1V0ALambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2609         fh2V0ALambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2610         if(bIsMCV0InJet)
2611         {
2612           fh2V0ALambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2613           Double_t valueEtaALInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2614           fh3V0ALambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaALInGen);
2615         }
2616       }
2617     }
2618   }
2619
2620   jetArraySel->Delete();
2621   delete jetArraySel;
2622   jetArrayPerp->Delete();
2623   delete jetArrayPerp;
2624   if(jetRnd)
2625     delete jetRnd;
2626   jetRnd = 0;
2627
2628   PostData(1, fOutputListStd);
2629   PostData(2, fOutputListQA);
2630   PostData(3, fOutputListCuts);
2631   PostData(4, fOutputListMC);
2632 //  if(fDebug>5) printf("TaskV0sInJetsEmcal: FillHistograms: End\n");
2633   return kFALSE;
2634 }
2635
2636 void AliAnalysisTaskV0sInJetsEmcal::FillQAHistogramV0(AliAODVertex* vtx, const AliAODv0* vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda)
2637 {
2638   if(!IsCandK0s && !IsCandLambda)
2639     return;
2640
2641 //  Double_t dMassK0s = vZero->MassK0Short();
2642 //  Double_t dMassLambda = vZero->MassLambda();
2643
2644   fh1QAV0Status[iIndexHisto]->Fill(vZero->GetOnFlyStatus());
2645
2646   AliAODTrack* trackNeg = (AliAODTrack*)vZero->GetDaughter(1); // negative track
2647   AliAODTrack* trackPos = (AliAODTrack*)vZero->GetDaughter(0); // positive track
2648
2649   Short_t fTotalCharge = 0;
2650   for(Int_t i = 0; i < 2; i++)
2651   {
2652     AliAODTrack* track = (AliAODTrack*)vZero->GetDaughter(i); // track
2653     // Tracks TPC OK
2654     fh1QAV0TPCRefit[iIndexHisto]->Fill(track->IsOn(AliAODTrack::kTPCrefit));
2655     Double_t nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
2656     fh1QAV0TPCRows[iIndexHisto]->Fill(nCrossedRowsTPC);
2657     Int_t findable = track->GetTPCNclsF();
2658     fh1QAV0TPCFindable[iIndexHisto]->Fill(findable);
2659     if(findable != 0)
2660     {
2661       fh1QAV0TPCRowsFind[iIndexHisto]->Fill(nCrossedRowsTPC / findable);
2662     }
2663     // Daughters: pseudo-rapidity cut
2664     fh1QAV0Eta[iIndexHisto]->Fill(track->Eta());
2665     if((nCrossedRowsTPC > (160. / (250. - 85.) * (255.*TMath::Abs(tan(track->Theta())) - 85.)) + 20.) && (track->Eta() < 0) && (track->Pt() > 0.15))
2666 //      if (IsCandK0s)
2667     {
2668       fh2QAV0EtaRows[iIndexHisto]->Fill(track->Eta(), nCrossedRowsTPC);
2669       fh2QAV0PtRows[iIndexHisto]->Fill(track->Pt(), nCrossedRowsTPC);
2670       fh2QAV0PhiRows[iIndexHisto]->Fill(track->Phi(), nCrossedRowsTPC);
2671       fh2QAV0NClRows[iIndexHisto]->Fill(findable, nCrossedRowsTPC);
2672       fh2QAV0EtaNCl[iIndexHisto]->Fill(track->Eta(), findable);
2673     }
2674
2675     // Daughters: transverse momentum cut
2676     fh1QAV0Pt[iIndexHisto]->Fill(track->Pt());
2677     fTotalCharge += track->Charge();
2678   }
2679   fh1QAV0Charge[iIndexHisto]->Fill(fTotalCharge);
2680
2681   // Daughters: Impact parameter of daughters to prim vtx
2682   fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaNegToPrimVertex()));
2683   fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaPosToPrimVertex()));
2684 //  fh2CutDCAVtx[iIndexHisto]->Fill(dMassK0s,TMath::Abs(vZero->DcaNegToPrimVertex()));
2685
2686   // Daughters: DCA
2687   fh1QAV0DCAV0[iIndexHisto]->Fill(vZero->DcaV0Daughters());
2688 //  fh2CutDCAV0[iIndexHisto]->Fill(dMassK0s,vZero->DcaV0Daughters());
2689
2690   // V0: Cosine of the pointing angle
2691   fh1QAV0Cos[iIndexHisto]->Fill(vZero->CosPointingAngle(vtx));
2692 //  fh2CutCos[iIndexHisto]->Fill(dMassK0s,vZero->CosPointingAngle(vtx));
2693
2694   // V0: Fiducial volume
2695   Double_t xyz[3];
2696   vZero->GetSecondaryVtx(xyz);
2697   Double_t r2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
2698   fh1QAV0R[iIndexHisto]->Fill(TMath::Sqrt(r2));
2699
2700   Double_t dAlpha = vZero->AlphaV0();
2701   Double_t dPtArm = vZero->PtArmV0();
2702
2703   if(IsCandK0s)
2704   {
2705     if(IsInPeakK0s)
2706     {
2707 //          fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2708 //          fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2709       fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2710       fh2QAV0PtPtK0sPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2711       fh2ArmPodK0s[iIndexHisto]->Fill(dAlpha, dPtArm);
2712     }
2713     fh2QAV0EtaEtaK0s[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2714     fh2QAV0PhiPhiK0s[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2715     fh1QAV0RapK0s[iIndexHisto]->Fill(vZero->RapK0Short());
2716   }
2717
2718   if(IsCandLambda)
2719   {
2720     if(IsInPeakLambda)
2721     {
2722 //          fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2723 //          fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2724       fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2725       fh2QAV0PtPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2726       fh2ArmPodLambda[iIndexHisto]->Fill(dAlpha, dPtArm);
2727     }
2728     fh2QAV0EtaEtaLambda[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2729     fh2QAV0PhiPhiLambda[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2730     fh1QAV0RapLambda[iIndexHisto]->Fill(vZero->RapLambda());
2731   }
2732
2733   fh2ArmPod[iIndexHisto]->Fill(dAlpha, dPtArm);
2734
2735 }
2736
2737 void AliAnalysisTaskV0sInJetsEmcal::FillCandidates(Double_t mK, Double_t mL, Double_t mAL, Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut/*cut index*/, Int_t iCent/*cent index*/)
2738 {
2739   if(isK)
2740   {
2741     fh1V0CounterCentK0s[iCent]->Fill(iCut);
2742     fh1V0InvMassK0sAll[iCut]->Fill(mK);
2743   }
2744   if(isL)
2745   {
2746     fh1V0CounterCentLambda[iCent]->Fill(iCut);
2747     fh1V0InvMassLambdaAll[iCut]->Fill(mL);
2748   }
2749   if(isAL)
2750   {
2751     fh1V0CounterCentALambda[iCent]->Fill(iCut);
2752     fh1V0InvMassALambdaAll[iCut]->Fill(mAL);
2753   }
2754 }
2755
2756 Bool_t AliAnalysisTaskV0sInJetsEmcal::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
2757 {
2758 // decides whether a particle is inside a jet cone
2759   if(!part1 || !part2)
2760     return kFALSE;
2761
2762   TVector3 vecMom2(part2->Px(), part2->Py(), part2->Pz());
2763   TVector3 vecMom1(part1->Px(), part1->Py(), part1->Pz());
2764   Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
2765   if(dR < dRMax) // momentum vectors of part1 and part2 are closer than dRMax
2766     return kTRUE;
2767   return kFALSE;
2768 }
2769
2770 Bool_t AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets(const TClonesArray* array, const AliVParticle* part, Double_t dDistance) const
2771 {
2772 // decides whether a cone overlaps with other jets
2773   if(!part)
2774   {
2775     if(fDebug > 0) printf("AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets: Error: No part\n");
2776     return kFALSE;
2777   }
2778   if(!array)
2779   {
2780     if(fDebug > 0) printf("AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets: Error: No array\n");
2781     return kFALSE;
2782   }
2783   Int_t iNJets = array->GetEntriesFast();
2784   if(iNJets <= 0)
2785   {
2786     if(fDebug > 2) printf("AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets: Warning: No jets\n");
2787     return kFALSE;
2788   }
2789   AliVParticle* jet = 0;
2790   for(Int_t iJet = 0; iJet < iNJets; iJet++)
2791   {
2792     jet = (AliVParticle*)array->At(iJet);
2793     if(!jet)
2794     {
2795       if(fDebug > 0) printf("AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets: Error: Failed to load jet %d/%d\n", iJet, iNJets);
2796       continue;
2797     }
2798     if(IsParticleInCone(part, jet, dDistance))
2799       return kTRUE;
2800   }
2801   return kFALSE;
2802 }
2803
2804 AliAODJet* AliAnalysisTaskV0sInJetsEmcal::GetRandomCone(const TClonesArray* array, Double_t dEtaConeMax, Double_t dDistance) const
2805 {
2806 // generate a random cone which does not overlap with selected jets
2807 //  printf("Generating random cone...\n");
2808   TLorentzVector vecCone;
2809   AliAODJet* part = 0;
2810   Double_t dEta, dPhi;
2811   Int_t iNTrialsMax = 10;
2812   Bool_t bStatus = kFALSE;
2813   for(Int_t iTry = 0; iTry < iNTrialsMax; iTry++)
2814   {
2815 //      printf("Try %d\n",iTry);
2816     dEta = dEtaConeMax * (2 * fRandom->Rndm() - 1.); // random eta in [-dEtaConeMax,+dEtaConeMax]
2817     dPhi = TMath::TwoPi() * fRandom->Rndm(); // random phi in [0,2*Pi]
2818     vecCone.SetPtEtaPhiM(1., dEta, dPhi, 0.);
2819     part = new AliAODJet(vecCone);
2820     if(!OverlapWithJets(array, part, dDistance))
2821     {
2822       bStatus = kTRUE;
2823 //          printf("Success\n");
2824       break;
2825     }
2826     else
2827       delete part;
2828   }
2829   if(!bStatus)
2830     part = 0;
2831   return part;
2832 }
2833
2834 AliEmcalJet* AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster(AliJetContainer* cont, Double_t dEtaConeMax) const
2835 {
2836 // sort kt clusters by pT/area and return the middle one, based on code in AliAnalysisTaskJetChem
2837   if(!cont)
2838   {
2839     if(fDebug > 0) printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Error: No container\n");
2840     return NULL;
2841   }
2842   Int_t iNClTot = cont->GetNJets(); // number of all clusters in the array
2843   Int_t iNCl = 0; // number of accepted clusters
2844
2845   // get list of densities
2846   std::vector<std::vector<Double_t> > vecListClusters; // vector that contains pairs [ index, density ]
2847 //  printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Loop over %d clusters.\n", iNClTot);
2848   for(Int_t ij = 0; ij < iNClTot; ij++)
2849   {
2850     AliEmcalJet* clusterBg = (AliEmcalJet*)(cont->GetAcceptJet(ij));
2851     if(!clusterBg)
2852       continue;
2853 //    printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Cluster %d/%d used as accepted cluster %d.\n", ij, iNClTot, int(vecListClusters.size()));
2854     Double_t dPtBg = clusterBg->Pt();
2855     Double_t dAreaBg = clusterBg->Area();
2856     Double_t dDensityBg = 0;
2857     if(dAreaBg > 0)
2858       dDensityBg = dPtBg / dAreaBg;
2859     std::vector<Double_t> vecCluster;
2860     vecCluster.push_back(ij);
2861     vecCluster.push_back(dDensityBg);
2862     vecListClusters.push_back(vecCluster);
2863   }
2864   iNCl = vecListClusters.size();
2865   if(iNCl < 3) // need at least 3 clusters (skipping 2 highest)
2866   {
2867 //    if(fDebug > 2) printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Warning: Too little clusters\n");
2868     return NULL;
2869   }
2870
2871 //  printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Original lists:\n");
2872 //  for(Int_t i = 0; i < iNCl; i++)
2873 //    printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
2874
2875   // sort list of indeces by density in descending order
2876   std::sort(vecListClusters.begin(), vecListClusters.end(), CompareClusters);
2877
2878 //  printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Sorted lists:\n");
2879 //  for(Int_t i = 0; i < iNCl; i++)
2880 //    printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
2881
2882   // get median cluster with median density
2883   AliEmcalJet* clusterMed = 0;
2884   Int_t iIndex = 0; // index of the median cluster in the sorted list
2885   Int_t iIndexMed = 0; // index of the median cluster in the original array
2886   if(TMath::Odd(iNCl))  // odd number of clusters
2887   {
2888     iIndex = (Int_t)(0.5 * (iNCl + 1)); // = (n - skip + 1)/2 + 1, skip = 2
2889 //    printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Odd, median index = %d/%d\n", iIndex, iNCl);
2890   }
2891   else // even number: picking randomly one of the two closest to median
2892   {
2893     Int_t iIndex1 = (Int_t)(0.5 * iNCl); // = (n - skip)/2 + 1, skip = 2
2894     Int_t iIndex2 = (Int_t)(0.5 * iNCl + 1); // = (n - skip)/2 + 1 + 1, skip = 2
2895     iIndex = ((fRandom->Rndm() > 0.5) ? iIndex1 : iIndex2);
2896 //    printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: Even, median index = %d or %d -> %d/%d\n", iIndex1, iIndex2, iIndex, iNCl);
2897   }
2898   iIndexMed = Int_t((vecListClusters[iIndex])[0]);
2899
2900 //  printf("AliAnalysisTaskV0sInJetsEmcal::GetMedianCluster: getting median cluster %d/%d ok, rho = %g\n", iIndexMed, iNClTot, (vecListClusters[iIndex])[1]);
2901   clusterMed = (AliEmcalJet*)(cont->GetAcceptJet(iIndexMed));
2902
2903   if(TMath::Abs(clusterMed->Eta()) > dEtaConeMax)
2904     return NULL;
2905
2906   return clusterMed;
2907 }
2908
2909 Double_t AliAnalysisTaskV0sInJetsEmcal::AreaCircSegment(Double_t dRadius, Double_t dDistance) const
2910 {
2911 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
2912   Double_t dEpsilon = 1e-2;
2913   Double_t dR = dRadius;
2914   Double_t dD = dDistance;
2915   if(TMath::Abs(dR) < dEpsilon)
2916   {
2917     if(fDebug > 0) printf("AliAnalysisTaskV0sInJetsEmcal::AreaCircSegment: Error: Too small radius: %f < %f\n", dR, dEpsilon);
2918     return 0.;
2919   }
2920   if(dD >= dR)
2921     return 0.;
2922   if(dD <= -dR)
2923     return TMath::Pi() * dR * dR;
2924   return dR * dR * TMath::ACos(dD / dR) - dD * TMath::Sqrt(dR * dR - dD * dD);
2925 }
2926
2927 Bool_t AliAnalysisTaskV0sInJetsEmcal::IsSelectedForJets(AliAODEvent* fAOD, Double_t dVtxZCut, Double_t dVtxR2Cut, Double_t dCentCutLo, Double_t dCentCutUp, Bool_t bCutDeltaZ, Double_t dDeltaZMax)
2928 {
2929 // event selection
2930   AliAODVertex* vertex = fAOD->GetPrimaryVertex();
2931   if(!vertex)
2932     return kFALSE;
2933   Int_t iNContribMin = 3;
2934   if(!fbIsPbPb)
2935     iNContribMin = 2;
2936   if(vertex->GetNContributors() < iNContribMin)
2937     return kFALSE;
2938   TString vtxTitle(vertex->GetTitle());
2939   if(vtxTitle.Contains("TPCVertex"))
2940     return kFALSE;
2941   Double_t zVertex = vertex->GetZ();
2942   if(TMath::Abs(zVertex) > dVtxZCut)
2943     return kFALSE;
2944   if(bCutDeltaZ)
2945   {
2946     AliAODVertex* vertexSPD = fAOD->GetPrimaryVertexSPD();
2947     if(!vertexSPD)
2948     {
2949 //          printf("IsSelectedForJets: Error: No SPD vertex\n");
2950       return kFALSE;
2951     }
2952     Double_t zVertexSPD = vertexSPD->GetZ();
2953     if(TMath::Abs(zVertex - zVertexSPD) > dDeltaZMax)
2954     {
2955 //          printf("IsSelectedForJets: Rejecting event due to delta z = %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
2956       return kFALSE;
2957     }
2958 //      printf("IsSelectedForJets: Event OK: %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
2959   }
2960   Double_t xVertex = vertex->GetX();
2961   Double_t yVertex = vertex->GetY();
2962   Double_t radiusSq = yVertex * yVertex + xVertex * xVertex;
2963   if(radiusSq > dVtxR2Cut)
2964     return kFALSE;
2965   Double_t centrality;
2966 //  centrality = fAOD->GetHeader()->GetCentrality();
2967   centrality = fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentile("V0M");
2968   if(fbIsPbPb)
2969   {
2970     if(centrality < 0)
2971       return kFALSE;
2972     if((dCentCutUp < 0) || (dCentCutLo < 0) || (dCentCutUp > 100) || (dCentCutLo > 100) || (dCentCutLo > dCentCutUp))
2973       return kFALSE;
2974     if((centrality < dCentCutLo) || (centrality > dCentCutUp))
2975       return kFALSE;
2976   }
2977   else
2978   {
2979     if(centrality != -1)
2980       return kFALSE;
2981   }
2982   return kTRUE;
2983 }
2984
2985 Int_t AliAnalysisTaskV0sInJetsEmcal::GetCentralityBinIndex(Double_t centrality)
2986 {
2987 // returns index of the centrality bin corresponding to the provided value of centrality
2988   if(centrality < 0 || centrality > fgkiCentBinRanges[fgkiNBinsCent - 1])
2989     return -1;
2990   for(Int_t i = 0; i < fgkiNBinsCent; i++)
2991   {
2992     if(centrality <= fgkiCentBinRanges[i])
2993       return i;
2994   }
2995   return -1;
2996 }
2997
2998 Int_t AliAnalysisTaskV0sInJetsEmcal::GetCentralityBinEdge(Int_t index)
2999 {
3000 // returns the upper edge of the centrality bin corresponding to the provided value of index
3001   if(index < 0 || index >= fgkiNBinsCent)
3002     return -1;
3003   return fgkiCentBinRanges[index];
3004 }
3005
3006 TString AliAnalysisTaskV0sInJetsEmcal::GetCentBinLabel(Int_t index)
3007 {
3008 // get string with centrality range for given bin
3009   TString lowerEdge = ((index == 0) ? "0" : Form("%d", GetCentralityBinEdge(index - 1)));
3010   TString upperEdge = Form("%d", GetCentralityBinEdge(index));
3011   return Form("%s-%s %%", lowerEdge.Data(), upperEdge.Data());
3012 }
3013
3014 Double_t AliAnalysisTaskV0sInJetsEmcal::MassPeakSigmaOld(Double_t pt, Int_t particle)
3015 {
3016 // estimation of the sigma of the invariant-mass peak as a function of pT and particle type
3017   switch(particle)
3018   {
3019     case 0: // K0S
3020       return 0.0044 + 0.0004 * (pt - 1.);
3021       break;
3022     case 1: // Lambda
3023       return 0.0023 + 0.00034 * (pt - 1.);
3024       break;
3025     default:
3026       return 0;
3027       break;
3028   }
3029 }
3030
3031 bool AliAnalysisTaskV0sInJetsEmcal::CompareClusters(const std::vector<Double_t> cluster1, const std::vector<Double_t> cluster2)
3032 {
3033   return (cluster1[1] > cluster2[1]);
3034 }