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