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