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