AliPhysics  b752f14 (b752f14)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskV0sInJetsEmcal.cxx
Go to the documentation of this file.
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 //-------------------------------------------------------------------------
17 // task for analysis of V0s (K0S, (anti-)Lambda) in charged jets
18 // fork of AliAnalysisTaskV0sInJets for the EMCal framework
19 // Author: Vit Kucera (vit.kucera@cern.ch)
20 //-------------------------------------------------------------------------
21 
22 #include <TChain.h>
23 #include <TH1D.h>
24 #include <TH2D.h>
25 #include <THnSparse.h>
26 #include <TClonesArray.h>
27 #include <TRandom3.h>
28 #include <TDatabasePDG.h>
29 #include <TPDGCode.h>
30 
31 #include "AliAnalysisManager.h"
32 #include "AliInputEventHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODMCHeader.h"
35 #include "AliLog.h"
36 #include "AliEventPoolManager.h"
37 #include "AliPIDResponse.h"
38 #include "AliAODTrack.h"
39 #include "AliAODMCParticle.h"
40 #include "AliMCEvent.h"
41 
42 #include "AliEmcalJet.h"
43 #include "AliJetContainer.h"
44 #include "AliParticleContainer.h"
45 //#include "AliClusterContainer.h"
46 //#include "AliPicoTrack.h"
47 
49 
51 
52 // upper edges of centrality bins
53 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[] = {10, 30, 50, 80}; // Alice Zimmermann
54 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[] = {10, 20, 40, 60, 80}; // Vit Kucera, initial binning
55 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[] = {5, 10, 20, 40, 60, 80}; // Iouri Belikov, LF analysis
56 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiCentBinRanges[] = {10}; // only central
57 // edges of centrality bins for event mixing
59 // edges of z_vtx bins for event mixing
60 Double_t AliAnalysisTaskV0sInJetsEmcal::fgkiZVtxMixBinRanges[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10};
61 // axis: pT of V0
66 // axis: pT of jets
68 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJet = sizeof(AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet) / sizeof(AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet[0]) - 1;
69 //const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJet] - (AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[0]) / 5.); // bin width 5 GeV/c
70 const Int_t AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[AliAnalysisTaskV0sInJetsEmcal::fgkiNBinsPtJet] - (AliAnalysisTaskV0sInJetsEmcal::fgkdBinsPtJet)[0]) / 10.); // bin width 10 GeV/c
71 // axis: K0S invariant mass
75 // axis: Lambda invariant mass
79 // delta phi range
80 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdDeltaPhiMin = -TMath::PiOver2(); // minimum delta-phi_V0-jet
81 const Double_t AliAnalysisTaskV0sInJetsEmcal::fgkdDeltaPhiMax = 3.*TMath::PiOver2(); // maximum delta-phi_V0-jet
82 
83 
84 // Default constructor
87  fAODIn(0),
88  fAODOut(0),
89  fEventMC(0),
90  fRandom(0),
91  fPoolMgr(0),
92  fOutputListStd(0),
93  fOutputListQA(0),
94  fOutputListCuts(0),
95  fOutputListMC(0),
96 
97  fbIsPbPb(0),
98  fbMCAnalysis(0),
99  fsGeneratorName(""),
100 
101  fdCutVertexZ(0),
102  fdCutVertexR2(0),
103  fdCutCentLow(0),
104  fdCutCentHigh(0),
105  fdCutDeltaZMax(0),
106  fdCentrality(0),
107 
108  fbCorrelations(0),
109  fiSizePool(0),
110  fiNJetsPerPool(0),
111  ffFractionMin(0),
112  fiNEventsMin(0),
113  fdDeltaEtaMax(0),
114 
115  fbTPCRefit(0),
116  fbRejectKinks(0),
117  fbFindableClusters(0),
118  fdCutNCrossedRowsTPCMin(0),
119  fdCutCrossedRowsOverFindMin(0),
120  fdCutCrossedRowsOverFindMax(0),
121  fdCutPtDaughterMin(0),
122  fdCutDCAToPrimVtxMin(0),
123  fdCutDCADaughtersMax(0),
124  fdCutEtaDaughterMax(0),
125  fdCutNSigmadEdxMax(0),
126  fdPtProtonPIDMax(0),
127  fbOnFly(0),
128  fdCutCPAKMin(0),
129  fdCutCPALMin(0),
130  fdCutRadiusDecayMin(0),
131  fdCutRadiusDecayMax(0),
132  fdCutEtaV0Max(0),
133  fdCutRapV0Max(0),
134  fdCutNTauKMax(0),
135  fdCutNTauLMax(0),
136  fbCutArmPod(0),
137  fbCutCross(0),
138 
139  fbJetSelection(0),
140  fdCutPtJetMin(0),
141  fdCutPtTrackJetMin(0),
142  fdCutAreaPercJetMin(0),
143  fdDistanceV0JetMax(0.4),
144  fiBgSubtraction(0),
145 
146  fbCompareTriggers(0),
147 
148  fJetsCont(0),
149  fJetsBgCont(0),
150  fTracksCont(0),
151 
152  fh1EventCounterCut(0),
153  fh1EventCent(0),
154  fh1EventCent2(0),
155  fh1EventCent2Jets(0),
156  fh1EventCent2NoJets(0),
157  fh2EventCentTracks(0),
158  fh1V0CandPerEvent(0),
159  fh1NRndConeCent(0),
160  fh1NMedConeCent(0),
161  fh1AreaExcluded(0),
162 
163  fh2CCK0s(0),
164  fh2CCLambda(0),
165  fh3CCMassCorrelBoth(0),
166  fh3CCMassCorrelKNotL(0),
167  fh3CCMassCorrelLNotK(0)
168 {
169  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
170  {
171  fh1QAV0Status[i] = 0;
172  fh1QAV0TPCRefit[i] = 0;
173  fh1QAV0TPCRows[i] = 0;
174  fh1QAV0TPCFindable[i] = 0;
175  fh2QAV0PtNCls[i] = 0;
176  fh2QAV0PtChi[i] = 0;
177  fh1QAV0TPCRowsFind[i] = 0;
178  fh1QAV0Eta[i] = 0;
179  fh2QAV0EtaRows[i] = 0;
180  fh2QAV0PtRows[i] = 0;
181  fh2QAV0PhiRows[i] = 0;
182  fh2QAV0NClRows[i] = 0;
183  fh2QAV0EtaNCl[i] = 0;
184 
185  fh2QAV0EtaPtK0sPeak[i] = 0;
186  fh2QAV0EtaEtaK0s[i] = 0;
187  fh2QAV0PhiPhiK0s[i] = 0;
188  fh1QAV0RapK0s[i] = 0;
189  fh2QAV0PtPtK0sPeak[i] = 0;
190  fh2ArmPodK0s[i] = 0;
191 
192  fh2QAV0EtaPtLambdaPeak[i] = 0;
193  fh2QAV0EtaEtaLambda[i] = 0;
194  fh2QAV0PhiPhiLambda[i] = 0;
195  fh1QAV0RapLambda[i] = 0;
196  fh2QAV0PtPtLambdaPeak[i] = 0;
197  fh2ArmPodLambda[i] = 0;
198 
200  fh2QAV0EtaEtaALambda[i] = 0;
201  fh2QAV0PhiPhiALambda[i] = 0;
202  fh1QAV0RapALambda[i] = 0;
203  fh2QAV0PtPtALambdaPeak[i] = 0;
204  fh2ArmPodALambda[i] = 0;
205 
206  fh2QAV0PhiPtK0sPeak[i] = 0;
207  fh2QAV0PhiPtLambdaPeak[i] = 0;
209  fh1QAV0Pt[i] = 0;
210  fh1QAV0Charge[i] = 0;
211  fh1QAV0DCAVtx[i] = 0;
212  fh1QAV0DCAV0[i] = 0;
213  fh1QAV0Cos[i] = 0;
214  fh1QAV0R[i] = 0;
215  fh1QACTau2D[i] = 0;
216  fh1QACTau3D[i] = 0;
217 
218  fh2ArmPod[i] = 0;
219 
220  /*
221  fh2CutTPCRowsK0s[i] = 0;
222  fh2CutTPCRowsLambda[i] = 0;
223  fh2CutPtPosK0s[i] = 0;
224  fh2CutPtNegK0s[i] = 0;
225  fh2CutPtPosLambda[i] = 0;
226  fh2CutPtNegLambda[i] = 0;
227  fh2CutDCAVtx[i] = 0;
228  fh2CutDCAV0[i] = 0;
229  fh2CutCos[i] = 0;
230  fh2CutR[i] = 0;
231  fh2CutEtaK0s[i] = 0;
232  fh2CutEtaLambda[i] = 0;
233  fh2CutRapK0s[i] = 0;
234  fh2CutRapLambda[i] = 0;
235  fh2CutCTauK0s[i] = 0;
236  fh2CutCTauLambda[i] = 0;
237  fh2CutPIDPosK0s[i] = 0;
238  fh2CutPIDNegK0s[i] = 0;
239  fh2CutPIDPosLambda[i] = 0;
240  fh2CutPIDNegLambda[i] = 0;
241 
242  fh2Tau3DVs2D[i] = 0;
243  */
244  }
245  for(Int_t i = 0; i < fgkiNCategV0; i++)
246  {
247  fh1V0InvMassK0sAll[i] = 0;
248  fh1V0InvMassLambdaAll[i] = 0;
249  fh1V0InvMassALambdaAll[i] = 0;
250  }
251  for(Int_t i = 0; i < fgkiNBinsCent; i++)
252  {
253  fh1EventCounterCutCent[i] = 0;
254  fh1V0CounterCentK0s[i] = 0;
255  fh1V0CounterCentLambda[i] = 0;
260  fhnPtDaughterK0s[i] = 0;
261  fhnPtDaughterLambda[i] = 0;
262  fhnPtDaughterALambda[i] = 0;
263  fh1V0InvMassK0sCent[i] = 0;
264  fh1V0InvMassLambdaCent[i] = 0;
266  fh1V0K0sPtMCGen[i] = 0;
267  fh2V0K0sPtMassMCRec[i] = 0;
268  fh1V0K0sPtMCRecFalse[i] = 0;
269  fh2V0K0sEtaPtMCGen[i] = 0;
270  fh3V0K0sEtaPtMassMCRec[i] = 0;
271  fh2V0K0sInJetPtMCGen[i] = 0;
275  fh2V0K0sMCResolMPt[i] = 0;
276  fh2V0K0sMCPtGenPtRec[i] = 0;
277  fh1V0LambdaPtMCGen[i] = 0;
278  fh2V0LambdaPtMassMCRec[i] = 0;
280  fh2V0LambdaEtaPtMCGen[i] = 0;
286  fh2V0LambdaMCResolMPt[i] = 0;
288  fhnV0LambdaInclMCFD[i] = 0;
289  fhnV0LambdaInJetsMCFD[i] = 0;
290  fhnV0LambdaBulkMCFD[i] = 0;
291  fh1V0XiPtMCGen[i] = 0;
292  fh1V0ALambdaPtMCGen[i] = 0;
295  fh2V0ALambdaEtaPtMCGen[i] = 0;
301  fh2V0ALambdaMCResolMPt[i] = 0;
303  fhnV0ALambdaInclMCFD[i] = 0;
304  fhnV0ALambdaInJetsMCFD[i] = 0;
305  fhnV0ALambdaBulkMCFD[i] = 0;
306  fh1V0AXiPtMCGen[i] = 0;
307 
308  // eta daughters
315 
316  // Inclusive
317  fhnV0InclusiveK0s[i] = 0;
318  fhnV0InclusiveLambda[i] = 0;
319  fhnV0InclusiveALambda[i] = 0;
320  // Cones
321  fhnV0InJetK0s[i] = 0;
322  fhnV0InPerpK0s[i] = 0;
323  fhnV0InRndK0s[i] = 0;
324  fhnV0InMedK0s[i] = 0;
325  fhnV0OutJetK0s[i] = 0;
326  fhnV0NoJetK0s[i] = 0;
327  fhnV0InJetLambda[i] = 0;
328  fhnV0InPerpLambda[i] = 0;
329  fhnV0InRndLambda[i] = 0;
330  fhnV0InMedLambda[i] = 0;
331  fhnV0OutJetLambda[i] = 0;
332  fhnV0NoJetLambda[i] = 0;
333  fhnV0InJetALambda[i] = 0;
334  fhnV0InPerpALambda[i] = 0;
335  fhnV0InRndALambda[i] = 0;
336  fhnV0InMedALambda[i] = 0;
337  fhnV0OutJetALambda[i] = 0;
338  fhnV0NoJetALambda[i] = 0;
339  // correlations
340  fhnV0CorrelSEK0s[i] = 0;
341  fhnV0CorrelMEK0s[i] = 0;
342  fhnV0CorrelSELambda[i] = 0;
343  fhnV0CorrelMELambda[i] = 0;
344 
345  fh2V0PtJetAngleK0s[i] = 0;
346  fh2V0PtJetAngleLambda[i] = 0;
347  fh2V0PtJetAngleALambda[i] = 0;
348 
349  fh1PtJet[i] = 0;
350  fh1EtaJet[i] = 0;
351  fh2EtaPtJet[i] = 0;
352  fh1PhiJet[i] = 0;
353  fh2PtJetPtTrackLeading[i] = 0;
354  fh2PtJetPtTrigger[i] = 0;
355  fh1PtTrigger[i] = 0;
356  fh1NJetPerEvent[i] = 0;
357  fh2EtaPhiRndCone[i] = 0;
358  fh2EtaPhiMedCone[i] = 0;
359  fh1DistanceJets[i] = 0;
360  fh1DistanceV0JetsK0s[i] = 0;
363 
364  fh1VtxZ[i] = 0;
365  fh1VtxZME[i] = 0;
366  fh2VtxXY[i] = 0;
367  }
368 }
369 
370 // Constructor
372  AliAnalysisTaskEmcalJet(name, kTRUE),
373  fAODIn(0),
374  fAODOut(0),
375  fEventMC(0),
376  fRandom(0),
377  fPoolMgr(0),
378  fOutputListStd(0),
379  fOutputListQA(0),
380  fOutputListCuts(0),
381  fOutputListMC(0),
382 
383  fbIsPbPb(0),
384  fbMCAnalysis(0),
385  fsGeneratorName(""),
386 
387  fdCutVertexZ(0),
388  fdCutVertexR2(0),
389  fdCutCentLow(0),
390  fdCutCentHigh(0),
391  fdCutDeltaZMax(0),
392  fdCentrality(0),
393 
394  fbCorrelations(0),
395  fiSizePool(0),
396  fiNJetsPerPool(0),
397  ffFractionMin(0),
398  fiNEventsMin(0),
399  fdDeltaEtaMax(0),
400 
401  fbTPCRefit(0),
402  fbRejectKinks(0),
403  fbFindableClusters(0),
404  fdCutNCrossedRowsTPCMin(0),
405  fdCutCrossedRowsOverFindMin(0),
406  fdCutCrossedRowsOverFindMax(0),
407  fdCutPtDaughterMin(0),
408  fdCutDCAToPrimVtxMin(0),
409  fdCutDCADaughtersMax(0),
410  fdCutEtaDaughterMax(0),
411  fdCutNSigmadEdxMax(0),
412  fdPtProtonPIDMax(0),
413  fbOnFly(0),
414  fdCutCPAKMin(0),
415  fdCutCPALMin(0),
416  fdCutRadiusDecayMin(0),
417  fdCutRadiusDecayMax(0),
418  fdCutEtaV0Max(0),
419  fdCutRapV0Max(0),
420  fdCutNTauKMax(0),
421  fdCutNTauLMax(0),
422  fbCutArmPod(0),
423  fbCutCross(0),
424 
425  fbJetSelection(0),
426  fdCutPtJetMin(0),
427  fdCutPtTrackJetMin(0),
428  fdCutAreaPercJetMin(0),
429  fdDistanceV0JetMax(0.4),
430  fiBgSubtraction(0),
431 
432  fbCompareTriggers(0),
433 
434  fJetsCont(0),
435  fJetsBgCont(0),
436  fTracksCont(0),
437 
438  fh1EventCounterCut(0),
439  fh1EventCent(0),
440  fh1EventCent2(0),
441  fh1EventCent2Jets(0),
442  fh1EventCent2NoJets(0),
443  fh2EventCentTracks(0),
444  fh1V0CandPerEvent(0),
445  fh1NRndConeCent(0),
446  fh1NMedConeCent(0),
447  fh1AreaExcluded(0),
448 
449  fh2CCK0s(0),
450  fh2CCLambda(0),
451  fh3CCMassCorrelBoth(0),
452  fh3CCMassCorrelKNotL(0),
453  fh3CCMassCorrelLNotK(0)
454 {
455  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
456  {
457  fh1QAV0Status[i] = 0;
458  fh1QAV0TPCRefit[i] = 0;
459  fh1QAV0TPCRows[i] = 0;
460  fh1QAV0TPCFindable[i] = 0;
461  fh2QAV0PtNCls[i] = 0;
462  fh2QAV0PtChi[i] = 0;
463  fh1QAV0TPCRowsFind[i] = 0;
464  fh1QAV0Eta[i] = 0;
465  fh2QAV0EtaRows[i] = 0;
466  fh2QAV0PtRows[i] = 0;
467  fh2QAV0PhiRows[i] = 0;
468  fh2QAV0NClRows[i] = 0;
469  fh2QAV0EtaNCl[i] = 0;
470 
471  fh2QAV0EtaPtK0sPeak[i] = 0;
472  fh2QAV0EtaEtaK0s[i] = 0;
473  fh2QAV0PhiPhiK0s[i] = 0;
474  fh1QAV0RapK0s[i] = 0;
475  fh2QAV0PtPtK0sPeak[i] = 0;
476  fh2ArmPodK0s[i] = 0;
477 
478  fh2QAV0EtaPtLambdaPeak[i] = 0;
479  fh2QAV0EtaEtaLambda[i] = 0;
480  fh2QAV0PhiPhiLambda[i] = 0;
481  fh1QAV0RapLambda[i] = 0;
482  fh2QAV0PtPtLambdaPeak[i] = 0;
483  fh2ArmPodLambda[i] = 0;
484 
486  fh2QAV0EtaEtaALambda[i] = 0;
487  fh2QAV0PhiPhiALambda[i] = 0;
488  fh1QAV0RapALambda[i] = 0;
489  fh2QAV0PtPtALambdaPeak[i] = 0;
490  fh2ArmPodALambda[i] = 0;
491 
492  fh2QAV0PhiPtK0sPeak[i] = 0;
493  fh2QAV0PhiPtLambdaPeak[i] = 0;
495  fh1QAV0Pt[i] = 0;
496  fh1QAV0Charge[i] = 0;
497  fh1QAV0DCAVtx[i] = 0;
498  fh1QAV0DCAV0[i] = 0;
499  fh1QAV0Cos[i] = 0;
500  fh1QAV0R[i] = 0;
501  fh1QACTau2D[i] = 0;
502  fh1QACTau3D[i] = 0;
503 
504  fh2ArmPod[i] = 0;
505 
506  /*
507  fh2CutTPCRowsK0s[i] = 0;
508  fh2CutTPCRowsLambda[i] = 0;
509  fh2CutPtPosK0s[i] = 0;
510  fh2CutPtNegK0s[i] = 0;
511  fh2CutPtPosLambda[i] = 0;
512  fh2CutPtNegLambda[i] = 0;
513  fh2CutDCAVtx[i] = 0;
514  fh2CutDCAV0[i] = 0;
515  fh2CutCos[i] = 0;
516  fh2CutR[i] = 0;
517  fh2CutEtaK0s[i] = 0;
518  fh2CutEtaLambda[i] = 0;
519  fh2CutRapK0s[i] = 0;
520  fh2CutRapLambda[i] = 0;
521  fh2CutCTauK0s[i] = 0;
522  fh2CutCTauLambda[i] = 0;
523  fh2CutPIDPosK0s[i] = 0;
524  fh2CutPIDNegK0s[i] = 0;
525  fh2CutPIDPosLambda[i] = 0;
526  fh2CutPIDNegLambda[i] = 0;
527 
528  fh2Tau3DVs2D[i] = 0;
529  */
530  }
531  for(Int_t i = 0; i < fgkiNCategV0; i++)
532  {
533  fh1V0InvMassK0sAll[i] = 0;
534  fh1V0InvMassLambdaAll[i] = 0;
535  fh1V0InvMassALambdaAll[i] = 0;
536  }
537  for(Int_t i = 0; i < fgkiNBinsCent; i++)
538  {
539  fh1EventCounterCutCent[i] = 0;
540  fh1V0CounterCentK0s[i] = 0;
541  fh1V0CounterCentLambda[i] = 0;
546  fhnPtDaughterK0s[i] = 0;
547  fhnPtDaughterLambda[i] = 0;
548  fhnPtDaughterALambda[i] = 0;
549  fh1V0InvMassK0sCent[i] = 0;
550  fh1V0InvMassLambdaCent[i] = 0;
552  fh1V0K0sPtMCGen[i] = 0;
553  fh2V0K0sPtMassMCRec[i] = 0;
554  fh1V0K0sPtMCRecFalse[i] = 0;
555  fh2V0K0sEtaPtMCGen[i] = 0;
556  fh3V0K0sEtaPtMassMCRec[i] = 0;
557  fh2V0K0sInJetPtMCGen[i] = 0;
561  fh2V0K0sMCResolMPt[i] = 0;
562  fh2V0K0sMCPtGenPtRec[i] = 0;
563  fh1V0LambdaPtMCGen[i] = 0;
564  fh2V0LambdaPtMassMCRec[i] = 0;
566  fh2V0LambdaEtaPtMCGen[i] = 0;
572  fh2V0LambdaMCResolMPt[i] = 0;
574  fhnV0LambdaInclMCFD[i] = 0;
575  fhnV0LambdaInJetsMCFD[i] = 0;
576  fhnV0LambdaBulkMCFD[i] = 0;
577  fh1V0XiPtMCGen[i] = 0;
578  fh1V0ALambdaPtMCGen[i] = 0;
581  fh2V0ALambdaEtaPtMCGen[i] = 0;
587  fh2V0ALambdaMCResolMPt[i] = 0;
589  fhnV0ALambdaInclMCFD[i] = 0;
590  fhnV0ALambdaInJetsMCFD[i] = 0;
591  fhnV0ALambdaBulkMCFD[i] = 0;
592  fh1V0AXiPtMCGen[i] = 0;
593 
594  // eta daughters
601 
602  // Inclusive
603  fhnV0InclusiveK0s[i] = 0;
604  fhnV0InclusiveLambda[i] = 0;
605  fhnV0InclusiveALambda[i] = 0;
606  // Cones
607  fhnV0InJetK0s[i] = 0;
608  fhnV0InPerpK0s[i] = 0;
609  fhnV0InRndK0s[i] = 0;
610  fhnV0InMedK0s[i] = 0;
611  fhnV0OutJetK0s[i] = 0;
612  fhnV0NoJetK0s[i] = 0;
613  fhnV0InJetLambda[i] = 0;
614  fhnV0InPerpLambda[i] = 0;
615  fhnV0InRndLambda[i] = 0;
616  fhnV0InMedLambda[i] = 0;
617  fhnV0OutJetLambda[i] = 0;
618  fhnV0NoJetLambda[i] = 0;
619  fhnV0InJetALambda[i] = 0;
620  fhnV0InPerpALambda[i] = 0;
621  fhnV0InRndALambda[i] = 0;
622  fhnV0InMedALambda[i] = 0;
623  fhnV0OutJetALambda[i] = 0;
624  fhnV0NoJetALambda[i] = 0;
625  // correlations
626  fhnV0CorrelSEK0s[i] = 0;
627  fhnV0CorrelMEK0s[i] = 0;
628  fhnV0CorrelSELambda[i] = 0;
629  fhnV0CorrelMELambda[i] = 0;
630 
631  fh2V0PtJetAngleK0s[i] = 0;
632  fh2V0PtJetAngleLambda[i] = 0;
633  fh2V0PtJetAngleALambda[i] = 0;
634 
635  fh1PtJet[i] = 0;
636  fh1EtaJet[i] = 0;
637  fh2EtaPtJet[i] = 0;
638  fh1PhiJet[i] = 0;
639  fh2PtJetPtTrackLeading[i] = 0;
640  fh2PtJetPtTrigger[i] = 0;
641  fh1PtTrigger[i] = 0;
642  fh1NJetPerEvent[i] = 0;
643  fh2EtaPhiRndCone[i] = 0;
644  fh2EtaPhiMedCone[i] = 0;
645  fh1DistanceJets[i] = 0;
646  fh1DistanceV0JetsK0s[i] = 0;
649 
650  fh1VtxZ[i] = 0;
651  fh1VtxZME[i] = 0;
652  fh2VtxXY[i] = 0;
653  }
654  // Define input and output slots here
655  // Input slot #0 works with a TChain
656  DefineInput(0, TChain::Class());
657  // Output slot #0 id reserved by the base class for AOD
658  // Output slot #1 writes into a TList container
659  DefineOutput(1, TList::Class());
660  DefineOutput(2, TList::Class());
661  DefineOutput(3, TList::Class());
662  DefineOutput(4, TList::Class());
663 }
664 
666 {
667  delete fRandom;
668  fRandom = 0;
669 }
670 
672 {
673  // Called once
674 
676 
677  if(fbJetSelection)
678  {
682  {
684  }
685  }
686 
687  // Initialise random-number generator
688  fRandom = new TRandom3(0);
689 
690  // Create event pool manager
691  if(fbCorrelations)
692  {
695  }
696 
697  // Create histograms
698 
699  fOutputListStd = new TList();
700  fOutputListStd->SetOwner();
701  fOutputListQA = new TList();
702  fOutputListQA->SetOwner();
703  fOutputListCuts = new TList();
704  fOutputListCuts->SetOwner();
705  fOutputListMC = new TList();
706  fOutputListMC->SetOwner();
707 
708  // event categories
709  const Int_t iNCategEvent = 12;
710  TString categEvent[iNCategEvent] = {
711  "coll. candid.",
712  "AOD OK",
713  "no SPD pile-up",
714  "PV contrib.",
715  "not TPC PV",
716  "|z| PV",
717  "|#Deltaz| SPD PV",
718  "r PV",
719  "cent.", //2
720  "with V0",
721  "with jets",
722  "jet selection"
723  };
724 
725  // labels for stages of V0 selection
726  TString categV0[fgkiNCategV0] = {
727  "all"/*0*/,
728  "mass range"/*1*/,
729  "rec. method"/*2*/,
730  "tracks TPC"/*3*/,
731  "track pt"/*4*/,
732  "DCA prim v"/*5*/,
733  "DCA daughters"/*6*/,
734  "CPA"/*7*/,
735  "volume"/*8*/,
736  "track #it{#eta}"/*9*/,
737  "V0 #it{y} & #it{#eta}"/*10*/,
738  "lifetime"/*11*/,
739  "PID"/*12*/,
740  "Arm.-Pod."/*13*/,
741  "cross-cont."/*14*/,
742  "inclusive"/*15*/,
743  "in jet event"/*16*/,
744  "in jet"/*17*/
745  };
746 
747  fh1EventCounterCut = new TH1D("fh1EventCounterCut", "Number of events after filtering;selection filter;counts", iNCategEvent, 0, iNCategEvent);
748  for(Int_t i = 0; i < iNCategEvent; i++)
749  fh1EventCounterCut->GetXaxis()->SetBinLabel(i + 1, categEvent[i].Data());
750  fh1EventCent2 = new TH1D("fh1EventCent2", "Number of events vs centrality;centrality;counts", 100, 0, 100);
751  fh1EventCent2Jets = new TH1D("fh1EventCent2Jets", "Number of sel.-jet events vs centrality;centrality;counts", 100, 0, 100);
752  fh1EventCent2NoJets = new TH1D("fh1EventCent2NoJets", "Number of no-jet events vs centrality;centrality;counts", 100, 0, 100);
753  fh2EventCentTracks = new TH2D("fh2EventCentTracks", "Number of tracks vs centrality;centrality;tracks;counts", 100, 0, 100, 150, 0, 15e3);
754  fh1EventCent = new TH1D("fh1EventCent", "Number of events in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
755  for(Int_t i = 0; i < fgkiNBinsCent; i++)
756  fh1EventCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
757  fh1NRndConeCent = new TH1D("fh1NRndConeCent", "Number of rnd. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
758  for(Int_t i = 0; i < fgkiNBinsCent; i++)
759  fh1NRndConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
760  fh1NMedConeCent = new TH1D("fh1NMedConeCent", "Number of med.-cl. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
761  for(Int_t i = 0; i < fgkiNBinsCent; i++)
762  fh1NMedConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
763  fh1AreaExcluded = new TH1D("fh1AreaExcluded", "Area of excluded cones in centrality bins;centrality;area", fgkiNBinsCent, 0, fgkiNBinsCent);
764  for(Int_t i = 0; i < fgkiNBinsCent; i++)
765  fh1AreaExcluded->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
775 
776  fh1V0CandPerEvent = new TH1D("fh1V0CandPerEvent", "Number of all V0 candidates per event;candidates;events", 300, 0, 30000);
778 
779  for(Int_t i = 0; i < fgkiNBinsCent; i++)
780  {
781  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);
782  for(Int_t j = 0; j < iNCategEvent; j++)
783  fh1EventCounterCutCent[i]->GetXaxis()->SetBinLabel(j + 1, categEvent[j].Data());
784  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, 200);
785  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, 200);
786  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, 200);
787  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);
788  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);
789  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);
790 
791  for(Int_t j = 0; j < fgkiNCategV0; j++)
792  {
793  fh1V0CounterCentK0s[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
794  fh1V0CounterCentLambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
795  fh1V0CounterCentALambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
796  }
804  }
805  // pt binning for V0 and jets
806  Int_t iNBinsPtV0 = fgkiNBinsPtV0Init;
807  Int_t iNBinsPtV0InJet = fgkiNBinsPtV0InitInJet;
808  Double_t dPtV0Min = fgkdBinsPtV0[0];
809  Double_t dPtV0Max = fgkdBinsPtV0[fgkiNBinsPtV0];
810  Int_t iNJetPtBins = fgkiNBinsPtJetInit;
811  Double_t dJetPtMin = fgkdBinsPtJet[0];
813 
814  fh2CCK0s = new TH2D("fh2CCK0s", "K0s candidates in Lambda peak", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
815  fh2CCLambda = new TH2D("fh2CCLambda", "Lambda candidates in K0s peak", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
816  Int_t binsCorrel[3] = {fgkiNBinsMassK0s, fgkiNBinsMassLambda, iNBinsPtV0};
817  Double_t xminCorrel[3] = {fgkdMassK0sMin, fgkdMassLambdaMin, dPtV0Min};
818  Double_t xmaxCorrel[3] = {fgkdMassK0sMax, fgkdMassLambdaMax, dPtV0Max};
819  fh3CCMassCorrelBoth = new THnSparseD("fh3CCMassCorrelBoth", "Mass correlation: K0S && Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
820  fh3CCMassCorrelKNotL = new THnSparseD("fh3CCMassCorrelKNotL", "Mass correlation: K0S, not Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
821  fh3CCMassCorrelLNotK = new THnSparseD("fh3CCMassCorrelLNotK", "Mass correlation: Lambda, not K0S;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
822  fOutputListQA->Add(fh2CCK0s);
827 
828  Double_t dStepEtaV0 = 0.05;
829  Double_t dRangeEtaV0Max = 0.8;
830  const Int_t iNBinsEtaV0 = 2 * Int_t(dRangeEtaV0Max / dStepEtaV0);
831 // printf("%s: %d\n", "iNBinsEtaV0", iNBinsEtaV0);
832  // inclusive
833  const Int_t iNDimIncl = 3;
834  Int_t binsKIncl[iNDimIncl] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
835  Double_t xminKIncl[iNDimIncl] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
836  Double_t xmaxKIncl[iNDimIncl] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
837  Int_t binsLIncl[iNDimIncl] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
838  Double_t xminLIncl[iNDimIncl] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
839  Double_t xmaxLIncl[iNDimIncl] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
840  // binning in jets
841  const Int_t iNDimInJC = 4;
842  Int_t binsKInJC[iNDimInJC] = {fgkiNBinsMassK0s, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins};
843  Double_t xminKInJC[iNDimInJC] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
844  Double_t xmaxKInJC[iNDimInJC] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
845  Int_t binsLInJC[iNDimInJC] = {fgkiNBinsMassLambda, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins};
846  Double_t xminLInJC[iNDimInJC] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
847  Double_t xmaxLInJC[iNDimInJC] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
848  // binning in V0-jet correlations
849  Int_t iNBinsDeltaPhi = 60;
850  const Int_t iNDimCorrel = 6;
851 // Double_t fdDeltaEtaMax = 2. * 0.7; // >= 2*eta_V0^max
852 // Int_t iNBinsDeltaEtaCorrel = 2 * Int_t(fdDeltaEtaMax / 0.1); // does not work if fdDeltaEtaMax is Float_t
853 // Int_t iNBinsDeltaEtaCorrel = 2 * Int_t(fdDeltaEtaMax / 0.7); // does not work if fdDeltaEtaMax is Float_t
854 // Int_t iNBinsDeltaEtaCorrel = 2 * Int_t(fdDeltaEtaMax / fdDeltaEtaMax); // does not work if fdDeltaEtaMax is Float_t
855  Int_t iNBinsDeltaEtaCorrel = 1;
856 // printf("%s: %d\n", "iNBinsDeltaEtaCorrel", iNBinsDeltaEtaCorrel);
857  Int_t binsKCorrel[iNDimCorrel] = {fgkiNBinsMassK0s / 2, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaPhi, iNBinsDeltaEtaCorrel};
858  Double_t xminKCorrel[iNDimCorrel] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, fgkdDeltaPhiMin, -fdDeltaEtaMax};
859  Double_t xmaxKCorrel[iNDimCorrel] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, fgkdDeltaPhiMax, fdDeltaEtaMax};
860  Int_t binsLCorrel[iNDimCorrel] = {fgkiNBinsMassLambda / 2, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaPhi, iNBinsDeltaEtaCorrel};
861  Double_t xminLCorrel[iNDimCorrel] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, fgkdDeltaPhiMin, -fdDeltaEtaMax};
862  Double_t xmaxLCorrel[iNDimCorrel] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, fgkdDeltaPhiMax, fdDeltaEtaMax};
863 // printf("binsKCorrel: %d %d %d %d %d %d\n", fgkiNBinsMassK0s / 2, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaPhi, iNBinsDeltaEtaCorrel);
864 
865  // binning eff inclusive vs eta-pT
866  Double_t dStepDeltaEta = 0.1;
867  Double_t dRangeDeltaEtaMax = 0.5;
868  const Int_t iNBinsDeltaEta = 2 * Int_t(dRangeDeltaEtaMax / dStepDeltaEta);
869 // printf("%s: %d\n", "iNBinsDeltaEta", iNBinsDeltaEta);
870  Int_t binsEtaK[3] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
871  Double_t xminEtaK[3] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
872  Double_t xmaxEtaK[3] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
873  Int_t binsEtaL[3] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
874  Double_t xminEtaL[3] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
875  Double_t xmaxEtaL[3] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
876  // binning eff in jets vs eta-pT
877  // associated
878  Int_t binsEtaKInRec[5] = {fgkiNBinsMassK0s, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
879  Double_t xminEtaKInRec[5] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
880  Double_t xmaxEtaKInRec[5] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
881  Int_t binsEtaLInRec[5] = {fgkiNBinsMassLambda, iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
882  Double_t xminEtaLInRec[5] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
883  Double_t xmaxEtaLInRec[5] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
884  // generated
885  Int_t binsEtaInGen[4] = {iNBinsPtV0InJet, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
886  Double_t xminEtaInGen[4] = {dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
887  Double_t xmaxEtaInGen[4] = {dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
888  // daughter eta: charge-etaD-ptD-etaV0-ptV0-ptJet
889  const Int_t iNDimEtaD = 6;
890  Int_t binsEtaDaughter[iNDimEtaD] = {2, 20, iNBinsPtV0, iNBinsEtaV0, iNBinsPtV0, iNJetPtBins};
891  Double_t xminEtaDaughter[iNDimEtaD] = {0, -1, dPtV0Min, -dRangeEtaV0Max, dPtV0Min, dJetPtMin};
892  Double_t xmaxEtaDaughter[iNDimEtaD] = {2, 1, dPtV0Max, dRangeEtaV0Max, dPtV0Max, dJetPtMax};
893  // daughter track pt: pos-neg-ptV0-ptJet-ptLead
894  const Int_t iNDimDaughter = 5;
895  Int_t binsDaughter[iNDimDaughter] = {80, 80, iNBinsPtV0, iNJetPtBins, 80};
896  Double_t xminDaughter[iNDimDaughter] = {0, 0, dPtV0Min, dJetPtMin, 0};
897  Double_t xmaxDaughter[iNDimDaughter] = {20, 20, dPtV0Max, dJetPtMax, 20};
898 
899  for(Int_t i = 0; i < fgkiNBinsCent; i++)
900  {
901  // Daughter track Pt vs pt
902  fhnPtDaughterK0s[i] = new THnSparseD(Form("fhnPtDaughterK0s_%d", i), Form("K0s: Daughter Pt vs Pt, cent: %s;neg #it{p}_{T} (GeV/#it{c});pos #it{p}_{T} (GeV/#it{c});p_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c});#it{p}_{T} leading track (GeV/#it{c});counts", GetCentBinLabel(i).Data()), iNDimDaughter, binsDaughter, xminDaughter, xmaxDaughter);
903  fhnPtDaughterLambda[i] = new THnSparseD(Form("fhnPtDaughterLambda_%d", i), Form("Lambda: Daughter Pt vs Pt, cent: %s;neg #it{p}_{T} (GeV/#it{c});pos #it{p}_{T} (GeV/#it{c});p_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c});#it{p}_{T} leading track (GeV/#it{c});counts", GetCentBinLabel(i).Data()), iNDimDaughter, binsDaughter, xminDaughter, xmaxDaughter);
904  fhnPtDaughterALambda[i] = new THnSparseD(Form("fhnPtDaughterALambda_%d", i), Form("ALambda: Daughter Pt vs Pt, cent: %s;neg #it{p}_{T} (GeV/#it{c});pos #it{p}_{T} (GeV/#it{c});p_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c});#it{p}_{T} leading track (GeV/#it{c});counts", GetCentBinLabel(i).Data()), iNDimDaughter, binsDaughter, xminDaughter, xmaxDaughter);
908 
909  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);
910  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);
911  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);
915  // Inclusive
916  fhnV0InclusiveK0s[i] = new THnSparseD(Form("fhnV0InclusiveK0s_C%d", i), "K0s: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};counts", iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
917  fhnV0InclusiveLambda[i] = new THnSparseD(Form("fhnV0InclusiveLambda_C%d", i), "Lambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
918  fhnV0InclusiveALambda[i] = new THnSparseD(Form("fhnV0InclusiveALambda_C%d", i), "ALambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
922  // In cones
923  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
924  fOutputListStd->Add(fhnV0InJetK0s[i]);
925  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
927  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
928  fOutputListStd->Add(fhnV0InRndK0s[i]);
929  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
930  fOutputListStd->Add(fhnV0InMedK0s[i]);
931  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
933  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
934  fOutputListStd->Add(fhnV0NoJetK0s[i]);
935  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
937  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
939  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
941  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
943  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
945  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
947  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
949  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{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
951  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
953  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
955  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
957  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});#it{#eta}_{V0}", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
959  // correlations
960  fhnV0CorrelSEK0s[i] = new THnSparseD(Form("fhnV0CorrelSEK0s_%d", i), Form("K0s: Correlations with jets in same events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c});#Delta#it{#phi}_{V0-jet};#Delta#it{#eta}_{V0-jet}", GetCentBinLabel(i).Data()), iNDimCorrel, binsKCorrel, xminKCorrel, xmaxKCorrel);
962  fhnV0CorrelMEK0s[i] = new THnSparseD(Form("fhnV0CorrelMEK0s_%d", i), Form("K0s: Correlations with jets in mixed events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c});#Delta#it{#phi}_{V0-jet};#Delta#it{#eta}_{V0-jet}", GetCentBinLabel(i).Data()), iNDimCorrel, binsKCorrel, xminKCorrel, xmaxKCorrel);
964  fhnV0CorrelSELambda[i] = new THnSparseD(Form("fhnV0CorrelSELambda_%d", i), Form("Lambda: Correlations with jets in same events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c});#Delta#it{#phi}_{V0-jet};#Delta#it{#eta}_{V0-jet}", GetCentBinLabel(i).Data()), iNDimCorrel, binsLCorrel, xminLCorrel, xmaxLCorrel);
966  fhnV0CorrelMELambda[i] = new THnSparseD(Form("fhnV0CorrelMELambda_%d", i), Form("Lambda: Correlations with jets in mixed events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c});#Delta#it{#phi}_{V0-jet};#Delta#it{#eta}_{V0-jet}", GetCentBinLabel(i).Data()), iNDimCorrel, binsLCorrel, xminLCorrel, xmaxLCorrel);
968 
969  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()), 2 * iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdDistanceV0JetMax + 0.1);
971  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()), 2 * iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdDistanceV0JetMax + 0.1);
973  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()), 2 * iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdDistanceV0JetMax + 0.1);
975 
976  // jet histograms
977  fh1PtJet[i] = new TH1D(Form("fh1PtJet_%d", i), Form("Jet pt spectrum, cent: %s;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), 2 * iNJetPtBins, dJetPtMin, dJetPtMax);
978  fOutputListStd->Add(fh1PtJet[i]);
979  fh1EtaJet[i] = new TH1D(Form("fh1EtaJet_%d", i), Form("Jet eta spectrum, cent: %s;#it{#eta} jet", GetCentBinLabel(i).Data()), 80, -1., 1.);
980  fOutputListStd->Add(fh1EtaJet[i]);
981  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., 2 * iNJetPtBins, dJetPtMin, dJetPtMax);
982  fOutputListStd->Add(fh2EtaPtJet[i]);
983  fh1PhiJet[i] = new TH1D(Form("fh1PhiJet_%d", i), Form("Jet phi spectrum, cent: %s;#it{#phi} jet", GetCentBinLabel(i).Data()), 90, 0., TMath::TwoPi());
984  fOutputListStd->Add(fh1PhiJet[i]);
985  fh2PtJetPtTrackLeading[i] = new TH2D(Form("fh2PtJetPtTrackLeading_%d", i), Form("jet pt vs leading track pt, cent: %s;#it{p}_{T}^{jet} (GeV/#it{c});#it{p}_{T} leading track (GeV/#it{c})", GetCentBinLabel(i).Data()), 4 * iNJetPtBins, dJetPtMin, dJetPtMax, 200, 0., 20);
987  fh2PtJetPtTrigger[i] = new TH2D(Form("fh2PtJetPtTrigger_%d", i), Form("jet pt vs trigger track pt, cent: %s;#it{p}_{T}^{jet} (GeV/#it{c});#it{p}_{T} trigger track (GeV/#it{c})", GetCentBinLabel(i).Data()), 4 * iNJetPtBins, dJetPtMin, dJetPtMax, 200, 0., 20);
989  fh1PtTrigger[i] = new TH1D(Form("fh1PtTrigger_%d", i), Form("trigger track pt, cent: %s;#it{p}_{T} trigger track (GeV/#it{c})", GetCentBinLabel(i).Data()), 200, 0., 20);
990  fOutputListStd->Add(fh1PtTrigger[i]);
991  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.);
993  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., 90, 0., TMath::TwoPi());
995  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., 90, 0., TMath::TwoPi());
997  fh1DistanceJets[i] = new TH1D(Form("fh1DistanceJets_%d", i), Form("Distance between jets in #eta-#phi, cent: %s;#it{D}", GetCentBinLabel(i).Data()), 40, 0., 4.);
999  fh1DistanceV0JetsK0s[i] = new TH1D(Form("fh1DistanceV0JetsK0s_%d", i), Form("Distance between V0 and the closest jet in #eta-#phi, cent: %s;#it{D}", GetCentBinLabel(i).Data()), 80, 0., 2.);
1000  fh1DistanceV0JetsLambda[i] = new TH1D(Form("fh1DistanceV0JetsLambda_%d", i), Form("Distance between V0 and the closest jet in #eta-#phi, cent: %s;#it{D}", GetCentBinLabel(i).Data()), 80, 0., 2.);
1001  fh1DistanceV0JetsALambda[i] = new TH1D(Form("fh1DistanceV0JetsALambda_%d", i), Form("Distance between V0 and the closest jet in #eta-#phi, cent: %s;#it{D}", GetCentBinLabel(i).Data()), 80, 0., 2.);
1005  // event histograms
1006  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);
1007  fOutputListQA->Add(fh1VtxZ[i]);
1008  fh1VtxZME[i] = new TH1D(Form("fh1VtxZME_%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);
1009  fOutputListQA->Add(fh1VtxZME[i]);
1010  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);
1011  fOutputListQA->Add(fh2VtxXY[i]);
1012  // fOutputListStd->Add([i]);
1013  if(fbMCAnalysis)
1014  {
1015  // inclusive pt
1016  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);
1017  fOutputListMC->Add(fh1V0K0sPtMCGen[i]);
1018  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);
1020  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);
1022  // inclusive pt-eta
1023  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);
1025  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);
1027  // in jet pt
1028  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()), iNBinsPtV0InJet, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
1030  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);
1032  // in jet pt-eta
1033  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);
1035  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);
1037 
1038  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);
1040  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);
1042 
1043  // inclusive pt
1044  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);
1046  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);
1048  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);
1050  // inclusive pt-eta
1051  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);
1053  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);
1055  // in jet pt
1056  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()), iNBinsPtV0InJet, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
1058  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);
1060  // in jet pt-eta
1061  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);
1063  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);
1065 
1066  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);
1068  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);
1070 
1071  // inclusive pt
1072  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);
1074  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);
1076  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);
1078  // inclusive pt-eta
1079  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);
1081  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);
1083  // in jet pt
1084  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()), iNBinsPtV0InJet, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
1086  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);
1088  // in jet pt-eta
1089  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);
1091  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);
1093 
1094  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);
1096  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);
1098 
1099  Int_t iNBinsPtXi = 80;
1100  Double_t dPtXiMin = 0;
1101  Double_t dPtXiMax = 8;
1102  const Int_t iNDimFD = 3;
1103  Int_t binsFD[iNDimFD] = {iNBinsPtV0, iNBinsPtXi, iNJetPtBins};
1104  Double_t xminFD[iNDimFD] = {dPtV0Min, dPtXiMin, dJetPtMin};
1105  Double_t xmaxFD[iNDimFD] = {dPtV0Max, dPtXiMax, dJetPtMax};
1106  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);
1108  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);
1110  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);
1112  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);
1113  fOutputListMC->Add(fh1V0XiPtMCGen[i]);
1114  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);
1116  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);
1118  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);
1120  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);
1121  fOutputListMC->Add(fh1V0AXiPtMCGen[i]);
1122 
1123  // daughter eta
1124  fhnV0K0sInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1125  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 gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1126  fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1127  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 gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1128  fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1129  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 gen daughter;pT gen daughter;eta gen V0;pT gen V0;pT rec jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1130 
1137  }
1138  }
1139 
1140  // QA Histograms
1141  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
1142  {
1143 // [i] = new TH1D(Form("%d",i),";;Counts",,,);
1144  fh1QAV0Status[i] = new TH1D(Form("fh1QAV0Status_%d", i), "QA: V0 status", 2, 0, 2);
1145  fh1QAV0TPCRefit[i] = new TH1D(Form("fh1QAV0TPCRefit_%d", i), "QA: TPC refit", 2, 0, 2);
1146  fh1QAV0TPCRows[i] = new TH1D(Form("fh1QAV0TPCRows_%d", i), "QA: TPC Rows", 160, 0, 160);
1147  fh1QAV0TPCFindable[i] = new TH1D(Form("fh1QAV0TPCFindable_%d", i), "QA: TPC Findable", 160, 0, 160);
1148  fh2QAV0PtNCls[i] = new TH2D(Form("fh2QAV0PtNCls_%d", i), "QA: Daughter Pt vs TPC clusters;pt;# TPC clusters", 100, 0, 10, 160, 0, 160);
1149  fh2QAV0PtChi[i] = new TH2D(Form("fh2QAV0PtChi_%d", i), "QA: Daughter Pt vs Chi2/ndf;pt;Chi2/ndf", 100, 0, 10, 100, 0, 100);
1150  fh1QAV0TPCRowsFind[i] = new TH1D(Form("fh1QAV0TPCRowsFind_%d", i), "QA: TPC Rows/Findable", 100, 0, 2);
1151  fh1QAV0Eta[i] = new TH1D(Form("fh1QAV0Eta_%d", i), "QA: Daughter Eta", 200, -2, 2);
1152  fh2QAV0EtaRows[i] = new TH2D(Form("fh2QAV0EtaRows_%d", i), "QA: Daughter Eta vs TPC rows;#eta;TPC rows", 200, -2, 2, 160, 0, 160);
1153  fh2QAV0PtRows[i] = new TH2D(Form("fh2QAV0PtRows_%d", i), "QA: Daughter Pt vs TPC rows;pt;TPC rows", 100, 0, 10, 160, 0, 160);
1154  fh2QAV0PhiRows[i] = new TH2D(Form("fh2QAV0PhiRows_%d", i), "QA: Daughter Phi vs TPC rows;#phi;TPC rows", 90, 0, TMath::TwoPi(), 160, 0, 160);
1155  fh2QAV0NClRows[i] = new TH2D(Form("fh2QAV0NClRows_%d", i), "QA: Daughter NCl vs TPC rows;findable clusters;TPC rows", 100, 0, 160, 160, 0, 160);
1156  fh2QAV0EtaNCl[i] = new TH2D(Form("fh2QAV0EtaNCl_%d", i), "QA: Daughter Eta vs NCl;#eta;findable clusters", 200, -2, 2, 160, 0, 160);
1157 
1158  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);
1159  fh2QAV0EtaEtaK0s[i] = new TH2D(Form("fh2QAV0EtaEtaK0s_%d", i), "QA: K0s: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1160  fh2QAV0PhiPhiK0s[i] = new TH2D(Form("fh2QAV0PhiPhiK0s_%d", i), "QA: K0s: Phi vs Phi Daughter", 90, 0, TMath::TwoPi(), 90, 0, TMath::TwoPi());
1161  fh1QAV0RapK0s[i] = new TH1D(Form("fh1QAV0RapK0s_%d", i), "QA: K0s: V0 Rapidity", 200, -2, 2);
1162  fh2QAV0PtPtK0sPeak[i] = new TH2D(Form("fh2QAV0PtPtK0sPeak_%d", i), "QA: K0s: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1163 
1164  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);
1165  fh2QAV0EtaEtaLambda[i] = new TH2D(Form("fh2QAV0EtaEtaLambda_%d", i), "QA: Lambda: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1166  fh2QAV0PhiPhiLambda[i] = new TH2D(Form("fh2QAV0PhiPhiLambda_%d", i), "QA: Lambda: Phi vs Phi Daughter", 90, 0, TMath::TwoPi(), 90, 0, TMath::TwoPi());
1167  fh1QAV0RapLambda[i] = new TH1D(Form("fh1QAV0RapLambda_%d", i), "QA: Lambda: V0 Rapidity", 200, -2, 2);
1168  fh2QAV0PtPtLambdaPeak[i] = new TH2D(Form("fh2QAV0PtPtLambdaPeak_%d", i), "QA: Lambda: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1169 
1170  fh2QAV0EtaPtALambdaPeak[i] = new TH2D(Form("fh2QAV0EtaPtALambdaPeak_%d", i), "QA: anti-Lambda: Daughter Eta vs V0 pt, peak;track eta;V0 pt", 200, -2, 2, iNBinsPtV0, dPtV0Min, dPtV0Max);
1171  fh2QAV0EtaEtaALambda[i] = new TH2D(Form("fh2QAV0EtaEtaALambda_%d", i), "QA: anti-Lambda: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1172  fh2QAV0PhiPhiALambda[i] = new TH2D(Form("fh2QAV0PhiPhiALambda_%d", i), "QA: anti-Lambda: Phi vs Phi Daughter", 90, 0, TMath::TwoPi(), 90, 0, TMath::TwoPi());
1173  fh1QAV0RapALambda[i] = new TH1D(Form("fh1QAV0RapALambda_%d", i), "QA: anti-Lambda: V0 Rapidity", 200, -2, 2);
1174  fh2QAV0PtPtALambdaPeak[i] = new TH2D(Form("fh2QAV0PtPtALambdaPeak_%d", i), "QA: anti-Lambda: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1175 
1176  fh2QAV0PhiPtK0sPeak[i] = new TH2D(Form("fh2QAV0PhiPtK0sPeak_%d", i), "QA: K0S: #phi-pt;#phi;pt", 90, 0, TMath::TwoPi(), iNBinsPtV0, dPtV0Min, dPtV0Max);
1177  fh2QAV0PhiPtLambdaPeak[i] = new TH2D(Form("fh2QAV0PhiPtLambdaPeak_%d", i), "QA: Lambda: #phi-pt;#phi;pt", 90, 0, TMath::TwoPi(), iNBinsPtV0, dPtV0Min, dPtV0Max);
1178  fh2QAV0PhiPtALambdaPeak[i] = new TH2D(Form("fh2QAV0PhiPtALambdaPeak_%d", i), "QA: anti-Lambda: #phi-pt;#phi;pt", 90, 0, TMath::TwoPi(), iNBinsPtV0, dPtV0Min, dPtV0Max);
1179  fh1QAV0Pt[i] = new TH1D(Form("fh1QAV0Pt_%d", i), "QA: Daughter Pt", 100, 0, 5);
1180  fh1QAV0Charge[i] = new TH1D(Form("fh1QAV0Charge_%d", i), "QA: V0 Charge", 3, -1, 2);
1181  fh1QAV0DCAVtx[i] = new TH1D(Form("fh1QAV0DCAVtx_%d", i), "QA: DCA daughters to primary vertex", 1000, 0, 10);
1182  fh1QAV0DCAV0[i] = new TH1D(Form("fh1QAV0DCAV0_%d", i), "QA: DCA daughters", 100, 0, 2);
1183  fh1QAV0Cos[i] = new TH1D(Form("fh1QAV0Cos_%d", i), "QA: CPA", 10000, 0.9, 1);
1184  fh1QAV0R[i] = new TH1D(Form("fh1QAV0R_%d", i), "QA: R", 1500, 0, 150);
1185  fh1QACTau2D[i] = new TH1D(Form("fh1QACTau2D_%d", i), "QA: K0s: c#tau 2D;mR/pt#tau", 100, 0, 10);
1186  fh1QACTau3D[i] = new TH1D(Form("fh1QACTau3D_%d", i), "QA: K0s: c#tau 3D;mL/p#tau", 100, 0, 10);
1187 
1188  fh2ArmPod[i] = new TH2D(Form("fh2ArmPod_%d", i), "Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1189  fh2ArmPodK0s[i] = new TH2D(Form("fh2ArmPodK0s_%d", i), "K0s: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1190  fh2ArmPodLambda[i] = new TH2D(Form("fh2ArmPodLambda_%d", i), "Lambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1191  fh2ArmPodALambda[i] = new TH2D(Form("fh2ArmPodALambda_%d", i), "ALambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1192 
1193  fOutputListQA->Add(fh1QAV0Status[i]);
1194  fOutputListQA->Add(fh1QAV0TPCRefit[i]);
1195  fOutputListQA->Add(fh1QAV0TPCRows[i]);
1197  fOutputListQA->Add(fh2QAV0PtNCls[i]);
1198  fOutputListQA->Add(fh2QAV0PtChi[i]);
1200  fOutputListQA->Add(fh1QAV0Eta[i]);
1201  fOutputListQA->Add(fh2QAV0EtaRows[i]);
1202  fOutputListQA->Add(fh2QAV0PtRows[i]);
1203  fOutputListQA->Add(fh2QAV0PhiRows[i]);
1204  fOutputListQA->Add(fh2QAV0NClRows[i]);
1205  fOutputListQA->Add(fh2QAV0EtaNCl[i]);
1206 
1210  fOutputListQA->Add(fh1QAV0RapK0s[i]);
1212 
1218 
1224 
1228  fOutputListQA->Add(fh1QAV0Pt[i]);
1229  fOutputListQA->Add(fh1QAV0Charge[i]);
1230  fOutputListQA->Add(fh1QAV0DCAVtx[i]);
1231  fOutputListQA->Add(fh1QAV0DCAV0[i]);
1232  fOutputListQA->Add(fh1QAV0Cos[i]);
1233  fOutputListQA->Add(fh1QAV0R[i]);
1234  fOutputListQA->Add(fh1QACTau2D[i]);
1235  fOutputListQA->Add(fh1QACTau3D[i]);
1236 
1237  fOutputListQA->Add(fh2ArmPod[i]);
1238  fOutputListQA->Add(fh2ArmPodK0s[i]);
1239  fOutputListQA->Add(fh2ArmPodLambda[i]);
1241 
1242  /*
1243  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);
1244  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);
1245  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);
1246  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);
1247  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);
1248  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);
1249  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);
1250  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);
1251  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);
1252  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);
1253  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);
1254  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);
1255  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);
1256  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);
1257  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);
1258  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);
1259  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);
1260  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);
1261  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);
1262  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);
1263  fh2Tau3DVs2D[i] = new TH2D(Form("fh2Tau3DVs2D_%d", i), "Decay 3D vs 2D;pt;3D/2D", 100, 0, 10, 200, 0.5, 1.5);
1264 
1265  fOutputListCuts->Add(fh2CutTPCRowsK0s[i]);
1266  fOutputListCuts->Add(fh2CutTPCRowsLambda[i]);
1267  fOutputListCuts->Add(fh2CutPtPosK0s[i]);
1268  fOutputListCuts->Add(fh2CutPtNegK0s[i]);
1269  fOutputListCuts->Add(fh2CutPtPosLambda[i]);
1270  fOutputListCuts->Add(fh2CutPtNegLambda[i]);
1271  fOutputListCuts->Add(fh2CutDCAVtx[i]);
1272  fOutputListCuts->Add(fh2CutDCAV0[i]);
1273  fOutputListCuts->Add(fh2CutCos[i]);
1274  fOutputListCuts->Add(fh2CutR[i]);
1275  fOutputListCuts->Add(fh2CutEtaK0s[i]);
1276  fOutputListCuts->Add(fh2CutEtaLambda[i]);
1277  fOutputListCuts->Add(fh2CutRapK0s[i]);
1278  fOutputListCuts->Add(fh2CutRapLambda[i]);
1279  fOutputListCuts->Add(fh2CutCTauK0s[i]);
1280  fOutputListCuts->Add(fh2CutCTauLambda[i]);
1281  fOutputListCuts->Add(fh2CutPIDPosK0s[i]);
1282  fOutputListCuts->Add(fh2CutPIDNegK0s[i]);
1283  fOutputListCuts->Add(fh2CutPIDPosLambda[i]);
1284  fOutputListCuts->Add(fh2CutPIDNegLambda[i]);
1285  fOutputListCuts->Add(fh2Tau3DVs2D[i]);
1286  */
1287  }
1288 
1289  for(Int_t i = 0; i < fgkiNCategV0; i++)
1290  {
1291  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);
1292  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);
1293  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);
1297  }
1298 
1299  for(Int_t i = 0; i < fOutputListStd->GetEntries(); ++i)
1300  {
1301  TH1* h1 = dynamic_cast<TH1*>(fOutputListStd->At(i));
1302  if(h1)
1303  {
1304  h1->Sumw2();
1305  continue;
1306  }
1307  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListStd->At(i));
1308  if(hn) hn->Sumw2();
1309  }
1310  for(Int_t i = 0; i < fOutputListQA->GetEntries(); ++i)
1311  {
1312  TH1* h1 = dynamic_cast<TH1*>(fOutputListQA->At(i));
1313  if(h1)
1314  {
1315  h1->Sumw2();
1316  continue;
1317  }
1318  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListQA->At(i));
1319  if(hn) hn->Sumw2();
1320  }
1321  for(Int_t i = 0; i < fOutputListCuts->GetEntries(); ++i)
1322  {
1323  TH1* h1 = dynamic_cast<TH1*>(fOutputListCuts->At(i));
1324  if(h1)
1325  {
1326  h1->Sumw2();
1327  continue;
1328  }
1329  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListCuts->At(i));
1330  if(hn) hn->Sumw2();
1331  }
1332  for(Int_t i = 0; i < fOutputListMC->GetEntries(); ++i)
1333  {
1334  TH1* h1 = dynamic_cast<TH1*>(fOutputListMC->At(i));
1335  if(h1)
1336  {
1337  h1->Sumw2();
1338  continue;
1339  }
1340  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListMC->At(i));
1341  if(hn) hn->Sumw2();
1342  }
1343 
1344  PostData(1, fOutputListStd);
1345  PostData(2, fOutputListQA);
1346  PostData(3, fOutputListCuts);
1347  PostData(4, fOutputListMC);
1348 }
1349 
1351 {
1353 
1354  if(fbJetSelection)
1355  {
1356  if(fJetsCont && fJetsCont->GetArray() == 0)
1357  fJetsCont = 0;
1358  if(fJetsBgCont && fJetsBgCont->GetArray() == 0)
1359  fJetsBgCont = 0;
1360  if(fbCompareTriggers)
1361  {
1362  if(fTracksCont && fTracksCont->GetArray() == 0)
1363  fTracksCont = 0;
1364  }
1365  }
1366 
1367  // Jet selection
1368  Double_t dCutEtaJetMax = fdCutEtaV0Max - fdDistanceV0JetMax; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1369  if(fbJetSelection && fJetsCont)
1370  {
1373  if(fdCutEtaV0Max > 0. && fdDistanceV0JetMax > 0.) fJetsCont->SetJetEtaLimits(-dCutEtaJetMax, dCutEtaJetMax);
1374  // Trigger selection
1376  {
1379  fTracksCont->SetParticleEtaLimits(-0.8, 0.8);
1380  }
1381  }
1382 
1383  printf("=======================================================\n");
1384  printf("%s: Configuration summary:\n", ClassName());
1385  printf("task name: %s\n", GetName());
1386  printf("-------------------------------------------------------\n");
1387  printf("collision system: %s\n", fbIsPbPb ? "Pb+Pb" : "p+p");
1388  printf("data type: %s\n", fbMCAnalysis ? "MC" : "real");
1389  if(fbMCAnalysis)
1390  printf("MC generator: %s\n", fsGeneratorName.Length() ? fsGeneratorName.Data() : "any");
1391  if(fbIsPbPb)
1392  printf("centrality range: %g-%g %%\n", fdCutCentLow, fdCutCentHigh);
1393  if(fdCutVertexZ > 0.) printf("max |z| of the prim vtx [cm]: %g\n", fdCutVertexZ);
1394  if(fdCutVertexR2 > 0.) printf("max r^2 of the prim vtx [cm^2]: %g\n", fdCutVertexR2);
1395  if(fdCutDeltaZMax > 0.) printf("max |Delta z| between nominal prim vtx and SPD vtx [cm]: %g\n", fdCutDeltaZMax);
1396  printf("-------------------------------------------------------\n");
1397  if(fbTPCRefit) printf("TPC refit for daughter tracks\n");
1398  if(fbRejectKinks) printf("reject kink-like production vertices of daughter tracks\n");
1399  if(fbFindableClusters) printf("require positive number of findable clusters\n");
1400  if(fdCutNCrossedRowsTPCMin > 0.) printf("min number of crossed TPC rows: %g\n", fdCutNCrossedRowsTPCMin);
1401  if(fdCutCrossedRowsOverFindMin > 0.) printf("min ratio crossed rows / findable clusters: %g\n", fdCutCrossedRowsOverFindMin);
1402  if(fdCutCrossedRowsOverFindMax > 0.) printf("max ratio crossed rows / findable clusters: %g\n", fdCutCrossedRowsOverFindMax);
1403  if(fdCutPtDaughterMin > 0.) printf("min pt of daughter tracks [GeV/c]: %g\n", fdCutPtDaughterMin);
1404  if(fdCutDCAToPrimVtxMin > 0.) printf("min DCA of daughters to the prim vtx [cm]: %g\n", fdCutDCAToPrimVtxMin);
1405  if(fdCutDCADaughtersMax > 0.) printf("max DCA between daughters [sigma of TPC tracking]: %g\n", fdCutDCADaughtersMax);
1406  if(fdCutEtaDaughterMax > 0.) printf("max |eta| of daughter tracks: %g\n", fdCutEtaDaughterMax);
1407  if(fdCutNSigmadEdxMax > 0. && (!fbIsPbPb || (fbIsPbPb && fdPtProtonPIDMax > 0.))) printf("max |Delta(dE/dx)| in the TPC [sigma dE/dx]: %g\n", fdCutNSigmadEdxMax);
1408  if(fdCutNSigmadEdxMax > 0. && fbIsPbPb && fdPtProtonPIDMax > 0.) printf("max pt of proton for applying PID cut [GeV/c]: %g\n", fdPtProtonPIDMax);
1409  printf("V0 reconstruction method: %s\n", fbOnFly ? "on-the-fly" : "offline");
1410  if(fdCutCPAKMin > 0.) printf("min CPA, K0S: %g\n", fdCutCPAKMin);
1411  if(fdCutCPALMin > 0.) printf("min CPA, (A)Lambda: %g\n", fdCutCPALMin);
1412  if(fdCutRadiusDecayMin > 0. && fdCutRadiusDecayMax > 0.) printf("R of the decay vertex [cm]: %g-%g\n", fdCutRadiusDecayMin, fdCutRadiusDecayMax);
1413  if(fdCutEtaV0Max > 0.) printf("max |eta| of V0: %g\n", fdCutEtaV0Max);
1414  if(fdCutRapV0Max > 0.) printf("max |y| of V0: %g\n", fdCutRapV0Max);
1415  if(fdCutNTauKMax > 0.) printf("max proper lifetime, K0S [tau]: %g\n", fdCutNTauKMax);
1416  if(fdCutNTauLMax > 0.) printf("max proper lifetime, (A)Lambda [tau]: %g\n", fdCutNTauLMax);
1417  if(fbCutArmPod) printf("Armenteros-Podolanski cut for K0S\n");
1418  if(fbCutCross) printf("cross-contamination cut\n");
1419  printf("-------------------------------------------------------\n");
1420  printf("analysis of V0s in jets: %s\n", fbJetSelection ? "yes" : "no");
1421  if(fbJetSelection)
1422  {
1423  printf("rho subtraction: ");
1424  if(fiBgSubtraction == 0) printf("none\n");
1425  else if(fiBgSubtraction == 1) printf("scalar\n");
1426  else if(fiBgSubtraction == 2) printf("vector\n");
1427  if(fdCutPtJetMin > 0.) printf("min jet pt [GeV/c]: %g\n", fdCutPtJetMin);
1428  if(fdCutPtTrackJetMin > 0.) printf("min pt of leading jet-track [GeV/c]: %g\n", fdCutPtTrackJetMin);
1429  if(fdCutAreaPercJetMin > 0.) printf("min area of jet [pi*R^2]: %g\n", fdCutAreaPercJetMin);
1430  if(fdCutEtaV0Max > 0. && fdDistanceV0JetMax > 0.) printf("max |eta| of jet: %g\n", dCutEtaJetMax);
1431  printf("max distance between V0 and jet axis: %g\n", fdDistanceV0JetMax);
1432  printf("angular correlations of V0s with jets: %s\n", fbCorrelations ? "yes" : "no");
1433  if(fbCorrelations) printf("max |delta-eta_V0-jet|: %g\n", fdDeltaEtaMax);
1434  printf("pt correlations of jets with trigger tracks: %s\n", fbCompareTriggers ? "yes" : "no");
1435  printf("-------------------------------------------------------\n");
1436  if(fJetsCont)
1437  {
1438  printf("Signal jet container parameters\n");
1439  printf("Jet R = %g\n", fJetsCont->GetJetRadius());
1440  fJetsCont->PrintCuts();
1441  }
1442  else
1443  printf("No signal jet container!\n");
1444  if(fJetsBgCont)
1445  {
1446  printf("Background jet container parameters\n");
1447  printf("Jet R = %g\n", fJetsBgCont->GetJetRadius());
1449  }
1450  else
1451  printf("No background jet container\n");
1452  }
1453  printf("=======================================================\n");
1454 }
1455 
1457 {
1458 // Run analysis code here, if needed. It will be executed before FillHistograms().
1459 
1460  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
1461 }
1462 
1464 {
1465  // Main loop, called for each event
1466  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Start");
1467 
1468  fh1EventCounterCut->Fill(0); // all available selected events (collision candidates)
1469 
1470  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); // input AOD
1471  if(!fAODIn)
1472  {
1473  AliError("No input AOD found!");
1474  return kFALSE;
1475  }
1476  if(fDebug > 1) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Loading AOD OK");
1477 
1478  TClonesArray* arrayMC = 0; // array particles in the MC event
1479  AliAODMCHeader* headerMC = 0; // MC header
1480  Int_t iNTracksMC = 0; // number of MC tracks
1481  Double_t dPrimVtxMCX = 0., dPrimVtxMCY = 0., dPrimVtxMCZ = 0.; // position of the MC primary vertex
1482 
1483  // Simulation info
1484  if(fbMCAnalysis)
1485  {
1486  fEventMC = MCEvent();
1487  arrayMC = (TClonesArray*)fAODIn->FindListObject(AliAODMCParticle::StdBranchName());
1488  if(!arrayMC || !fEventMC)
1489  {
1490  AliError("No MC array/event found!");
1491  return kFALSE;
1492  }
1493  if(fDebug > 1) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "MC array found");
1494  iNTracksMC = arrayMC->GetEntriesFast();
1495  if(fDebug > 2) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("There are %d MC tracks in this event", iNTracksMC));
1496  headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
1497  if(!headerMC)
1498  {
1499  AliError("No MC header found!");
1500  return kFALSE;
1501  }
1502  // get position of the MC primary vertex
1503  dPrimVtxMCX = headerMC->GetVtxX();
1504  dPrimVtxMCY = headerMC->GetVtxY();
1505  dPrimVtxMCZ = headerMC->GetVtxZ();
1506  }
1507 
1508  // PID Response Task object
1509  AliPIDResponse* fPIDResponse = 0;
1510  if(fdCutNSigmadEdxMax > 0. && (!fbIsPbPb || (fbIsPbPb && fdPtProtonPIDMax > 0.)))
1511  {
1512  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1513  AliInputEventHandler* inputHandler = (AliInputEventHandler*)mgr->GetInputEventHandler();
1514  fPIDResponse = inputHandler->GetPIDResponse();
1515  if(!fPIDResponse)
1516  {
1517  AliError("No PID response object found!");
1518  return kFALSE;
1519  }
1520  }
1521 
1522  // AOD files are OK
1523  fh1EventCounterCut->Fill(1);
1524 
1525  // Event selection
1527  {
1528  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Event rejected");
1529  return kFALSE;
1530  }
1531  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Event accepted: cent. %g", fdCentrality));
1532  if(!fbIsPbPb)
1533  fdCentrality = 0.; // select the first bin for p+p data
1534  Int_t iCentIndex = GetCentralityBinIndex(fdCentrality); // get index of centrality bin
1535  if(iCentIndex < 0)
1536  {
1537  AliError(Form("Event is out of histogram range: cent: %g!", fdCentrality));
1538  return kFALSE;
1539  }
1540  fh1EventCounterCut->Fill(8); // selected events (centrality OK)
1541  fh1EventCounterCutCent[iCentIndex]->Fill(8);
1542 
1543  UInt_t iNTracks = fAODIn->GetNumberOfTracks(); // get number of tracks in event
1544  if(fDebug > 2) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("There are %d tracks in this event", iNTracks));
1545 
1546 // Double_t dMagField = fAODIn->GetMagneticField();
1547 // if(fDebug > 2) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Magnetic field: %g", dMagField));
1548 
1549  Int_t iNV0s = fAODIn->GetNumberOfV0s(); // get the number of V0 candidates
1550  if(!iNV0s)
1551  {
1552  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "No V0s found in event");
1553  }
1554 
1555  //===== Event is OK for the analysis =====
1556  fh1EventCent->Fill(iCentIndex);
1557  fh1EventCent2->Fill(fdCentrality);
1558  fh2EventCentTracks->Fill(fdCentrality, iNTracks);
1559 
1560  AliEventPool* pool = 0;
1561  Bool_t bPoolReady = kFALSE; // status of pool
1562  Double_t dZVtxME = -100; // z_vtx for events used in mixed events
1563  if(fbCorrelations)
1564  {
1565  AliAODVertex* vertex = fAODIn->GetPrimaryVertex();
1566  dZVtxME = vertex->GetZ();
1567  pool = fPoolMgr->GetEventPool(fdCentrality, dZVtxME);
1568  if(!pool)
1569  {
1570  AliError(Form("No pool found for centrality = %g, z_vtx = %g!", fdCentrality, dZVtxME));
1571  return kFALSE;
1572  }
1573  bPoolReady = pool->IsReady();
1574  if(fDebug > 2)
1575  {
1576  // pool->SetDebug(1);
1577  // pool->PrintInfo();
1578  printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Pool %d-%d: events in pool = %d, jets in pool = %d", pool->MultBinIndex(), pool->ZvtxBinIndex(), pool->GetCurrentNEvents(), pool->NTracksInPool()));
1579  }
1580  }
1581 
1582  if(iNV0s)
1583  {
1584  fh1EventCounterCut->Fill(9); // events with V0s
1585  fh1EventCounterCutCent[iCentIndex]->Fill(9);
1586  }
1587 
1588  AliAODv0* v0 = 0; // pointer to V0 candidates
1589  TVector3 vecV0Momentum; // 3D vector of V0 momentum
1590  Double_t dMassV0K0s = 0; // invariant mass of the K0s candidate
1591  Double_t dMassV0Lambda = 0; // invariant mass of the Lambda candidate
1592  Double_t dMassV0ALambda = 0; // invariant mass of the Lambda candidate
1593  Int_t iNV0CandTot = 0; // counter of all V0 candidates at the beginning
1594  Int_t iNV0CandK0s = 0; // counter of K0s candidates at the end
1595  Int_t iNV0CandLambda = 0; // counter of Lambda candidates at the end
1596  Int_t iNV0CandALambda = 0; // counter of Lambda candidates at the end
1597 
1598  Bool_t bPrintCuts = 0; // print out which cuts are applied
1599  Bool_t bPrintJetSelection = 0; // print out which jets are selected
1600 
1601  // Other cuts
1602  Double_t dNSigmaMassMax = 3.; // [sigma m] max difference between candidate mass and real particle mass (used only for mass peak method of signal extraction)
1603  Double_t dDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
1604 
1605  // Mean lifetime
1606  Double_t dCTauK0s = 2.6844; // [cm] c*tau of K0S
1607  Double_t dCTauLambda = 7.89; // [cm] c*tau of Lambda
1608 
1609  // particle masses from PDG
1610  Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1611  Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1612 
1613  // PDG codes of used particles
1614  Int_t iPdgCodePion = 211;
1615  Int_t iPdgCodeProton = 2212;
1616  Int_t iPdgCodeK0s = 310;
1617  Int_t iPdgCodeLambda = 3122;
1618 
1619  Double_t dCutEtaJetMax = fdCutEtaV0Max - fdDistanceV0JetMax; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1620  Double_t dRadiusExcludeCone = 2 * fdDistanceV0JetMax; // radius of cones around jets excluded for V0 outside jets
1621  Bool_t bLeadingJetOnly = 0; // consider only leading jets
1622 
1623  Int_t iNJet = 0; // number of reconstructed jets in the input
1624  TClonesArray* arrayJetSel = new TClonesArray("AliAODJet", 0); // object where the selected jets are copied as AliAODJet
1625  Int_t iNJetSel = 0; // number of selected reconstructed jets
1626  TClonesArray* arrayJetPerp = new TClonesArray("AliAODJet", 0); // object where the perp. cones are stored
1627  Int_t iNJetPerp = 0; // number of perpendicular cones
1628  TObjArray* arrayMixedEventAdd = new TObjArray; // array of jets from this event to be added to the pool
1629  arrayMixedEventAdd->SetOwner(kTRUE);
1630 
1631  AliAODJet* jet = 0; // pointer to a jet
1632  AliAODJet* jetPerp = 0; // pointer to a perp. cone
1633  AliAODJet* jetRnd = 0; // pointer to a rand. cone
1634  AliEmcalJet* jetMed = 0; // pointer to a median cluster
1635  TVector3 vecJetMomentum; // 3D vector of jet momentum
1636  TVector3 vecJetMomentumPair; // 3D vector of another jet momentum for calculating distance between jets
1637  Bool_t bJetEventGood = kTRUE; // indicator of good jet events
1638  Double_t dRho = 0; // average bg pt density
1639  TLorentzVector vecJetSel; // 4-momentum of selected jet
1640  TLorentzVector vecPerpPlus; // 4-momentum of perpendicular cone plus
1641  TLorentzVector vecPerpMinus; // 4-momentum of perpendicular cone minus
1642  Double_t dPtJet = 0; // pt of jet associated with V0
1643  Double_t dPtJetTrackLeading = 0; // pt of leading track of jet associated with V0
1644 
1645  if(fbJetSelection) // analysis of V0s in jets is switched on
1646  {
1647  if(!fJetsCont)
1648  {
1649  AliError("No signal jet container!");
1650  bJetEventGood = kFALSE;
1651  }
1652  if(bJetEventGood)
1653  iNJet = fJetsCont->GetNJets();
1654  if(bJetEventGood && !iNJet) // check whether there are some jets
1655  {
1656  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "No jets in array");
1657  bJetEventGood = kFALSE;
1658  }
1659  if(fbIsPbPb && bJetEventGood && !fJetsBgCont)
1660  {
1661  AliWarning("No bg jet container!");
1662 // bJetEventGood = kFALSE;
1663  }
1664  }
1665  else // no in-jet analysis
1666  bJetEventGood = kFALSE;
1667 
1668  // select good jets and copy them to another array
1669  if(bJetEventGood)
1670  {
1671  if(fiBgSubtraction)
1672  {
1673  dRho = fJetsCont->GetRhoVal();
1674  if(fDebug > 4) printf("%s::%s: %s\n", ClassName(), __func__, Form("Loaded rho value: %g", dRho));
1675  }
1676  if(bLeadingJetOnly)
1677  iNJet = 1; // only leading jets
1678  if(fDebug > 2) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Jet selection for %d jets", iNJet));
1679  for(Int_t iJet = 0; iJet < iNJet; iJet++)
1680  {
1681  AliEmcalJet* jetSel = (AliEmcalJet*)(fJetsCont->GetAcceptJet(iJet)); // load a jet in the list
1682  if(bLeadingJetOnly)
1683  jetSel = fJetsCont->GetLeadingJet();
1684  if(!jetSel)
1685  {
1686  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Jet %d not accepted in container", iJet));
1687  continue;
1688  }
1689 
1690  Double_t dPtJetCorr = jetSel->Pt(); // raw pt
1691  Double_t dEtaJetCorr = jetSel->Eta(); // raw eta
1692  Double_t dPhiJetCorr = jetSel->Phi(); // raw phi
1693  Double_t dPtTrackJet = fJetsCont->GetLeadingHadronPt(jetSel); // pt of leading track
1694 
1695  if(fiBgSubtraction == 1) // scalar subtraction: correct pt only
1696  {
1697  dPtJetCorr = jetSel->PtSub(dRho);
1698  }
1699  else if(fiBgSubtraction == 2) // vector subtraction: correct momentum vector
1700  {
1701  vecJetSel = jetSel->SubtractRhoVect(dRho);
1702  dPtJetCorr = vecJetSel.Pt();
1703  dEtaJetCorr = vecJetSel.Eta();
1704  dPhiJetCorr = TVector2::Phi_0_2pi(vecJetSel.Phi());
1705  }
1706 
1707  if(bPrintJetSelection)
1708  if(fDebug > 4) printf("jet: i = %d, pT = %g, eta = %g, phi = %g, pt lead tr = %g ", iJet, dPtJetCorr, dEtaJetCorr, dPhiJetCorr, dPtTrackJet);
1709  if(fdCutPtJetMin > 0. && dPtJetCorr < fdCutPtJetMin) // selection of high-pt jets, needs to be applied on the pt after bg subtraction
1710  {
1711  if(bPrintJetSelection)
1712  if(fDebug > 4) printf("rejected (pt)\n");
1713  continue;
1714  }
1715  if(bPrintJetSelection)
1716  if(fDebug > 4) printf("accepted\n");
1717  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Jet %d with pt %g passed selection", iJet, dPtJetCorr));
1718 
1719  vecJetSel.SetPtEtaPhiM(dPtJetCorr, dEtaJetCorr, dPhiJetCorr, 0.);
1720  vecPerpPlus = vecJetSel;
1721  vecPerpMinus = vecJetSel;
1722  vecPerpPlus.RotateZ(TMath::Pi() / 2.); // rotate vector by +90 deg around z
1723  vecPerpMinus.RotateZ(-TMath::Pi() / 2.); // rotate vector by -90 deg around z
1724  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Adding perp. cones number %d, %d, pt %g", iNJetPerp, iNJetPerp + 1, vecPerpPlus.Pt()));
1725  new((*arrayJetPerp)[iNJetPerp++]) AliAODJet(vecPerpPlus); // write perp. cone to the array
1726  new((*arrayJetPerp)[iNJetPerp++]) AliAODJet(vecPerpMinus); // write perp. cone to the array
1727  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Adding jet number %d", iNJetSel));
1728  new((*arrayJetSel)[iNJetSel++]) AliAODJet(vecJetSel); // copy selected jet to the array
1729  ((AliAODJet*)arrayJetSel->At(iNJetSel - 1))->SetPtLeading(dPtTrackJet);
1730  fh2PtJetPtTrackLeading[iCentIndex]->Fill(dPtJetCorr, dPtTrackJet); // pt_jet vs pt of leading jet track
1731  if(fbCorrelations)
1732  arrayMixedEventAdd->Add(new TLorentzVector(vecJetSel)); // copy selected jet to the list of new jets for event mixing
1733  }
1734  if(fDebug > 3) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Added jets: %d", iNJetSel));
1735  iNJetSel = arrayJetSel->GetEntriesFast();
1736  if(fDebug > 2) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Selected jets in array: %d", iNJetSel));
1737  // fill jet spectra
1738  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1739  {
1740  jet = (AliAODJet*)arrayJetSel->At(iJet); // load a jet in the list
1741  fh1PtJet[iCentIndex]->Fill(jet->Pt()); // pt spectrum of selected jets
1742  fh1EtaJet[iCentIndex]->Fill(jet->Eta()); // eta spectrum of selected jets
1743  fh2EtaPtJet[iCentIndex]->Fill(jet->Eta(), jet->Pt()); // eta-pT spectrum of selected jets
1744  fh1PhiJet[iCentIndex]->Fill(jet->Phi()); // phi spectrum of selected jets
1745  Double_t dAreaExcluded = TMath::Pi() * dRadiusExcludeCone * dRadiusExcludeCone; // area of the cone
1746  dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, fdCutEtaV0Max - jet->Eta()); // positive eta overhang
1747  dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, fdCutEtaV0Max + jet->Eta()); // negative eta overhang
1748  fh1AreaExcluded->Fill(iCentIndex, dAreaExcluded);
1749  // calculate distances between all pairs of selected jets in the event
1750  vecJetMomentum.SetXYZ(jet->Px(), jet->Py(), jet->Pz()); // set the vector of jet momentum
1751  for(Int_t iJetPair = iJet + 1; iJetPair < iNJetSel; iJetPair++)
1752  {
1753  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Calculating distance for jet pair: %d-%d/%d\n", iJet, iJetPair, iNJetSel));
1754  jet = (AliAODJet*)arrayJetSel->At(iJetPair); // load another jet in the list
1755  vecJetMomentumPair.SetXYZ(jet->Px(), jet->Py(), jet->Pz()); // set the vector of the second jet momentum
1756  fh1DistanceJets[iCentIndex]->Fill(vecJetMomentum.DeltaR(vecJetMomentumPair));
1757  }
1758  }
1759  jet = 0; // just to be sure
1760  }
1761  fh1NJetPerEvent[iCentIndex]->Fill(iNJetSel);
1762 
1763  if(bJetEventGood) // there should be some reconstructed jets
1764  {
1765  fh1EventCounterCut->Fill(10); // events with jet(s)
1766  fh1EventCounterCutCent[iCentIndex]->Fill(10); // events with jet(s)
1767  if(iNJetSel)
1768  {
1769  fh1EventCounterCut->Fill(11); // events with selected jets
1770  fh1EventCounterCutCent[iCentIndex]->Fill(11);
1771  }
1772  }
1773  if(iNJetSel)
1775  else
1777 
1778  if(iNJetSel)
1779  {
1780  jetRnd = GetRandomCone(arrayJetSel, dCutEtaJetMax, 2 * fdDistanceV0JetMax);
1781  if(jetRnd)
1782  {
1783  fh1NRndConeCent->Fill(iCentIndex);
1784  fh2EtaPhiRndCone[iCentIndex]->Fill(jetRnd->Eta(), jetRnd->Phi());
1785  }
1786  if(fJetsBgCont)
1787  {
1788  jetMed = GetMedianCluster(fJetsBgCont, dCutEtaJetMax);
1789  if(jetMed)
1790  {
1791  fh1NMedConeCent->Fill(iCentIndex);
1792  fh2EtaPhiMedCone[iCentIndex]->Fill(jetMed->Eta(), jetMed->Phi());
1793  }
1794  }
1795  }
1796 
1797  if(fbJetSelection && fbCompareTriggers) // Correlations of pt_jet with pt_trigger-track
1798  {
1799  if(!fTracksCont)
1800  AliError("No track container!");
1801  else
1802  {
1803  AliAODJet* jetTrig = 0;
1804  fTracksCont->ResetCurrentID();
1805  AliVTrack* track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
1806  if(track) // there are some accepted trigger tracks
1807  {
1808  while(track) // loop over selected tracks
1809  {
1810  fh1PtTrigger[iCentIndex]->Fill(track->Pt());
1811  if(iNJetSel) // there are some accepted jets
1812  {
1813  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1814  {
1815  jetTrig = (AliAODJet*)arrayJetSel->At(iJet);
1816  fh2PtJetPtTrigger[iCentIndex]->Fill(jetTrig->Pt(), track->Pt());
1817  }
1818  }
1819  else // there are no accepted jets
1820  {
1821  fh2PtJetPtTrigger[iCentIndex]->Fill(0., track->Pt());
1822  }
1823  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()); // load next accepted trigger track
1824  }
1825  }
1826  else // there are no accepted trigger tracks
1827  {
1828  if(iNJetSel) // there are some accepted jets
1829  {
1830  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1831  {
1832  jetTrig = (AliAODJet*)arrayJetSel->At(iJet);
1833  fh2PtJetPtTrigger[iCentIndex]->Fill(jetTrig->Pt(), 0.);
1834  }
1835  }
1836  else // there are no accepted jets
1837  fh2PtJetPtTrigger[iCentIndex]->Fill(0., 0.);
1838  }
1839  }
1840  }
1841 
1842  // Loading primary vertex info
1843  AliAODVertex* primVtx = fAODIn->GetPrimaryVertex(); // get the primary vertex
1844  Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
1845  primVtx->GetXYZ(dPrimVtxPos);
1846  fh1VtxZ[iCentIndex]->Fill(dPrimVtxPos[2]);
1847  fh2VtxXY[iCentIndex]->Fill(dPrimVtxPos[0], dPrimVtxPos[1]);
1848 
1849  //===== Start of loop over V0 candidates =====
1850  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Start of V0 loop");
1851  for(Int_t iV0 = 0; iV0 < iNV0s; iV0++)
1852  {
1853  v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1854  if(!v0)
1855  continue;
1856 
1857  iNV0CandTot++;
1858 
1859  // Initialization of status indicators
1860  Bool_t bIsCandidateK0s = kTRUE; // candidate for K0s
1861  Bool_t bIsCandidateLambda = kTRUE; // candidate for Lambda
1862  Bool_t bIsCandidateALambda = kTRUE; // candidate for anti-Lambda
1863  Bool_t bIsInPeakK0s = kFALSE; // candidate within the K0s mass peak
1864  Bool_t bIsInPeakLambda = kFALSE; // candidate within the Lambda mass peak
1865  Bool_t bIsInPeakALambda = kFALSE; // candidate within the anti-Lambda mass peak
1866  Bool_t bIsInConeJet = kFALSE; // candidate within the jet cones
1867  Bool_t bIsInConePerp = kFALSE; // candidate within a perpendicular cone
1868  Bool_t bIsInConeRnd = kFALSE; // candidate within the random cone
1869  Bool_t bIsInConeMed = kFALSE; // candidate within the median-cluster cone
1870  Bool_t bIsOutsideCones = kFALSE; // candidate outside excluded cones
1871 
1872  // Invariant mass calculation
1873  dMassV0K0s = v0->MassK0Short();
1874  dMassV0Lambda = v0->MassLambda();
1875  dMassV0ALambda = v0->MassAntiLambda();
1876 
1877  Int_t iCutIndex = 0; // indicator of current selection step
1878  // 0
1879  // All V0 candidates
1880  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1881  iCutIndex++;
1882 
1883  Double_t dPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
1884  vecV0Momentum = TVector3(v0->Px(), v0->Py(), v0->Pz()); // set the vector of V0 momentum
1885 
1886  // Sigma of the mass peak window
1887  Double_t dMassPeakWindowK0s = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 0);
1888  Double_t dMassPeakWindowLambda = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 1);
1889  if(!fbIsPbPb) // p-p
1890  {
1891  dMassPeakWindowK0s = 0.010; // LF p-p
1892  dMassPeakWindowLambda = 0.005; // LF p-p
1893  }
1894 
1895  // Invariant mass peak selection
1896  if(TMath::Abs(dMassV0K0s - dMassPDGK0s) < dMassPeakWindowK0s)
1897  bIsInPeakK0s = kTRUE;
1898  if(TMath::Abs(dMassV0Lambda - dMassPDGLambda) < dMassPeakWindowLambda)
1899  bIsInPeakLambda = kTRUE;
1900  if(TMath::Abs(dMassV0ALambda - dMassPDGLambda) < dMassPeakWindowLambda)
1901  bIsInPeakALambda = kTRUE;
1902 
1903  // QA histograms before cuts
1904  FillQAHistogramV0(primVtx, v0, 0, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, bIsInPeakK0s, bIsInPeakLambda, bIsInPeakALambda);
1905 
1906  // Skip candidates outside the histogram range
1907  if((dMassV0K0s < fgkdMassK0sMin) || (dMassV0K0s >= fgkdMassK0sMax))
1908  bIsCandidateK0s = kFALSE;
1909  if((dMassV0Lambda < fgkdMassLambdaMin) || (dMassV0Lambda >= fgkdMassLambdaMax))
1910  bIsCandidateLambda = kFALSE;
1911  if((dMassV0ALambda < fgkdMassLambdaMin) || (dMassV0ALambda >= fgkdMassLambdaMax))
1912  bIsCandidateALambda = kFALSE;
1913  if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1914  continue;
1915 
1916  // Retrieving all relevant properties of the V0 candidate
1917  Bool_t bOnFlyStatus = v0->GetOnFlyStatus(); // online (on fly) reconstructed vs offline reconstructed
1918  const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1919  const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
1920  Double_t dPtDaughterPos = trackPos->Pt(); // transverse momentum of a daughter track calculated as if primary, != v0->PtProng(0)
1921  Double_t dPtDaughterNeg = trackNeg->Pt(); // != v0->PtProng(1)
1922  Double_t dNRowsPos = trackPos->GetTPCClusterInfo(2, 1); // crossed TPC pad rows of a daughter track
1923  Double_t dNRowsNeg = trackNeg->GetTPCClusterInfo(2, 1);
1924  Double_t dFindablePos = Double_t(trackPos->GetTPCNclsF()); // Findable clusters
1925  Double_t dFindableNeg = Double_t(trackNeg->GetTPCNclsF());
1926  Double_t dDCAToPrimVtxPos = TMath::Abs(v0->DcaPosToPrimVertex()); // dca of a daughter to the primary vertex
1927  Double_t dDCAToPrimVtxNeg = TMath::Abs(v0->DcaNegToPrimVertex());
1928  Double_t dDCADaughters = v0->DcaV0Daughters(); // dca between daughters
1929  Double_t dCPA = v0->CosPointingAngle(primVtx); // cosine of the pointing angle
1930  Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
1931 // Double_t dSecVtxPos[3] = {v0->DecayVertexV0X(),v0->DecayVertexV0Y(),v0->DecayVertexV0Z()}; // V0 vertex position
1932  v0->GetSecondaryVtx(dSecVtxPos);
1933  Double_t dRadiusDecay = TMath::Sqrt(dSecVtxPos[0] * dSecVtxPos[0] + dSecVtxPos[1] * dSecVtxPos[1]); // distance of the V0 vertex from the z-axis
1934  Double_t dEtaDaughterPos = trackPos->Eta(); // pseudorapidity of a daughter track calculated as if primary, != v0->EtaProng(0)
1935  Double_t dEtaDaughterNeg = trackNeg->Eta(); // != v0->EtaProng(1);
1936  Double_t dRapK0s = v0->RapK0Short(); // rapidity calculated for K0s assumption
1937  Double_t dRapLambda = v0->RapLambda(); // rapidity calculated for Lambda assumption
1938  Double_t dEtaV0 = v0->Eta(); // V0 pseudorapidity
1939  Double_t dPhiV0 = v0->Phi(); // V0 azimuth
1940  Double_t dDecayPath[3];
1941  for(Int_t iPos = 0; iPos < 3; iPos++)
1942  dDecayPath[iPos] = dSecVtxPos[iPos] - dPrimVtxPos[iPos]; // vector of the V0 path
1943  Double_t dDecLen = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]); // path length L
1944  Double_t dDecLen2D = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]); // transverse path length R
1945  Double_t dLOverP = dDecLen / v0->P(); // L/p
1946  Double_t dROverPt = dDecLen2D / dPtV0; // R/pT
1947  Double_t dMLOverPK0s = dMassPDGK0s * dLOverP; // m*L/p = c*(proper lifetime)
1948 // Double_t dMLOverPLambda = dMassPDGLambda*dLOverP; // m*L/p
1949  Double_t dMROverPtK0s = dMassPDGK0s * dROverPt; // m*R/pT
1950  Double_t dMROverPtLambda = dMassPDGLambda * dROverPt; // m*R/pT
1951  Double_t dNSigmaPosPion = (fPIDResponse ? TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kPion)) : 0.); // difference between measured and expected signal of the dE/dx in the TPC
1952  Double_t dNSigmaPosProton = (fPIDResponse ? TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kProton)) : 0.);
1953  Double_t dNSigmaNegPion = (fPIDResponse ? TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kPion)) : 0.);
1954  Double_t dNSigmaNegProton = (fPIDResponse ? TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kProton)) : 0.);
1955  Double_t dAlpha = v0->AlphaV0(); // Armenteros-Podolanski alpha
1956  Double_t dPtArm = v0->PtArmV0(); // Armenteros-Podolanski pT
1957  AliAODVertex* prodVtxDaughterPos = (AliAODVertex*)(trackPos->GetProdVertex()); // production vertex of the positive daughter track
1958  Char_t cTypeVtxProdPos = prodVtxDaughterPos->GetType(); // type of the production vertex
1959  AliAODVertex* prodVtxDaughterNeg = (AliAODVertex*)(trackNeg->GetProdVertex()); // production vertex of the negative daughter track
1960  Char_t cTypeVtxProdNeg = prodVtxDaughterNeg->GetType(); // type of the production vertex
1961 
1962 // fh2Tau3DVs2D[0]->Fill(dPtV0, dLOverP / dROverPt);
1963 
1964  // Cut vs mass histograms before cuts
1965  /*
1966  if(bIsCandidateK0s)
1967  {
1968  fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsPos);
1969  fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsNeg);
1970  fh2CutPtPosK0s[0]->Fill(dMassV0K0s, dPtDaughterPos);
1971  fh2CutPtNegK0s[0]->Fill(dMassV0K0s, dPtDaughterNeg);
1972  fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
1973  fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
1974  fh2CutDCAV0[0]->Fill(dMassV0K0s, dDCADaughters);
1975  fh2CutCos[0]->Fill(dMassV0K0s, dCPA);
1976  fh2CutR[0]->Fill(dMassV0K0s, dRadiusDecay);
1977  fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterPos);
1978  fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterNeg);
1979  fh2CutRapK0s[0]->Fill(dMassV0K0s, dRapK0s);
1980  fh2CutCTauK0s[0]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
1981  fh2CutPIDPosK0s[0]->Fill(dMassV0K0s, dNSigmaPosPion);
1982  fh2CutPIDNegK0s[0]->Fill(dMassV0K0s, dNSigmaNegPion);
1983  }
1984  if(bIsCandidateLambda)
1985  {
1986  fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsPos);
1987  fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsNeg);
1988  fh2CutPtPosLambda[0]->Fill(dMassV0Lambda, dPtDaughterPos);
1989  fh2CutPtNegLambda[0]->Fill(dMassV0Lambda, dPtDaughterNeg);
1990  fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterPos);
1991  fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterNeg);
1992  fh2CutRapLambda[0]->Fill(dMassV0Lambda, dRapLambda);
1993  fh2CutCTauLambda[0]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
1994  fh2CutPIDPosLambda[0]->Fill(dMassV0Lambda, dNSigmaPosProton);
1995  fh2CutPIDNegLambda[0]->Fill(dMassV0Lambda, dNSigmaNegPion);
1996  }
1997  */
1998 
1999  //===== Start of reconstruction cutting =====
2000 
2001  // 1
2002  // All V0 candidates
2003  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2004  iCutIndex++;
2005 
2006  // Start of global cuts
2007  // 2
2008  // Reconstruction method
2009  if(bPrintCuts) printf("Rec: Applying cut: Reconstruction method: %s\n", (fbOnFly ? "on-the-fly" : "offline"));
2010  if(bOnFlyStatus != fbOnFly)
2011  continue;
2012  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2013  iCutIndex++;
2014 
2015  // 3
2016  // Tracks TPC OK
2017  if(bPrintCuts) printf("Rec: Applying cut: Correct charge of daughters\n");
2018  if(!trackNeg || !trackPos)
2019  continue;
2020  if(trackNeg->Charge() == trackPos->Charge()) // daughters have different charge?
2021  continue;
2022  if(trackNeg->Charge() != -1) // daughters have expected charge?
2023  continue;
2024  if(trackPos->Charge() != 1) // daughters have expected charge?
2025  continue;
2026 
2027  if(fbTPCRefit)
2028  {
2029  if(bPrintCuts) printf("Rec: Applying cut: TPC refit\n");
2030  if(!trackNeg->IsOn(AliAODTrack::kTPCrefit)) // TPC refit is ON?
2031  continue;
2032  if(!trackPos->IsOn(AliAODTrack::kTPCrefit))
2033  continue;
2034  }
2035 
2036  if(fbRejectKinks)
2037  {
2038  if(bPrintCuts) printf("Rec: Applying cut: Type of production vertex of daughter: No kinks\n");
2039  if(cTypeVtxProdNeg == AliAODVertex::kKink) // kink daughter rejection
2040  continue;
2041  if(cTypeVtxProdPos == AliAODVertex::kKink)
2042  continue;
2043  }
2044 
2045  if(fbFindableClusters)
2046  {
2047  if(bPrintCuts) printf("Rec: Applying cut: Positive number of findable clusters\n");
2048  if(dFindableNeg <= 0.)
2049  continue;
2050  if(dFindablePos <= 0.)
2051  continue;
2052  }
2053 
2054  if(fdCutNCrossedRowsTPCMin > 0.)
2055  {
2056  if(bPrintCuts) printf("Rec: Applying cut: Number of TPC rows >= %g\n", fdCutNCrossedRowsTPCMin);
2057  if(dNRowsNeg < fdCutNCrossedRowsTPCMin) // Crossed TPC padrows
2058  continue;
2059  if(dNRowsPos < fdCutNCrossedRowsTPCMin)
2060  continue;
2061  }
2062 
2064  {
2065  if(bPrintCuts) printf("Rec: Applying cut: rows/findable >= %g\n", fdCutCrossedRowsOverFindMin);
2066  if(dNRowsNeg / dFindableNeg < fdCutCrossedRowsOverFindMin)
2067  continue;
2068  if(dNRowsPos / dFindablePos < fdCutCrossedRowsOverFindMin)
2069  continue;
2070  }
2071 
2073  {
2074  if(bPrintCuts) printf("Rec: Applying cut: rows/findable <= %g\n", fdCutCrossedRowsOverFindMax);
2075  if(dNRowsNeg / dFindableNeg > fdCutCrossedRowsOverFindMax)
2076  continue;
2077  if(dNRowsPos / dFindablePos > fdCutCrossedRowsOverFindMax)
2078  continue;
2079  }
2080 
2081  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2082  iCutIndex++;
2083 
2084  // 4
2085  // Daughters: transverse momentum cut
2086  if(fdCutPtDaughterMin > 0.)
2087  {
2088  if(bPrintCuts) printf("Rec: Applying cut: Daughter pt >= %g\n", fdCutPtDaughterMin);
2089  if((dPtDaughterNeg < fdCutPtDaughterMin) || (dPtDaughterPos < fdCutPtDaughterMin))
2090  continue;
2091  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2092  }
2093  iCutIndex++;
2094 
2095  // 5
2096  // Daughters: Impact parameter of daughters to prim vtx
2097  if(fdCutDCAToPrimVtxMin > 0.)
2098  {
2099  if(bPrintCuts) printf("Rec: Applying cut: Daughter DCA to prim vtx >= %g\n", fdCutDCAToPrimVtxMin);
2100  if((dDCAToPrimVtxNeg < fdCutDCAToPrimVtxMin) || (dDCAToPrimVtxPos < fdCutDCAToPrimVtxMin))
2101  continue;
2102  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2103  }
2104  iCutIndex++;
2105 
2106  // 6
2107  // Daughters: DCA
2108  if(fdCutDCADaughtersMax > 0.)
2109  {
2110  if(bPrintCuts) printf("Rec: Applying cut: DCA between daughters <= %g\n", fdCutDCADaughtersMax);
2111  if(dDCADaughters > fdCutDCADaughtersMax)
2112  continue;
2113  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2114  }
2115  iCutIndex++;
2116 
2117  // 7
2118  // V0: Cosine of the pointing angle
2119  if(fdCutCPAKMin > 0.)
2120  {
2121  if(bPrintCuts) printf("Rec: Applying cut: CPA >= %g (K)\n", fdCutCPAKMin);
2122  if(dCPA < fdCutCPAKMin)
2123  bIsCandidateK0s = kFALSE;
2124  }
2125  if(fdCutCPALMin > 0.)
2126  {
2127  if(bPrintCuts) printf("Rec: Applying cut: CPA >= %g (L, AL)\n", fdCutCPALMin);
2128  if(dCPA < fdCutCPALMin)
2129  {
2130  bIsCandidateLambda = kFALSE;
2131  bIsCandidateALambda = kFALSE;
2132  }
2133  }
2134  if(fdCutCPAKMin > 0. || fdCutCPALMin > 0.)
2135  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2136  iCutIndex++;
2137 
2138  // 8
2139  // V0: Fiducial volume
2140  if(fdCutRadiusDecayMin > 0. && fdCutRadiusDecayMax > 0.)
2141  {
2142  if(bPrintCuts) printf("Rec: Applying cut: Decay radius >= %g, <= %g\n", fdCutRadiusDecayMin, fdCutRadiusDecayMax);
2143  if((dRadiusDecay < fdCutRadiusDecayMin) || (dRadiusDecay > fdCutRadiusDecayMax))
2144  continue;
2145  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2146  }
2147  iCutIndex++;
2148 
2149  // 9
2150  // Daughters: pseudorapidity cut
2151  if(fdCutEtaDaughterMax > 0.)
2152  {
2153  if(bPrintCuts) printf("Rec: Applying cut: Daughter |eta| < %g\n", fdCutEtaDaughterMax);
2154  if((TMath::Abs(dEtaDaughterNeg) > fdCutEtaDaughterMax) || (TMath::Abs(dEtaDaughterPos) > fdCutEtaDaughterMax))
2155  continue;
2156  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2157  }
2158  iCutIndex++;
2159  // End of global cuts
2160 
2161  // Start of particle-dependent cuts
2162  // 10
2163  // V0: pseudorapidity cut & rapidity cut
2164  if(fdCutEtaV0Max > 0.)
2165  {
2166  if(bPrintCuts) printf("Rec: Applying cut: V0 |eta| < %g\n", fdCutEtaV0Max);
2167  if(TMath::Abs(dEtaV0) > fdCutEtaV0Max)
2168  {
2169  bIsCandidateK0s = kFALSE;
2170  bIsCandidateLambda = kFALSE;
2171  bIsCandidateALambda = kFALSE;
2172  }
2173  }
2174  if(fdCutRapV0Max > 0.)
2175  {
2176  if(bPrintCuts) printf("Rec: Applying cut: V0 |y| < %g\n", fdCutRapV0Max);
2177  if(TMath::Abs(dRapK0s) > fdCutRapV0Max)
2178  bIsCandidateK0s = kFALSE;
2179  if(TMath::Abs(dRapLambda) > fdCutRapV0Max)
2180  {
2181  bIsCandidateLambda = kFALSE;
2182  bIsCandidateALambda = kFALSE;
2183  }
2184  }
2185  if(fdCutEtaV0Max > 0. || fdCutRapV0Max > 0.)
2186  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2187  iCutIndex++;
2188 
2189  // 11
2190  // Lifetime cut
2191  if(fdCutNTauKMax > 0.)
2192  {
2193  if(bPrintCuts) printf("Rec: Applying cut: Proper lifetime < %g (K)\n", fdCutNTauKMax);
2194  if(dMROverPtK0s > fdCutNTauKMax * dCTauK0s)
2195  bIsCandidateK0s = kFALSE;
2196  }
2197  if(fdCutNTauLMax > 0.)
2198  {
2199  if(bPrintCuts) printf("Rec: Applying cut: Proper lifetime < %g (L, AL)\n", fdCutNTauLMax);
2200  if(dMROverPtLambda > fdCutNTauLMax * dCTauLambda)
2201  {
2202  bIsCandidateLambda = kFALSE;
2203  bIsCandidateALambda = kFALSE;
2204  }
2205  }
2206  if(fdCutNTauKMax > 0. || fdCutNTauLMax > 0.)
2207  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2208  iCutIndex++;
2209 
2210  // 12
2211  // Daughter PID
2212  if(fdCutNSigmadEdxMax > 0.)
2213  {
2214  if(fbIsPbPb && fdPtProtonPIDMax > 0.) // Pb-Pb
2215  {
2216  if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (proton below %g GeV/c) < %g\n", fdPtProtonPIDMax, fdCutNSigmadEdxMax);
2217  if((dPtDaughterPos < fdPtProtonPIDMax) && (dNSigmaPosProton > fdCutNSigmadEdxMax)) // p+
2218  bIsCandidateLambda = kFALSE;
2219  if((dPtDaughterNeg < fdPtProtonPIDMax) && (dNSigmaNegProton > fdCutNSigmadEdxMax)) // p-
2220  bIsCandidateALambda = kFALSE;
2221  }
2222  else // p-p
2223  {
2224  if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (both daughters): < %g\n", fdCutNSigmadEdxMax);
2225  if(dNSigmaPosPion > fdCutNSigmadEdxMax || dNSigmaNegPion > fdCutNSigmadEdxMax) // pi+, pi-
2226  bIsCandidateK0s = kFALSE;
2227  if(dNSigmaPosProton > fdCutNSigmadEdxMax || dNSigmaNegPion > fdCutNSigmadEdxMax) // p+, pi-
2228  bIsCandidateLambda = kFALSE;
2229  if(dNSigmaNegProton > fdCutNSigmadEdxMax || dNSigmaPosPion > fdCutNSigmadEdxMax) // p-, pi+
2230  bIsCandidateALambda = kFALSE;
2231  }
2232  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2233  }
2234  iCutIndex++;
2235 
2236  Double_t valueCorrel[3] = {dMassV0K0s, dMassV0Lambda, dPtV0};
2237  if(bIsCandidateK0s && bIsCandidateLambda)
2238  fh3CCMassCorrelBoth->Fill(valueCorrel); // correlation of mass distribution of candidates selected as both K0s and Lambda
2239  if(bIsCandidateK0s && !bIsCandidateLambda)
2240  fh3CCMassCorrelKNotL->Fill(valueCorrel); // correlation of mass distribution of candidates selected as K0s and not Lambda
2241  if(!bIsCandidateK0s && bIsCandidateLambda)
2242  fh3CCMassCorrelLNotK->Fill(valueCorrel); // correlation of mass distribution of candidates selected as not K0s and Lambda
2243 
2244  // 13
2245  // Armenteros-Podolanski cut
2246  if(fbCutArmPod)
2247  {
2248  if(bPrintCuts) printf("Rec: Applying cut: Armenteros-Podolanski (K0S) pT > %g * |alpha|\n", 0.2);
2249  if(dPtArm < TMath::Abs(0.2 * dAlpha))
2250  bIsCandidateK0s = kFALSE;
2251  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2252  }
2253  iCutIndex++;
2254 
2255  // 14
2256  // Cross-contamination
2257  /*
2258  if(bIsInPeakK0s)
2259  {
2260  if(bIsCandidateLambda) // Lambda candidates in K0s peak, excluded from Lambda candidates by CC cut
2261  fh2CCLambda->Fill(dMassV0Lambda, dPtV0);
2262  }
2263  if(bIsInPeakLambda)
2264  {
2265  if(bIsCandidateK0s) // K0s candidates in Lambda peak, excluded from K0s candidates by CC cut
2266  fh2CCK0s->Fill(dMassV0K0s, dPtV0);
2267  }
2268  */
2269  if(fbCutCross)
2270  {
2271  if(bIsInPeakK0s)
2272  {
2273  bIsCandidateLambda = kFALSE;
2274  bIsCandidateALambda = kFALSE;
2275  }
2276  if(bIsInPeakLambda)
2277  {
2278  bIsCandidateK0s = kFALSE;
2279  }
2280  if(bIsInPeakALambda)
2281  {
2282  bIsCandidateK0s = kFALSE;
2283  }
2284  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2285  }
2286  iCutIndex++;
2287  // End of particle-dependent cuts
2288 
2289  //===== End of reconstruction cutting =====
2290 
2291  if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
2292  continue;
2293 
2294  // Selection of V0s in jet cones, perpendicular cones, random cones, outside cones
2295  if(iNJetSel && (bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda))
2296  {
2297  Double_t dDMin, dD = 0; // minimal / current value of V0-jet Distance (used for estimation of closest jet to V0)
2298  dDMin = 20.;
2299  for(Int_t iJet = 0; iJet < iNJetSel; iJet++) // could be included in loop beneath
2300  {
2301  // finding the closest jet to the v0
2302  jet = (AliAODJet*)arrayJetSel->At(iJet); // load a jet in the list
2303  if(!jet)
2304  continue;
2305  dD = GetD(v0, jet);
2306  if(dD < dDMin)
2307  dDMin = dD;
2308  }
2309  if(bIsCandidateK0s)
2310  fh1DistanceV0JetsK0s[iCentIndex]->Fill(dDMin);
2311  if(bIsCandidateLambda)
2312  fh1DistanceV0JetsLambda[iCentIndex]->Fill(dDMin);
2313  if(bIsCandidateALambda)
2314  fh1DistanceV0JetsALambda[iCentIndex]->Fill(dDMin);
2315 
2316  // Selection of V0s in jet cones
2317  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for V0 %d %d in %d jet cones", bIsCandidateK0s, bIsCandidateLambda, iNJetSel));
2318  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2319  {
2320  jet = (AliAODJet*)arrayJetSel->At(iJet); // load a jet in the list
2321  if(!jet)
2322  continue;
2323  vecJetMomentum.SetXYZ(jet->Px(), jet->Py(), jet->Pz()); // set the vector of jet momentum
2324  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Checking if V0 %d %d in jet cone %d", bIsCandidateK0s, bIsCandidateLambda, iJet));
2325  if(IsParticleInCone(v0, jet, fdDistanceV0JetMax)) // If good jet in event, find out whether V0 is in that jet
2326  {
2327  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("V0 %d %d found in jet cone %d", bIsCandidateK0s, bIsCandidateLambda, iJet));
2328  bIsInConeJet = kTRUE;
2329  dPtJetTrackLeading = jet->GetPtLeading();
2330  dPtJet = jet->Pt();
2331  break;
2332  }
2333  }
2334  // Selection of V0s in perp. cones
2335  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for V0 %d %d in %d perp. cones", bIsCandidateK0s, bIsCandidateLambda, iNJetPerp));
2336  for(Int_t iJet = 0; iJet < iNJetPerp; iJet++)
2337  {
2338  jetPerp = (AliAODJet*)arrayJetPerp->At(iJet); // load a jet in the list
2339  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Checking if V0 %d %d in perp. cone %d", bIsCandidateK0s, bIsCandidateLambda, iJet));
2340  if(IsParticleInCone(v0, jetPerp, fdDistanceV0JetMax)) // V0 in perp. cone
2341  {
2342  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("V0 %d %d found in perp. cone %d", bIsCandidateK0s, bIsCandidateLambda, iJet));
2343  bIsInConePerp = kTRUE;
2344  break;
2345  }
2346  }
2347  // Selection of V0s in random cones
2348  if(jetRnd)
2349  {
2350  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for V0 %d %d in the rnd. cone", bIsCandidateK0s, bIsCandidateLambda));
2351  if(IsParticleInCone(v0, jetRnd, fdDistanceV0JetMax)) // V0 in rnd. cone?
2352  {
2353  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("V0 %d %d found in the rnd. cone", bIsCandidateK0s, bIsCandidateLambda));
2354  bIsInConeRnd = kTRUE;
2355  }
2356  }
2357  // Selection of V0s in median-cluster cones
2358  if(jetMed)
2359  {
2360  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for V0 %d %d in the med. cone", bIsCandidateK0s, bIsCandidateLambda));
2361  if(IsParticleInCone(v0, jetMed, fdDistanceV0JetMax)) // V0 in med. cone?
2362  {
2363  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("V0 %d %d found in the med. cone", bIsCandidateK0s, bIsCandidateLambda));
2364  bIsInConeMed = kTRUE;
2365  }
2366  }
2367  // Selection of V0s outside jet cones
2368  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for V0 %d %d outside jet cones", bIsCandidateK0s, bIsCandidateLambda));
2369  if(!OverlapWithJets(arrayJetSel, v0, dRadiusExcludeCone)) // V0 oustide jet cones
2370  {
2371  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("V0 %d %d found outside jet cones", bIsCandidateK0s, bIsCandidateLambda));
2372  bIsOutsideCones = kTRUE;
2373  }
2374  }
2375 
2376  // QA histograms after cuts
2377  FillQAHistogramV0(primVtx, v0, 1, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, bIsInPeakK0s, bIsInPeakLambda, bIsInPeakALambda);
2378  // Cut vs mass histograms after cuts
2379  /*
2380  if(bIsCandidateK0s)
2381  {
2382  fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsPos);
2383  fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsNeg);
2384  fh2CutPtPosK0s[1]->Fill(dMassV0K0s, dPtDaughterPos);
2385  fh2CutPtNegK0s[1]->Fill(dMassV0K0s, dPtDaughterNeg);
2386  fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
2387  fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
2388  fh2CutDCAV0[1]->Fill(dMassV0K0s, dDCADaughters);
2389  fh2CutCos[1]->Fill(dMassV0K0s, dCPA);
2390  fh2CutR[1]->Fill(dMassV0K0s, dRadiusDecay);
2391  fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterPos);
2392  fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterNeg);
2393  fh2CutRapK0s[1]->Fill(dMassV0K0s, dRapK0s);
2394  fh2CutCTauK0s[1]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
2395  fh2CutPIDPosK0s[1]->Fill(dMassV0K0s, dNSigmaPosPion);
2396  fh2CutPIDNegK0s[1]->Fill(dMassV0K0s, dNSigmaNegPion);
2397  }
2398  if(bIsCandidateLambda)
2399  {
2400  fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsPos);
2401  fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsNeg);
2402  fh2CutPtPosLambda[1]->Fill(dMassV0Lambda, dPtDaughterPos);
2403  fh2CutPtNegLambda[1]->Fill(dMassV0Lambda, dPtDaughterNeg);
2404  fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterPos);
2405  fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterNeg);
2406  fh2CutRapLambda[1]->Fill(dMassV0Lambda, dRapLambda);
2407  fh2CutCTauLambda[1]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
2408  fh2CutPIDPosLambda[1]->Fill(dMassV0Lambda, dNSigmaPosProton);
2409  fh2CutPIDNegLambda[1]->Fill(dMassV0Lambda, dNSigmaNegPion);
2410  }
2411  */
2412 
2413  //===== Start of filling V0 spectra =====
2414 
2415  Double_t dAngle = TMath::Pi(); // angle between V0 momentum and jet momentum
2416  if(bIsInConeJet)
2417  {
2418  dAngle = vecV0Momentum.Angle(vecJetMomentum);
2419  }
2420 
2421  // iCutIndex = 15
2422  if(bIsCandidateK0s)
2423  {
2424  // 15 K0s candidates after cuts
2425 // printf("K0S: i = %d, m = %g, pT = %g, eta = %g, phi = %g\n",iV0,dMassV0K0s,dPtV0,dEtaV0,dPhiV0);
2426  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex, iCentIndex);
2427  Double_t valueKIncl[3] = {dMassV0K0s, dPtV0, dEtaV0};
2428  fhnV0InclusiveK0s[iCentIndex]->Fill(valueKIncl);
2429  fh1V0InvMassK0sCent[iCentIndex]->Fill(dMassV0K0s);
2430 
2431  fh1QACTau2D[1]->Fill(dMROverPtK0s / dCTauK0s);
2432  fh1QACTau3D[1]->Fill(dMLOverPK0s / dCTauK0s);
2433 // fh2Tau3DVs2D[1]->Fill(dPtV0, dLOverP / dROverPt);
2434 
2435  if(iNJetSel)
2436  {
2437  // 16 K0s in jet events
2438  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 1, iCentIndex);
2439  }
2440  if(bIsInConeJet)
2441  {
2442  // 17 K0s in jets
2443  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 2, iCentIndex);
2444  Double_t valueKInJC[4] = {dMassV0K0s, dPtV0, dEtaV0, jet->Pt()};
2445  fhnV0InJetK0s[iCentIndex]->Fill(valueKInJC);
2446  fh2V0PtJetAngleK0s[iCentIndex]->Fill(jet->Pt(), dAngle);
2447  Double_t valuesDaughter[5] = {dPtDaughterPos, dPtDaughterNeg, dPtV0, dPtJet, dPtJetTrackLeading};
2448  fhnPtDaughterK0s[iCentIndex]->Fill(valuesDaughter);
2449  }
2450  if(bIsOutsideCones)
2451  {
2452  Double_t valueKOutJC[3] = {dMassV0K0s, dPtV0, dEtaV0};
2453  fhnV0OutJetK0s[iCentIndex]->Fill(valueKOutJC);
2454  }
2455  if(bIsInConePerp)
2456  {
2457  Double_t valueKInPC[4] = {dMassV0K0s, dPtV0, dEtaV0, jetPerp->Pt()};
2458  fhnV0InPerpK0s[iCentIndex]->Fill(valueKInPC);
2459  }
2460  if(bIsInConeRnd)
2461  {
2462  Double_t valueKInRnd[3] = {dMassV0K0s, dPtV0, dEtaV0};
2463  fhnV0InRndK0s[iCentIndex]->Fill(valueKInRnd);
2464  }
2465  if(bIsInConeMed)
2466  {
2467  Double_t valueKInMed[3] = {dMassV0K0s, dPtV0, dEtaV0};
2468  fhnV0InMedK0s[iCentIndex]->Fill(valueKInMed);
2469  }
2470  if(!iNJetSel)
2471  {
2472  Double_t valueKNoJet[3] = {dMassV0K0s, dPtV0, dEtaV0};
2473  fhnV0NoJetK0s[iCentIndex]->Fill(valueKNoJet);
2474  }
2475  iNV0CandK0s++;
2476  }
2477  if(bIsCandidateLambda)
2478  {
2479  // 15 Lambda candidates after cuts
2480 // printf("La: i = %d, m = %g, pT = %g, eta = %g, phi = %g\n",iV0,dMassV0Lambda,dPtV0,dEtaV0,dPhiV0);
2481  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex, iCentIndex);
2482  Double_t valueLIncl[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2483  fhnV0InclusiveLambda[iCentIndex]->Fill(valueLIncl);
2484  fh1V0InvMassLambdaCent[iCentIndex]->Fill(dMassV0Lambda);
2485  if(iNJetSel)
2486  {
2487  // 16 Lambda in jet events
2488  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 1, iCentIndex);
2489  }
2490  if(bIsInConeJet)
2491  {
2492  // 17 Lambda in jets
2493  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 2, iCentIndex);
2494  Double_t valueLInJC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jet->Pt()};
2495  fhnV0InJetLambda[iCentIndex]->Fill(valueLInJC);
2496  fh2V0PtJetAngleLambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2497  Double_t valuesDaughter[5] = {dPtDaughterPos, dPtDaughterNeg, dPtV0, dPtJet, dPtJetTrackLeading};
2498  fhnPtDaughterLambda[iCentIndex]->Fill(valuesDaughter);
2499  }
2500  if(bIsOutsideCones)
2501  {
2502  Double_t valueLOutJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2503  fhnV0OutJetLambda[iCentIndex]->Fill(valueLOutJet);
2504  }
2505  if(bIsInConePerp)
2506  {
2507  Double_t valueLInPC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jetPerp->Pt()};
2508  fhnV0InPerpLambda[iCentIndex]->Fill(valueLInPC);
2509  }
2510  if(bIsInConeRnd)
2511  {
2512  Double_t valueLInRnd[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2513  fhnV0InRndLambda[iCentIndex]->Fill(valueLInRnd);
2514  }
2515  if(bIsInConeMed)
2516  {
2517  Double_t valueLInMed[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2518  fhnV0InMedLambda[iCentIndex]->Fill(valueLInMed);
2519  }
2520  if(!iNJetSel)
2521  {
2522  Double_t valueLNoJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2523  fhnV0NoJetLambda[iCentIndex]->Fill(valueLNoJet);
2524  }
2525  iNV0CandLambda++;
2526  }
2527  if(bIsCandidateALambda)
2528  {
2529  // 15 ALambda candidates after cuts
2530 // printf("AL: i = %d, m = %g, pT = %g, eta = %g, phi = %g\n",iV0,dMassV0ALambda,dPtV0,dEtaV0,dPhiV0);
2531  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex, iCentIndex);
2532  Double_t valueALIncl[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2533  fhnV0InclusiveALambda[iCentIndex]->Fill(valueALIncl);
2534  fh1V0InvMassALambdaCent[iCentIndex]->Fill(dMassV0ALambda);
2535  if(iNJetSel)
2536  {
2537  // 16 ALambda in jet events
2538  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 1, iCentIndex);
2539  }
2540  if(bIsInConeJet)
2541  {
2542  // 17 ALambda in jets
2543  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 2, iCentIndex);
2544  Double_t valueLInJC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jet->Pt()};
2545  fhnV0InJetALambda[iCentIndex]->Fill(valueLInJC);
2546  fh2V0PtJetAngleALambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2547  Double_t valuesDaughter[5] = {dPtDaughterPos, dPtDaughterNeg, dPtV0, dPtJet, dPtJetTrackLeading};
2548  fhnPtDaughterALambda[iCentIndex]->Fill(valuesDaughter);
2549  }
2550  if(bIsOutsideCones)
2551  {
2552  Double_t valueALOutJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2553  fhnV0OutJetALambda[iCentIndex]->Fill(valueALOutJet);
2554  }
2555  if(bIsInConePerp)
2556  {
2557  Double_t valueLInPC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jetPerp->Pt()};
2558  fhnV0InPerpALambda[iCentIndex]->Fill(valueLInPC);
2559  }
2560  if(bIsInConeRnd)
2561  {
2562  Double_t valueALInRnd[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2563  fhnV0InRndALambda[iCentIndex]->Fill(valueALInRnd);
2564  }
2565  if(bIsInConeMed)
2566  {
2567  Double_t valueALInMed[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2568  fhnV0InMedALambda[iCentIndex]->Fill(valueALInMed);
2569  }
2570  if(!iNJetSel)
2571  {
2572  Double_t valueALNoJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2573  fhnV0NoJetALambda[iCentIndex]->Fill(valueALNoJet);
2574  }
2575  iNV0CandALambda++;
2576  }
2577  // V0-jet correlations
2578  if(fbCorrelations && iNJetSel)
2579  {
2580  Double_t dDPhi, dDEta;
2581  // Fill V0-jet correlations in same events
2582  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2583  {
2584  AliAODJet* jetCorrel = (AliAODJet*)arrayJetSel->At(iJet); // load a jet in the list
2585  dDPhi = GetNormalPhi(dPhiV0 - jetCorrel->Phi());
2586  dDEta = dEtaV0 - jetCorrel->Eta();
2587  if(TMath::Abs(dDEta) > fdDeltaEtaMax)
2588  continue;
2589  if(bIsCandidateK0s)
2590  {
2591  Double_t valueKCorrel[6] = {dMassV0K0s, dPtV0, dEtaV0, jetCorrel->Pt(), dDPhi, dDEta};
2592  fhnV0CorrelSEK0s[iCentIndex]->Fill(valueKCorrel);
2593  }
2594  if(bIsCandidateLambda)
2595  {
2596  Double_t valueLCorrel[6] = {dMassV0Lambda, dPtV0, dEtaV0, jetCorrel->Pt(), dDPhi, dDEta};
2597  fhnV0CorrelSELambda[iCentIndex]->Fill(valueLCorrel);
2598  }
2599  }
2600  // Fill V0-jet correlations in mixed events
2601  if(bPoolReady)
2602  {
2603  for(Int_t iMix = 0; iMix < pool->GetCurrentNEvents(); iMix++)
2604  {
2605  TObjArray* arrayMixedEvent = pool->GetEvent(iMix);
2606  if(!arrayMixedEvent)
2607  continue;
2608  for(Int_t iJet = 0; iJet < arrayMixedEvent->GetEntriesFast(); iJet++)
2609  {
2610  TLorentzVector* jetMixed = (TLorentzVector*)arrayMixedEvent->At(iJet);
2611  if(!jetMixed)
2612  continue;
2613  dDPhi = GetNormalPhi(dPhiV0 - jetMixed->Phi());
2614  dDEta = dEtaV0 - jetMixed->Eta();
2615  if(TMath::Abs(dDEta) > fdDeltaEtaMax)
2616  continue;
2617  if(bIsCandidateK0s)
2618  {
2619  Double_t valueKCorrel[6] = {dMassV0K0s, dPtV0, dEtaV0, jetMixed->Pt(), dDPhi, dDEta};
2620  fhnV0CorrelMEK0s[iCentIndex]->Fill(valueKCorrel);
2621  }
2622  if(bIsCandidateLambda)
2623  {
2624  Double_t valueLCorrel[6] = {dMassV0Lambda, dPtV0, dEtaV0, jetMixed->Pt(), dDPhi, dDEta};
2625  fhnV0CorrelMELambda[iCentIndex]->Fill(valueLCorrel);
2626  }
2627  }
2628  }
2629  }
2630  }
2631  //===== End of filling V0 spectra =====
2632 
2633 
2634  //===== Association of reconstructed V0 candidates with MC particles =====
2635  if(fbMCAnalysis)
2636  {
2637  // Associate selected candidates only
2638 // if ( !(bIsCandidateK0s && bIsInPeakK0s) && !(bIsCandidateLambda && bIsInPeakLambda) ) // signal candidates
2639  if(!(bIsCandidateK0s) && !(bIsCandidateLambda) && !(bIsCandidateALambda)) // chosen candidates with any mass
2640  continue;
2641 
2642  // Get MC labels of reconstructed daughter tracks
2643  Int_t iLabelPos = TMath::Abs(trackPos->GetLabel());
2644  Int_t iLabelNeg = TMath::Abs(trackNeg->GetLabel());
2645 
2646  // Make sure MC daughters are in the array range
2647  if((iLabelNeg < 0) || (iLabelNeg >= iNTracksMC) || (iLabelPos < 0) || (iLabelPos >= iNTracksMC))
2648  continue;
2649 
2650  // Get MC particles corresponding to reconstructed daughter tracks
2651  AliAODMCParticle* particleMCDaughterNeg = (AliAODMCParticle*)arrayMC->At(iLabelNeg);
2652  AliAODMCParticle* particleMCDaughterPos = (AliAODMCParticle*)arrayMC->At(iLabelPos);
2653  if(!particleMCDaughterNeg || !particleMCDaughterPos)
2654  continue;
2655 
2656  // Make sure MC daughter particles are not physical primary
2657  if((particleMCDaughterNeg->IsPhysicalPrimary()) || (particleMCDaughterPos->IsPhysicalPrimary()))
2658  continue;
2659 
2660  // Get identities of MC daughter particles
2661  Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
2662  Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
2663 
2664  // Get index of the mother particle for each MC daughter particle
2665  Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
2666  Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
2667 
2668  if((iIndexMotherNeg < 0) || (iIndexMotherNeg >= iNTracksMC) || (iIndexMotherPos < 0) || (iIndexMotherPos >= iNTracksMC))
2669  continue;
2670 
2671  // Check whether MC daughter particles have the same mother
2672  if(iIndexMotherNeg != iIndexMotherPos)
2673  continue;
2674 
2675  // Get the MC mother particle of both MC daughter particles
2676  AliAODMCParticle* particleMCMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherPos);
2677  if(!particleMCMother)
2678  continue;
2679 
2680  // Get identity of the MC mother particle
2681  Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
2682 
2683  // Skip not interesting particles
2684  if((iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda))
2685  continue;
2686 
2687  // Check identity of the MC mother particle and the decay channel
2688  // Is MC mother particle K0S?
2689  Bool_t bV0MCIsK0s = ((iPdgCodeMother == iPdgCodeK0s) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2690  // Is MC mother particle Lambda?
2691  Bool_t bV0MCIsLambda = ((iPdgCodeMother == +iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodeProton) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2692  // Is MC mother particle anti-Lambda?
2693  Bool_t bV0MCIsALambda = ((iPdgCodeMother == -iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodeProton));
2694 
2695  Double_t dPtV0Gen = particleMCMother->Pt();
2696  Double_t dRapV0Gen = particleMCMother->Y();
2697  Double_t dEtaV0Gen = particleMCMother->Eta();
2698 // Double_t dPhiV0Gen = particleMCMother->Phi();
2699 
2700  // V0 pseudorapidity cut applied on generated particles
2701  if(fdCutEtaV0Max > 0.)
2702  {
2703  if(bPrintCuts) printf("Rec->Gen: Applying cut: V0 |eta|: < %g\n", fdCutEtaV0Max);
2704  if((TMath::Abs(dEtaV0Gen) > fdCutEtaV0Max))
2705  continue;
2706  }
2707  // V0 rapidity cut applied on generated particles
2708  if(fdCutRapV0Max > 0.)
2709  {
2710  if(bPrintCuts) printf("Rec->Gen: Applying cut: V0 |y|: < %g\n", fdCutRapV0Max);
2711  if((TMath::Abs(dRapV0Gen) > fdCutRapV0Max))
2712  continue;
2713  }
2714 
2715  // Select only particles from a specific generator
2716  if(!IsFromGoodGenerator(iIndexMotherPos))
2717  continue;
2718 
2719  // Is MC mother particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2720 // Bool_t bV0MCIsPrimary = particleMCMother->IsPhysicalPrimary();
2721  // Get the MC mother particle of the MC mother particle
2722  Int_t iIndexMotherOfMother = particleMCMother->GetMother();
2723  AliAODMCParticle* particleMCMotherOfMother = 0;
2724  if(iIndexMotherOfMother >= 0)
2725  particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2726  // Get identity of the MC mother particle of the MC mother particle if it exists
2727  Int_t iPdgCodeMotherOfMother = 0;
2728  if(particleMCMotherOfMother)
2729  iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2730  // 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*-)
2731 // Bool_t bV0MCComesFromSigma = kFALSE; // Is MC mother particle daughter of a Sigma?
2732 // if ( (particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && ( (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) ) )
2733 // bV0MCComesFromSigma = kTRUE;
2734  // Should MC mother particle be considered as primary when it is Lambda?
2735 // Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2736  // Check if the MC mother particle of the MC mother particle is a Xi (3322 - Xi0, 3312 - Xi-)
2737  Bool_t bV0MCComesFromXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == 3322) || (iPdgCodeMotherOfMother == 3312))); // Is MC mother particle daughter of a Xi?
2738  Bool_t bV0MCComesFromAXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == -3322) || (iPdgCodeMotherOfMother == -3312))); // Is MC mother particle daughter of a anti-Xi?
2739 
2740  // Get the distance between production point of the MC mother particle and the primary vertex
2741  Double_t dx = dPrimVtxMCX - particleMCMother->Xv();
2742  Double_t dy = dPrimVtxMCY - particleMCMother->Yv();
2743  Double_t dz = dPrimVtxMCZ - particleMCMother->Zv();
2744  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2745  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2746 
2747  // phi, eta resolution for V0-reconstruction
2748 // Double_t dResolutionV0Eta = particleMCMother->Eta()-v0->Eta();
2749 // Double_t dResolutionV0Phi = particleMCMother->Phi()-v0->Phi();
2750 
2751  // K0s
2752 // if (bIsCandidateK0s && bIsInPeakK0s) // selected candidates in peak
2753  if(bIsCandidateK0s) // selected candidates with any mass
2754  {
2755 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2756  if(bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2757  {
2758  fh2V0K0sPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0K0s);
2759  Double_t valueEtaK[3] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen};
2760  fh3V0K0sEtaPtMassMCRec[iCentIndex]->Fill(valueEtaK);
2761 
2762  Double_t valueEtaDKNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2763  fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKNeg);
2764  Double_t valueEtaDKPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2765  fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKPos);
2766 
2767  fh2V0K0sMCResolMPt[iCentIndex]->Fill(dMassV0K0s - dMassPDGK0s, dPtV0);
2768  fh2V0K0sMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2769  if(bIsInConeJet) // true V0 associated to a candidate in jet
2770  {
2771  Double_t valueKInJCMC[4] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2772  fh3V0K0sInJetPtMassMCRec[iCentIndex]->Fill(valueKInJCMC);
2773  Double_t valueEtaKIn[5] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2774  fh4V0K0sInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaKIn);
2775 
2776  Double_t valueEtaDKJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2777  fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCNeg);
2778  Double_t valueEtaDKJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2779  fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCPos);
2780  }
2781  }
2782  if(bV0MCIsK0s && !bV0MCIsPrimaryDist) // not primary K0s
2783  {
2784  fh1V0K0sPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2785  }
2786  }
2787  // Lambda
2788 // if (bIsCandidateLambda && bIsInPeakLambda) // selected candidates in peak
2789  if(bIsCandidateLambda) // selected candidates with any mass
2790  {
2791 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2792  if(bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2793  {
2794  fh2V0LambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0Lambda);
2795  Double_t valueEtaL[3] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen};
2796  fh3V0LambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaL);
2797 
2798  Double_t valueEtaDLNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2799  fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLNeg);
2800  Double_t valueEtaDLPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2801  fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLPos);
2802 
2803  fh2V0LambdaMCResolMPt[iCentIndex]->Fill(dMassV0Lambda - dMassPDGLambda, dPtV0);
2804  fh2V0LambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2805  if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2806  {
2807  Double_t valueLInJCMC[4] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2808  fh3V0LambdaInJetPtMassMCRec[iCentIndex]->Fill(valueLInJCMC);
2809  Double_t valueEtaLIn[5] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2810  fh4V0LambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaLIn);
2811 
2812  Double_t valueEtaDLJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2813  fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCNeg);
2814  Double_t valueEtaDLJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2815  fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCPos);
2816  }
2817  }
2818  // Fill the feed-down histograms
2819  if(bV0MCIsLambda && bV0MCComesFromXi)
2820  {
2821  Double_t valueFDLIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2822  fhnV0LambdaInclMCFD[iCentIndex]->Fill(valueFDLIncl);
2823  if(bIsInConeRnd)
2824  {
2825  fhnV0LambdaBulkMCFD[iCentIndex]->Fill(valueFDLIncl);
2826  }
2827  if(bIsInConeJet)
2828  {
2829  Double_t valueFDLInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2830  fhnV0LambdaInJetsMCFD[iCentIndex]->Fill(valueFDLInJets);
2831  }
2832  }
2833  if(bV0MCIsLambda && !bV0MCIsPrimaryDist && !bV0MCComesFromXi) // not primary Lambda
2834  {
2835  fh1V0LambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2836  }
2837  }
2838  // anti-Lambda
2839 // if (bIsCandidateALambda && bIsInPeakALambda) // selected candidates in peak
2840  if(bIsCandidateALambda) // selected candidates with any mass
2841  {
2842 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2843  if(bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2844  {
2845  fh2V0ALambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0ALambda);
2846  Double_t valueEtaAL[3] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen};
2847  fh3V0ALambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaAL);
2848 
2849  Double_t valueEtaDALNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2850  fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALNeg);
2851  Double_t valueEtaDALPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2852  fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALPos);
2853 
2854  fh2V0ALambdaMCResolMPt[iCentIndex]->Fill(dMassV0ALambda - dMassPDGLambda, dPtV0);
2855  fh2V0ALambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2856  if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2857  {
2858  Double_t valueALInJCMC[4] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2859  fh3V0ALambdaInJetPtMassMCRec[iCentIndex]->Fill(valueALInJCMC);
2860  Double_t valueEtaALIn[5] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2861  fh4V0ALambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaALIn);
2862 
2863  Double_t valueEtaDALJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2864  fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCNeg);
2865  Double_t valueEtaDALJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2866  fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCPos);
2867  }
2868  }
2869  // Fill the feed-down histograms
2870  if(bV0MCIsALambda && bV0MCComesFromAXi)
2871  {
2872  Double_t valueFDALIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2873  fhnV0ALambdaInclMCFD[iCentIndex]->Fill(valueFDALIncl);
2874  if(bIsInConeRnd)
2875  {
2876  fhnV0ALambdaBulkMCFD[iCentIndex]->Fill(valueFDALIncl);
2877  }
2878  if(bIsInConeJet)
2879  {
2880  Double_t valueFDALInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2881  fhnV0ALambdaInJetsMCFD[iCentIndex]->Fill(valueFDALInJets);
2882  }
2883  }
2884  if(bV0MCIsALambda && !bV0MCIsPrimaryDist && !bV0MCComesFromAXi) // not primary anti-Lambda
2885  {
2886  fh1V0ALambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2887  }
2888  }
2889  }
2890  //===== End Association of reconstructed V0 candidates with MC particles =====
2891  }
2892  //===== End of V0 loop =====
2893  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "End of V0 loop");
2894 
2895  if(fbCorrelations && iNJetSel)
2896  {
2897  if(bPoolReady)
2898  fh1VtxZME[iCentIndex]->Fill(dZVtxME); // fill z_vtx if event was used for event mixing
2899  pool->UpdatePool(arrayMixedEventAdd); // update the pool with jets from this event
2900  }
2901 
2902  fh1V0CandPerEvent->Fill(iNV0CandTot);
2903  fh1V0CandPerEventCentK0s[iCentIndex]->Fill(iNV0CandK0s);
2904  fh1V0CandPerEventCentLambda[iCentIndex]->Fill(iNV0CandLambda);
2905  fh1V0CandPerEventCentALambda[iCentIndex]->Fill(iNV0CandALambda);
2906 
2907  // Spectra of generated particles
2908  if(fbMCAnalysis)
2909  {
2910  for(Int_t iPartMC = 0; iPartMC < iNTracksMC; iPartMC++)
2911  {
2912  // Get MC particle
2913  AliAODMCParticle* particleMC = (AliAODMCParticle*)arrayMC->At(iPartMC);
2914  if(!particleMC)
2915  continue;
2916 
2917  // Get identity of MC particle
2918  Int_t iPdgCodeParticleMC = particleMC->GetPdgCode();
2919  // Fill Xi spectrum (3322 - Xi0, 3312 - Xi-)
2920 // if ( (iPdgCodeParticleMC==3322) || (iPdgCodeParticleMC==3312) )
2921  if((iPdgCodeParticleMC == 3312) && (TMath::Abs(particleMC->Y()) < 0.5) && IsFromGoodGenerator(iPartMC))
2922  {
2923  fh1V0XiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2924  }
2925  else if((iPdgCodeParticleMC == -3312) && (TMath::Abs(particleMC->Y()) < 0.5) && IsFromGoodGenerator(iPartMC))
2926  {
2927  fh1V0AXiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2928  }
2929  // Skip not interesting particles
2930  if((iPdgCodeParticleMC != iPdgCodeK0s) && (TMath::Abs(iPdgCodeParticleMC) != iPdgCodeLambda))
2931  continue;
2932 
2933  // Check identity of the MC V0 particle
2934  // Is MC V0 particle K0S?
2935  Bool_t bV0MCIsK0s = (iPdgCodeParticleMC == iPdgCodeK0s);
2936  // Is MC V0 particle Lambda?
2937  Bool_t bV0MCIsLambda = (iPdgCodeParticleMC == +iPdgCodeLambda);
2938  // Is MC V0 particle anti-Lambda?
2939  Bool_t bV0MCIsALambda = (iPdgCodeParticleMC == -iPdgCodeLambda);
2940 
2941  Double_t dPtV0Gen = particleMC->Pt();
2942  Double_t dRapV0Gen = particleMC->Y();
2943  Double_t dEtaV0Gen = particleMC->Eta();
2944 
2945  // V0 pseudorapidity cut
2946  if(fdCutEtaV0Max > 0.)
2947  {
2948  if(bPrintCuts) printf("Gen: Applying cut: V0 |eta|: < %g\n", fdCutEtaV0Max);
2949  if((TMath::Abs(dEtaV0Gen) > fdCutEtaV0Max))
2950  continue;
2951  }
2952  // V0 rapidity cut
2953  if(fdCutRapV0Max > 0.)
2954  {
2955  if(bPrintCuts) printf("Gen: Applying cut: V0 |y|: < %g\n", fdCutRapV0Max);
2956  if((TMath::Abs(dRapV0Gen) > fdCutRapV0Max))
2957  continue;
2958  }
2959  /*
2960  // Is MC V0 particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2961  Bool_t bV0MCIsPrimary = particleMC->IsPhysicalPrimary();
2962 
2963  // Get the MC mother particle of the MC V0 particle
2964  Int_t iIndexMotherOfMother = particleMC->GetMother();
2965  AliAODMCParticle* particleMCMotherOfMother = 0;
2966  if (iIndexMotherOfMother >= 0)
2967  particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2968  // Get identity of the MC mother particle of the MC V0 particle if it exists
2969  Int_t iPdgCodeMotherOfMother = 0;
2970  if (particleMCMotherOfMother)
2971  iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2972  // Check if the MC mother particle is a physical primary Sigma
2973  Bool_t bV0MCComesFromSigma = kFALSE;
2974  if ((particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) )
2975  bV0MCComesFromSigma = kTRUE;
2976  // Should the MC V0 particle be considered as primary when it is Lambda?
2977  Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2978  */
2979 
2980  // Get the distance between the production point of the MC V0 particle and the primary vertex
2981  Double_t dx = dPrimVtxMCX - particleMC->Xv();
2982  Double_t dy = dPrimVtxMCY - particleMC->Yv();
2983  Double_t dz = dPrimVtxMCZ - particleMC->Zv();
2984  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2985  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2986 
2987  // Select only primary-like MC V0 particles
2988  if(!bV0MCIsPrimaryDist)
2989  continue;
2990 
2991  // Select only particles from a specific generator
2992  if(!IsFromGoodGenerator(iPartMC))
2993  continue;
2994 
2995  // Check whether the MC V0 particle is in a MC jet
2996  AliAODJet* jetMC = 0;
2997  Bool_t bIsMCV0InJet = kFALSE;
2998  if(iNJetSel)
2999  {
3000  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Searching for gen V0 in %d MC jets", iNJetSel));
3001  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
3002  {
3003  jetMC = (AliAODJet*)arrayJetSel->At(iJet); // load a jet in the list
3004  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Checking if gen V0 in MC jet %d", iJet));
3005  if(IsParticleInCone(particleMC, jetMC, fdDistanceV0JetMax)) // If good jet in event, find out whether V0 is in that jet
3006  {
3007  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("gen V0 found in MC jet %d", iJet));
3008  bIsMCV0InJet = kTRUE;
3009  break;
3010  }
3011  }
3012  }
3013 
3014  // K0s
3015 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
3016  if(bV0MCIsK0s) // well reconstructed candidates
3017  {
3018  fh1V0K0sPtMCGen[iCentIndex]->Fill(dPtV0Gen);
3019  fh2V0K0sEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
3020  if(bIsMCV0InJet)
3021  {
3022  fh2V0K0sInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
3023  Double_t valueEtaKInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
3024  fh3V0K0sInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaKInGen);
3025  }
3026  }
3027  // Lambda
3028 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
3029  if(bV0MCIsLambda) // well reconstructed candidates
3030  {
3031  fh1V0LambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
3032  fh2V0LambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
3033  if(bIsMCV0InJet)
3034  {
3035  fh2V0LambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
3036  Double_t valueEtaLInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
3037  fh3V0LambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaLInGen);
3038  }
3039  }
3040  // anti-Lambda
3041 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
3042  if(bV0MCIsALambda) // well reconstructed candidates
3043  {
3044  fh1V0ALambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
3045  fh2V0ALambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
3046  if(bIsMCV0InJet)
3047  {
3048  fh2V0ALambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
3049  Double_t valueEtaALInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
3050  fh3V0ALambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaALInGen);
3051  }
3052  }
3053  }
3054  }
3055 
3056  arrayJetSel->Delete();
3057  delete arrayJetSel;
3058  arrayJetPerp->Delete();
3059  delete arrayJetPerp;
3060  if(jetRnd)
3061  delete jetRnd;
3062  jetRnd = 0;
3063 
3064  PostData(1, fOutputListStd);
3065  PostData(2, fOutputListQA);
3066  PostData(3, fOutputListCuts);
3067  PostData(4, fOutputListMC);
3068 
3069  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "End");
3070 
3071  return kFALSE; // Must be false to avoid calling PostData from AliAnalysisTaskEmcal. Otherwise, slot 1 is not stored.
3072 }
3073 
3074 void AliAnalysisTaskV0sInJetsEmcal::FillQAHistogramV0(AliAODVertex* vtx, const AliAODv0* vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsCandALambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda, Bool_t IsInPeakALambda)
3075 {
3076  if(!IsCandK0s && !IsCandLambda && !IsCandALambda)
3077  return;
3078 
3079 // Double_t dMassK0s = vZero->MassK0Short();
3080 // Double_t dMassLambda = vZero->MassLambda();
3081 
3082  fh1QAV0Status[iIndexHisto]->Fill(vZero->GetOnFlyStatus());
3083 
3084  AliAODTrack* trackNeg = (AliAODTrack*)vZero->GetDaughter(1); // negative track
3085  AliAODTrack* trackPos = (AliAODTrack*)vZero->GetDaughter(0); // positive track
3086 
3087  Short_t fTotalCharge = 0;
3088  for(Int_t i = 0; i < 2; i++)
3089  {
3090  AliAODTrack* track = (AliAODTrack*)vZero->GetDaughter(i); // track
3091  // Tracks TPC OK
3092  fh1QAV0TPCRefit[iIndexHisto]->Fill(track->IsOn(AliAODTrack::kTPCrefit));
3093  Double_t nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
3094  fh1QAV0TPCRows[iIndexHisto]->Fill(nCrossedRowsTPC);
3095  Int_t findable = track->GetTPCNclsF();
3096  fh1QAV0TPCFindable[iIndexHisto]->Fill(findable);
3097  if(findable != 0)
3098  {
3099  fh1QAV0TPCRowsFind[iIndexHisto]->Fill(nCrossedRowsTPC / findable);
3100  }
3101  // Daughters: pseudo-rapidity cut
3102  fh1QAV0Eta[iIndexHisto]->Fill(track->Eta());
3103 // if((nCrossedRowsTPC > (160. / (250. - 85.) * (255.*TMath::Abs(tan(track->Theta())) - 85.)) + 20.) && (track->Eta() < 0) && (track->Pt() > 0.15))
3104 // if (IsCandK0s)
3105 // {
3106  fh2QAV0EtaRows[iIndexHisto]->Fill(track->Eta(), nCrossedRowsTPC);
3107  fh2QAV0PtRows[iIndexHisto]->Fill(track->Pt(), nCrossedRowsTPC);
3108  fh2QAV0PhiRows[iIndexHisto]->Fill(track->Phi(), nCrossedRowsTPC);
3109  fh2QAV0NClRows[iIndexHisto]->Fill(findable, nCrossedRowsTPC);
3110  fh2QAV0EtaNCl[iIndexHisto]->Fill(track->Eta(), findable);
3111  fh2QAV0PtNCls[iIndexHisto]->Fill(track->Pt(), track->GetTPCNcls());
3112  fh2QAV0PtChi[iIndexHisto]->Fill(track->Pt(), track->Chi2perNDF());
3113 // }
3114 
3115  // Daughters: transverse momentum cut
3116  fh1QAV0Pt[iIndexHisto]->Fill(track->Pt());
3117  fTotalCharge += track->Charge();
3118  }
3119  fh1QAV0Charge[iIndexHisto]->Fill(fTotalCharge);
3120 
3121  // Daughters: Impact parameter of daughters to prim vtx
3122  fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaNegToPrimVertex()));
3123  fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaPosToPrimVertex()));
3124 // fh2CutDCAVtx[iIndexHisto]->Fill(dMassK0s,TMath::Abs(vZero->DcaNegToPrimVertex()));
3125 
3126  // Daughters: DCA
3127  fh1QAV0DCAV0[iIndexHisto]->Fill(vZero->DcaV0Daughters());
3128 // fh2CutDCAV0[iIndexHisto]->Fill(dMassK0s,vZero->DcaV0Daughters());
3129 
3130  // V0: Cosine of the pointing angle
3131  fh1QAV0Cos[iIndexHisto]->Fill(vZero->CosPointingAngle(vtx));
3132 // fh2CutCos[iIndexHisto]->Fill(dMassK0s,vZero->CosPointingAngle(vtx));
3133 
3134  // V0: Fiducial volume
3135  Double_t xyz[3];
3136  vZero->GetSecondaryVtx(xyz);
3137  Double_t r2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
3138  fh1QAV0R[iIndexHisto]->Fill(TMath::Sqrt(r2));
3139 
3140  Double_t dAlpha = vZero->AlphaV0();
3141  Double_t dPtArm = vZero->PtArmV0();
3142 
3143  if(IsCandK0s)
3144  {
3145  if(IsInPeakK0s)
3146  {
3147 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
3148 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
3149  fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
3150  fh2QAV0PtPtK0sPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
3151  fh2ArmPodK0s[iIndexHisto]->Fill(dAlpha, dPtArm);
3152  fh2QAV0PhiPtK0sPeak[iIndexHisto]->Fill(vZero->Phi(), vZero->Pt());
3153  }
3154  fh2QAV0EtaEtaK0s[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
3155  fh2QAV0PhiPhiK0s[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
3156  fh1QAV0RapK0s[iIndexHisto]->Fill(vZero->RapK0Short());
3157  }
3158 
3159  if(IsCandLambda)
3160  {
3161  if(IsInPeakLambda)
3162  {
3163 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
3164 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
3165  fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
3166  fh2QAV0PtPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
3167  fh2ArmPodLambda[iIndexHisto]->Fill(dAlpha, dPtArm);
3168  fh2QAV0PhiPtLambdaPeak[iIndexHisto]->Fill(vZero->Phi(), vZero->Pt());
3169  }
3170  fh2QAV0EtaEtaLambda[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
3171  fh2QAV0PhiPhiLambda[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
3172  fh1QAV0RapLambda[iIndexHisto]->Fill(vZero->RapLambda());
3173  }
3174 
3175  if(IsCandALambda)
3176  {
3177  if(IsInPeakALambda)
3178  {
3179 // fh2QAV0EtaPtALambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
3180 // fh2QAV0EtaPtALambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
3181  fh2QAV0EtaPtALambdaPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
3182  fh2QAV0PtPtALambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
3183  fh2ArmPodALambda[iIndexHisto]->Fill(dAlpha, dPtArm);
3184  fh2QAV0PhiPtALambdaPeak[iIndexHisto]->Fill(vZero->Phi(), vZero->Pt());
3185  }
3186  fh2QAV0EtaEtaALambda[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
3187  fh2QAV0PhiPhiALambda[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
3188  fh1QAV0RapALambda[iIndexHisto]->Fill(vZero->RapLambda());
3189  }
3190 
3191  fh2ArmPod[iIndexHisto]->Fill(dAlpha, dPtArm);
3192 }
3193 
3194 void AliAnalysisTaskV0sInJetsEmcal::FillCandidates(Double_t mK, Double_t mL, Double_t mAL, Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut/*cut index*/, Int_t iCent/*cent index*/)
3195 {
3196  if(isK)
3197  {
3198  fh1V0CounterCentK0s[iCent]->Fill(iCut);
3199  fh1V0InvMassK0sAll[iCut]->Fill(mK);
3200  }
3201  if(isL)
3202  {
3203  fh1V0CounterCentLambda[iCent]->Fill(iCut);
3204  fh1V0InvMassLambdaAll[iCut]->Fill(mL);
3205  }
3206  if(isAL)
3207  {
3208  fh1V0CounterCentALambda[iCent]->Fill(iCut);
3209  fh1V0InvMassALambdaAll[iCut]->Fill(mAL);
3210  }
3211 }
3212 
3213 Bool_t AliAnalysisTaskV0sInJetsEmcal::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
3214 {
3215 // decides whether a particle is inside a jet cone
3216  if(!part1 || !part2)
3217  return kFALSE;
3218 
3219  TVector3 vecMom2(part2->Px(), part2->Py(), part2->Pz());
3220  TVector3 vecMom1(part1->Px(), part1->Py(), part1->Pz());
3221  Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
3222  if(dR < dRMax) // momentum vectors of part1 and part2 are closer than dRMax
3223  return kTRUE;
3224  return kFALSE;
3225 }
3226 
3227 Bool_t AliAnalysisTaskV0sInJetsEmcal::OverlapWithJets(const TClonesArray* array, const AliVParticle* part, Double_t dDistance) const
3228 {
3229 // decides whether a cone overlaps with other jets
3230  if(!part)
3231  {
3232  AliError("No particle!");
3233  return kFALSE;
3234  }
3235  if(!array)
3236  {
3237  AliError("No jet array!");
3238  return kFALSE;
3239  }
3240  Int_t iNJets = array->GetEntriesFast();
3241  if(iNJets <= 0)
3242  {
3243  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Warning: No jets");
3244  return kFALSE;
3245  }
3246  AliVParticle* jet = 0;
3247  for(Int_t iJet = 0; iJet < iNJets; iJet++)
3248  {
3249  jet = (AliVParticle*)array->At(iJet);
3250  if(!jet)
3251  {
3252  AliError(Form("Failed to load jet %d/%d!", iJet, iNJets));
3253  continue;
3254  }
3255  if(IsParticleInCone(part, jet, dDistance))
3256  return kTRUE;
3257  }
3258  return kFALSE;
3259 }
3260 
3261 AliAODJet* AliAnalysisTaskV0sInJetsEmcal::GetRandomCone(const TClonesArray* array, Double_t dEtaConeMax, Double_t dDistance) const
3262 {
3263 // generate a random cone which does not overlap with selected jets
3264  if(fDebug > 3) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Generating random cone...");
3265  TLorentzVector vecCone;
3266  AliAODJet* part = 0;
3267  Double_t dEta, dPhi;
3268  Int_t iNTrialsMax = 10;
3269  Bool_t bStatus = kFALSE;
3270  for(Int_t iTry = 0; iTry < iNTrialsMax; iTry++)
3271  {
3272  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Try %d", iTry));
3273  dEta = dEtaConeMax * (2 * fRandom->Rndm() - 1.); // random eta in [-dEtaConeMax,+dEtaConeMax]
3274  dPhi = TMath::TwoPi() * fRandom->Rndm(); // random phi in [0,2*Pi]
3275  vecCone.SetPtEtaPhiM(1., dEta, dPhi, 0.);
3276  part = new AliAODJet(vecCone);
3277  if(!OverlapWithJets(array, part, dDistance))
3278  {
3279  bStatus = kTRUE;
3280  if(fDebug > 1) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Random cone successfully generated");
3281  break;
3282  }
3283  else
3284  delete part;
3285  }
3286  if(!bStatus)
3287  {
3288  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Failed to find a random cone");
3289  part = 0;
3290  }
3291  return part;
3292 }
3293 
3295 {
3296 // sort kt clusters by pT/area and return the middle one, based on code in AliAnalysisTaskJetChem
3297  if(!cont)
3298  {
3299  AliError("No jet container!");
3300  return NULL;
3301  }
3302  Int_t iNClTot = cont->GetNJets(); // number of all clusters in the array
3303  Int_t iNCl = 0; // number of accepted clusters
3304 
3305  // get list of densities
3306  std::vector<std::vector<Double_t> > vecListClusters; // vector that contains pairs [ index, density ]
3307  if(fDebug > 3) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Loop over %d clusters.", iNClTot));
3308  for(Int_t ij = 0; ij < iNClTot; ij++)
3309  {
3310  AliEmcalJet* clusterBg = (AliEmcalJet*)(cont->GetAcceptJet(ij));
3311  if(!clusterBg)
3312  continue;
3313  if(fDebug > 4) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Cluster %d/%d used as accepted cluster %d.", ij, iNClTot, int(vecListClusters.size())));
3314  Double_t dPtBg = clusterBg->Pt();
3315  Double_t dAreaBg = clusterBg->Area();
3316  Double_t dDensityBg = 0;
3317  if(dAreaBg > 0)
3318  dDensityBg = dPtBg / dAreaBg;
3319  std::vector<Double_t> vecCluster;
3320  vecCluster.push_back(ij);
3321  vecCluster.push_back(dDensityBg);
3322  vecListClusters.push_back(vecCluster);
3323  }
3324  iNCl = vecListClusters.size();
3325  if(iNCl < 3) // need at least 3 clusters (skipping 2 highest)
3326  {
3327  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Warning: Too little clusters");
3328  return NULL;
3329  }
3330 
3331 // printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Original lists:");
3332 // for(Int_t i = 0; i < iNCl; i++)
3333 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3334 
3335  // sort list of indeces by density in descending order
3336  std::sort(vecListClusters.begin(), vecListClusters.end(), CompareClusters);
3337 
3338 // printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Sorted lists:");
3339 // for(Int_t i = 0; i < iNCl; i++)
3340 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3341 
3342  // get median cluster with median density
3343  AliEmcalJet* clusterMed = 0;
3344  Int_t iIndex = 0; // index of the median cluster in the sorted list
3345  Int_t iIndexMed = 0; // index of the median cluster in the original array
3346  if(TMath::Odd(iNCl)) // odd number of clusters
3347  {
3348  iIndex = (Int_t)(0.5 * (iNCl + 1)); // = (n - skip + 1)/2 + 1, skip = 2
3349 // printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Odd, median index = %d/%d", iIndex, iNCl));
3350  }
3351  else // even number: picking randomly one of the two closest to median
3352  {
3353  Int_t iIndex1 = (Int_t)(0.5 * iNCl); // = (n - skip)/2 + 1, skip = 2
3354  Int_t iIndex2 = (Int_t)(0.5 * iNCl + 1); // = (n - skip)/2 + 1 + 1, skip = 2
3355  iIndex = ((fRandom->Rndm() > 0.5) ? iIndex1 : iIndex2);
3356 // printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Even, median index = %d or %d -> %d/%d", iIndex1, iIndex2, iIndex, iNCl));
3357  }
3358  iIndexMed = Int_t((vecListClusters[iIndex])[0]);
3359 
3360  if(fDebug > 1) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Getting median cluster %d/%d ok, rho = %g", iIndexMed, iNClTot, (vecListClusters[iIndex])[1]));
3361  clusterMed = (AliEmcalJet*)(cont->GetAcceptJet(iIndexMed));
3362 
3363  if(TMath::Abs(clusterMed->Eta()) > dEtaConeMax)
3364  {
3365  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Chosen median cluster is out of window |eta| < %g, eta = %g", dEtaConeMax, clusterMed->Eta()));
3366  return NULL;
3367  }
3368 
3369  return clusterMed;
3370 }
3371 
3373 {
3374 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
3375  Double_t dEpsilon = 1e-2;
3376  Double_t dR = dRadius;
3377  Double_t dD = dDistance;
3378  if(TMath::Abs(dR) < dEpsilon)
3379  {
3380  AliError(Form("Too small radius: %g < %g!", dR, dEpsilon));
3381  return 0.;
3382  }
3383  if(dD >= dR)
3384  return 0.;
3385  if(dD <= -dR)
3386  return TMath::Pi() * dR * dR;
3387  return dR * dR * TMath::ACos(dD / dR) - dD * TMath::Sqrt(dR * dR - dD * dD);
3388 }
3389 
3390 Double_t AliAnalysisTaskV0sInJetsEmcal::GetD(const AliVParticle* part1, const AliVParticle* part2) const
3391 {
3392 // return D between two particles
3393  if(!part1 || !part2)
3394  return -1;
3395 
3396  TVector3 vecMom2(part2->Px(), part2->Py(), part2->Pz());
3397  TVector3 vecMom1(part1->Px(), part1->Py(), part1->Pz());
3398  Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
3399  return dR;
3400 }
3401 
3403 {
3404 // event selection
3405  if(!fbIsPbPb)
3406  {
3407  if(fAOD->IsPileupFromSPD())
3408  {
3409  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "SPD pile-up");
3410  return kFALSE;
3411  }
3412  fh1EventCounterCut->Fill(2); // not pile-up from SPD
3413  }
3414  AliAODVertex* vertex = fAOD->GetPrimaryVertex();
3415  if(!vertex)
3416  {
3417  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "No vertex");
3418  return kFALSE;
3419  }
3420  Int_t iNContribMin = 3;
3421  if(vertex->GetNContributors() < iNContribMin)
3422  {
3423  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Not enough contributors, %d", vertex->GetNContributors()));
3424  return kFALSE;
3425  }
3426  fh1EventCounterCut->Fill(3); // enough contributors
3427  TString vtxTitle(vertex->GetTitle());
3428  if(vtxTitle.Contains("TPCVertex"))
3429  {
3430  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "TPC vertex");
3431  return kFALSE;
3432  }
3433  fh1EventCounterCut->Fill(4); // not TPC vertex only
3434  Double_t zVertex = vertex->GetZ();
3435  if(TMath::Abs(zVertex) > dVtxZCut)
3436  {
3437  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Cut on z, %g", zVertex));
3438  return kFALSE;
3439  }
3440  fh1EventCounterCut->Fill(5); // PV z coordinate within range
3441  if(dDeltaZMax > 0.) // cut on |delta z| between SPD vertex and nominal primary vertex
3442  {
3443  AliAODVertex* vertexSPD = fAOD->GetPrimaryVertexSPD();
3444  if(!vertexSPD)
3445  {
3446  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "No SPD vertex");
3447  return kFALSE;
3448  }
3449  Double_t zVertexSPD = vertexSPD->GetZ();
3450  if(TMath::Abs(zVertex - zVertexSPD) > dDeltaZMax)
3451  {
3452  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Cut on Delta z = %g - %g = %g", zVertex, zVertexSPD, zVertex - zVertexSPD));
3453  return kFALSE;
3454  }
3455  fh1EventCounterCut->Fill(6); // delta z within range
3456  }
3457  Double_t xVertex = vertex->GetX();
3458  Double_t yVertex = vertex->GetY();
3459  Double_t radiusSq = yVertex * yVertex + xVertex * xVertex;
3460  if(radiusSq > dVtxR2Cut)
3461  {
3462  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Cut on r, %g", radiusSq));
3463  return kFALSE;
3464  }
3465  fh1EventCounterCut->Fill(7); // radius within range
3466  if(fbIsPbPb)
3467  {
3468  fdCentrality = ((AliVAODHeader*)fAOD->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
3469  if(fdCentrality < 0)
3470  {
3471  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Negative centrality");
3472  return kFALSE;
3473  }
3474  if((dCentCutUp < 0) || (dCentCutLo < 0) || (dCentCutUp > 100) || (dCentCutLo > 100) || (dCentCutLo > dCentCutUp))
3475  {
3476  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, "Wrong centrality limits");
3477  return kFALSE;
3478  }
3479  if((fdCentrality < dCentCutLo) || (fdCentrality > dCentCutUp))
3480  {
3481  if(fDebug > 0) printf("%s %s::%s: %s\n", GetName(), ClassName(), __func__, Form("Centrality cut, %g", fdCentrality));
3482  return kFALSE;
3483  }
3484  }
3485  return kTRUE;
3486 }
3487 
3489 {
3490 // returns index of the centrality bin corresponding to the provided value of centrality
3492  return -1;
3493  for(Int_t i = 0; i < fgkiNBinsCent; i++)
3494  {
3495  if(centrality <= fgkiCentBinRanges[i])
3496  return i;
3497  }
3498  return -1;
3499 }
3500 
3502 {
3503 // returns the upper edge of the centrality bin corresponding to the provided value of index
3504  if(index < 0 || index >= fgkiNBinsCent)
3505  return -1;
3506  return fgkiCentBinRanges[index];
3507 }
3508 
3510 {
3511 // get string with centrality range for given bin
3512  TString lowerEdge = ((index == 0) ? "0" : Form("%d", GetCentralityBinEdge(index - 1)));
3513  TString upperEdge = Form("%d", GetCentralityBinEdge(index));
3514  return Form("%s-%s %%", lowerEdge.Data(), upperEdge.Data());
3515 }
3516 
3518 {
3519 // estimation of the sigma of the invariant-mass peak as a function of pT and particle type
3520  switch(particle)
3521  {
3522  case 0: // K0S
3523  return 0.0044 + 0.0004 * (pt - 1.);
3524  break;
3525  case 1: // Lambda
3526  return 0.0023 + 0.00034 * (pt - 1.);
3527  break;
3528  default:
3529  return 0;
3530  break;
3531  }
3532 }
3533 
3534 bool AliAnalysisTaskV0sInJetsEmcal::CompareClusters(const std::vector<Double_t> cluster1, const std::vector<Double_t> cluster2)
3535 {
3536  return (cluster1[1] > cluster2[1]);
3537 }
3538 
3540 {
3541  if(!fEventMC)
3542  {
3543  AliError("No MC event!");
3544  return kFALSE;
3545  }
3546  if(fsGeneratorName.Length())
3547  {
3548  TString sGenName = "";
3549  Bool_t bCocktailOK = fEventMC->GetCocktailGenerator(index, sGenName);
3550  if(!bCocktailOK || !sGenName.Contains(fsGeneratorName.Data()))
3551  return kFALSE;
3552  }
3553  return kTRUE;
3554 }
THnSparse * fhnV0InJetK0s[fgkiNBinsCent]
V0 inclusive, in a centrality bin, m_V0; pt_V0; eta_V0.
THnSparse * fhnV0InRndALambda[fgkiNBinsCent]
THnSparse * fh3V0LambdaInJetPtMassMCRec[fgkiNBinsCent]
TH1D * fh1DistanceV0JetsALambda[fgkiNBinsCent]
distance in eta-phi between V0 and the closest jet
TH2D * fh2QAV0PtChi[fgkiNQAIndeces]
pt vs TPC clusters
Double_t Area() const
Definition: AliEmcalJet.h:123
void SetParticlePtCut(Double_t cut)
AliEmcalJet * GetMedianCluster(AliJetContainer *cont, Double_t dEtaConeMax) const
TH2D * fh2QAV0PhiPtLambdaPeak[fgkiNQAIndeces]
K0S candidate in peak: azimuth; pt.
Double_t GetRhoVal() const
TH1D * fh1V0InvMassALambdaAll[fgkiNCategV0]
number of ALambda candidates after various cuts
double Double_t
Definition: External.C:58
TList * fOutputListMC
Output list for checking cuts.
THnSparse * fhnPtDaughterK0s[fgkiNBinsCent]
V0 in no-jet events, in a centrality bin, m_V0; pt_V0; eta_V0.
TH1D * fh1QAV0TPCRows[fgkiNQAIndeces]
TPC refit on vs off.
TH2D * fh2QAV0NClRows[fgkiNQAIndeces]
azimuth vs TPC rows
THnSparse * fh3V0K0sInJetPtMassMCRec[fgkiNBinsCent]
pt spectrum of generated K0s in jet
TH2D * fh2PtJetPtTrackLeading[fgkiNBinsCent]
jet phi
AliJetContainer * GetJetContainer(Int_t i=0) const
THnSparse * fhnV0ALambdaBulkMCFD[fgkiNBinsCent]
THnSparse * fh3V0K0sInJetEtaPtMCGen[fgkiNBinsCent]
mass-pt spectrum of successfully reconstructed K0s in jet
Double_t Eta() const
Definition: AliEmcalJet.h:114
THnSparse * fh3CCMassCorrelKNotL
mass correlation of candidates
TH1D * fh1V0K0sPtMCGen[fgkiNBinsCent]
pt jet vs angle V0-jet, in centrality bins
Double_t Phi() const
Definition: AliEmcalJet.h:110
void SetPtBiasJetTrack(Float_t b)
Container with name, TClonesArray and cuts for particles.
THnSparse * fhnV0InMedLambda[fgkiNBinsCent]
TH2D * fh2V0LambdaMCResolMPt[fgkiNBinsCent]
V0 in jets, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
THnSparse * fh3CCMassCorrelBoth
Lambda candidates in K0s peak.
TH1D * fh1QACTau2D[fgkiNQAIndeces]
radial distance between prim vtx and decay vertex
centrality
char Char_t
Definition: External.C:18
TH1D * fh1EventCounterCutCent[fgkiNBinsCent]
number of events for different selection steps
Double_t GetD(const AliVParticle *part1, const AliVParticle *part2) const
void ExecOnce()
Perform steps needed to initialize the analysis.
TH1D * fh1V0InvMassK0sCent[fgkiNBinsCent]
number of K0s candidates per event, in centrality bins
THnSparse * fhnV0OutJetLambda[fgkiNBinsCent]
THnSparse * fhnV0InPerpK0s[fgkiNBinsCent]
V0 in jet cones, in a centrality bin, m_V0; pt_V0; eta_V0; pt_jet.
TH2D * fh2VtxXY[fgkiNBinsCent]
z coordinate of the primary vertex for events used in mixed events
TH1D * fh1V0CandPerEventCentK0s[fgkiNBinsCent]
Armenteros-Podolanski.
TH2D * fh2V0K0sPtMassMCRec[fgkiNBinsCent]
pt spectrum of all generated K0s in event
TH1D * fh1EtaJet[fgkiNBinsCent]
pt spectra of jets for normalisation of in-jet V0 spectra
TH2D * fh2V0PtJetAngleALambda[fgkiNBinsCent]
pt correlations, in a centrality bin, pt_neg-daughter;pt_pos-daughter;pt_V0;pt_jet;pt_leading-track ...
TH1D * fh1V0CounterCentK0s[fgkiNBinsCent]
pseudorapidity vs clusters
void SetPercAreaCut(Float_t p)
THnSparse * fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
THnSparse * fh3CCMassCorrelLNotK
mass correlation of candidates
THnSparse * fhnV0ALambdaInJetsMCFD[fgkiNBinsCent]
TH2D * fh2V0LambdaInJetPtMCGen[fgkiNBinsCent]
V0 inclusive, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
TH1D * fh1VtxZME[fgkiNBinsCent]
z coordinate of the primary vertex
TList * fOutputListQA
Output list for standard analysis results.
THnSparse * fh3V0ALambdaEtaPtMassMCRec[fgkiNBinsCent]
THnSparse * fhnV0InJetLambda[fgkiNBinsCent]
TH2D * fh2CCK0s
Armenteros-Podolanski.
THnSparse * fhnV0NoJetK0s[fgkiNBinsCent]
V0 outside jet cones, in a centrality bin, m_V0; pt_V0; eta_V0.
TH2D * fh2QAV0PhiRows[fgkiNQAIndeces]
pt vs TPC rows
TH1D * fh1QAV0TPCRefit[fgkiNQAIndeces]
online vs offline reconstructed V0 candidates
THnSparse * fhnV0LambdaInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH2D * fh2V0K0sEtaPtMCGen[fgkiNBinsCent]
pt spectrum of false reconstructed K0s in event
AliMCEvent * fEventMC
Output AOD event.
THnSparse * fhnV0InRndK0s[fgkiNBinsCent]
V0 in perpendicular cones, in a centrality bin, m_V0; pt_V0; eta_V0; pt_jet.
TH2D * fh2EtaPhiRndCone[fgkiNBinsCent]
number of generated random cones in centrality bins
TH1D * fh1V0InvMassK0sAll[fgkiNCategV0]
number of K0s candidates after various cuts
THnSparse * fhnV0K0sInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
mass-eta-pt spectrum of successfully reconstructed K0s in jet
static Double_t fgkiCentMixBinRanges[fgkiNBinsCentMix+1]
THnSparse * fhnV0InPerpLambda[fgkiNBinsCent]
AliAODJet * GetRandomCone(const TClonesArray *array, Double_t dEtaConeMax, Double_t dDistance) const
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH1D * fh1AreaExcluded
median-cluster cone eta-phi
TH1D * fh1QAV0RapK0s[fgkiNQAIndeces]
daughters azimuth vs azimuth
AliEmcalJet * GetLeadingJet(const char *opt="")
Bool_t fbIsPbPb
Output list for MC related results.
static const Int_t fgkiNCategV0
distance in eta-phi between V0 and the closest jet
AliEventPoolManager * fPoolMgr
random-number generator
ClassImp(AliAnalysisTaskV0sInJetsEmcal) const Int_t AliAnalysisTaskV0sInJetsEmcal
TH1D * fh1VtxZ[fgkiNBinsCent]
number of tracks vs centrality
THnSparse * fhnV0CorrelSELambda[fgkiNBinsCent]
pt correlations, in a centrality bin, pt_neg-daughter;pt_pos-daughter;pt_V0;pt_jet;pt_leading-track ...
TH1D * fh1V0K0sPtMCRecFalse[fgkiNBinsCent]
pt-mass spectrum of successfully reconstructed K0s in event
int Int_t
Definition: External.C:63
TH1D * fh1DistanceV0JetsK0s[fgkiNBinsCent]
distance in eta-phi between jets within events
TH1D * fh1V0CandPerEvent
xy coordinates of the primary vertex
TH2D * fh2ArmPodK0s[fgkiNQAIndeces]
daughters pt vs pt, in mass peak
TH1D * fh1EventCent
number of events for different selection steps and different centralities
unsigned int UInt_t
Definition: External.C:33
THnSparse * fhnV0InclusiveLambda[fgkiNBinsCent]
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
static const Int_t fgkiCentBinRanges[fgkiNBinsCent]
THnSparse * fh3V0ALambdaInJetPtMassMCRec[fgkiNBinsCent]
THnSparse * fhnV0LambdaInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH2D * fh2ArmPod[fgkiNQAIndeces]
lifetime calculated in xyz
TH2D * fh2V0K0sMCResolMPt[fgkiNBinsCent]
V0 in jets, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
TH1D * fh1DistanceV0JetsLambda[fgkiNBinsCent]
distance in eta-phi between V0 and the closest jet
TH1D * fh1QAV0DCAVtx[fgkiNQAIndeces]
charge
TH1D * fh1NJetPerEvent[fgkiNBinsCent]
pt of trigger track
void SetFilterHybridTracks(Bool_t f)
THnSparse * fhnV0InRndLambda[fgkiNBinsCent]
Int_t GetNJets() const
THnSparse * fh3V0LambdaInJetEtaPtMCGen[fgkiNBinsCent]
TH2D * fh2EtaPtJet[fgkiNBinsCent]
jet eta
Definition: External.C:228
Definition: External.C:212
TH1D * fh1QAV0Charge[fgkiNQAIndeces]
pt daughter
THnSparse * fh3V0LambdaEtaPtMassMCRec[fgkiNBinsCent]
TH1D * fh1PhiJet[fgkiNBinsCent]
jet eta-pT
TLorentzVector SubtractRhoVect(Double_t rho, Bool_t save=kFALSE)
THnSparse * fh4V0LambdaInJetEtaPtMassMCRec[fgkiNBinsCent]
void FillQAHistogramV0(AliAODVertex *vtx, const AliAODv0 *vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsCandALambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda, Bool_t IsInPeakALambda)
Double_t PtSub() const
Definition: AliEmcalJet.h:151
THnSparse * fhnV0NoJetALambda[fgkiNBinsCent]
TH1D * fh1DistanceJets[fgkiNBinsCent]
area of excluded cones for outside-cones V0s
TH1D * fh1QAV0DCAV0[fgkiNQAIndeces]
DCA of daughters to prim vtx.
THnSparse * fh4V0K0sInJetEtaPtMassMCRec[fgkiNBinsCent]
eta-pt spectrum of generated K0s in jet
AliJetContainer * fJetsBgCont
Signal Jets.
TH1D * fh1EventCent2Jets
number of events for different centralities
TH1D * fh1NRndConeCent
number of jets per event
TList * fOutputListCuts
Output list for quality assurance.
TH2D * fh2V0ALambdaMCResolMPt[fgkiNBinsCent]
V0 in jets, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
THnSparse * fhnV0LambdaBulkMCFD[fgkiNBinsCent]
TH2D * fh2V0K0sMCPtGenPtRec[fgkiNBinsCent]
K0s mass resolution vs pt.
THnSparse * fhnV0InPerpALambda[fgkiNBinsCent]
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
TH1D * fh1V0CounterCentLambda[fgkiNBinsCent]
K0s generated pt vs reconstructed pt.
TH1D * fh1EventCent2NoJets
number of events for different centralities
THnSparse * fhnV0InMedALambda[fgkiNBinsCent]
THnSparse * fhnV0OutJetALambda[fgkiNBinsCent]
THnSparse * fhnV0CorrelMEK0s[fgkiNBinsCent]
V0-jet phi,eta correlations in same events, in a centrality bin, m_V0; pt_V0; eta_V0; pt_jet; delta-p...
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TH2D * fh2V0ALambdaInJetPtMCGen[fgkiNBinsCent]
V0 inclusive, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
THnSparse * fhnV0InclusiveALambda[fgkiNBinsCent]
short Short_t
Definition: External.C:23
THnSparse * fhnV0InclusiveK0s[fgkiNBinsCent]
V0 invariant mass, in centrality bins.
TH2D * fh2QAV0EtaPtK0sPeak[fgkiNQAIndeces]
V0 invariant mass for each selection steps.
Double_t Pt() const
Definition: AliEmcalJet.h:102
TH1D * fh1QAV0TPCRowsFind[fgkiNQAIndeces]
pt vs Chi2/ndf
TH1D * fh1PtJet[fgkiNBinsCent]
number of V0 cand per event
TH2D * fh2PtJetPtTrigger[fgkiNBinsCent]
pt_jet; pt of leading jet track
THnSparse * fhnV0CorrelSEK0s[fgkiNBinsCent]
pt correlations, in a centrality bin, pt_neg-daughter;pt_pos-daughter;pt_V0;pt_jet;pt_leading-track ...
TH2D * fh2QAV0EtaNCl[fgkiNQAIndeces]
clusters vs TPC rows
TH1D * fh1V0InvMassLambdaAll[fgkiNCategV0]
number of Lambda candidates after various cuts
TH1D * fh1PtTrigger[fgkiNBinsCent]
pt_jet; pt of trigger track
static Double_t fgkiZVtxMixBinRanges[fgkiNBinsZVtxMix+1]
void SetParticleEtaLimits(Double_t min, Double_t max)
Double_t MassPeakSigmaOld(Double_t pt, Int_t particle)
Float_t GetJetRadius() const
TH1D * fh1EventCent2
number of events for different centralities
Int_t GetCentralityBinIndex(Double_t centrality)
THnSparse * fhnV0K0sInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
eta-pt-mass spectrum of successfully reconstructed K0s in event
TH2D * fh2QAV0PtRows[fgkiNQAIndeces]
pseudorapidity vs TPC rows
TH1D * fh1QACTau3D[fgkiNQAIndeces]
lifetime calculated in xy
THnSparse * fh3V0ALambdaInJetEtaPtMCGen[fgkiNBinsCent]
AliTrackContainer * fTracksCont
Background Jets.
AliEmcalJet * GetAcceptJet(Int_t i) const
THnSparse * fhnV0NoJetLambda[fgkiNBinsCent]
TH2D * fh2QAV0PtPtK0sPeak[fgkiNQAIndeces]
V0 rapidity.
TH2D * fh2V0PtJetAngleK0s[fgkiNBinsCent]
V0-jet phi,eta correlations in mixed events, in a centrality bin, m_V0; pt_V0; eta_V0; pt_jet; delta-...
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH2D * fh2EventCentTracks
number of events for different centralities
Bool_t IsSelectedForJets(AliAODEvent *fAOD, Double_t dVtxZCut, Double_t dVtxR2Cut, Double_t dCentCutLo, Double_t dCentCutUp, Double_t dDeltaZMax=-1)
TH2D * fh2QAV0EtaRows[fgkiNQAIndeces]
pseudorapidity
TH2D * fh2EtaPhiMedCone[fgkiNBinsCent]
number of found median-cluster cones in centrality bins
THnSparse * fhnV0InMedK0s[fgkiNBinsCent]
V0 in random cones, in a centrality bin, m_V0; pt_V0; eta_V0.
THnSparse * fhnV0CorrelMELambda[fgkiNBinsCent]
THnSparse * fhnV0ALambdaInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH2D * fh2CCLambda
K0s candidates in Lambda peak.
void FillCandidates(Double_t mK, Double_t mL, Double_t mAL, Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut, Int_t iCent)
void UserCreateOutputObjects()
Main initialization function on the worker.
TH2D * fh2QAV0EtaEtaK0s[fgkiNQAIndeces]
daughters pseudorapidity vs V0 pt, in mass peak
Bool_t FillHistograms()
Function filling histograms.
Bool_t OverlapWithJets(const TClonesArray *array, const AliVParticle *cone, Double_t dDistance) const
virtual AliVParticle * GetNextAcceptParticle()
THnSparse * fhnV0ALambdaInclMCFD[fgkiNBinsCent]
bool Bool_t
Definition: External.C:53
void SetJetEtaLimits(Float_t min, Float_t max)
THnSparse * fhnV0InJetALambda[fgkiNBinsCent]
THnSparse * fh4V0ALambdaInJetEtaPtMassMCRec[fgkiNBinsCent]
TH2D * fh2QAV0PhiPhiK0s[fgkiNQAIndeces]
daughters pseudorapidity vs pseudorapidity
TH1D * fh1QAV0Cos[fgkiNQAIndeces]
DCA between daughters.
THnSparse * fhnPtDaughterALambda[fgkiNBinsCent]
Double_t AreaCircSegment(Double_t dRadius, Double_t dDistance) const
TH2D * fh2QAV0PtNCls[fgkiNQAIndeces]
findable clusters
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Bool_t IsParticleInCone(const AliVParticle *part1, const AliVParticle *part2, Double_t dRMax) const
THnSparse * fhnV0LambdaInJetsMCFD[fgkiNBinsCent]
TH2D * fh2QAV0PhiPtALambdaPeak[fgkiNQAIndeces]
Lambda candidate in peak: azimuth; pt.
TH1D * fh1QAV0Pt[fgkiNQAIndeces]
anti-Lambda candidate in peak: azimuth; pt
static bool CompareClusters(const std::vector< Double_t > cluster1, const std::vector< Double_t > cluster2)
TH1D * fh1QAV0R[fgkiNQAIndeces]
cosine of pointing angle (CPA)
THnSparse * fh3V0K0sEtaPtMassMCRec[fgkiNBinsCent]
eta-pt spectrum of all generated K0s in event
TH2D * fh2V0K0sInJetPtMCGen[fgkiNBinsCent]
V0 inclusive, reconstructed: charge_daughter; eta_daughter; pt_daughter; eta_V0; pt_V0; pt_jet...
THnSparse * fhnPtDaughterLambda[fgkiNBinsCent]
THnSparse * fhnV0LambdaInclMCFD[fgkiNBinsCent]
TH1D * fh1QAV0Eta[fgkiNQAIndeces]
ratio rows/clusters
THnSparse * fhnV0OutJetK0s[fgkiNBinsCent]
V0 in medium cones, in a centrality bin, m_V0; pt_V0; eta_V0.
TH1D * fh1QAV0TPCFindable[fgkiNQAIndeces]
crossed TPC pad rows