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