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