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