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