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