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