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