AliPhysics  5eaf189 (5eaf189)
AliAnalysisTaskV0sInJets.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 // task for analysis of V0s (K0S, (anti-)Lambda) in charged jets
17 // Author: Vit Kucera (vit.kucera@cern.ch)
18 
19 #include "TChain.h"
20 #include "TTree.h"
21 #include "TH1D.h"
22 #include "TH2D.h"
23 #include "THnSparse.h"
24 #include "TCanvas.h"
25 
26 #include "AliAnalysisTask.h"
27 #include "AliAnalysisManager.h"
28 
29 #include "AliESDEvent.h"
30 #include "AliAODEvent.h"
31 #include "AliAODTrack.h"
32 #include <TDatabasePDG.h>
33 #include <TPDGCode.h>
34 #include "AliPIDResponse.h"
35 #include "AliInputEventHandler.h"
36 #include "AliAODMCHeader.h"
37 #include "AliAODMCParticle.h"
38 #include "TClonesArray.h"
39 //#include "AliEventInfoObject.cxx"
40 //#include "AliV0Object.cxx"
41 //#include "AliJetObject.cxx"
42 #include "TRandom3.h"
43 
45 
47 
48 // upper edges of centrality bins
49 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10, 30, 50, 80}; // Alice Zimmermann
50 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10, 20, 40, 60, 80}; // Vit Kucera, initial binning
51 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {5, 10, 20, 40, 60, 80}; // Iouri Belikov, LF analysis
52 const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10}; // only central
53 
54 // axis: pT of V0
58 // axis: pT of jets
60 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJet = sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet) / sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet[0]) - 1;
61 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[AliAnalysisTaskV0sInJets::fgkiNBinsPtJet] - (AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[0]) / 5.); // bin width 5 GeV/c
62 // axis: K0S invariant mass
66 // axis: Lambda invariant mass
70 
71 
72 // Default constructor
75  fAODIn(0),
76  fAODOut(0),
77  fOutputListStd(0),
78  fOutputListQA(0),
79  fOutputListCuts(0),
80  fOutputListMC(0),
81 // ftreeOut(0),
82 
83  fiAODAnalysis(1),
84  fbIsPbPb(1),
85 
86  fdCutDCAToPrimVtxMin(0.1),
87  fdCutDCADaughtersMax(1.),
88  fdCutNSigmadEdxMax(3),
89  fdCutCPAMin(0.998),
90  fdCutNTauMax(5),
91 
92  fsJetBranchName(0),
93  fsJetBgBranchName(0),
94  fdCutPtJetMin(0),
95  fdCutPtTrackMin(5),
96  fdRadiusJet(0.4),
97  fdRadiusJetBg(0.4),
98  fbJetSelection(0),
99  fbMCAnalysis(0),
100 // fbTreeOutput(0),
101  fRandom(0),
102 
103  fdCutVertexZ(10),
104  fdCutVertexR2(1),
105  fdCutCentLow(0),
106  fdCutCentHigh(80),
107 
108 /*
109  fBranchV0Rec(0),
110  fBranchV0Gen(0),
111  fBranchJet(0),
112  fEventInfo(0),
113 */
114  fdCentrality(0),
115  fh1EventCounterCut(0),
116  fh1EventCent(0),
117  fh1EventCent2(0),
118  fh1EventCent2Jets(0),
119  fh1EventCent2NoJets(0),
120  fh2EventCentTracks(0),
121  fh1V0CandPerEvent(0),
122  fh1NRndConeCent(0),
123  fh1NMedConeCent(0),
124  fh1AreaExcluded(0),
125 
126  fh2CCK0s(0),
127  fh2CCLambda(0),
128  fh3CCMassCorrelBoth(0),
129  fh3CCMassCorrelKNotL(0),
130  fh3CCMassCorrelLNotK(0)
131 {
132  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
133  {
134  fh1QAV0Status[i] = 0;
135  fh1QAV0TPCRefit[i] = 0;
136  fh1QAV0TPCRows[i] = 0;
137  fh1QAV0TPCFindable[i] = 0;
138  fh1QAV0TPCRowsFind[i] = 0;
139  fh1QAV0Eta[i] = 0;
140  fh2QAV0EtaRows[i] = 0;
141  fh2QAV0PtRows[i] = 0;
142  fh2QAV0PhiRows[i] = 0;
143  fh2QAV0NClRows[i] = 0;
144  fh2QAV0EtaNCl[i] = 0;
145 
146  fh2QAV0EtaPtK0sPeak[i] = 0;
147  fh2QAV0EtaEtaK0s[i] = 0;
148  fh2QAV0PhiPhiK0s[i] = 0;
149  fh1QAV0RapK0s[i] = 0;
150  fh2QAV0PtPtK0sPeak[i] = 0;
151  fh2ArmPodK0s[i] = 0;
152 
153  fh2QAV0EtaPtLambdaPeak[i] = 0;
154  fh2QAV0EtaEtaLambda[i] = 0;
155  fh2QAV0PhiPhiLambda[i] = 0;
156  fh1QAV0RapLambda[i] = 0;
157  fh2QAV0PtPtLambdaPeak[i] = 0;
158  fh2ArmPodLambda[i] = 0;
159 
161  fh2QAV0EtaEtaALambda[i] = 0;
162  fh2QAV0PhiPhiALambda[i] = 0;
163  fh1QAV0RapALambda[i] = 0;
164  fh2QAV0PtPtALambdaPeak[i] = 0;
165  fh2ArmPodALambda[i] = 0;
166 
167  fh1QAV0Pt[i] = 0;
168  fh1QAV0Charge[i] = 0;
169  fh1QAV0DCAVtx[i] = 0;
170  fh1QAV0DCAV0[i] = 0;
171  fh1QAV0Cos[i] = 0;
172  fh1QAV0R[i] = 0;
173  fh1QACTau2D[i] = 0;
174  fh1QACTau3D[i] = 0;
175 
176  fh2ArmPod[i] = 0;
177 
178  /*
179  fh2CutTPCRowsK0s[i] = 0;
180  fh2CutTPCRowsLambda[i] = 0;
181  fh2CutPtPosK0s[i] = 0;
182  fh2CutPtNegK0s[i] = 0;
183  fh2CutPtPosLambda[i] = 0;
184  fh2CutPtNegLambda[i] = 0;
185  fh2CutDCAVtx[i] = 0;
186  fh2CutDCAV0[i] = 0;
187  fh2CutCos[i] = 0;
188  fh2CutR[i] = 0;
189  fh2CutEtaK0s[i] = 0;
190  fh2CutEtaLambda[i] = 0;
191  fh2CutRapK0s[i] = 0;
192  fh2CutRapLambda[i] = 0;
193  fh2CutCTauK0s[i] = 0;
194  fh2CutCTauLambda[i] = 0;
195  fh2CutPIDPosK0s[i] = 0;
196  fh2CutPIDNegK0s[i] = 0;
197  fh2CutPIDPosLambda[i] = 0;
198  fh2CutPIDNegLambda[i] = 0;
199 
200  fh2Tau3DVs2D[i] = 0;
201  */
202  }
203  for(Int_t i = 0; i < fgkiNCategV0; i++)
204  {
205  fh1V0InvMassK0sAll[i] = 0;
206  fh1V0InvMassLambdaAll[i] = 0;
207  fh1V0InvMassALambdaAll[i] = 0;
208  }
209  for(Int_t i = 0; i < fgkiNBinsCent; i++)
210  {
211  fh1EventCounterCutCent[i] = 0;
212  fh1V0CounterCentK0s[i] = 0;
213  fh1V0CounterCentLambda[i] = 0;
218  fh1V0InvMassK0sCent[i] = 0;
219  fh1V0InvMassLambdaCent[i] = 0;
221  fh1V0K0sPtMCGen[i] = 0;
222  fh2V0K0sPtMassMCRec[i] = 0;
223  fh1V0K0sPtMCRecFalse[i] = 0;
224  fh2V0K0sEtaPtMCGen[i] = 0;
225  fh3V0K0sEtaPtMassMCRec[i] = 0;
226  fh2V0K0sInJetPtMCGen[i] = 0;
230  fh2V0K0sMCResolMPt[i] = 0;
231  fh2V0K0sMCPtGenPtRec[i] = 0;
232  fh1V0LambdaPtMCGen[i] = 0;
233  fh2V0LambdaPtMassMCRec[i] = 0;
235  fh2V0LambdaEtaPtMCGen[i] = 0;
241  fh2V0LambdaMCResolMPt[i] = 0;
243  fhnV0LambdaInclMCFD[i] = 0;
244  fhnV0LambdaInJetsMCFD[i] = 0;
245  fhnV0LambdaBulkMCFD[i] = 0;
246  fh1V0XiPtMCGen[i] = 0;
247  fh1V0ALambdaPt[i] = 0;
248  fh1V0ALambdaPtMCGen[i] = 0;
249  fh1V0ALambdaPtMCRec[i] = 0;
252  fh2V0ALambdaEtaPtMCGen[i] = 0;
259  fh2V0ALambdaMCResolMPt[i] = 0;
261  fhnV0ALambdaInclMCFD[i] = 0;
262  fhnV0ALambdaInJetsMCFD[i] = 0;
263  fhnV0ALambdaBulkMCFD[i] = 0;
264  fh1V0AXiPtMCGen[i] = 0;
265 
266  // eta daughters
267 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
269 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
271 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
273 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
275 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
277 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
279 
280  // Inclusive
281  fhnV0InclusiveK0s[i] = 0;
282  fhnV0InclusiveLambda[i] = 0;
283  fhnV0InclusiveALambda[i] = 0;
284  // Cones
285  fhnV0InJetK0s[i] = 0;
286  fhnV0InPerpK0s[i] = 0;
287  fhnV0InRndK0s[i] = 0;
288  fhnV0InMedK0s[i] = 0;
289  fhnV0OutJetK0s[i] = 0;
290  fhnV0NoJetK0s[i] = 0;
291  fhnV0InJetLambda[i] = 0;
292  fhnV0InPerpLambda[i] = 0;
293  fhnV0InRndLambda[i] = 0;
294  fhnV0InMedLambda[i] = 0;
295  fhnV0OutJetLambda[i] = 0;
296  fhnV0NoJetLambda[i] = 0;
297  fhnV0InJetALambda[i] = 0;
298  fhnV0InPerpALambda[i] = 0;
299  fhnV0InRndALambda[i] = 0;
300  fhnV0InMedALambda[i] = 0;
301  fhnV0OutJetALambda[i] = 0;
302  fhnV0NoJetALambda[i] = 0;
303 
304  fh2V0PtJetAngleK0s[i] = 0;
305  fh2V0PtJetAngleLambda[i] = 0;
306  fh2V0PtJetAngleALambda[i] = 0;
307  fh1DCAInK0s[i] = 0;
308  fh1DCAInLambda[i] = 0;
309  fh1DCAInALambda[i] = 0;
310  fh1DCAOutK0s[i] = 0;
311  fh1DCAOutLambda[i] = 0;
312  fh1DCAOutALambda[i] = 0;
313 // fh1DeltaZK0s[i] = 0;
314 // fh1DeltaZLambda[i] = 0;
315 // fh1DeltaZALambda[i] = 0;
316 
317  fh1PtJet[i] = 0;
318  fh1EtaJet[i] = 0;
319  fh2EtaPtJet[i] = 0;
320  fh1PhiJet[i] = 0;
321  fh1NJetPerEvent[i] = 0;
322  fh2EtaPhiRndCone[i] = 0;
323  fh2EtaPhiMedCone[i] = 0;
324 
325  fh1VtxZ[i] = 0;
326  fh2VtxXY[i] = 0;
327  }
328 }
329 
330 // Constructor
332  AliAnalysisTaskSE(name),
333  fAODIn(0),
334  fAODOut(0),
335  fOutputListStd(0),
336  fOutputListQA(0),
337  fOutputListCuts(0),
338  fOutputListMC(0),
339 // ftreeOut(0),
340 
341  fiAODAnalysis(1),
342  fbIsPbPb(1),
343 
347  fdCutCPAMin(0.998),
348  fdCutNTauMax(5),
349 
350  fsJetBranchName(0),
352  fdCutPtJetMin(0),
353  fdCutPtTrackMin(5),
354  fdRadiusJet(0.4),
355  fdRadiusJetBg(0.4),
356  fbJetSelection(0),
357  fbMCAnalysis(0),
358 // fbTreeOutput(0),
359  fRandom(0),
360 
361  fdCutVertexZ(10),
362  fdCutVertexR2(1),
363  fdCutCentLow(0),
364  fdCutCentHigh(80),
365 /*
366  fBranchV0Rec(0),
367  fBranchV0Gen(0),
368  fBranchJet(0),
369  fEventInfo(0),
370 */
371  fdCentrality(0),
373  fh1EventCent(0),
374  fh1EventCent2(0),
379  fh1NRndConeCent(0),
380  fh1NMedConeCent(0),
381  fh1AreaExcluded(0),
382 
383  fh2CCK0s(0),
384  fh2CCLambda(0),
388 {
389  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
390  {
391  fh1QAV0Status[i] = 0;
392  fh1QAV0TPCRefit[i] = 0;
393  fh1QAV0TPCRows[i] = 0;
394  fh1QAV0TPCFindable[i] = 0;
395  fh1QAV0TPCRowsFind[i] = 0;
396  fh1QAV0Eta[i] = 0;
397  fh2QAV0EtaRows[i] = 0;
398  fh2QAV0PtRows[i] = 0;
399  fh2QAV0PhiRows[i] = 0;
400  fh2QAV0NClRows[i] = 0;
401  fh2QAV0EtaNCl[i] = 0;
402 
403  fh2QAV0EtaPtK0sPeak[i] = 0;
404  fh2QAV0EtaEtaK0s[i] = 0;
405  fh2QAV0PhiPhiK0s[i] = 0;
406  fh1QAV0RapK0s[i] = 0;
407  fh2QAV0PtPtK0sPeak[i] = 0;
408  fh2ArmPodK0s[i] = 0;
409 
410  fh2QAV0EtaPtLambdaPeak[i] = 0;
411  fh2QAV0EtaEtaLambda[i] = 0;
412  fh2QAV0PhiPhiLambda[i] = 0;
413  fh1QAV0RapLambda[i] = 0;
414  fh2QAV0PtPtLambdaPeak[i] = 0;
415  fh2ArmPodLambda[i] = 0;
416 
418  fh2QAV0EtaEtaALambda[i] = 0;
419  fh2QAV0PhiPhiALambda[i] = 0;
420  fh1QAV0RapALambda[i] = 0;
421  fh2QAV0PtPtALambdaPeak[i] = 0;
422  fh2ArmPodALambda[i] = 0;
423 
424  fh1QAV0Pt[i] = 0;
425  fh1QAV0Charge[i] = 0;
426  fh1QAV0DCAVtx[i] = 0;
427  fh1QAV0DCAV0[i] = 0;
428  fh1QAV0Cos[i] = 0;
429  fh1QAV0R[i] = 0;
430  fh1QACTau2D[i] = 0;
431  fh1QACTau3D[i] = 0;
432 
433  fh2ArmPod[i] = 0;
434 
435  /*
436  fh2CutTPCRowsK0s[i] = 0;
437  fh2CutTPCRowsLambda[i] = 0;
438  fh2CutPtPosK0s[i] = 0;
439  fh2CutPtNegK0s[i] = 0;
440  fh2CutPtPosLambda[i] = 0;
441  fh2CutPtNegLambda[i] = 0;
442  fh2CutDCAVtx[i] = 0;
443  fh2CutDCAV0[i] = 0;
444  fh2CutCos[i] = 0;
445  fh2CutR[i] = 0;
446  fh2CutEtaK0s[i] = 0;
447  fh2CutEtaLambda[i] = 0;
448  fh2CutRapK0s[i] = 0;
449  fh2CutRapLambda[i] = 0;
450  fh2CutCTauK0s[i] = 0;
451  fh2CutCTauLambda[i] = 0;
452  fh2CutPIDPosK0s[i] = 0;
453  fh2CutPIDNegK0s[i] = 0;
454  fh2CutPIDPosLambda[i] = 0;
455  fh2CutPIDNegLambda[i] = 0;
456 
457  fh2Tau3DVs2D[i] = 0;
458  */
459  }
460  for(Int_t i = 0; i < fgkiNCategV0; i++)
461  {
462  fh1V0InvMassK0sAll[i] = 0;
463  fh1V0InvMassLambdaAll[i] = 0;
464  fh1V0InvMassALambdaAll[i] = 0;
465  }
466  for(Int_t i = 0; i < fgkiNBinsCent; i++)
467  {
468  fh1EventCounterCutCent[i] = 0;
469  fh1V0CounterCentK0s[i] = 0;
470  fh1V0CounterCentLambda[i] = 0;
475  fh1V0InvMassK0sCent[i] = 0;
476  fh1V0InvMassLambdaCent[i] = 0;
478  fh1V0K0sPtMCGen[i] = 0;
479  fh2V0K0sPtMassMCRec[i] = 0;
480  fh1V0K0sPtMCRecFalse[i] = 0;
481  fh2V0K0sEtaPtMCGen[i] = 0;
482  fh3V0K0sEtaPtMassMCRec[i] = 0;
483  fh2V0K0sInJetPtMCGen[i] = 0;
487  fh2V0K0sMCResolMPt[i] = 0;
488  fh2V0K0sMCPtGenPtRec[i] = 0;
489  fh1V0LambdaPtMCGen[i] = 0;
490  fh2V0LambdaPtMassMCRec[i] = 0;
492  fh2V0LambdaEtaPtMCGen[i] = 0;
498  fh2V0LambdaMCResolMPt[i] = 0;
500  fhnV0LambdaInclMCFD[i] = 0;
501  fhnV0LambdaInJetsMCFD[i] = 0;
502  fhnV0LambdaBulkMCFD[i] = 0;
503  fh1V0XiPtMCGen[i] = 0;
504  fh1V0ALambdaPt[i] = 0;
505  fh1V0ALambdaPtMCGen[i] = 0;
506  fh1V0ALambdaPtMCRec[i] = 0;
509  fh2V0ALambdaEtaPtMCGen[i] = 0;
516  fh2V0ALambdaMCResolMPt[i] = 0;
518  fhnV0ALambdaInclMCFD[i] = 0;
519  fhnV0ALambdaInJetsMCFD[i] = 0;
520  fhnV0ALambdaBulkMCFD[i] = 0;
521  fh1V0AXiPtMCGen[i] = 0;
522 
523  // eta daughters
524 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
526 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
528 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
530 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
532 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
534 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
536 
537  // Inclusive
538  fhnV0InclusiveK0s[i] = 0;
539  fhnV0InclusiveLambda[i] = 0;
540  fhnV0InclusiveALambda[i] = 0;
541  // Cones
542  fhnV0InJetK0s[i] = 0;
543  fhnV0InPerpK0s[i] = 0;
544  fhnV0InRndK0s[i] = 0;
545  fhnV0InMedK0s[i] = 0;
546  fhnV0OutJetK0s[i] = 0;
547  fhnV0NoJetK0s[i] = 0;
548  fhnV0InJetLambda[i] = 0;
549  fhnV0InPerpLambda[i] = 0;
550  fhnV0InRndLambda[i] = 0;
551  fhnV0InMedLambda[i] = 0;
552  fhnV0OutJetLambda[i] = 0;
553  fhnV0NoJetLambda[i] = 0;
554  fhnV0InJetALambda[i] = 0;
555  fhnV0InPerpALambda[i] = 0;
556  fhnV0InRndALambda[i] = 0;
557  fhnV0InMedALambda[i] = 0;
558  fhnV0OutJetALambda[i] = 0;
559  fhnV0NoJetALambda[i] = 0;
560 
561  fh2V0PtJetAngleK0s[i] = 0;
562  fh2V0PtJetAngleLambda[i] = 0;
563  fh2V0PtJetAngleALambda[i] = 0;
564  fh1DCAInK0s[i] = 0;
565  fh1DCAInLambda[i] = 0;
566  fh1DCAInALambda[i] = 0;
567  fh1DCAOutK0s[i] = 0;
568  fh1DCAOutLambda[i] = 0;
569  fh1DCAOutALambda[i] = 0;
570 // fh1DeltaZK0s[i] = 0;
571 // fh1DeltaZLambda[i] = 0;
572 // fh1DeltaZALambda[i] = 0;
573 
574  fh1PtJet[i] = 0;
575  fh1EtaJet[i] = 0;
576  fh2EtaPtJet[i] = 0;
577  fh1PhiJet[i] = 0;
578  fh1NJetPerEvent[i] = 0;
579  fh2EtaPhiRndCone[i] = 0;
580  fh2EtaPhiMedCone[i] = 0;
581 
582  fh1VtxZ[i] = 0;
583  fh2VtxXY[i] = 0;
584  }
585  // Define input and output slots here
586  // Input slot #0 works with a TChain
587  DefineInput(0, TChain::Class());
588  // Output slot #0 id reserved by the base class for AOD
589  // Output slot #1 writes into a TList container
590  DefineOutput(1, TList::Class());
591  DefineOutput(2, TList::Class());
592  DefineOutput(3, TList::Class());
593  DefineOutput(4, TList::Class());
594  DefineOutput(5, TTree::Class());
595 }
596 
598 {
599 /*
600  if (fBranchV0Rec)
601  fBranchV0Rec->Delete();
602  delete fBranchV0Rec;
603  fBranchV0Rec = 0;
604  if (fBranchV0Gen)
605  fBranchV0Gen->Delete();
606  delete fBranchV0Gen;
607  fBranchV0Gen = 0;
608  if (fBranchJet)
609  fBranchJet->Delete();
610  delete fBranchJet;
611  fBranchJet = 0;
612  if (fEventInfo)
613  fEventInfo->Delete();
614  delete fEventInfo;
615  fEventInfo = 0;
616 */
617  delete fRandom;
618  fRandom = 0;
619 }
620 
622 {
623  // Create histograms
624  // Called once
625 
626  fRandom = new TRandom3(0);
627 
628 /*
629  if (!fBranchV0Rec && fbTreeOutput)
630  {
631 // fBranchV0Rec = new TClonesArray("AliAODv0",0);
632  fBranchV0Rec = new TClonesArray("AliV0Object",0);
633  fBranchV0Rec->SetName("branch_V0Rec");
634  }
635  if (!fBranchV0Gen && fbTreeOutput)
636  {
637  fBranchV0Gen = new TClonesArray("AliAODMCParticle",0);
638  fBranchV0Gen->SetName("branch_V0Gen");
639  }
640  if (!fBranchJet && fbTreeOutput)
641  {
642 // fBranchJet = new TClonesArray("AliAODJet",0);
643  fBranchJet = new TClonesArray("AliJetObject",0);
644  fBranchJet->SetName("branch_Jet");
645  }
646  if (!fEventInfo && fbTreeOutput)
647  {
648  fEventInfo = new AliEventInfoObject();
649  fEventInfo->SetName("eventInfo");
650  }
651  Int_t dSizeBuffer = 32000; // default 32000
652  if (fbTreeOutput)
653  {
654  ftreeOut = new TTree("treeV0","Tree V0");
655  ftreeOut->Branch("branch_V0Rec",&fBranchV0Rec,dSizeBuffer,2);
656  ftreeOut->Branch("branch_V0Gen",&fBranchV0Gen,dSizeBuffer,2);
657  ftreeOut->Branch("branch_Jet",&fBranchJet,dSizeBuffer,2);
658  ftreeOut->Branch("eventInfo",&fEventInfo,dSizeBuffer,2);
659  }
660 */
661 
662  fOutputListStd = new TList();
663  fOutputListStd->SetOwner();
664  fOutputListQA = new TList();
665  fOutputListQA->SetOwner();
666  fOutputListCuts = new TList();
667  fOutputListCuts->SetOwner();
668  fOutputListMC = new TList();
669  fOutputListMC->SetOwner();
670 
671  // event categories
672  const Int_t iNCategEvent = 6;
673  TString categEvent[iNCategEvent] = {"coll. candid.", "AOD OK", "vtx & cent", "with V0", "with jets", "jet selection"};
674  // labels for stages of V0 selection
675  TString categV0[fgkiNCategV0] = {"all"/*0*/, "mass range"/*1*/, "rec. method"/*2*/, "tracks TPC"/*3*/, "track pt"/*4*/, "DCA prim v"/*5*/, "DCA daughters"/*6*/, "CPA"/*7*/, "volume"/*8*/, "track #it{#eta}"/*9*/, "V0 #it{y} & #it{#eta}"/*10*/, "lifetime"/*11*/, "PID"/*12*/, "Arm.-Pod."/*13*/, "inclusive"/*14*/, "in jet event"/*15*/, "in jet"/*16*/};
676 
677  fh1EventCounterCut = new TH1D("fh1EventCounterCut", "Number of events after filtering;selection filter;counts", iNCategEvent, 0, iNCategEvent);
678  for(Int_t i = 0; i < iNCategEvent; i++)
679  fh1EventCounterCut->GetXaxis()->SetBinLabel(i + 1, categEvent[i].Data());
680  fh1EventCent2 = new TH1D("fh1EventCent2", "Number of events vs centrality;centrality;counts", 100, 0, 100);
681  fh1EventCent2Jets = new TH1D("fh1EventCent2Jets", "Number of sel.-jet events vs centrality;centrality;counts", 100, 0, 100);
682  fh1EventCent2NoJets = new TH1D("fh1EventCent2NoJets", "Number of no-jet events vs centrality;centrality;counts", 100, 0, 100);
683  fh2EventCentTracks = new TH2D("fh2EventCentTracks", "Number of tracks vs centrality;centrality;tracks;counts", 100, 0, 100, 150, 0, 15e3);
684  fh1EventCent = new TH1D("fh1EventCent", "Number of events in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
685  for(Int_t i = 0; i < fgkiNBinsCent; i++)
686  fh1EventCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
687  fh1NRndConeCent = new TH1D("fh1NRndConeCent", "Number of rnd. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
688  for(Int_t i = 0; i < fgkiNBinsCent; i++)
689  fh1NRndConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
690  fh1NMedConeCent = new TH1D("fh1NMedConeCent", "Number of med.-cl. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
691  for(Int_t i = 0; i < fgkiNBinsCent; i++)
692  fh1NMedConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
693  fh1AreaExcluded = new TH1D("fh1AreaExcluded", "Area of excluded cones in centrality bins;centrality;area", fgkiNBinsCent, 0, fgkiNBinsCent);
694  for(Int_t i = 0; i < fgkiNBinsCent; i++)
695  fh1AreaExcluded->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
705 
706  fh1V0CandPerEvent = new TH1D("fh1V0CandPerEvent", "Number of all V0 candidates per event;candidates;events", 1000, 0, 1000);
708 
709  for(Int_t i = 0; i < fgkiNBinsCent; i++)
710  {
711  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);
712  for(Int_t j = 0; j < iNCategEvent; j++)
713  fh1EventCounterCutCent[i]->GetXaxis()->SetBinLabel(j + 1, categEvent[j].Data());
714  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, 100);
715  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, 100);
716  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, 100);
717  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);
718  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);
719  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);
720  for(Int_t j = 0; j < fgkiNCategV0; j++)
721  {
722  fh1V0CounterCentK0s[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
723  fh1V0CounterCentLambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
724  fh1V0CounterCentALambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
725  }
733  }
734  // pt binning for V0 and jets
735  Int_t iNBinsPtV0 = fgkiNBinsPtV0Init;
736  Double_t dPtV0Min = fgkdBinsPtV0[0];
737  Double_t dPtV0Max = fgkdBinsPtV0[fgkiNBinsPtV0];
738  Int_t iNJetPtBins = fgkiNBinsPtJetInit;
739  Double_t dJetPtMin = fgkdBinsPtJet[0];
741 
742  fh2CCK0s = new TH2D("fh2CCK0s", "K0s candidates in Lambda peak", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
743  fh2CCLambda = new TH2D("fh2CCLambda", "Lambda candidates in K0s peak", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
744  Int_t binsCorrel[3] = {fgkiNBinsMassK0s, fgkiNBinsMassLambda, iNBinsPtV0};
745  Double_t xminCorrel[3] = {fgkdMassK0sMin, fgkdMassLambdaMin, dPtV0Min};
746  Double_t xmaxCorrel[3] = {fgkdMassK0sMax, fgkdMassLambdaMax, dPtV0Max};
747 // Int_t binsCorrel[3] = {200, 200, iNBinsPtV0};
748 // Double_t xminCorrel[3] = {0, 0, dPtV0Min};
749 // Double_t xmaxCorrel[3] = {2, 2, dPtV0Max};
750  fh3CCMassCorrelBoth = new THnSparseD("fh3CCMassCorrelBoth", "Mass correlation: K0S && Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
751  fh3CCMassCorrelKNotL = new THnSparseD("fh3CCMassCorrelKNotL", "Mass correlation: K0S, not Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
752  fh3CCMassCorrelLNotK = new THnSparseD("fh3CCMassCorrelLNotK", "Mass correlation: Lambda, not K0S;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
753  fOutputListQA->Add(fh2CCK0s);
758 
759  Double_t dStepEtaV0 = 0.025;
760  Double_t dRangeEtaV0Max = 0.8;
761  const Int_t iNBinsEtaV0 = 2 * Int_t(dRangeEtaV0Max / dStepEtaV0);
762  // inclusive
763  const Int_t iNDimIncl = 3;
764  Int_t binsKIncl[iNDimIncl] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
765  Double_t xminKIncl[iNDimIncl] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
766  Double_t xmaxKIncl[iNDimIncl] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
767  Int_t binsLIncl[iNDimIncl] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
768  Double_t xminLIncl[iNDimIncl] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
769  Double_t xmaxLIncl[iNDimIncl] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
770  // binning in jets
771  const Int_t iNDimInJC = 4;
772  Int_t binsKInJC[iNDimInJC] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
773  Double_t xminKInJC[iNDimInJC] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
774  Double_t xmaxKInJC[iNDimInJC] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
775  Int_t binsLInJC[iNDimInJC] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
776  Double_t xminLInJC[iNDimInJC] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
777  Double_t xmaxLInJC[iNDimInJC] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
778 
779  // binning eff inclusive vs eta-pT
780  Double_t dStepDeltaEta = 0.1;
781  Double_t dRangeDeltaEtaMax = 0.5;
782  const Int_t iNBinsDeltaEta = 2 * Int_t(dRangeDeltaEtaMax / dStepDeltaEta);
783  Int_t binsEtaK[3] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
784  Double_t xminEtaK[3] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
785  Double_t xmaxEtaK[3] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
786  Int_t binsEtaL[3] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
787  Double_t xminEtaL[3] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
788  Double_t xmaxEtaL[3] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
789  // binning eff in jets vs eta-pT
790  // associated
791  Int_t binsEtaKInRec[5] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
792  Double_t xminEtaKInRec[5] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
793  Double_t xmaxEtaKInRec[5] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
794  Int_t binsEtaLInRec[5] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
795  Double_t xminEtaLInRec[5] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
796  Double_t xmaxEtaLInRec[5] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
797  // generated
798  Int_t binsEtaInGen[4] = {iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
799  Double_t xminEtaInGen[4] = {dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
800  Double_t xmaxEtaInGen[4] = {dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
801  // daughter eta: charge-etaD-ptD-etaV0-ptV0-ptJet
802  const Int_t iNDimEtaD = 6;
803  Int_t binsEtaDaughter[iNDimEtaD] = {2, 20, iNBinsPtV0, iNBinsEtaV0, iNBinsPtV0, iNJetPtBins};
804  Double_t xminEtaDaughter[iNDimEtaD] = {0, -1, dPtV0Min, -dRangeEtaV0Max, dPtV0Min, dJetPtMin};
805  Double_t xmaxEtaDaughter[iNDimEtaD] = {2, 1, dPtV0Max, dRangeEtaV0Max, dPtV0Max, dJetPtMax};
806 
807  for(Int_t i = 0; i < fgkiNBinsCent; i++)
808  {
809  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);
810  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);
811  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);
815  // Inclusive
816  fhnV0InclusiveK0s[i] = new THnSparseD(Form("fhnV0InclusiveK0s_C%d", i), "K0s: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
817  fhnV0InclusiveLambda[i] = new THnSparseD(Form("fhnV0InclusiveLambda_C%d", i), "Lambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
818  fhnV0InclusiveALambda[i] = new THnSparseD(Form("fhnV0InclusiveALambda_C%d", i), "ALambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
822  // In cones
823  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
824  fOutputListStd->Add(fhnV0InJetK0s[i]);
825  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
827  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
828  fOutputListStd->Add(fhnV0InRndK0s[i]);
829  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
830  fOutputListStd->Add(fhnV0InMedK0s[i]);
831  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
833  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
834  fOutputListStd->Add(fhnV0NoJetK0s[i]);
835  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
837  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
839  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
841  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
843  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
845  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
847  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
849  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{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
851  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
853  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
855  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
857  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})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
859 
860  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()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
862  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()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
864  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()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
866 
867  fh1DCAInK0s[i] = new TH1D(Form("fh1DCAInK0s_%d", i), Form("K0s in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
868  fOutputListQA->Add(fh1DCAInK0s[i]);
869  fh1DCAInLambda[i] = new TH1D(Form("fh1DCAInLambda_%d", i), Form("Lambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
870  fOutputListQA->Add(fh1DCAInLambda[i]);
871  fh1DCAInALambda[i] = new TH1D(Form("fh1DCAInALambda_%d", i), Form("ALambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
873 
874  fh1DCAOutK0s[i] = new TH1D(Form("fh1DCAOutK0s_%d", i), Form("K0s outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
875  fOutputListQA->Add(fh1DCAOutK0s[i]);
876  fh1DCAOutLambda[i] = new TH1D(Form("fh1DCAOutLambda_%d", i), Form("Lambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
878  fh1DCAOutALambda[i] = new TH1D(Form("fh1DCAOutALambda_%d", i), Form("ALambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
880 
881  /*
882  fh1DeltaZK0s[i] = new TH1D(Form("fh1DeltaZK0s_%d", i), Form("K0s: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
883  fOutputListQA->Add(fh1DeltaZK0s[i]);
884  fh1DeltaZLambda[i] = new TH1D(Form("fh1DeltaZLambda_%d", i), Form("Lambda: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
885  fOutputListQA->Add(fh1DeltaZLambda[i]);
886  fh1DeltaZALambda[i] = new TH1D(Form("fh1DeltaZALambda_%d", i), Form("ALambda: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
887  fOutputListQA->Add(fh1DeltaZALambda[i]);
888  */
889 
890  // jet histograms
891  fh1PtJet[i] = new TH1D(Form("fh1PtJet_%d", i), Form("Jet pt spectrum, cent: %s;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax);
892  fOutputListStd->Add(fh1PtJet[i]);
893  fh1EtaJet[i] = new TH1D(Form("fh1EtaJet_%d", i), Form("Jet eta spectrum, cent: %s;#it{#eta} jet", GetCentBinLabel(i).Data()), 80, -1., 1.);
894  fOutputListStd->Add(fh1EtaJet[i]);
895  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., iNJetPtBins, dJetPtMin, dJetPtMax);
896  fOutputListStd->Add(fh2EtaPtJet[i]);
897  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., 100, 0., TMath::TwoPi());
899  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., 100, 0., TMath::TwoPi());
901  fh1PhiJet[i] = new TH1D(Form("fh1PhiJet_%d", i), Form("Jet phi spectrum, cent: %s;#it{#phi} jet", GetCentBinLabel(i).Data()), 100, 0., TMath::TwoPi());
902  fOutputListStd->Add(fh1PhiJet[i]);
903  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.);
905  // event histograms
906  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);
907  fOutputListQA->Add(fh1VtxZ[i]);
908  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);
909  fOutputListQA->Add(fh2VtxXY[i]);
910  // fOutputListStd->Add([i]);
911  if(fbMCAnalysis)
912  {
913  // inclusive pt
914  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);
916  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);
918  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);
920  // inclusive pt-eta
921  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);
923  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);
925  // in jet pt
926  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()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
928  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);
930  // in jet pt-eta
931  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);
933  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);
935 
936  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);
938  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);
940 
941  // inclusive pt
942  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);
944  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);
946  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);
948  // inclusive pt-eta
949  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);
951  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);
953  // in jet pt
954  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()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
956  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);
958  // in jet pt-eta
959  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);
961  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);
963 
964  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);
966  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);
968 
969  // inclusive pt
970  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);
972  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);
974  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);
976  // inclusive pt-eta
977  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);
979  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);
981  // in jet pt
982  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()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
984  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);
986  // in jet pt-eta
987  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);
989  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);
991 
992  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);
994  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);
996 
997  Int_t iNBinsPtXi = 80;
998  Double_t dPtXiMin = 0;
999  Double_t dPtXiMax = 8;
1000  const Int_t iNDimFD = 3;
1001  Int_t binsFD[iNDimFD] = {iNBinsPtV0, iNBinsPtXi, iNJetPtBins};
1002  Double_t xminFD[iNDimFD] = {dPtV0Min, dPtXiMin, dJetPtMin};
1003  Double_t xmaxFD[iNDimFD] = {dPtV0Max, dPtXiMax, dJetPtMax};
1004  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);
1006  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);
1008  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);
1010  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);
1011  fOutputListMC->Add(fh1V0XiPtMCGen[i]);
1012  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);
1014  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);
1016  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);
1018  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);
1019  fOutputListMC->Add(fh1V0AXiPtMCGen[i]);
1020 
1021  // daughter eta
1022 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1023  fhnV0K0sInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1024 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1025  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 daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1026 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1027  fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1028 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1029  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 daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1030 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1031  fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1032 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1033  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 daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1034 
1035 // fOutputListMC->Add(fhnV0K0sInclDaughterEtaPtPtMCGen[i]);
1037 // fOutputListMC->Add(fhnV0K0sInJetsDaughterEtaPtPtMCGen[i]);
1039 // fOutputListMC->Add(fhnV0LambdaInclDaughterEtaPtPtMCGen[i]);
1041 // fOutputListMC->Add(fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i]);
1043 // fOutputListMC->Add(fhnV0ALambdaInclDaughterEtaPtPtMCGen[i]);
1045 // fOutputListMC->Add(fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i]);
1047  }
1048  }
1049 
1050  // QA Histograms
1051  for(Int_t i = 0; i < fgkiNQAIndeces; i++)
1052  {
1053 // [i] = new TH1D(Form("%d",i),";;Counts",,,);
1054  fh1QAV0Status[i] = new TH1D(Form("fh1QAV0Status_%d", i), "QA: V0 status", 2, 0, 2);
1055  fh1QAV0TPCRefit[i] = new TH1D(Form("fh1QAV0TPCRefit_%d", i), "QA: TPC refit", 2, 0, 2);
1056  fh1QAV0TPCRows[i] = new TH1D(Form("fh1QAV0TPCRows_%d", i), "QA: TPC Rows", 160, 0, 160);
1057  fh1QAV0TPCFindable[i] = new TH1D(Form("fh1QAV0TPCFindable_%d", i), "QA: TPC Findable", 160, 0, 160);
1058  fh1QAV0TPCRowsFind[i] = new TH1D(Form("fh1QAV0TPCRowsFind_%d", i), "QA: TPC Rows/Findable", 100, 0, 2);
1059  fh1QAV0Eta[i] = new TH1D(Form("fh1QAV0Eta_%d", i), "QA: Daughter Eta", 200, -2, 2);
1060  fh2QAV0EtaRows[i] = new TH2D(Form("fh2QAV0EtaRows_%d", i), "QA: Daughter Eta vs TPC rows;#eta;TPC rows", 200, -2, 2, 160, 0, 160);
1061  fh2QAV0PtRows[i] = new TH2D(Form("fh2QAV0PtRows_%d", i), "QA: Daughter Pt vs TPC rows;pt;TPC rows", 100, 0, 10, 160, 0, 160);
1062  fh2QAV0PhiRows[i] = new TH2D(Form("fh2QAV0PhiRows_%d", i), "QA: Daughter Phi vs TPC rows;#phi;TPC rows", 100, 0, TMath::TwoPi(), 160, 0, 160);
1063  fh2QAV0NClRows[i] = new TH2D(Form("fh2QAV0NClRows_%d", i), "QA: Daughter NCl vs TPC rows;findable clusters;TPC rows", 100, 0, 160, 160, 0, 160);
1064  fh2QAV0EtaNCl[i] = new TH2D(Form("fh2QAV0EtaNCl_%d", i), "QA: Daughter Eta vs NCl;#eta;findable clusters", 200, -2, 2, 160, 0, 160);
1065 
1066  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);
1067  fh2QAV0EtaEtaK0s[i] = new TH2D(Form("fh2QAV0EtaEtaK0s_%d", i), "QA: K0s: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1068  fh2QAV0PhiPhiK0s[i] = new TH2D(Form("fh2QAV0PhiPhiK0s_%d", i), "QA: K0s: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1069  fh1QAV0RapK0s[i] = new TH1D(Form("fh1QAV0RapK0s_%d", i), "QA: K0s: V0 Rapidity", 200, -2, 2);
1070  fh2QAV0PtPtK0sPeak[i] = new TH2D(Form("fh2QAV0PtPtK0sPeak_%d", i), "QA: K0s: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1071 
1072  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);
1073  fh2QAV0EtaEtaLambda[i] = new TH2D(Form("fh2QAV0EtaEtaLambda_%d", i), "QA: Lambda: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1074  fh2QAV0PhiPhiLambda[i] = new TH2D(Form("fh2QAV0PhiPhiLambda_%d", i), "QA: Lambda: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1075  fh1QAV0RapLambda[i] = new TH1D(Form("fh1QAV0RapLambda_%d", i), "QA: Lambda: V0 Rapidity", 200, -2, 2);
1076  fh2QAV0PtPtLambdaPeak[i] = new TH2D(Form("fh2QAV0PtPtLambdaPeak_%d", i), "QA: Lambda: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1077 
1078  fh1QAV0Pt[i] = new TH1D(Form("fh1QAV0Pt_%d", i), "QA: Daughter Pt", 100, 0, 5);
1079  fh1QAV0Charge[i] = new TH1D(Form("fh1QAV0Charge_%d", i), "QA: V0 Charge", 3, -1, 2);
1080  fh1QAV0DCAVtx[i] = new TH1D(Form("fh1QAV0DCAVtx_%d", i), "QA: DCA daughters to primary vertex", 100, 0, 10);
1081  fh1QAV0DCAV0[i] = new TH1D(Form("fh1QAV0DCAV0_%d", i), "QA: DCA daughters", 100, 0, 2);
1082  fh1QAV0Cos[i] = new TH1D(Form("fh1QAV0Cos_%d", i), "QA: CPA", 10000, 0.9, 1);
1083  fh1QAV0R[i] = new TH1D(Form("fh1QAV0R_%d", i), "QA: R", 1500, 0, 150);
1084  fh1QACTau2D[i] = new TH1D(Form("fh1QACTau2D_%d", i), "QA: K0s: c#tau 2D;mR/pt#tau", 100, 0, 10);
1085  fh1QACTau3D[i] = new TH1D(Form("fh1QACTau3D_%d", i), "QA: K0s: c#tau 3D;mL/p#tau", 100, 0, 10);
1086 
1087  fh2ArmPod[i] = new TH2D(Form("fh2ArmPod_%d", i), "Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1088  fh2ArmPodK0s[i] = new TH2D(Form("fh2ArmPodK0s_%d", i), "K0s: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1089  fh2ArmPodLambda[i] = new TH2D(Form("fh2ArmPodLambda_%d", i), "Lambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1090  fh2ArmPodALambda[i] = new TH2D(Form("fh2ArmPodALambda_%d", i), "ALambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1091 
1092  fOutputListQA->Add(fh1QAV0Status[i]);
1093  fOutputListQA->Add(fh1QAV0TPCRefit[i]);
1094  fOutputListQA->Add(fh1QAV0TPCRows[i]);
1097  fOutputListQA->Add(fh1QAV0Eta[i]);
1098  fOutputListQA->Add(fh2QAV0EtaRows[i]);
1099  fOutputListQA->Add(fh2QAV0PtRows[i]);
1100  fOutputListQA->Add(fh2QAV0PhiRows[i]);
1101  fOutputListQA->Add(fh2QAV0NClRows[i]);
1102  fOutputListQA->Add(fh2QAV0EtaNCl[i]);
1103 
1107  fOutputListQA->Add(fh1QAV0RapK0s[i]);
1109 
1115 
1116  fOutputListQA->Add(fh1QAV0Pt[i]);
1117  fOutputListQA->Add(fh1QAV0Charge[i]);
1118  fOutputListQA->Add(fh1QAV0DCAVtx[i]);
1119  fOutputListQA->Add(fh1QAV0DCAV0[i]);
1120  fOutputListQA->Add(fh1QAV0Cos[i]);
1121  fOutputListQA->Add(fh1QAV0R[i]);
1122  fOutputListQA->Add(fh1QACTau2D[i]);
1123  fOutputListQA->Add(fh1QACTau3D[i]);
1124 
1125  fOutputListQA->Add(fh2ArmPod[i]);
1126  fOutputListQA->Add(fh2ArmPodK0s[i]);
1127  fOutputListQA->Add(fh2ArmPodLambda[i]);
1129 
1130  /*
1131  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);
1132  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);
1133  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);
1134  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);
1135  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);
1136  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);
1137  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);
1138  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);
1139  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);
1140  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);
1141  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);
1142  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);
1143  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);
1144  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);
1145  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);
1146  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);
1147  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);
1148  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);
1149  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);
1150  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);
1151  fh2Tau3DVs2D[i] = new TH2D(Form("fh2Tau3DVs2D_%d", i), "Decay 3D vs 2D;pt;3D/2D", 100, 0, 10, 200, 0.5, 1.5);
1152 
1153  fOutputListCuts->Add(fh2CutTPCRowsK0s[i]);
1154  fOutputListCuts->Add(fh2CutTPCRowsLambda[i]);
1155  fOutputListCuts->Add(fh2CutPtPosK0s[i]);
1156  fOutputListCuts->Add(fh2CutPtNegK0s[i]);
1157  fOutputListCuts->Add(fh2CutPtPosLambda[i]);
1158  fOutputListCuts->Add(fh2CutPtNegLambda[i]);
1159  fOutputListCuts->Add(fh2CutDCAVtx[i]);
1160  fOutputListCuts->Add(fh2CutDCAV0[i]);
1161  fOutputListCuts->Add(fh2CutCos[i]);
1162  fOutputListCuts->Add(fh2CutR[i]);
1163  fOutputListCuts->Add(fh2CutEtaK0s[i]);
1164  fOutputListCuts->Add(fh2CutEtaLambda[i]);
1165  fOutputListCuts->Add(fh2CutRapK0s[i]);
1166  fOutputListCuts->Add(fh2CutRapLambda[i]);
1167  fOutputListCuts->Add(fh2CutCTauK0s[i]);
1168  fOutputListCuts->Add(fh2CutCTauLambda[i]);
1169  fOutputListCuts->Add(fh2CutPIDPosK0s[i]);
1170  fOutputListCuts->Add(fh2CutPIDNegK0s[i]);
1171  fOutputListCuts->Add(fh2CutPIDPosLambda[i]);
1172  fOutputListCuts->Add(fh2CutPIDNegLambda[i]);
1173  fOutputListCuts->Add(fh2Tau3DVs2D[i]);
1174  */
1175  }
1176 
1177  for(Int_t i = 0; i < fgkiNCategV0; i++)
1178  {
1179  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);
1180  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);
1181  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);
1185  }
1186 
1187  for(Int_t i = 0; i < fOutputListStd->GetEntries(); ++i)
1188  {
1189  TH1* h1 = dynamic_cast<TH1*>(fOutputListStd->At(i));
1190  if(h1)
1191  {
1192  h1->Sumw2();
1193  continue;
1194  }
1195  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListStd->At(i));
1196  if(hn) hn->Sumw2();
1197  }
1198  for(Int_t i = 0; i < fOutputListQA->GetEntries(); ++i)
1199  {
1200  TH1* h1 = dynamic_cast<TH1*>(fOutputListQA->At(i));
1201  if(h1)
1202  {
1203  h1->Sumw2();
1204  continue;
1205  }
1206  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListQA->At(i));
1207  if(hn) hn->Sumw2();
1208  }
1209  for(Int_t i = 0; i < fOutputListCuts->GetEntries(); ++i)
1210  {
1211  TH1* h1 = dynamic_cast<TH1*>(fOutputListCuts->At(i));
1212  if(h1)
1213  {
1214  h1->Sumw2();
1215  continue;
1216  }
1217  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListCuts->At(i));
1218  if(hn) hn->Sumw2();
1219  }
1220  for(Int_t i = 0; i < fOutputListMC->GetEntries(); ++i)
1221  {
1222  TH1* h1 = dynamic_cast<TH1*>(fOutputListMC->At(i));
1223  if(h1)
1224  {
1225  h1->Sumw2();
1226  continue;
1227  }
1228  THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListMC->At(i));
1229  if(hn) hn->Sumw2();
1230  }
1231 
1232  PostData(1, fOutputListStd);
1233  PostData(2, fOutputListQA);
1234  PostData(3, fOutputListCuts);
1235  PostData(4, fOutputListMC);
1236 // if (fbTreeOutput)
1237 // PostData(5,ftreeOut);
1238 }
1239 
1241 {
1242  // Main loop, called for each event
1243  if(fDebug > 5) printf("TaskV0sInJets: UserExec: Start\n");
1244 /*
1245  // reset branches for each event
1246  if (fBranchV0Rec)
1247  fBranchV0Rec->Clear();
1248  if (fBranchV0Gen)
1249  fBranchV0Gen->Clear();
1250  if (fBranchJet)
1251  fBranchJet->Clear();
1252  if (fEventInfo)
1253  fEventInfo->Reset();
1254 */
1255  if(!fiAODAnalysis)
1256  return;
1257 
1258  if(fDebug > 2) printf("TaskV0sInJets: AOD analysis\n");
1259  fh1EventCounterCut->Fill(0); // all available selected events (collision candidates)
1260 
1261  if(fDebug > 5) printf("TaskV0sInJets: UserExec: Loading AOD\n");
1262  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); // input AOD
1263  fAODOut = AODEvent(); // output AOD
1264  if(!fAODOut)
1265  {
1266  if(fDebug > 0) printf("TaskV0sInJets: No output AOD found\n");
1267  return;
1268  }
1269  if(!fAODIn)
1270  {
1271  if(fDebug > 0) printf("TaskV0sInJets: No input AOD found\n");
1272  return;
1273  }
1274  if(fDebug > 5) printf("TaskV0sInJets: UserExec: Loading AOD OK\n");
1275 
1276  TClonesArray* arrayMC = 0; // array particles in the MC event
1277  AliAODMCHeader* headerMC = 0; // MC header
1278  Int_t iNTracksMC = 0; // number of MC tracks
1279  Double_t dPrimVtxMCX = 0., dPrimVtxMCY = 0., dPrimVtxMCZ = 0.; // position of the MC primary vertex
1280 
1281  // Simulation info
1282  if(fbMCAnalysis)
1283  {
1284  arrayMC = (TClonesArray*)fAODIn->FindListObject(AliAODMCParticle::StdBranchName());
1285  if(!arrayMC)
1286  {
1287  if(fDebug > 0) printf("TaskV0sInJets: No MC array found\n");
1288  return;
1289  }
1290  if(fDebug > 5) printf("TaskV0sInJets: MC array found\n");
1291  iNTracksMC = arrayMC->GetEntriesFast();
1292  if(fDebug > 5) printf("TaskV0sInJets: There are %d MC tracks in this event\n", iNTracksMC);
1293 // if (!iNTracksMC)
1294 // return;
1295  headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
1296  if(!headerMC)
1297  {
1298  if(fDebug > 0) printf("TaskV0sInJets: No MC header found\n");
1299  return;
1300  }
1301  // get position of the MC primary vertex
1302  dPrimVtxMCX = headerMC->GetVtxX();
1303  dPrimVtxMCY = headerMC->GetVtxY();
1304  dPrimVtxMCZ = headerMC->GetVtxZ();
1305  }
1306 
1307  // PID Response Task object
1308  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1309  AliInputEventHandler* inputHandler = (AliInputEventHandler*)mgr->GetInputEventHandler();
1310  AliPIDResponse* fPIDResponse = inputHandler->GetPIDResponse();
1311  if(!fPIDResponse)
1312  {
1313  if(fDebug > 0) printf("TaskV0sInJets: No PID response object found\n");
1314  return;
1315  }
1316 
1317  // AOD files are OK
1318  fh1EventCounterCut->Fill(1);
1319 
1320  // Event selection
1321  if(!IsSelectedForJets(fAODIn, fdCutVertexZ, fdCutVertexR2, fdCutCentLow, fdCutCentHigh, 1, 0.1)) // cut on |delta z| in 2011 data between SPD vertex and nominal primary vertex
1322 // if (!IsSelectedForJets(fAODIn,fdCutVertexZ,fdCutVertexR2,fdCutCentLow,fdCutCentHigh)) // no need for cutting in 2010 data
1323  {
1324  if(fDebug > 5) printf("TaskV0sInJets: Event rejected\n");
1325  return;
1326  }
1327 
1328 // fdCentrality = ((AliVAODHeader*)fAODIn->GetHeader())->GetCentrality(); // event centrality
1329  fdCentrality = ((AliVAODHeader*)fAODIn->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M"); // event centrality
1330  if(!fbIsPbPb)
1331  fdCentrality = 0.;
1332  Int_t iCentIndex = GetCentralityBinIndex(fdCentrality); // get index of centrality bin
1333  if(iCentIndex < 0)
1334  {
1335  if(fDebug > 5) printf("TaskV0sInJets: Event is out of histogram range\n");
1336  return;
1337  }
1338  fh1EventCounterCut->Fill(2); // selected events (vertex, centrality)
1339  fh1EventCounterCutCent[iCentIndex]->Fill(2);
1340 
1341  UInt_t iNTracks = fAODIn->GetNumberOfTracks(); // get number of tracks in event
1342  if(fDebug > 5) printf("TaskV0sInJets: There are %d tracks in this event\n", iNTracks);
1343 // if (!iNTracks)
1344 // return;
1345 
1346  Int_t iNV0s = fAODIn->GetNumberOfV0s(); // get the number of V0 candidates
1347  if(!iNV0s)
1348  {
1349  if(fDebug > 2) printf("TaskV0sInJets: No V0s found in event\n");
1350 // return;
1351  }
1352 
1353  //===== Event is OK for the analysis =====
1354  fh1EventCent->Fill(iCentIndex);
1355  fh1EventCent2->Fill(fdCentrality);
1356  fh2EventCentTracks->Fill(fdCentrality, iNTracks);
1357 
1358 // if (fbTreeOutput)
1359 // fEventInfo->SetAODEvent(fAODIn);
1360 // printf("V0sInJets: EventInfo: Centrality: %f\n",fEventInfo->GetCentrality());
1361 
1362  if(iNV0s)
1363  {
1364  fh1EventCounterCut->Fill(3); // events with V0s
1365  fh1EventCounterCutCent[iCentIndex]->Fill(3);
1366  }
1367 
1368 // Int_t iNV0SelV0Rec = 0;
1369 // Int_t iNV0SelV0Gen = 0;
1370 
1371  AliAODv0* v0 = 0; // pointer to V0 candidates
1372 // AliV0Object* objectV0 = 0;
1373  TVector3 vecV0Momentum; // 3D vector of V0 momentum
1374  Double_t dMassV0K0s = 0; // invariant mass of the K0s candidate
1375  Double_t dMassV0Lambda = 0; // invariant mass of the Lambda candidate
1376  Double_t dMassV0ALambda = 0; // invariant mass of the Lambda candidate
1377  Int_t iNV0CandTot = 0; // counter of all V0 candidates at the beginning
1378  Int_t iNV0CandK0s = 0; // counter of K0s candidates at the end
1379  Int_t iNV0CandLambda = 0; // counter of Lambda candidates at the end
1380  Int_t iNV0CandALambda = 0; // counter of Lambda candidates at the end
1381 
1382  Bool_t bUseOldCuts = 0; // old reconstruction cuts
1383  Bool_t bUseAliceCuts = 0; // cuts used by Alice Zimmermann
1384  Bool_t bUseIouriCuts = 0; // cuts used by Iouri
1385  Bool_t bPrintCuts = 0; // print out which cuts are applied
1386  Bool_t bPrintJetSelection = 0; // print out which jets are selected
1387 
1388  // Values of V0 reconstruction cuts:
1389  // Daughter tracks
1390  Int_t iRefit = AliAODTrack::kTPCrefit; // TPC refit for daughter tracks
1391  Double_t dDCAToPrimVtxMin = fdCutDCAToPrimVtxMin; // 0.1; // [cm] min DCA of daughters to the prim vtx
1392  Double_t dDCADaughtersMax = fdCutDCADaughtersMax; // 1.; // [sigma of TPC tracking] max DCA between daughters
1393  Double_t dEtaDaughterMax = 0.8; // max |pseudorapidity| of daughter tracks
1394  Double_t dNSigmadEdxMax = fdCutNSigmadEdxMax;// 3.; // [sigma dE/dx] max difference between measured and expected signal of dE/dx in the TPC
1395  Double_t dPtProtonPIDMax = 1.; // [GeV/c] maxium pT of proton for applying PID cut
1396  // V0 candidate
1397  Bool_t bOnFly = 0; // on-the-fly (yes) or offline (no) reconstructed
1398  Double_t dCPAMin = fdCutCPAMin;// 0.998; // min cosine of the pointing angle
1399  Double_t dRadiusDecayMin = 5.; // [cm] min radial distance of the decay vertex
1400  Double_t dRadiusDecayMax = 100.; // [cm] max radial distance of the decay vertex
1401  Double_t dEtaMax = 0.7; // max |pseudorapidity| of V0
1402  Double_t dNTauMax = fdCutNTauMax; // 5.0; // [tau] max proper lifetime in multiples of the mean lifetime
1403 
1404  // Old cuts Start
1405  Double_t dNCrossedRowsTPCMin = 70.; // min number of crossed TPC rows (turned off)
1406 // Double_t dCrossedRowsOverFindMin = 0.8; // min ratio crossed rows / findable clusters (turned off)
1407 // Double_t dCrossedRowsOverFindMax = 1e3; // max ratio crossed rows / findable clusters (turned off)
1408  Double_t dPtDaughterMin = 0.150; // [GeV/c] min transverse momentum of daughter tracks (turned off)
1409  Double_t dRapMax = 0.75; // max |rapidity| of V0 (turned off)
1410  // Old cuts End
1411 
1412  // Other cuts
1413  Double_t dNSigmaMassMax = 3.; // [sigma m] max difference between candidate mass and real particle mass (used only for mass peak method of signal extraction)
1414  Double_t dDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
1415 
1416  // Selection of active cuts
1417  Bool_t bCutEtaDaughter = 1; // daughter pseudorapidity
1418  Bool_t bCutRapV0 = 0; // V0 rapidity
1419  Bool_t bCutEtaV0 = 1; // V0 pseudorapidity
1420  Bool_t bCutTau = 1; // V0 lifetime
1421  Bool_t bCutPid = 1; // PID (TPC dE/dx)
1422  Bool_t bCutArmPod = 1; // Armenteros-Podolanski for K0S
1423 // Bool_t bCutCross = 0; // cross contamination
1424 
1425  if(bUseOldCuts)
1426  {
1427  bCutRapV0 = 1;
1428  dEtaMax = 0.75;
1429  dNTauMax = 3.0;
1430  }
1431  else if(bUseAliceCuts)
1432  {
1433 // bOnFly = 1;
1434  dEtaMax = 0.75;
1435  dNTauMax = 5.0;
1436  }
1437  else if(bUseIouriCuts)
1438  {
1439  bCutRapV0 = 1;
1440  bCutEtaV0 = 0;
1441  dNTauMax = 3.0;
1442  dRapMax = 0.5;
1443  }
1444 
1445  Double_t dCTauK0s = 2.6844; // [cm] c tau of K0S
1446  Double_t dCTauLambda = 7.89; // [cm] c tau of Lambda
1447 
1448  // Load PDG values of particle masses
1449  Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1450  Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1451 
1452  // PDG codes of used particles
1453  Int_t iPdgCodePion = 211;
1454  Int_t iPdgCodeProton = 2212;
1455  Int_t iPdgCodeK0s = 310;
1456  Int_t iPdgCodeLambda = 3122;
1457 
1458  // Jet selection: fdCutPtJetMin, fdCutPtTrackMin
1459  Double_t dJetEtaWindow = dEtaMax - fdRadiusJet; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1460  Double_t dCutJetAreaMin = 0.6 * TMath::Pi() * fdRadiusJet * fdRadiusJet; // minimum jet area
1461  Double_t dRadiusExcludeCone = 2 * fdRadiusJet; // radius of cones around jets excluded for V0 outside jets
1462  Bool_t bLeadingJetOnly = 0;
1463 
1464  if(bUseAliceCuts)
1465  {
1466  fdCutPtJetMin = 5;
1467  fdCutPtTrackMin = 5;
1468  dCutJetAreaMin = 0;
1469  bLeadingJetOnly = 0;
1470  }
1471 
1472 // Int_t iNJetAll = 0; // number of reconstructed jets in fBranchJet
1473 // iNJetAll = fBranchJet->GetEntriesFast(); // number of reconstructed jets
1474  TClonesArray* jetArray = 0; // object where the input jets are stored
1475  TClonesArray* jetArrayBg = 0; // object where the kt clusters are stored
1476  Int_t iNJet = 0; // number of reconstructed jets in the input
1477  TClonesArray* jetArraySel = new TClonesArray("AliAODJet", 0); // object where the selected jets are copied
1478  Int_t iNJetSel = 0; // number of selected reconstructed jets
1479 // iNJetSel = jetArraySel->GetEntriesFast(); // number of selected reconstructed jets
1480  TClonesArray* jetArrayPerp = new TClonesArray("AliAODJet", 0); // object where the perp. cones are stored
1481  Int_t iNJetPerp = 0; // number of perpendicular cones
1482 
1483  AliAODJet* jet = 0; // pointer to a jet
1484 // AliJetObject* objectJet = 0;
1485  AliAODJet* jetPerp = 0; // pointer to a perp. cone
1486  AliAODJet* jetRnd = 0; // pointer to a rand. cone
1487  AliAODJet* jetMed = 0; // pointer to a median cluster
1488  TVector3 vecJetMomentum; // 3D vector of jet momentum
1489 // TVector3 vecPerpConeMomentum; // 3D vector of perpendicular cone momentum
1490 // TVector3 vecRndConeMomentum; // 3D vector of random cone momentum
1491  Bool_t bJetEventGood = kTRUE; // indicator of good jet events
1492 
1493 // printf("iNJetAll, iNJetSel: %d %d\n",iNJetAll,iNJetSel);
1494 
1495  if(fbJetSelection) // analysis of V0s in jets is switched on
1496  {
1497  jetArray = (TClonesArray*)(fAODOut->FindListObject(fsJetBranchName.Data())); // find object with jets in the output AOD
1498  if(!jetArray)
1499  {
1500  if(fDebug > 0) printf("TaskV0sInJets: No array of name: %s\n", fsJetBranchName.Data());
1501  bJetEventGood = kFALSE;
1502  }
1503  if(bJetEventGood)
1504  iNJet = jetArray->GetEntriesFast();
1505  if(bJetEventGood && !iNJet) // check whether there are some jets
1506  {
1507  if(fDebug > 2) printf("TaskV0sInJets: No jets in array\n");
1508  bJetEventGood = kFALSE;
1509  }
1510  if(bJetEventGood)
1511  {
1512 // printf("TaskV0sInJets: Loading bg array of name: %s\n",fsJetBgBranchName.Data());
1513  jetArrayBg = (TClonesArray*)(fAODOut->FindListObject(fsJetBgBranchName.Data())); // find object with jets in the output AOD
1514  if(!jetArrayBg)
1515  {
1516  if(fDebug > 0) printf("TaskV0sInJets: No bg array of name: %s\n", fsJetBgBranchName.Data());
1517 // bJetEventGood = kFALSE;
1518  }
1519  }
1520  }
1521  else // no in-jet analysis
1522  bJetEventGood = kFALSE;
1523 
1524  // select good jets and copy them to another array
1525  if(bJetEventGood)
1526  {
1527  if(bLeadingJetOnly)
1528  iNJet = 1; // only leading jets
1529  if(fDebug > 5) printf("TaskV0sInJets: Jet selection for %d jets\n", iNJet);
1530  for(Int_t iJet = 0; iJet < iNJet; iJet++)
1531  {
1532  AliAODJet* jetSel = (AliAODJet*)jetArray->At(iJet); // load a jet in the list
1533  if(!jetSel)
1534  {
1535  if(fDebug > 0) printf("TaskV0sInJets: Cannot load jet %d\n", iJet);
1536  continue;
1537  }
1538  if(bPrintJetSelection)
1539  if(fDebug > 7) printf("jet: i = %d, pT = %f, eta = %f, phi = %f, pt lead tr = %f ", iJet, jetSel->Pt(), jetSel->Eta(), jetSel->Phi(), jetSel->GetPtLeading());
1540 // printf("TaskV0sInJets: Checking pt > %.2f for jet %d with pt %.2f\n",fdCutPtJetMin,iJet,jetSel->Pt());
1541  if(jetSel->Pt() < fdCutPtJetMin) // selection of high-pt jets
1542  {
1543  if(bPrintJetSelection)
1544  if(fDebug > 7) printf("rejected (pt)\n");
1545  continue;
1546  }
1547 // printf("TaskV0sInJets: Checking |eta| < %.2f for jet %d with |eta| %.2f\n",dJetEtaWindow,iJet,TMath::Abs(jetSel->Eta()));
1548  if(TMath::Abs(jetSel->Eta()) > dJetEtaWindow) // selection of jets in the chosen pseudorapidity range
1549  {
1550  if(bPrintJetSelection)
1551  if(fDebug > 7) printf("rejected (eta)\n");
1552  continue;
1553  }
1554  if(!bUseOldCuts)
1555  {
1556  if(jetSel->EffectiveAreaCharged() < dCutJetAreaMin)
1557  continue;
1558  }
1559  Int_t iNTracksInJet = 0;
1560  Double_t dPtLeadTrack = 0; // pt of the leading track
1561 // Int_t iLeadTrack = 0;
1562  iNTracksInJet = jetSel->GetRefTracks()->GetEntriesFast(); // number od tracks that constitute the jet
1563 // printf("TaskV0sInJets: Searching for leading track from %d tracks in jet %d\n",iNTracksInJet,iJet);
1564  if(fdCutPtTrackMin > 0) // a positive min leading track pt is set
1565  {
1566  for(Int_t j = 0; j < iNTracksInJet; j++) // find the track with the highest pt
1567  {
1568  AliAODTrack* track = (AliAODTrack*)jetSel->GetTrack(j); // is this the leading track?
1569  if(!track)
1570  continue;
1571 // printf("TaskV0sInJets: %d: %.2f\n",j,track->Pt());
1572  if(track->Pt() > dPtLeadTrack)
1573  {
1574  dPtLeadTrack = track->Pt();
1575 // iLeadTrack = j;
1576  }
1577  }
1578 // printf("Leading track pT: my: %f, ali: %f\n",dPtLeadTrack,jetSel->GetPtLeading());
1579 // printf("TaskV0sInJets: Checking leading track pt > %.2f for pt %.2f of track %d in jet %d\n",fdCutPtTrackMin,dPtLeadTrack,iLeadTrack,iJet);
1580  if(dPtLeadTrack < fdCutPtTrackMin) // selection of high-pt jet-track events
1581  {
1582  if(bPrintJetSelection)
1583  if(fDebug > 7) printf("rejected (track pt)\n");
1584  continue;
1585  }
1586  }
1587  if(bPrintJetSelection)
1588  if(fDebug > 7) printf("accepted\n");
1589  if(fDebug > 5) printf("TaskV0sInJets: Jet %d with pt %.2f passed selection\n", iJet, jetSel->Pt());
1590 /*
1591  if (fbTreeOutput)
1592  {
1593 // new ((*fBranchJet)[iNJetAll++]) AliAODJet(*((AliAODJet*)jetSel));
1594  objectJet = new ((*fBranchJet)[iNJetAll++]) AliJetObject(jetSel); // copy selected jet to the array
1595 // objectJet->SetPtLeadingTrack(dPtLeadTrack);
1596  objectJet->SetRadius(fdRadiusJet);
1597  objectJet = 0;
1598  }
1599 */
1600  TLorentzVector vecPerpPlus(*(jetSel->MomentumVector()));
1601  vecPerpPlus.RotateZ(TMath::Pi() / 2.); // rotate vector by 90 deg around z
1602  TLorentzVector vecPerpMinus(*(jetSel->MomentumVector()));
1603  vecPerpMinus.RotateZ(-TMath::Pi() / 2.); // rotate vector by -90 deg around z
1604 // AliAODJet jetTmp = AliAODJet(vecPerp);
1605  if(fDebug > 5) printf("TaskV0sInJets: Adding perp. cones number %d, %d\n", iNJetPerp, iNJetPerp + 1);
1606 // printf("TaskV0sInJets: Adding perp. cone number %d: pT %f, phi %f, eta %f, pT %f, phi %f, eta %f\n",iNJetSel,vecPerp.Pt(),vecPerp.Phi(),vecPerp.Eta(),jetTmp.Pt(),jetTmp.Phi(),jetTmp.Eta());
1607  new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpPlus); // write perp. cone to the array
1608  new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpMinus); // write perp. cone to the array
1609  if(fDebug > 5) printf("TaskV0sInJets: Adding jet number %d\n", iNJetSel);
1610 // printf("TaskV0sInJets: Adding jet number %d: pT %f, phi %f, eta %f\n",iNJetSel,jetSel->Pt(),jetSel->Phi(),jetSel->Eta());
1611  new((*jetArraySel)[iNJetSel++]) AliAODJet(*((AliAODJet*)jetSel)); // copy selected jet to the array
1612  }
1613  if(fDebug > 5) printf("TaskV0sInJets: Added jets: %d\n", iNJetSel);
1614  iNJetSel = jetArraySel->GetEntriesFast();
1615  if(fDebug > 2) printf("TaskV0sInJets: Selected jets in array: %d\n", iNJetSel);
1616  fh1NJetPerEvent[iCentIndex]->Fill(iNJetSel);
1617  // fill jet spectra
1618  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1619  {
1620  jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
1621  fh1PtJet[iCentIndex]->Fill(jet->Pt()); // pt spectrum of selected jets
1622  fh1EtaJet[iCentIndex]->Fill(jet->Eta()); // eta spectrum of selected jets
1623  fh2EtaPtJet[iCentIndex]->Fill(jet->Eta(), jet->Pt()); // eta-pT spectrum of selected jets
1624  fh1PhiJet[iCentIndex]->Fill(jet->Phi()); // phi spectrum of selected jets
1625  Double_t dAreaExcluded = TMath::Pi() * dRadiusExcludeCone * dRadiusExcludeCone; // area of the cone
1626  dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax - jet->Eta()); // positive eta overhang
1627  dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax + jet->Eta()); // negative eta overhang
1628  fh1AreaExcluded->Fill(iCentIndex, dAreaExcluded);
1629  }
1630  jet = 0;
1631  }
1632 
1633  if(bJetEventGood) // there should be some reconstructed jets
1634  {
1635  fh1EventCounterCut->Fill(4); // events with jet(s)
1636  fh1EventCounterCutCent[iCentIndex]->Fill(4); // events with jet(s)
1637  if(iNJetSel)
1638  {
1639  fh1EventCounterCut->Fill(5); // events with selected jets
1640  fh1EventCounterCutCent[iCentIndex]->Fill(5);
1641  }
1642  }
1643  if(iNJetSel)
1645  else
1647 
1648  if(iNJetSel)
1649  {
1650  jetRnd = GetRandomCone(jetArraySel, dJetEtaWindow, 2 * fdRadiusJet);
1651  if(jetRnd)
1652  {
1653  fh1NRndConeCent->Fill(iCentIndex);
1654  fh2EtaPhiRndCone[iCentIndex]->Fill(jetRnd->Eta(), jetRnd->Phi());
1655  }
1656  jetMed = GetMedianCluster(jetArrayBg, dJetEtaWindow);
1657  if(jetMed)
1658  {
1659  fh1NMedConeCent->Fill(iCentIndex);
1660  fh2EtaPhiMedCone[iCentIndex]->Fill(jetMed->Eta(), jetMed->Phi());
1661  }
1662  }
1663 
1664  // Loading primary vertex info
1665  AliAODVertex* primVtx = fAODIn->GetPrimaryVertex(); // get the primary vertex
1666  Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
1667  primVtx->GetXYZ(dPrimVtxPos);
1668  fh1VtxZ[iCentIndex]->Fill(dPrimVtxPos[2]);
1669  fh2VtxXY[iCentIndex]->Fill(dPrimVtxPos[0], dPrimVtxPos[1]);
1670 
1671  //===== Start of loop over V0 candidates =====
1672  if(fDebug > 2) printf("TaskV0sInJets: Start of V0 loop\n");
1673  for(Int_t iV0 = 0; iV0 < iNV0s; iV0++)
1674  {
1675  v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1676  if(!v0)
1677  continue;
1678 
1679  iNV0CandTot++;
1680 
1681  // Initialization of status indicators
1682  Bool_t bIsCandidateK0s = kTRUE; // candidate for K0s
1683  Bool_t bIsCandidateLambda = kTRUE; // candidate for Lambda
1684  Bool_t bIsCandidateALambda = kTRUE; // candidate for Lambda
1685  Bool_t bIsInPeakK0s = kFALSE; // candidate within the K0s mass peak
1686  Bool_t bIsInPeakLambda = kFALSE; // candidate within the Lambda mass peak
1687  Bool_t bIsInConeJet = kFALSE; // candidate within the jet cones
1688  Bool_t bIsInConePerp = kFALSE; // candidate within the perpendicular cone
1689  Bool_t bIsInConeRnd = kFALSE; // candidate within the random cone
1690  Bool_t bIsInConeMed = kFALSE; // candidate within the median-cluster cone
1691  Bool_t bIsOutsideCones = kFALSE; // candidate outside excluded cones
1692 
1693  // Invariant mass calculation
1694  dMassV0K0s = v0->MassK0Short();
1695  dMassV0Lambda = v0->MassLambda();
1696  dMassV0ALambda = v0->MassAntiLambda();
1697 
1698  Int_t iCutIndex = 0; // indicator of current selection step
1699  // 0
1700  // All V0 candidates
1701  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1702  iCutIndex++;
1703 
1704  // Skip candidates outside the histogram range
1705  if((dMassV0K0s < fgkdMassK0sMin) || (dMassV0K0s >= fgkdMassK0sMax))
1706  bIsCandidateK0s = kFALSE;
1707  if((dMassV0Lambda < fgkdMassLambdaMin) || (dMassV0Lambda >= fgkdMassLambdaMax))
1708  bIsCandidateLambda = kFALSE;
1709  if((dMassV0ALambda < fgkdMassLambdaMin) || (dMassV0ALambda >= fgkdMassLambdaMax))
1710  bIsCandidateALambda = kFALSE;
1711  if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1712  continue;
1713 
1714  Double_t dPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
1715  vecV0Momentum = TVector3(v0->Px(), v0->Py(), v0->Pz()); // set the vector of V0 momentum
1716 
1717  // Sigma of the mass peak window
1718  Double_t dMassPeakWindowK0s = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 0);
1719  Double_t dMassPeakWindowLambda = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 1);
1720 // Double_t dMassPeakWindowK0s = dNSigmaMassMax*MassPeakSigma(iCentIndex,dPtV0,0);
1721 // Double_t dMassPeakWindowLambda = dNSigmaMassMax*MassPeakSigma(iCentIndex,dPtV0,1);
1722 
1723  // Invariant mass peak selection
1724  if(TMath::Abs(dMassV0K0s - dMassPDGK0s) < dMassPeakWindowK0s)
1725  bIsInPeakK0s = kTRUE;
1726  if(TMath::Abs(dMassV0Lambda - dMassPDGLambda) < dMassPeakWindowLambda)
1727  bIsInPeakLambda = kTRUE;
1728 
1729  // Retrieving all relevant properties of the V0 candidate
1730  Bool_t bOnFlyStatus = v0->GetOnFlyStatus(); // online (on fly) reconstructed vs offline reconstructed
1731  const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1732  const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
1733  Double_t dPtDaughterPos = trackPos->Pt(); // transverse momentum of a daughter track
1734  Double_t dPtDaughterNeg = trackNeg->Pt();
1735  Double_t dNRowsPos = trackPos->GetTPCClusterInfo(2, 1); // crossed TPC pad rows of a daughter track
1736  Double_t dNRowsNeg = trackNeg->GetTPCClusterInfo(2, 1);
1737  Double_t dDCAToPrimVtxPos = TMath::Abs(v0->DcaPosToPrimVertex()); // dca of a daughter to the primary vertex
1738  Double_t dDCAToPrimVtxNeg = TMath::Abs(v0->DcaNegToPrimVertex());
1739  Double_t dDCADaughters = v0->DcaV0Daughters(); // dca between daughters
1740  Double_t dCPA = v0->CosPointingAngle(primVtx); // cosine of the pointing angle
1741  Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
1742 // Double_t dSecVtxPos[3] = {v0->DecayVertexV0X(),v0->DecayVertexV0Y(),v0->DecayVertexV0Z()}; // V0 vertex position
1743  v0->GetSecondaryVtx(dSecVtxPos);
1744  Double_t dRadiusDecay = TMath::Sqrt(dSecVtxPos[0] * dSecVtxPos[0] + dSecVtxPos[1] * dSecVtxPos[1]); // distance of the V0 vertex from the z-axis
1745  Double_t dEtaDaughterNeg = trackNeg->Eta(); // = v0->EtaProng(1), pseudorapidity of a daughter track
1746  Double_t dEtaDaughterPos = trackPos->Eta(); // = v0->EtaProng(0)
1747  Double_t dRapK0s = v0->RapK0Short(); // rapidity calculated for K0s assumption
1748  Double_t dRapLambda = v0->RapLambda(); // rapidity calculated for Lambda assumption
1749  Double_t dEtaV0 = v0->Eta(); // V0 pseudorapidity
1750 // Double_t dPhiV0 = v0->Phi(); // V0 pseudorapidity
1751  Double_t dDecayPath[3];
1752  for(Int_t iPos = 0; iPos < 3; iPos++)
1753  dDecayPath[iPos] = dSecVtxPos[iPos] - dPrimVtxPos[iPos]; // vector of the V0 path
1754  Double_t dDecLen = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]); // path length L
1755  Double_t dDecLen2D = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]); // transverse path length R
1756  Double_t dLOverP = dDecLen / v0->P(); // L/p
1757  Double_t dROverPt = dDecLen2D / dPtV0; // R/pT
1758  Double_t dMLOverPK0s = dMassPDGK0s * dLOverP; // m*L/p = c*(proper lifetime)
1759 // Double_t dMLOverPLambda = dMassPDGLambda*dLOverP; // m*L/p
1760  Double_t dMROverPtK0s = dMassPDGK0s * dROverPt; // m*R/pT
1761  Double_t dMROverPtLambda = dMassPDGLambda * dROverPt; // m*R/pT
1762  Double_t dNSigmaPosPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kPion)); // difference between measured and expected signal of the dE/dx in the TPC
1763  Double_t dNSigmaPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kProton));
1764  Double_t dNSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kPion));
1765  Double_t dNSigmaNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kProton));
1766  Double_t dAlpha = v0->AlphaV0(); // Armenteros-Podolanski alpha
1767  Double_t dPtArm = v0->PtArmV0(); // Armenteros-Podolanski pT
1768  AliAODVertex* prodVtxDaughterPos = (AliAODVertex*)(trackPos->GetProdVertex()); // production vertex of the positive daughter track
1769  Char_t cTypeVtxProdPos = prodVtxDaughterPos->GetType(); // type of the production vertex
1770  AliAODVertex* prodVtxDaughterNeg = (AliAODVertex*)(trackNeg->GetProdVertex()); // production vertex of the negative daughter track
1771  Char_t cTypeVtxProdNeg = prodVtxDaughterNeg->GetType(); // type of the production vertex
1772 
1773 // fh2Tau3DVs2D[0]->Fill(dPtV0, dLOverP / dROverPt);
1774 
1775  // QA histograms before cuts
1776  FillQAHistogramV0(primVtx, v0, 0, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
1777  // Cut vs mass histograms before cuts
1778  /*
1779  if(bIsCandidateK0s)
1780  {
1781  fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsPos);
1782  fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsNeg);
1783  fh2CutPtPosK0s[0]->Fill(dMassV0K0s, dPtDaughterPos);
1784  fh2CutPtNegK0s[0]->Fill(dMassV0K0s, dPtDaughterNeg);
1785  fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
1786  fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
1787  fh2CutDCAV0[0]->Fill(dMassV0K0s, dDCADaughters);
1788  fh2CutCos[0]->Fill(dMassV0K0s, dCPA);
1789  fh2CutR[0]->Fill(dMassV0K0s, dRadiusDecay);
1790  fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterPos);
1791  fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterNeg);
1792  fh2CutRapK0s[0]->Fill(dMassV0K0s, dRapK0s);
1793  fh2CutCTauK0s[0]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
1794  fh2CutPIDPosK0s[0]->Fill(dMassV0K0s, dNSigmaPosPion);
1795  fh2CutPIDNegK0s[0]->Fill(dMassV0K0s, dNSigmaNegPion);
1796  }
1797  if(bIsCandidateLambda)
1798  {
1799  fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsPos);
1800  fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsNeg);
1801  fh2CutPtPosLambda[0]->Fill(dMassV0Lambda, dPtDaughterPos);
1802  fh2CutPtNegLambda[0]->Fill(dMassV0Lambda, dPtDaughterNeg);
1803  fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterPos);
1804  fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterNeg);
1805  fh2CutRapLambda[0]->Fill(dMassV0Lambda, dRapLambda);
1806  fh2CutCTauLambda[0]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
1807  fh2CutPIDPosLambda[0]->Fill(dMassV0Lambda, dNSigmaPosProton);
1808  fh2CutPIDNegLambda[0]->Fill(dMassV0Lambda, dNSigmaNegPion);
1809  }
1810  */
1811 
1812  //===== Start of reconstruction cutting =====
1813 
1814  // 1
1815  // All V0 candidates
1816  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1817  iCutIndex++;
1818 
1819  // Start of global cuts
1820  // 2
1821  // Reconstruction method
1822  if(bPrintCuts) printf("Rec: Applying cut: Reconstruction method: on-the-fly? %s\n", (bOnFly ? "yes" : "no"));
1823  if(bOnFlyStatus != bOnFly)
1824  continue;
1825  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1826  iCutIndex++;
1827 
1828  // 3
1829  // Tracks TPC OK
1830  if(bPrintCuts) printf("Rec: Applying cut: Correct charge of daughters\n");
1831  if(!trackNeg || !trackPos)
1832  continue;
1833  if(trackNeg->Charge() == trackPos->Charge()) // daughters have different charge?
1834  continue;
1835  if(trackNeg->Charge() != -1) // daughters have expected charge?
1836  continue;
1837  if(trackPos->Charge() != 1) // daughters have expected charge?
1838  continue;
1839 
1840  if(bPrintCuts) printf("Rec: Applying cut: TPC refit: %d\n", iRefit);
1841  if(!trackNeg->IsOn(iRefit)) // TPC refit is ON?
1842  continue;
1843  if(bPrintCuts) printf("Rec: Applying cut: Type of production vertex of daughter: Not %d\n", AliAODVertex::kKink);
1844  if(cTypeVtxProdNeg == AliAODVertex::kKink) // kink daughter rejection
1845  continue;
1846  // Old cuts Start
1847  if(bUseOldCuts)
1848  {
1849  if(bPrintCuts) printf("Rec: Applying cut: Number of TPC rows: > %f\n", dNCrossedRowsTPCMin);
1850  if(dNRowsNeg < dNCrossedRowsTPCMin) // Crossed TPC padrows
1851  continue;
1852 // Int_t findable = trackNeg->GetTPCNclsF(); // Findable clusters
1853 // if (findable <= 0)
1854 // continue;
1855 // if (dNRowsNeg/findable < dCrossedRowsOverFindMin)
1856 // continue;
1857 // if (dNRowsNeg/findable > dCrossedRowsOverFindMax)
1858 // continue;
1859  }
1860  // Old cuts End
1861 
1862  if(!trackPos->IsOn(iRefit))
1863  continue;
1864  if(cTypeVtxProdPos == AliAODVertex::kKink) // kink daughter rejection
1865  continue;
1866  // Old cuts Start
1867  if(bUseOldCuts)
1868  {
1869  if(dNRowsPos < dNCrossedRowsTPCMin)
1870  continue;
1871 // findable = trackPos->GetTPCNclsF();
1872 // if (findable <= 0)
1873 // continue;
1874 // if (dNRowsPos/findable < dCrossedRowsOverFindMin)
1875 // continue;
1876 // if (dNRowsPos/findable > dCrossedRowsOverFindMax)
1877 // continue;
1878  }
1879  // Old cuts End
1880 
1881  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1882  iCutIndex++;
1883 
1884  // 4
1885  // Daughters: transverse momentum cut
1886  if(bUseOldCuts)
1887  {
1888  if(bPrintCuts) printf("Rec: Applying cut: Daughter pT: > %f\n", dPtDaughterMin);
1889  if((dPtDaughterNeg < dPtDaughterMin) || (dPtDaughterPos < dPtDaughterMin))
1890  continue;
1891  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1892  }
1893  iCutIndex++;
1894 
1895  // 5
1896  // Daughters: Impact parameter of daughters to prim vtx
1897  if(bPrintCuts) printf("Rec: Applying cut: Daughter DCA to prim vtx: > %f\n", dDCAToPrimVtxMin);
1898  if((dDCAToPrimVtxNeg < dDCAToPrimVtxMin) || (dDCAToPrimVtxPos < dDCAToPrimVtxMin))
1899  continue;
1900  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1901  iCutIndex++;
1902 
1903  // 6
1904  // Daughters: DCA
1905  if(bPrintCuts) printf("Rec: Applying cut: DCA between daughters: < %f\n", dDCADaughtersMax);
1906  if(dDCADaughters > dDCADaughtersMax)
1907  continue;
1908  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1909  iCutIndex++;
1910 
1911  // 7
1912  // V0: Cosine of the pointing angle
1913  if(bPrintCuts) printf("Rec: Applying cut: CPA: > %f\n", dCPAMin);
1914  if(dCPA < dCPAMin)
1915  continue;
1916  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1917  iCutIndex++;
1918 
1919  // 8
1920  // V0: Fiducial volume
1921  if(bPrintCuts) printf("Rec: Applying cut: Decay radius: > %f, < %f\n", dRadiusDecayMin, dRadiusDecayMax);
1922  if((dRadiusDecay < dRadiusDecayMin) || (dRadiusDecay > dRadiusDecayMax))
1923  continue;
1924  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1925  iCutIndex++;
1926 
1927  // 9
1928  // Daughters: pseudorapidity cut
1929  if(bCutEtaDaughter)
1930  {
1931  if(bPrintCuts) printf("Rec: Applying cut: Daughter |eta|: < %f\n", dEtaDaughterMax);
1932  if((TMath::Abs(dEtaDaughterNeg) > dEtaDaughterMax) || (TMath::Abs(dEtaDaughterPos) > dEtaDaughterMax))
1933  continue;
1934  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1935  }
1936  iCutIndex++;
1937  // End of global cuts
1938 
1939  // Start of particle-dependent cuts
1940  // 10
1941  // V0: rapidity cut & pseudorapidity cut
1942  if(bCutRapV0)
1943  {
1944  if(bPrintCuts) printf("Rec: Applying cut: V0 |y|: < %f\n", dRapMax);
1945  if(TMath::Abs(dRapK0s) > dRapMax)
1946  bIsCandidateK0s = kFALSE;
1947  if(TMath::Abs(dRapLambda) > dRapMax)
1948  {
1949  bIsCandidateLambda = kFALSE;
1950  bIsCandidateALambda = kFALSE;
1951  }
1952  }
1953  if(bCutEtaV0)
1954  {
1955  if(bPrintCuts) printf("Rec: Applying cut: V0 |eta|: < %f\n", dEtaMax);
1956  if(TMath::Abs(dEtaV0) > dEtaMax)
1957  {
1958  bIsCandidateK0s = kFALSE;
1959  bIsCandidateLambda = kFALSE;
1960  bIsCandidateALambda = kFALSE;
1961  }
1962  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1963  }
1964  iCutIndex++;
1965 
1966  // 11
1967  // Lifetime cut
1968  if(bCutTau)
1969  {
1970  if(bPrintCuts) printf("Rec: Applying cut: Proper lifetime: < %f\n", dNTauMax);
1971  if(dMROverPtK0s > dNTauMax * dCTauK0s)
1972  bIsCandidateK0s = kFALSE;
1973  if(dMROverPtLambda > dNTauMax * dCTauLambda)
1974  {
1975  bIsCandidateLambda = kFALSE;
1976  bIsCandidateALambda = kFALSE;
1977  }
1978  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1979  }
1980  iCutIndex++;
1981 
1982  // 12
1983  // Daughter PID
1984  if(bCutPid)
1985  {
1986  if(bUseOldCuts)
1987  {
1988  if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (both daughters): < %f\n", dNSigmadEdxMax);
1989  if(dNSigmaPosPion > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax) // pi+, pi-
1990  bIsCandidateK0s = kFALSE;
1991  if(dNSigmaPosProton > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax) // p+, pi-
1992  bIsCandidateLambda = kFALSE;
1993  if(dNSigmaNegProton > dNSigmadEdxMax || dNSigmaPosPion > dNSigmadEdxMax) // p-, pi+
1994  bIsCandidateALambda = kFALSE;
1995  }
1996  else
1997  {
1998  if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (proton below %f GeV/c): < %f\n", dPtProtonPIDMax, dNSigmadEdxMax);
1999  if((dPtDaughterPos < dPtProtonPIDMax) && (dNSigmaPosProton > dNSigmadEdxMax)) // p+
2000  bIsCandidateLambda = kFALSE;
2001  if((dPtDaughterNeg < dPtProtonPIDMax) && (dNSigmaNegProton > dNSigmadEdxMax)) // p-
2002  bIsCandidateALambda = kFALSE;
2003  }
2004  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2005  }
2006  iCutIndex++;
2007 
2008  Double_t valueCorrel[3] = {dMassV0K0s, dMassV0Lambda, dPtV0};
2009  if(bIsCandidateK0s && bIsCandidateLambda)
2010  fh3CCMassCorrelBoth->Fill(valueCorrel); // correlation of mass distribution of candidates selected as both K0s and Lambda
2011  if(bIsCandidateK0s && !bIsCandidateLambda)
2012  fh3CCMassCorrelKNotL->Fill(valueCorrel); // correlation of mass distribution of candidates selected as K0s and not Lambda
2013  if(!bIsCandidateK0s && bIsCandidateLambda)
2014  fh3CCMassCorrelLNotK->Fill(valueCorrel); // correlation of mass distribution of candidates selected as not K0s and Lambda
2015 
2016  // 13
2017  // Armenteros-Podolanski cut
2018  if(bCutArmPod)
2019  {
2020  if(bPrintCuts) printf("Rec: Applying cut: Armenteros-Podolanski (K0S): pT > %f * |alpha|\n", 0.2);
2021  if(dPtArm < TMath::Abs(0.2 * dAlpha))
2022  bIsCandidateK0s = kFALSE;
2023 // if(dPtArm < 0.025)
2024 // {
2025 // bIsCandidateLambda = kFALSE;
2026 // bIsCandidateALambda = kFALSE;
2027 // }
2028  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2029  }
2030  iCutIndex++;
2031 
2032  // Cross contamination
2033  if(bIsInPeakK0s)
2034  {
2035  if(bIsCandidateLambda) // Lambda candidates in K0s peak, excluded from Lambda candidates by CC cut
2036  fh2CCLambda->Fill(dMassV0Lambda, dPtV0);
2037  }
2038  if(bIsInPeakLambda)
2039  {
2040  if(bIsCandidateK0s) // K0s candidates in Lambda peak, excluded from K0s candidates by CC cut
2041  fh2CCK0s->Fill(dMassV0K0s, dPtV0);
2042  }
2043 // if (bCutCross)
2044 // {
2045 // if (bIsInPeakK0s)
2046 // bIsCandidateLambda = kFALSE;
2047 // if (bIsInPeakLambda)
2048 // bIsCandidateK0s = kFALSE;
2049 // FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2050 // }
2051 // iCutIndex++;
2052 
2053  // End of particle-dependent cuts
2054 
2055  //===== End of reconstruction cutting =====
2056 
2057  if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
2058  continue;
2059 
2060 /*
2061  if(fDebug>5) printf("TaskV0sInJets: Adding selected V0 to branch\n");
2062  // Add selected candidates to the output tree branch
2063  if ((bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda) && fbTreeOutput)
2064  {
2065  objectV0 = new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliV0Object(v0,primVtx);
2066 // new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliAODv0(*((AliAODv0*)v0));
2067  objectV0->SetIsCandidateK0S(bIsCandidateK0s);
2068  objectV0->SetIsCandidateLambda(bIsCandidateLambda);
2069  objectV0->SetIsCandidateALambda(bIsCandidateALambda);
2070  objectV0->SetNSigmaPosProton(dNSigmaPosProton);
2071  objectV0->SetNSigmaNegProton(dNSigmaNegProton);
2072  }
2073 */
2074 
2075  // Selection of V0s in jet cones, perpendicular cones, random cones, outside cones
2076  if(bJetEventGood && iNJetSel && (bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda))
2077  {
2078  // Selection of V0s in jet cones
2079  if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in %d jet cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetSel);
2080  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2081  {
2082  jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2083  vecJetMomentum = TVector3(jet->Px(), jet->Py(), jet->Pz()); // set the vector of jet momentum
2084  if(fDebug > 5) printf("TaskV0sInJets: Checking if V0 %d %d in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2085  if(IsParticleInCone(v0, jet, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2086  {
2087  if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2088  bIsInConeJet = kTRUE;
2089  break;
2090  }
2091  }
2092  // Selection of V0s in perp. cones
2093  if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in %d perp. cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetSel);
2094  for(Int_t iJet = 0; iJet < iNJetPerp; iJet++)
2095  {
2096  jetPerp = (AliAODJet*)jetArrayPerp->At(iJet); // load a jet in the list
2097  if(fDebug > 5) printf("TaskV0sInJets: Checking if V0 %d %d in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2098  if(IsParticleInCone(v0, jetPerp, fdRadiusJet)) // V0 in perp. cone
2099  {
2100  if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2101  bIsInConePerp = kTRUE;
2102  break;
2103  }
2104  }
2105  // Selection of V0s in random cones
2106  if(jetRnd)
2107  {
2108  if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2109  if(IsParticleInCone(v0, jetRnd, fdRadiusJet)) // V0 in rnd. cone?
2110  {
2111  if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2112  bIsInConeRnd = kTRUE;
2113  }
2114  }
2115  // Selection of V0s in median-cluster cones
2116  if(jetMed)
2117  {
2118  if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2119  if(IsParticleInCone(v0, jetMed, fdRadiusJet)) // V0 in med. cone?
2120  {
2121  if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2122  bIsInConeMed = kTRUE;
2123  }
2124  }
2125  // Selection of V0s outside jet cones
2126  if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2127  if(!OverlapWithJets(jetArraySel, v0, dRadiusExcludeCone)) // V0 oustide jet cones
2128  {
2129  if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2130  bIsOutsideCones = kTRUE;
2131  }
2132  }
2133 
2134  // QA histograms after cuts
2135  FillQAHistogramV0(primVtx, v0, 1, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
2136  // Cut vs mass histograms after cuts
2137  /*
2138  if(bIsCandidateK0s)
2139  {
2140  fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsPos);
2141  fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsNeg);
2142  fh2CutPtPosK0s[1]->Fill(dMassV0K0s, dPtDaughterPos);
2143  fh2CutPtNegK0s[1]->Fill(dMassV0K0s, dPtDaughterNeg);
2144  fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
2145  fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
2146  fh2CutDCAV0[1]->Fill(dMassV0K0s, dDCADaughters);
2147  fh2CutCos[1]->Fill(dMassV0K0s, dCPA);
2148  fh2CutR[1]->Fill(dMassV0K0s, dRadiusDecay);
2149  fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterPos);
2150  fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterNeg);
2151  fh2CutRapK0s[1]->Fill(dMassV0K0s, dRapK0s);
2152  fh2CutCTauK0s[1]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
2153  fh2CutPIDPosK0s[1]->Fill(dMassV0K0s, dNSigmaPosPion);
2154  fh2CutPIDNegK0s[1]->Fill(dMassV0K0s, dNSigmaNegPion);
2155  fh1DeltaZK0s[iCentIndex]->Fill(dDecayPath[2]);
2156  }
2157  if(bIsCandidateLambda)
2158  {
2159  fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsPos);
2160  fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsNeg);
2161  fh2CutPtPosLambda[1]->Fill(dMassV0Lambda, dPtDaughterPos);
2162  fh2CutPtNegLambda[1]->Fill(dMassV0Lambda, dPtDaughterNeg);
2163  fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterPos);
2164  fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterNeg);
2165  fh2CutRapLambda[1]->Fill(dMassV0Lambda, dRapLambda);
2166  fh2CutCTauLambda[1]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
2167  fh2CutPIDPosLambda[1]->Fill(dMassV0Lambda, dNSigmaPosProton);
2168  fh2CutPIDNegLambda[1]->Fill(dMassV0Lambda, dNSigmaNegPion);
2169  fh1DeltaZLambda[iCentIndex]->Fill(dDecayPath[2]);
2170  }
2171  */
2172 
2173  //===== Start of filling V0 spectra =====
2174 
2175  Double_t dAngle = TMath::Pi(); // angle between V0 momentum and jet momentum
2176  if(bIsInConeJet)
2177  {
2178  dAngle = vecV0Momentum.Angle(vecJetMomentum);
2179  }
2180 
2181  // iCutIndex = 14
2182  if(bIsCandidateK0s)
2183  {
2184  // 14 K0s candidates after cuts
2185 // printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0K0s,dPtV0,dEtaV0,dPhiV0);
2186  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex, iCentIndex);
2187  Double_t valueKIncl[3] = {dMassV0K0s, dPtV0, dEtaV0};
2188  fhnV0InclusiveK0s[iCentIndex]->Fill(valueKIncl);
2189  fh1V0InvMassK0sCent[iCentIndex]->Fill(dMassV0K0s);
2190 
2191  fh1QACTau2D[1]->Fill(dMROverPtK0s / dCTauK0s);
2192  fh1QACTau3D[1]->Fill(dMLOverPK0s / dCTauK0s);
2193 // fh2Tau3DVs2D[1]->Fill(dPtV0, dLOverP / dROverPt);
2194 
2195  if(iNJetSel)
2196  {
2197  // 15 K0s in jet events
2198  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 1, iCentIndex);
2199  }
2200  if(bIsInConeJet)
2201  {
2202  // 16 K0s in jets
2203  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 2, iCentIndex);
2204  Double_t valueKInJC[4] = {dMassV0K0s, dPtV0, dEtaV0, jet->Pt()};
2205  fhnV0InJetK0s[iCentIndex]->Fill(valueKInJC);
2206  fh2V0PtJetAngleK0s[iCentIndex]->Fill(jet->Pt(), dAngle);
2207  }
2208  if(bIsOutsideCones)
2209  {
2210  Double_t valueKOutJC[3] = {dMassV0K0s, dPtV0, dEtaV0};
2211  fhnV0OutJetK0s[iCentIndex]->Fill(valueKOutJC);
2212  }
2213  if(bIsInConePerp)
2214  {
2215  Double_t valueKInPC[4] = {dMassV0K0s, dPtV0, dEtaV0, jetPerp->Pt()};
2216  fhnV0InPerpK0s[iCentIndex]->Fill(valueKInPC);
2217  }
2218  if(bIsInConeRnd)
2219  {
2220  Double_t valueKInRnd[3] = {dMassV0K0s, dPtV0, dEtaV0};
2221  fhnV0InRndK0s[iCentIndex]->Fill(valueKInRnd);
2222  }
2223  if(bIsInConeMed)
2224  {
2225  Double_t valueKInMed[3] = {dMassV0K0s, dPtV0, dEtaV0};
2226  fhnV0InMedK0s[iCentIndex]->Fill(valueKInMed);
2227  }
2228  if(!iNJetSel)
2229  {
2230  Double_t valueKNoJet[3] = {dMassV0K0s, dPtV0, dEtaV0};
2231  fhnV0NoJetK0s[iCentIndex]->Fill(valueKNoJet);
2232  }
2233  iNV0CandK0s++;
2234  }
2235  if(bIsCandidateLambda)
2236  {
2237  // 14 Lambda candidates after cuts
2238 // printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0Lambda,dPtV0,dEtaV0,dPhiV0);
2239  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex, iCentIndex);
2240  Double_t valueLIncl[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2241  fhnV0InclusiveLambda[iCentIndex]->Fill(valueLIncl);
2242  fh1V0InvMassLambdaCent[iCentIndex]->Fill(dMassV0Lambda);
2243  if(iNJetSel)
2244  {
2245  // 15 Lambda in jet events
2246  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 1, iCentIndex);
2247  }
2248  if(bIsInConeJet)
2249  {
2250  // 16 Lambda in jets
2251  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 2, iCentIndex);
2252  Double_t valueLInJC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jet->Pt()};
2253  fhnV0InJetLambda[iCentIndex]->Fill(valueLInJC);
2254  fh2V0PtJetAngleLambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2255  }
2256  if(bIsOutsideCones)
2257  {
2258  Double_t valueLOutJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2259  fhnV0OutJetLambda[iCentIndex]->Fill(valueLOutJet);
2260  }
2261  if(bIsInConePerp)
2262  {
2263  Double_t valueLInPC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jetPerp->Pt()};
2264  fhnV0InPerpLambda[iCentIndex]->Fill(valueLInPC);
2265  }
2266  if(bIsInConeRnd)
2267  {
2268  Double_t valueLInRnd[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2269  fhnV0InRndLambda[iCentIndex]->Fill(valueLInRnd);
2270  }
2271  if(bIsInConeMed)
2272  {
2273  Double_t valueLInMed[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2274  fhnV0InMedLambda[iCentIndex]->Fill(valueLInMed);
2275  }
2276  if(!iNJetSel)
2277  {
2278  Double_t valueLNoJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2279  fhnV0NoJetLambda[iCentIndex]->Fill(valueLNoJet);
2280  }
2281  iNV0CandLambda++;
2282  }
2283  if(bIsCandidateALambda)
2284  {
2285  // 14 ALambda candidates after cuts
2286 // printf("AL: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0ALambda,dPtV0,dEtaV0,dPhiV0);
2287  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex, iCentIndex);
2288  Double_t valueALIncl[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2289  fhnV0InclusiveALambda[iCentIndex]->Fill(valueALIncl);
2290  fh1V0InvMassALambdaCent[iCentIndex]->Fill(dMassV0ALambda);
2291  if(iNJetSel)
2292  {
2293  // 15 ALambda in jet events
2294  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 1, iCentIndex);
2295  }
2296  if(bIsInConeJet)
2297  {
2298  // 16 ALambda in jets
2299  FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 2, iCentIndex);
2300  Double_t valueLInJC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jet->Pt()};
2301  fhnV0InJetALambda[iCentIndex]->Fill(valueLInJC);
2302  fh2V0PtJetAngleALambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2303  }
2304  if(bIsOutsideCones)
2305  {
2306  Double_t valueALOutJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2307  fhnV0OutJetALambda[iCentIndex]->Fill(valueALOutJet);
2308  }
2309  if(bIsInConePerp)
2310  {
2311  Double_t valueLInPC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jetPerp->Pt()};
2312  fhnV0InPerpALambda[iCentIndex]->Fill(valueLInPC);
2313  }
2314  if(bIsInConeRnd)
2315  {
2316  Double_t valueALInRnd[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2317  fhnV0InRndALambda[iCentIndex]->Fill(valueALInRnd);
2318  }
2319  if(bIsInConeMed)
2320  {
2321  Double_t valueALInMed[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2322  fhnV0InMedALambda[iCentIndex]->Fill(valueALInMed);
2323  }
2324  if(!iNJetSel)
2325  {
2326  Double_t valueALNoJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2327  fhnV0NoJetALambda[iCentIndex]->Fill(valueALNoJet);
2328  }
2329  iNV0CandALambda++;
2330  }
2331  //===== End of filling V0 spectra =====
2332 
2333 
2334  //===== Association of reconstructed V0 candidates with MC particles =====
2335  if(fbMCAnalysis)
2336  {
2337  // Associate selected candidates only
2338 // if ( !(bIsCandidateK0s && bIsInPeakK0s) && !(bIsCandidateLambda && bIsInPeakLambda) ) // signal candidates
2339  if(!(bIsCandidateK0s) && !(bIsCandidateLambda) && !(bIsCandidateALambda)) // chosen candidates with any mass
2340  continue;
2341 
2342  // Get MC labels of reconstructed daughter tracks
2343  Int_t iLabelPos = TMath::Abs(trackPos->GetLabel());
2344  Int_t iLabelNeg = TMath::Abs(trackNeg->GetLabel());
2345 
2346  // Make sure MC daughters are in the array range
2347  if((iLabelNeg < 0) || (iLabelNeg >= iNTracksMC) || (iLabelPos < 0) || (iLabelPos >= iNTracksMC))
2348  continue;
2349 
2350  // Get MC particles corresponding to reconstructed daughter tracks
2351  AliAODMCParticle* particleMCDaughterNeg = (AliAODMCParticle*)arrayMC->At(iLabelNeg);
2352  AliAODMCParticle* particleMCDaughterPos = (AliAODMCParticle*)arrayMC->At(iLabelPos);
2353  if(!particleMCDaughterNeg || !particleMCDaughterPos)
2354  continue;
2355 
2356  // Make sure MC daughter particles are not physical primary
2357  if((particleMCDaughterNeg->IsPhysicalPrimary()) || (particleMCDaughterPos->IsPhysicalPrimary()))
2358  continue;
2359 
2360  // Get identities of MC daughter particles
2361  Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
2362  Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
2363 
2364  // Get index of the mother particle for each MC daughter particle
2365  Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
2366  Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
2367 
2368  if((iIndexMotherNeg < 0) || (iIndexMotherNeg >= iNTracksMC) || (iIndexMotherPos < 0) || (iIndexMotherPos >= iNTracksMC))
2369  continue;
2370 
2371  // Check whether MC daughter particles have the same mother
2372  if(iIndexMotherNeg != iIndexMotherPos)
2373  continue;
2374 
2375  // Get the MC mother particle of both MC daughter particles
2376  AliAODMCParticle* particleMCMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherPos);
2377  if(!particleMCMother)
2378  continue;
2379 
2380  // Get identity of the MC mother particle
2381  Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
2382 
2383  // Skip not interesting particles
2384  if((iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda))
2385  continue;
2386 
2387  // Check identity of the MC mother particle and the decay channel
2388  // Is MC mother particle K0S?
2389  Bool_t bV0MCIsK0s = ((iPdgCodeMother == iPdgCodeK0s) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2390  // Is MC mother particle Lambda?
2391  Bool_t bV0MCIsLambda = ((iPdgCodeMother == +iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodeProton) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2392  // Is MC mother particle anti-Lambda?
2393  Bool_t bV0MCIsALambda = ((iPdgCodeMother == -iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodeProton));
2394 
2395  Double_t dPtV0Gen = particleMCMother->Pt();
2396 // Double_t dRapV0MC = particleMCMother->Y();
2397  Double_t dEtaV0Gen = particleMCMother->Eta();
2398 // Double_t dPhiV0Gen = particleMCMother->Phi();
2399 
2400  // Is MC mother particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2401 // Bool_t bV0MCIsPrimary = particleMCMother->IsPhysicalPrimary();
2402  // Get the MC mother particle of the MC mother particle
2403  Int_t iIndexMotherOfMother = particleMCMother->GetMother();
2404  AliAODMCParticle* particleMCMotherOfMother = 0;
2405  if(iIndexMotherOfMother >= 0)
2406  particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2407  // Get identity of the MC mother particle of the MC mother particle if it exists
2408  Int_t iPdgCodeMotherOfMother = 0;
2409  if(particleMCMotherOfMother)
2410  iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2411  // 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*-)
2412 // Bool_t bV0MCComesFromSigma = kFALSE; // Is MC mother particle daughter of a Sigma?
2413 // if ( (particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && ( (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) ) )
2414 // bV0MCComesFromSigma = kTRUE;
2415  // Should MC mother particle be considered as primary when it is Lambda?
2416 // Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2417  // Check if the MC mother particle of the MC mother particle is a Xi (3322 - Xi0, 3312 - Xi-)
2418  Bool_t bV0MCComesFromXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == 3322) || (iPdgCodeMotherOfMother == 3312))); // Is MC mother particle daughter of a Xi?
2419  Bool_t bV0MCComesFromAXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == -3322) || (iPdgCodeMotherOfMother == -3312))); // Is MC mother particle daughter of a anti-Xi?
2420 
2421  // Get the distance between production point of the MC mother particle and the primary vertex
2422  Double_t dx = dPrimVtxMCX - particleMCMother->Xv();
2423  Double_t dy = dPrimVtxMCY - particleMCMother->Yv();
2424  Double_t dz = dPrimVtxMCZ - particleMCMother->Zv();
2425  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2426  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2427 
2428  // phi, eta resolution for V0-reconstruction
2429 // Double_t dResolutionV0Eta = particleMCMother->Eta()-v0->Eta();
2430 // Double_t dResolutionV0Phi = particleMCMother->Phi()-v0->Phi();
2431 
2432 /*
2433  if (fbTreeOutput)
2434  {
2435  objectV0->SetPtTrue(dPtV0Gen);
2436  objectV0->SetEtaTrue(dEtaV0Gen);
2437  objectV0->SetPhiTrue(dPhiV0Gen);
2438  objectV0->SetPDGCode(iPdgCodeMother);
2439  objectV0->SetPDGCodeMother(iPdgCodeMotherOfMother);
2440  }
2441 */
2442 
2443  // K0s
2444 // if (bIsCandidateK0s && bIsInPeakK0s) // selected candidates in peak
2445  if(bIsCandidateK0s) // selected candidates with any mass
2446  {
2447 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2448  if(bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2449  {
2450 // if (fbTreeOutput)
2451 // objectV0->SetOrigin(1);
2452  fh2V0K0sPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0K0s);
2453  Double_t valueEtaK[3] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen};
2454  fh3V0K0sEtaPtMassMCRec[iCentIndex]->Fill(valueEtaK);
2455 
2456  Double_t valueEtaDKNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2457  fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKNeg);
2458  Double_t valueEtaDKPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2459  fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKPos);
2460 
2461  fh2V0K0sMCResolMPt[iCentIndex]->Fill(dMassV0K0s - dMassPDGK0s, dPtV0);
2462  fh2V0K0sMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2463  if(bIsInConeJet) // true V0 associated to a candidate in jet
2464  {
2465  Double_t valueKInJCMC[4] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2466  fh3V0K0sInJetPtMassMCRec[iCentIndex]->Fill(valueKInJCMC);
2467  Double_t valueEtaKIn[5] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2468  fh4V0K0sInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaKIn);
2469 
2470  Double_t valueEtaDKJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2471  fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCNeg);
2472  Double_t valueEtaDKJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2473  fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCPos);
2474  }
2475  }
2476  if(bV0MCIsK0s && !bV0MCIsPrimaryDist) // not primary K0s
2477  {
2478 // if (fbTreeOutput)
2479 // objectV0->SetOrigin(-1);
2480  fh1V0K0sPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2481  }
2482  }
2483  // Lambda
2484 // if (bIsCandidateLambda && bIsInPeakLambda) // selected candidates in peak
2485  if(bIsCandidateLambda) // selected candidates with any mass
2486  {
2487 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2488  if(bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2489  {
2490 // if (fbTreeOutput)
2491 // objectV0->SetOrigin(1);
2492  fh2V0LambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0Lambda);
2493  Double_t valueEtaL[3] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen};
2494  fh3V0LambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaL);
2495 
2496  Double_t valueEtaDLNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2497  fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLNeg);
2498  Double_t valueEtaDLPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2499  fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLPos);
2500 
2501  fh2V0LambdaMCResolMPt[iCentIndex]->Fill(dMassV0Lambda - dMassPDGLambda, dPtV0);
2502  fh2V0LambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2503  if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2504  {
2505  Double_t valueLInJCMC[4] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2506  fh3V0LambdaInJetPtMassMCRec[iCentIndex]->Fill(valueLInJCMC);
2507  Double_t valueEtaLIn[5] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2508  fh4V0LambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaLIn);
2509 
2510  Double_t valueEtaDLJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2511  fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCNeg);
2512  Double_t valueEtaDLJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2513  fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCPos);
2514  }
2515  }
2516  // Fill the feed-down histograms
2517  if(bV0MCIsLambda && bV0MCComesFromXi)
2518  {
2519 // if (fbTreeOutput)
2520 // objectV0->SetOrigin(2);
2521  Double_t valueFDLIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2522  fhnV0LambdaInclMCFD[iCentIndex]->Fill(valueFDLIncl);
2523  if(bIsInConeRnd)
2524  {
2525  fhnV0LambdaBulkMCFD[iCentIndex]->Fill(valueFDLIncl);
2526  }
2527  if(bIsInConeJet)
2528  {
2529  Double_t valueFDLInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2530  fhnV0LambdaInJetsMCFD[iCentIndex]->Fill(valueFDLInJets);
2531  }
2532  }
2533  if(bV0MCIsLambda && !bV0MCIsPrimaryDist && !bV0MCComesFromXi) // not primary Lambda
2534  {
2535 // if (fbTreeOutput)
2536 // objectV0->SetOrigin(-1);
2537  fh1V0LambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2538  }
2539  }
2540  // anti-Lambda
2541 // if (bIsCandidateALambda && bIsInPeakALambda) // selected candidates in peak
2542  if(bIsCandidateALambda) // selected candidates with any mass
2543  {
2544 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2545  if(bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2546  {
2547 // if (fbTreeOutput)
2548 // objectV0->SetOrigin(1);
2549  fh2V0ALambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0ALambda);
2550  Double_t valueEtaAL[3] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen};
2551  fh3V0ALambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaAL);
2552 
2553  Double_t valueEtaDALNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2554  fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALNeg);
2555  Double_t valueEtaDALPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2556  fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALPos);
2557 
2558  fh2V0ALambdaMCResolMPt[iCentIndex]->Fill(dMassV0ALambda - dMassPDGLambda, dPtV0);
2559  fh2V0ALambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2560  if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2561  {
2562  Double_t valueALInJCMC[4] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2563  fh3V0ALambdaInJetPtMassMCRec[iCentIndex]->Fill(valueALInJCMC);
2564  Double_t valueEtaALIn[5] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2565  fh4V0ALambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaALIn);
2566 
2567  Double_t valueEtaDALJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2568  fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCNeg);
2569  Double_t valueEtaDALJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2570  fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCPos);
2571  }
2572  }
2573  // Fill the feed-down histograms
2574  if(bV0MCIsALambda && bV0MCComesFromAXi)
2575  {
2576 // if (fbTreeOutput)
2577 // objectV0->SetOrigin(2);
2578  Double_t valueFDALIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2579  fhnV0ALambdaInclMCFD[iCentIndex]->Fill(valueFDALIncl);
2580  if(bIsInConeRnd)
2581  {
2582  fhnV0ALambdaBulkMCFD[iCentIndex]->Fill(valueFDALIncl);
2583  }
2584  if(bIsInConeJet)
2585  {
2586  Double_t valueFDALInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2587  fhnV0ALambdaInJetsMCFD[iCentIndex]->Fill(valueFDALInJets);
2588  }
2589  }
2590  if(bV0MCIsALambda && !bV0MCIsPrimaryDist && !bV0MCComesFromAXi) // not primary anti-Lambda
2591  {
2592 // if (fbTreeOutput)
2593 // objectV0->SetOrigin(-1);
2594  fh1V0ALambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2595  }
2596  }
2597  }
2598  //===== End Association of reconstructed V0 candidates with MC particles =====
2599  }
2600  //===== End of V0 loop =====
2601 
2602  fh1V0CandPerEvent->Fill(iNV0CandTot);
2603  fh1V0CandPerEventCentK0s[iCentIndex]->Fill(iNV0CandK0s);
2604  fh1V0CandPerEventCentLambda[iCentIndex]->Fill(iNV0CandLambda);
2605  fh1V0CandPerEventCentALambda[iCentIndex]->Fill(iNV0CandALambda);
2606 
2607  if(fDebug > 2) printf("TaskV0sInJets: End of V0 loop\n");
2608 
2609  // Spectra of generated particles
2610  if(fbMCAnalysis)
2611  {
2612  for(Int_t iPartMC = 0; iPartMC < iNTracksMC; iPartMC++)
2613  {
2614  // Get MC particle
2615  AliAODMCParticle* particleMC = (AliAODMCParticle*)arrayMC->At(iPartMC);
2616  if(!particleMC)
2617  continue;
2618 
2619  // Get identity of MC particle
2620  Int_t iPdgCodeParticleMC = particleMC->GetPdgCode();
2621  // Fill Xi spectrum (3322 - Xi0, 3312 - Xi-)
2622 // if ( (iPdgCodeParticleMC==3322) || (iPdgCodeParticleMC==3312) )
2623  if((iPdgCodeParticleMC == 3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2624  {
2625 // if (fbTreeOutput)
2626 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2627  fh1V0XiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2628  }
2629  if((iPdgCodeParticleMC == -3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2630  {
2631 // if (fbTreeOutput)
2632 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2633  fh1V0AXiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2634  }
2635  // Skip not interesting particles
2636  if((iPdgCodeParticleMC != iPdgCodeK0s) && (TMath::Abs(iPdgCodeParticleMC) != iPdgCodeLambda))
2637  continue;
2638 
2639  // Check identity of the MC V0 particle
2640  // Is MC V0 particle K0S?
2641  Bool_t bV0MCIsK0s = (iPdgCodeParticleMC == iPdgCodeK0s);
2642  // Is MC V0 particle Lambda?
2643  Bool_t bV0MCIsLambda = (iPdgCodeParticleMC == +iPdgCodeLambda);
2644  // Is MC V0 particle anti-Lambda?
2645  Bool_t bV0MCIsALambda = (iPdgCodeParticleMC == -iPdgCodeLambda);
2646 
2647  Double_t dPtV0Gen = particleMC->Pt();
2648  Double_t dRapV0Gen = particleMC->Y();
2649  Double_t dEtaV0Gen = particleMC->Eta();
2650 
2651  // V0 rapidity cut
2652  if(bCutRapV0)
2653  {
2654  if(bPrintCuts) printf("Gen: Applying cut: V0 |y|: < %f\n", dRapMax);
2655  if((TMath::Abs(dRapV0Gen) > dRapMax))
2656  continue;
2657  }
2658  // V0 pseudorapidity cut
2659  if(bCutEtaV0)
2660  {
2661  if(bPrintCuts) printf("Gen: Applying cut: V0 |eta|: < %f\n", dEtaMax);
2662  if((TMath::Abs(dEtaV0Gen) > dEtaMax))
2663  continue;
2664  }
2665 /*
2666  // Is MC V0 particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2667  Bool_t bV0MCIsPrimary = particleMC->IsPhysicalPrimary();
2668 
2669  // Get the MC mother particle of the MC V0 particle
2670  Int_t iIndexMotherOfMother = particleMC->GetMother();
2671  AliAODMCParticle* particleMCMotherOfMother = 0;
2672  if (iIndexMotherOfMother >= 0)
2673  particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2674  // Get identity of the MC mother particle of the MC V0 particle if it exists
2675  Int_t iPdgCodeMotherOfMother = 0;
2676  if (particleMCMotherOfMother)
2677  iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2678  // Check if the MC mother particle is a physical primary Sigma
2679  Bool_t bV0MCComesFromSigma = kFALSE;
2680  if ((particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) )
2681  bV0MCComesFromSigma = kTRUE;
2682  // Should the MC V0 particle be considered as primary when it is Lambda?
2683  Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2684 */
2685  // Reject non primary particles
2686 // if (!bV0MCIsPrimaryLambda)
2687 // continue;
2688 
2689  // Get the distance between the production point of the MC V0 particle and the primary vertex
2690  Double_t dx = dPrimVtxMCX - particleMC->Xv();
2691  Double_t dy = dPrimVtxMCY - particleMC->Yv();
2692  Double_t dz = dPrimVtxMCZ - particleMC->Zv();
2693  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2694  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2695 
2696  // Check whether the MC V0 particle is in a MC jet
2697  AliAODJet* jetMC = 0;
2698  Bool_t bIsMCV0InJet = kFALSE;
2699  if(iNJetSel)
2700  {
2701  if(fDebug > 5) printf("TaskV0sInJets: Searching for gen V0 in %d MC jets\n", iNJetSel);
2702  for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2703  {
2704  jetMC = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2705  if(fDebug > 5) printf("TaskV0sInJets: Checking if gen V0 in MC jet %d\n", iJet);
2706  if(IsParticleInCone(particleMC, jetMC, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2707  {
2708  if(fDebug > 5) printf("TaskV0sInJets: gen V0 found in MC jet %d\n", iJet);
2709  bIsMCV0InJet = kTRUE;
2710  break;
2711  }
2712  }
2713  }
2714 
2715  // Select only primary-like MC V0 particles
2716  // K0s
2717 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2718  if(bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2719  {
2720 // if (fbTreeOutput)
2721 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2722  fh1V0K0sPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2723  fh2V0K0sEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2724  if(bIsMCV0InJet)
2725  {
2726  fh2V0K0sInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2727  Double_t valueEtaKInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2728  fh3V0K0sInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaKInGen);
2729  }
2730  }
2731  // Lambda
2732 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2733  if(bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2734  {
2735 // if (fbTreeOutput)
2736 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2737  fh1V0LambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2738  fh2V0LambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2739  if(bIsMCV0InJet)
2740  {
2741  fh2V0LambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2742  Double_t valueEtaLInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2743  fh3V0LambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaLInGen);
2744  }
2745  }
2746  // anti-Lambda
2747 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2748  if(bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2749  {
2750 // if (fbTreeOutput)
2751 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2752  fh1V0ALambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2753  fh2V0ALambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2754  if(bIsMCV0InJet)
2755  {
2756  fh2V0ALambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2757  Double_t valueEtaALInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2758  fh3V0ALambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaALInGen);
2759  }
2760  }
2761  }
2762  }
2763 
2764 // if (fbTreeOutput)
2765 // ftreeOut->Fill();
2766 
2767  jetArraySel->Delete();
2768  delete jetArraySel;
2769  jetArrayPerp->Delete();
2770  delete jetArrayPerp;
2771  if(jetRnd)
2772  delete jetRnd;
2773  jetRnd = 0;
2774 
2775  PostData(1, fOutputListStd);
2776  PostData(2, fOutputListQA);
2777  PostData(3, fOutputListCuts);
2778  PostData(4, fOutputListMC);
2779 // if (fbTreeOutput)
2780 // PostData(5,ftreeOut);
2781 // if(fDebug>5) printf("TaskV0sInJets: UserExec: End\n");
2782 }
2783 
2784 void AliAnalysisTaskV0sInJets::FillQAHistogramV0(AliAODVertex* vtx, const AliAODv0* vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda)
2785 {
2786  if(!IsCandK0s && !IsCandLambda)
2787  return;
2788 
2789 // Double_t dMassK0s = vZero->MassK0Short();
2790 // Double_t dMassLambda = vZero->MassLambda();
2791 
2792  fh1QAV0Status[iIndexHisto]->Fill(vZero->GetOnFlyStatus());
2793 
2794  AliAODTrack* trackNeg = (AliAODTrack*)vZero->GetDaughter(1); // negative track
2795  AliAODTrack* trackPos = (AliAODTrack*)vZero->GetDaughter(0); // positive track
2796 
2797  Short_t fTotalCharge = 0;
2798  for(Int_t i = 0; i < 2; i++)
2799  {
2800  AliAODTrack* track = (AliAODTrack*)vZero->GetDaughter(i); // track
2801  // Tracks TPC OK
2802  fh1QAV0TPCRefit[iIndexHisto]->Fill(track->IsOn(AliAODTrack::kTPCrefit));
2803  Double_t nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
2804  fh1QAV0TPCRows[iIndexHisto]->Fill(nCrossedRowsTPC);
2805  Int_t findable = track->GetTPCNclsF();
2806  fh1QAV0TPCFindable[iIndexHisto]->Fill(findable);
2807  if(findable != 0)
2808  {
2809  fh1QAV0TPCRowsFind[iIndexHisto]->Fill(nCrossedRowsTPC / findable);
2810  }
2811  // Daughters: pseudo-rapidity cut
2812  fh1QAV0Eta[iIndexHisto]->Fill(track->Eta());
2813  if((nCrossedRowsTPC > (160. / (250. - 85.) * (255.*TMath::Abs(tan(track->Theta())) - 85.)) + 20.) && (track->Eta() < 0) && (track->Pt() > 0.15))
2814 // if (IsCandK0s)
2815  {
2816  fh2QAV0EtaRows[iIndexHisto]->Fill(track->Eta(), nCrossedRowsTPC);
2817  fh2QAV0PtRows[iIndexHisto]->Fill(track->Pt(), nCrossedRowsTPC);
2818  fh2QAV0PhiRows[iIndexHisto]->Fill(track->Phi(), nCrossedRowsTPC);
2819  fh2QAV0NClRows[iIndexHisto]->Fill(findable, nCrossedRowsTPC);
2820  fh2QAV0EtaNCl[iIndexHisto]->Fill(track->Eta(), findable);
2821  }
2822 
2823  // Daughters: transverse momentum cut
2824  fh1QAV0Pt[iIndexHisto]->Fill(track->Pt());
2825  fTotalCharge += track->Charge();
2826  }
2827  fh1QAV0Charge[iIndexHisto]->Fill(fTotalCharge);
2828 
2829  // Daughters: Impact parameter of daughters to prim vtx
2830  fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaNegToPrimVertex()));
2831  fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaPosToPrimVertex()));
2832 // fh2CutDCAVtx[iIndexHisto]->Fill(dMassK0s,TMath::Abs(vZero->DcaNegToPrimVertex()));
2833 
2834  // Daughters: DCA
2835  fh1QAV0DCAV0[iIndexHisto]->Fill(vZero->DcaV0Daughters());
2836 // fh2CutDCAV0[iIndexHisto]->Fill(dMassK0s,vZero->DcaV0Daughters());
2837 
2838  // V0: Cosine of the pointing angle
2839  fh1QAV0Cos[iIndexHisto]->Fill(vZero->CosPointingAngle(vtx));
2840 // fh2CutCos[iIndexHisto]->Fill(dMassK0s,vZero->CosPointingAngle(vtx));
2841 
2842  // V0: Fiducial volume
2843  Double_t xyz[3];
2844  vZero->GetSecondaryVtx(xyz);
2845  Double_t r2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
2846  fh1QAV0R[iIndexHisto]->Fill(TMath::Sqrt(r2));
2847 
2848  Double_t dAlpha = vZero->AlphaV0();
2849  Double_t dPtArm = vZero->PtArmV0();
2850 
2851  if(IsCandK0s)
2852  {
2853  if(IsInPeakK0s)
2854  {
2855 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2856 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2857  fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2858  fh2QAV0PtPtK0sPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2859  fh2ArmPodK0s[iIndexHisto]->Fill(dAlpha, dPtArm);
2860  }
2861  fh2QAV0EtaEtaK0s[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2862  fh2QAV0PhiPhiK0s[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2863  fh1QAV0RapK0s[iIndexHisto]->Fill(vZero->RapK0Short());
2864  }
2865 
2866  if(IsCandLambda)
2867  {
2868  if(IsInPeakLambda)
2869  {
2870 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2871 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2872  fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2873  fh2QAV0PtPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2874  fh2ArmPodLambda[iIndexHisto]->Fill(dAlpha, dPtArm);
2875  }
2876  fh2QAV0EtaEtaLambda[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2877  fh2QAV0PhiPhiLambda[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2878  fh1QAV0RapLambda[iIndexHisto]->Fill(vZero->RapLambda());
2879  }
2880 
2881  fh2ArmPod[iIndexHisto]->Fill(dAlpha, dPtArm);
2882 
2883 }
2884 
2885 void AliAnalysisTaskV0sInJets::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*/)
2886 {
2887  if(isK)
2888  {
2889  fh1V0CounterCentK0s[iCent]->Fill(iCut);
2890  fh1V0InvMassK0sAll[iCut]->Fill(mK);
2891  }
2892  if(isL)
2893  {
2894  fh1V0CounterCentLambda[iCent]->Fill(iCut);
2895  fh1V0InvMassLambdaAll[iCut]->Fill(mL);
2896  }
2897  if(isAL)
2898  {
2899  fh1V0CounterCentALambda[iCent]->Fill(iCut);
2900  fh1V0InvMassALambdaAll[iCut]->Fill(mAL);
2901  }
2902 }
2903 
2904 Bool_t AliAnalysisTaskV0sInJets::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
2905 {
2906 // decides whether a particle is inside a jet cone
2907  if(!part1 || !part2)
2908  return kFALSE;
2909 
2910  TVector3 vecMom2(part2->Px(), part2->Py(), part2->Pz());
2911  TVector3 vecMom1(part1->Px(), part1->Py(), part1->Pz());
2912  Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
2913  if(dR < dRMax) // momentum vectors of part1 and part2 are closer than dRMax
2914  return kTRUE;
2915  return kFALSE;
2916 }
2917 
2918 Bool_t AliAnalysisTaskV0sInJets::OverlapWithJets(const TClonesArray* array, const AliVParticle* part, Double_t dDistance) const
2919 {
2920 // decides whether a cone overlaps with other jets
2921  if(!part)
2922  {
2923  if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No part\n");
2924  return kFALSE;
2925  }
2926  if(!array)
2927  {
2928  if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No array\n");
2929  return kFALSE;
2930  }
2931  Int_t iNJets = array->GetEntriesFast();
2932  if(iNJets <= 0)
2933  {
2934  if(fDebug > 2) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Warning: No jets\n");
2935  return kFALSE;
2936  }
2937  AliVParticle* jet = 0;
2938  for(Int_t iJet = 0; iJet < iNJets; iJet++)
2939  {
2940  jet = (AliVParticle*)array->At(iJet);
2941  if(!jet)
2942  {
2943  if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: Failed to load jet %d/%d\n", iJet, iNJets);
2944  continue;
2945  }
2946  if(IsParticleInCone(part, jet, dDistance))
2947  return kTRUE;
2948  }
2949  return kFALSE;
2950 }
2951 
2952 AliAODJet* AliAnalysisTaskV0sInJets::GetRandomCone(const TClonesArray* array, Double_t dEtaConeMax, Double_t dDistance) const
2953 {
2954 // generate a random cone which does not overlap with selected jets
2955 // printf("Generating random cone...\n");
2956  TLorentzVector vecCone;
2957  AliAODJet* part = 0;
2958  Double_t dEta, dPhi;
2959  Int_t iNTrialsMax = 10;
2960  Bool_t bStatus = kFALSE;
2961  for(Int_t iTry = 0; iTry < iNTrialsMax; iTry++)
2962  {
2963 // printf("Try %d\n",iTry);
2964  dEta = dEtaConeMax * (2 * fRandom->Rndm() - 1.); // random eta in [-dEtaConeMax,+dEtaConeMax]
2965  dPhi = TMath::TwoPi() * fRandom->Rndm(); // random phi in [0,2*Pi]
2966  vecCone.SetPtEtaPhiM(1., dEta, dPhi, 0.);
2967  part = new AliAODJet(vecCone);
2968  if(!OverlapWithJets(array, part, dDistance))
2969  {
2970  bStatus = kTRUE;
2971 // printf("Success\n");
2972  break;
2973  }
2974  else
2975  delete part;
2976  }
2977  if(!bStatus)
2978  part = 0;
2979  return part;
2980 }
2981 
2982 AliAODJet* AliAnalysisTaskV0sInJets::GetMedianCluster(const TClonesArray* array, Double_t dEtaConeMax) const
2983 {
2984 // sort kt clusters by pT/area and return the middle one, based on code in AliAnalysisTaskJetChem
2985  if(!array)
2986  {
2987  if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Error: No array\n");
2988  return NULL;
2989  }
2990  Int_t iNClTot = array->GetEntriesFast(); // number of all clusters in the array
2991  Int_t iNCl = 0; // number of accepted clusters
2992 
2993  // get list of densities
2994  std::vector<std::vector<Double_t> > vecListClusters; // vector that contains pairs [ index, density ]
2995 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Loop over %d clusters.\n", iNClTot);
2996  for(Int_t ij = 0; ij < iNClTot; ij++)
2997  {
2998  AliAODJet* clusterBg = (AliAODJet*)(array->At(ij));
2999  if(!clusterBg)
3000  {
3001 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: cluster %d/%d not ok\n", ij, iNClTot);
3002  return NULL;
3003  }
3004  if(TMath::Abs(clusterBg->Eta()) > 0.9 - fdRadiusJetBg)
3005  continue;
3006 // fh2EtaPhiMedCone[0]->Fill(clusterBg->Eta(), clusterBg->Phi());
3007 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Cluster %d/%d used as accepted cluster %d.\n", ij, iNClTot, int(vecListClusters.size()));
3008  Double_t dPtBg = clusterBg->Pt();
3009  Double_t dAreaBg = clusterBg->EffectiveAreaCharged();
3010  Double_t dDensityBg = 0;
3011  if(dAreaBg > 0)
3012  dDensityBg = dPtBg / dAreaBg;
3013  std::vector<Double_t> vecCluster;
3014  vecCluster.push_back(ij);
3015  vecCluster.push_back(dDensityBg);
3016  vecListClusters.push_back(vecCluster);
3017  }
3018  iNCl = vecListClusters.size();
3019  if(iNCl < 3) // need at least 3 clusters (skipping 2 highest)
3020  {
3021 // if(fDebug > 2) printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Warning: Too little clusters\n");
3022  return NULL;
3023  }
3024 
3025 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Original lists:\n");
3026 // for(Int_t i = 0; i < iNCl; i++)
3027 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3028 
3029  // sort list of indeces by density in descending order
3030  std::sort(vecListClusters.begin(), vecListClusters.end(), CompareClusters);
3031 
3032 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Sorted lists:\n");
3033 // for(Int_t i = 0; i < iNCl; i++)
3034 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3035 
3036  // get median cluster with median density
3037  AliAODJet* clusterMed = 0;
3038  Int_t iIndex = 0; // index of the median cluster in the sorted list
3039  Int_t iIndexMed = 0; // index of the median cluster in the original array
3040  if(TMath::Odd(iNCl)) // odd number of clusters
3041  {
3042  iIndex = (Int_t)(0.5 * (iNCl + 1)); // = (n - skip + 1)/2 + 1, skip = 2
3043 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Odd, median index = %d/%d\n", iIndex, iNCl);
3044  }
3045  else // even number: picking randomly one of the two closest to median
3046  {
3047  Int_t iIndex1 = (Int_t)(0.5 * iNCl); // = (n - skip)/2 + 1, skip = 2
3048  Int_t iIndex2 = (Int_t)(0.5 * iNCl + 1); // = (n - skip)/2 + 1 + 1, skip = 2
3049  iIndex = ((fRandom->Rndm() > 0.5) ? iIndex1 : iIndex2);
3050 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Even, median index = %d or %d -> %d/%d\n", iIndex1, iIndex2, iIndex, iNCl);
3051  }
3052  iIndexMed = Int_t((vecListClusters[iIndex])[0]);
3053 
3054 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: getting median cluster %d/%d ok, rho = %g\n", iIndexMed, iNClTot, (vecListClusters[iIndex])[1]);
3055  clusterMed = (AliAODJet*)(array->At(iIndexMed));
3056 
3057  if(TMath::Abs(clusterMed->Eta()) > dEtaConeMax)
3058  return NULL;
3059 
3060  return clusterMed;
3061 }
3062 
3064 {
3065 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
3066  Double_t dEpsilon = 1e-2;
3067  Double_t dR = dRadius;
3068  Double_t dD = dDistance;
3069  if(TMath::Abs(dR) < dEpsilon)
3070  {
3071  if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::AreaCircSegment: Error: Too small radius: %f < %f\n", dR, dEpsilon);
3072  return 0.;
3073  }
3074  if(dD >= dR)
3075  return 0.;
3076  if(dD <= -dR)
3077  return TMath::Pi() * dR * dR;
3078  return dR * dR * TMath::ACos(dD / dR) - dD * TMath::Sqrt(dR * dR - dD * dD);
3079 }
3080 
3081 Bool_t AliAnalysisTaskV0sInJets::IsSelectedForJets(AliAODEvent* fAOD, Double_t dVtxZCut, Double_t dVtxR2Cut, Double_t dCentCutLo, Double_t dCentCutUp, Bool_t bCutDeltaZ, Double_t dDeltaZMax)
3082 {
3083 // event selection
3084  AliAODVertex* vertex = fAOD->GetPrimaryVertex();
3085  if(!vertex)
3086  return kFALSE;
3087  Int_t iNContribMin = 3;
3088  if(!fbIsPbPb)
3089  iNContribMin = 2;
3090  if(vertex->GetNContributors() < iNContribMin)
3091  return kFALSE;
3092  TString vtxTitle(vertex->GetTitle());
3093  if(vtxTitle.Contains("TPCVertex"))
3094  return kFALSE;
3095  Double_t zVertex = vertex->GetZ();
3096  if(TMath::Abs(zVertex) > dVtxZCut)
3097  return kFALSE;
3098  if(bCutDeltaZ)
3099  {
3100  AliAODVertex* vertexSPD = fAOD->GetPrimaryVertexSPD();
3101  if(!vertexSPD)
3102  {
3103 // printf("IsSelectedForJets: Error: No SPD vertex\n");
3104  return kFALSE;
3105  }
3106  Double_t zVertexSPD = vertexSPD->GetZ();
3107  if(TMath::Abs(zVertex - zVertexSPD) > dDeltaZMax)
3108  {
3109 // printf("IsSelectedForJets: Rejecting event due to delta z = %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
3110  return kFALSE;
3111  }
3112 // printf("IsSelectedForJets: Event OK: %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
3113  }
3114  Double_t xVertex = vertex->GetX();
3115  Double_t yVertex = vertex->GetY();
3116  Double_t radiusSq = yVertex * yVertex + xVertex * xVertex;
3117  if(radiusSq > dVtxR2Cut)
3118  return kFALSE;
3120 // centrality = fAOD->GetHeader()->GetCentrality();
3121  centrality = ((AliVAODHeader*)fAOD->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
3122  if(fbIsPbPb)
3123  {
3124  if(centrality < 0)
3125  return kFALSE;
3126  if((dCentCutUp < 0) || (dCentCutLo < 0) || (dCentCutUp > 100) || (dCentCutLo > 100) || (dCentCutLo > dCentCutUp))
3127  return kFALSE;
3128  if((centrality < dCentCutLo) || (centrality > dCentCutUp))
3129  return kFALSE;
3130  }
3131  else
3132  {
3133  if(centrality != -1)
3134  return kFALSE;
3135  }
3136  return kTRUE;
3137 }
3138 
3140 {
3141 // returns index of the centrality bin corresponding to the provided value of centrality
3143  return -1;
3144  for(Int_t i = 0; i < fgkiNBinsCent; i++)
3145  {
3146  if(centrality <= fgkiCentBinRanges[i])
3147  return i;
3148  }
3149  return -1;
3150 }
3151 
3153 {
3154 // returns the upper edge of the centrality bin corresponding to the provided value of index
3155  if(index < 0 || index >= fgkiNBinsCent)
3156  return -1;
3157  return fgkiCentBinRanges[index];
3158 }
3159 
3161 {
3162 // get string with centrality range for given bin
3163  TString lowerEdge = ((index == 0) ? "0" : Form("%d", GetCentralityBinEdge(index - 1)));
3164  TString upperEdge = Form("%d", GetCentralityBinEdge(index));
3165  return Form("%s-%s %%", lowerEdge.Data(), upperEdge.Data());
3166 }
3167 
3169 {
3170 // estimation of the sigma of the invariant-mass peak as a function of pT and particle type
3171  switch(particle)
3172  {
3173  case 0: // K0S
3174  return 0.0044 + 0.0004 * (pt - 1.);
3175  break;
3176  case 1: // Lambda
3177  return 0.0023 + 0.00034 * (pt - 1.);
3178  break;
3179  default:
3180  return 0;
3181  break;
3182  }
3183 }
3184 
3185 bool AliAnalysisTaskV0sInJets::CompareClusters(const std::vector<Double_t> cluster1, const std::vector<Double_t> cluster2)
3186 {
3187  return (cluster1[1] > cluster2[1]);
3188 }
TH1D * fh1V0ALambdaPtMCRecFalse[fgkiNBinsCent]
Bool_t IsParticleInCone(const AliVParticle *part1, const AliVParticle *part2, Double_t dRMax) const
THnSparseD * fhnV0ALambdaInJetsMCFD[fgkiNBinsCent]
TH1D * fh1V0CandPerEventCentALambda[fgkiNBinsCent]
TList * fOutputListMC
Output list for checking cuts.
THnSparse * fh3V0ALambdaEtaPtMassMCRec[fgkiNBinsCent]
TH1D * fh1V0AXiPtMCGen[fgkiNBinsCent]
TH1D * fh1QAV0R[fgkiNQAIndeces]
cosine of pointing angle (CPA)
Int_t GetCentralityBinIndex(Double_t centrality)
THnSparse * fh3V0K0sInJetPtMassMCRec[fgkiNBinsCent]
pt spectrum of generated K0s in jet
THnSparse * fhnV0OutJetK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt vs jet pt, in centrality bins.
THnSparse * fhnV0ALambdaInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH1D * fh1NMedConeCent
random cone eta-pT
TH1D * fh1DCAOutALambda[fgkiNBinsCent]
double Double_t
Definition: External.C:58
THnSparse * fhnV0InPerpK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt vs jet pt, in centrality bins.
TH2D * fh2V0LambdaPtMassMCRec[fgkiNBinsCent]
TH2D * fh2QAV0EtaPtALambdaPeak[fgkiNQAIndeces]
THnSparse * fhnV0InJetK0s[fgkiNBinsCent]
V0 inv mass vs pt before and after cuts, in centrality bins.
static const Int_t fgkiNBinsPtV0Init
THnSparse * fhnV0K0sInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
mass-eta-pt spectrum of successfully reconstructed K0s in jet
THnSparse * fhnV0InRndALambda[fgkiNBinsCent]
THnSparse * fh3CCMassCorrelKNotL
mass correlation of candidates
TH1D * fh1V0K0sPtMCGen[fgkiNBinsCent]
DCA between daughters of V0 outside jets, in centrality bins.
TH2D * fh2QAV0PtRows[fgkiNQAIndeces]
pseudorapidity vs TPC rows
AliAODEvent * fAODOut
Input AOD event.
TH1D * fh1V0CounterCentALambda[fgkiNBinsCent]
TH1D * fh1DCAInLambda[fgkiNBinsCent]
TH1D * fh1V0LambdaPtMCRecFalse[fgkiNBinsCent]
THnSparse * fhnV0InMedLambda[fgkiNBinsCent]
TH1D * fh1PhiJet[fgkiNBinsCent]
jet eta-pT
TList * fOutputListStd
Output AOD event.
static const Double_t fgkdBinsPtJet[2]
virtual void UserExec(Option_t *option)
TH2D * fh2ArmPodALambda[fgkiNQAIndeces]
TH1D * fh1V0CandPerEventCentLambda[fgkiNBinsCent]
THnSparse * fh3V0LambdaEtaPtMassMCRec[fgkiNBinsCent]
centrality
char Char_t
Definition: External.C:18
TList * fOutputListQA
Output list for standard analysis results.
static bool CompareClusters(const std::vector< Double_t > cluster1, const std::vector< Double_t > cluster2)
TH1D * fh1NJetPerEvent[fgkiNBinsCent]
jet phi
THnSparse * fh3V0K0sInJetEtaPtMCGen[fgkiNBinsCent]
mass-pt spectrum of successfully reconstructed K0s in jet
static const Double_t fgkdMassK0sMax
Bool_t OverlapWithJets(const TClonesArray *array, const AliVParticle *cone, Double_t dDistance) const
TH1D * fh1QAV0Pt[fgkiNQAIndeces]
TH1D * fh1V0ALambdaPtMCRec[fgkiNBinsCent]
THnSparseD * fhnV0ALambdaBulkMCFD[fgkiNBinsCent]
THnSparse * fh4V0ALambdaInJetEtaPtMassMCRec[fgkiNBinsCent]
static const Int_t fgkiNBinsMassLambda
static const Double_t fgkdMassLambdaMax
TH1D * fh1V0InvMassLambdaCent[fgkiNBinsCent]
TH1D * fh1QAV0TPCFindable[fgkiNQAIndeces]
crossed TPC pad rows
THnSparseD * fhnV0LambdaInclMCFD[fgkiNBinsCent]
TH2D * fh2CCLambda
K0s candidates in Lambda peak.
TH1D * fh1V0InvMassK0sAll[fgkiNCategV0]
number of K0s candidates after various cuts
TH2D * fh2V0PtJetAngleLambda[fgkiNBinsCent]
AliAODJet * GetRandomCone(const TClonesArray *array, Double_t dEtaConeMax, Double_t dDistance) const
TH1D * fh1EventCent
number of events for different selection steps and different centralities
void FillCandidates(Double_t mK, Double_t mL, Double_t mAL, Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut, Int_t iCent)
TH2D * fh2V0ALambdaInJetPtMCRec[fgkiNBinsCent]
TH1D * fh1QAV0DCAVtx[fgkiNQAIndeces]
charge
TH2D * fh2ArmPodLambda[fgkiNQAIndeces]
Double_t fdCutVertexZ
random-number generator
THnSparse * fhnV0InRndK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt vs jet pt, in centrality bins.
THnSparseD * fhnV0LambdaInJetsMCFD[fgkiNBinsCent]
THnSparse * fh3V0LambdaInJetEtaPtMCGen[fgkiNBinsCent]
TH1D * fh1QAV0Status[fgkiNQAIndeces]
TH1D * fh1DCAOutLambda[fgkiNBinsCent]
THnSparse * fh3CCMassCorrelLNotK
mass correlation of candidates
TH1D * fh1V0ALambdaPtMCGen[fgkiNBinsCent]
THnSparse * fhnV0InJetLambda[fgkiNBinsCent]
THnSparse * fhnV0NoJetALambda[fgkiNBinsCent]
TH1D * fh1VtxZ[fgkiNBinsCent]
number of tracks vs centrality
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
TH1D * fh1V0InvMassLambdaAll[fgkiNCategV0]
number of Lambda candidates after various cuts
TH2D * fh2V0ALambdaInJetPtMCGen[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
TH2D * fh2QAV0PtPtK0sPeak[fgkiNQAIndeces]
V0 rapidity.
Double_t AreaCircSegment(Double_t dRadius, Double_t dDistance) const
TH1D * fh1PtJet[fgkiNBinsCent]
number of V0 cand per event
TH1D * fh1V0CandPerEventCentK0s[fgkiNBinsCent]
Armenteros-Podolanski.
TH2D * fh2QAV0NClRows[fgkiNQAIndeces]
azimuth vs TPC rows
TH2D * fh2V0K0sMCPtGenPtRec[fgkiNBinsCent]
K0s mass resolution vs pt.
TH2D * fh2V0K0sPtMassMCRec[fgkiNBinsCent]
pt spectrum of all generated K0s in event
THnSparse * fhnV0OutJetALambda[fgkiNBinsCent]
THnSparse * fh4V0LambdaInJetEtaPtMassMCRec[fgkiNBinsCent]
TH1D * fh1DCAInK0s[fgkiNBinsCent]
pt jet vs angle V0-jet, in centrality bins
static const Int_t fgkiNCategV0
area of excluded cones for outside-cones V0s
TH2D * fh2V0K0sInJetPtMCGen[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
THnSparse * fhnV0InclusiveK0s[fgkiNBinsCent]
V0 invariant mass, in centrality bins.
THnSparse * fhnV0InclusiveLambda[fgkiNBinsCent]
TH2D * fh2EtaPhiRndCone[fgkiNBinsCent]
number of generated random cones in centrality bins
TH2D * fh2V0ALambdaPtMassMCRec[fgkiNBinsCent]
TH1D * fh1V0ALambdaPt[fgkiNBinsCent]
int Int_t
Definition: External.C:63
TH1D * fh1V0XiPtMCGen[fgkiNBinsCent]
THnSparse * fhnV0InMedALambda[fgkiNBinsCent]
unsigned int UInt_t
Definition: External.C:33
TH1D * fh1QAV0RapALambda[fgkiNQAIndeces]
THnSparse * fh3CCMassCorrelBoth
Lambda candidates in K0s peak.
TH1D * fh1QACTau3D[fgkiNQAIndeces]
lifetime calculated in xy
TH1D * fh1NRndConeCent
number of jets per event
TH2D * fh2QAV0PtPtALambdaPeak[fgkiNQAIndeces]
THnSparse * fhnV0InPerpALambda[fgkiNBinsCent]
TH2D * fh2ArmPod[fgkiNQAIndeces]
lifetime calculated in xyz
TH2D * fh2V0K0sMCResolMPt[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
Definition: External.C:228
Definition: External.C:212
THnSparse * fhnV0LambdaInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH2D * fh2QAV0EtaEtaK0s[fgkiNQAIndeces]
daughters pseudorapidity vs V0 pt, in mass peak
THnSparse * fhnV0InclusiveALambda[fgkiNBinsCent]
THnSparse * fhnV0InPerpLambda[fgkiNBinsCent]
THnSparse * fhnV0OutJetLambda[fgkiNBinsCent]
TH1D * fh1QACTau2D[fgkiNQAIndeces]
radial distance between prim vtx and decay vertex
TH2D * fh2QAV0PtPtLambdaPeak[fgkiNQAIndeces]
TH2D * fh2EtaPtJet[fgkiNBinsCent]
jet eta
TH1D * fh1AreaExcluded
median-cluster cone eta-phi
TH2D * fh2QAV0PhiRows[fgkiNQAIndeces]
pt vs TPC rows
Int_t fiAODAnalysis
Output list for MC related results.
THnSparse * fhnV0NoJetK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt, in centrality bins.
THnSparse * fh4V0K0sInJetEtaPtMassMCRec[fgkiNBinsCent]
eta-pt spectrum of generated K0s in jet
TH2D * fh2V0LambdaEtaPtMCGen[fgkiNBinsCent]
THnSparse * fhnV0LambdaInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
static const Int_t fgkiCentBinRanges[fgkiNBinsCent]
TH2D * fh2QAV0EtaEtaALambda[fgkiNQAIndeces]
THnSparse * fhnV0NoJetLambda[fgkiNBinsCent]
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)
TH1D * fh1QAV0DCAV0[fgkiNQAIndeces]
DCA of daughters to prim vtx.
short Short_t
Definition: External.C:23
THnSparseD * fhnV0LambdaBulkMCFD[fgkiNBinsCent]
TH2D * fh2V0LambdaInJetPtMCGen[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
TH1D * fh1QAV0RapLambda[fgkiNQAIndeces]
TH2D * fh2QAV0PhiPhiLambda[fgkiNQAIndeces]
TH1D * fh1QAV0Cos[fgkiNQAIndeces]
DCA between daughters.
TH2D * fh2V0PtJetAngleK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt, in centrality bins.
TH2D * fh2EtaPhiMedCone[fgkiNBinsCent]
number of found median-cluster cones in centrality bins
TH2D * fh2QAV0EtaEtaLambda[fgkiNQAIndeces]
THnSparse * fhnV0InJetALambda[fgkiNBinsCent]
TH2D * fh2V0PtJetAngleALambda[fgkiNBinsCent]
TH1D * fh1V0CandPerEvent
xy coordinates of the primary vertex
TH1D * fh1DCAOutK0s[fgkiNBinsCent]
DCA between daughters of V0 inside jets, in centrality bins.
TH2D * fh2QAV0PhiPhiALambda[fgkiNQAIndeces]
static const Double_t fgkdBinsPtV0[2]
TH1D * fh1V0LambdaPtMCGen[fgkiNBinsCent]
TH1D * fh1EventCent2Jets
number of events for different centralities
Double_t MassPeakSigmaOld(Double_t pt, Int_t particle)
THnSparse * fhnV0InMedK0s[fgkiNBinsCent]
V0 invariant mass vs V0 pt vs jet pt, in centrality bins.
TH1D * fh1EtaJet[fgkiNBinsCent]
pt spectra of jets for normalisation of in-jet V0 spectra
TH1D * fh1QAV0TPCRowsFind[fgkiNQAIndeces]
findable clusters
TH1D * fh1EventCent2NoJets
number of events for different centralities
TH1D * fh1V0InvMassALambdaAll[fgkiNCategV0]
number of ALambda candidates after various cuts
TH1D * fh1QAV0RapK0s[fgkiNQAIndeces]
daughters azimuth vs azimuth
TH1D * fh1V0InvMassK0sCent[fgkiNBinsCent]
number of K0s candidates per event, in centrality bins
THnSparse * fh3V0LambdaInJetPtMassMCRec[fgkiNBinsCent]
TH2D * fh2CCK0s
Armenteros-Podolanski.
static const Double_t fgkdMassK0sMin
TH1D * fh1V0K0sPtMCRecFalse[fgkiNBinsCent]
pt-mass spectrum of successfully reconstructed K0s in event
TList * fOutputListCuts
Output list for quality assurance.
TH2D * fh2V0ALambdaMCResolMPt[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
TH2D * fh2ArmPodK0s[fgkiNQAIndeces]
daughters pt vs pt, in mass peak
void FillQAHistogramV0(AliAODVertex *vtx, const AliAODv0 *vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda)
const char Option_t
Definition: External.C:48
TH2D * fh2QAV0EtaPtK0sPeak[fgkiNQAIndeces]
V0 invariant mass, selection steps.
THnSparse * fhnV0K0sInclDaughterEtaPtPtMCRec[fgkiNBinsCent]
eta-pt-mass spectrum of successfully reconstructed K0s in event
TH1D * fh1V0CounterCentLambda[fgkiNBinsCent]
K0s generated pt vs reconstructed pt.
THnSparseD * fhnV0ALambdaInclMCFD[fgkiNBinsCent]
TH2D * fh2VtxXY[fgkiNBinsCent]
z coordinate of the primary vertex
THnSparse * fh3V0ALambdaInJetEtaPtMCGen[fgkiNBinsCent]
bool Bool_t
Definition: External.C:53
TH2D * fh2V0ALambdaEtaPtMCGen[fgkiNBinsCent]
THnSparse * fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[fgkiNBinsCent]
TH1D * fh1QAV0Eta[fgkiNQAIndeces]
ratio rows/clusters
TH1D * fh1V0CounterCentK0s[fgkiNBinsCent]
pseudorapidity vs clusters
AliAODJet * GetMedianCluster(const TClonesArray *array, Double_t dEtaConeMax) const
TH2D * fh2V0LambdaMCPtGenPtRec[fgkiNBinsCent]
TH1D * fh1V0InvMassALambdaCent[fgkiNBinsCent]
THnSparse * fhnV0InRndLambda[fgkiNBinsCent]
TH2D * fh2V0ALambdaMCPtGenPtRec[fgkiNBinsCent]
TH1D * fh1QAV0Charge[fgkiNQAIndeces]
pt
TH2D * fh2QAV0EtaNCl[fgkiNQAIndeces]
clusters vs TPC rows
TH1D * fh1QAV0TPCRefit[fgkiNQAIndeces]
online vs offline reconstructed V0 candidates
TH1D * fh1EventCounterCutCent[fgkiNBinsCent]
number of events for different selection steps
TH1D * fh1QAV0TPCRows[fgkiNQAIndeces]
TPC refit on vs off.
Bool_t IsSelectedForJets(AliAODEvent *fAOD, Double_t dVtxZCut, Double_t dVtxR2Cut, Double_t dCentCutLo, Double_t dCentCutUp, Bool_t bCutDeltaZ=kFALSE, Double_t dDeltaZMax=100.)
static const Double_t fgkdMassLambdaMin
THnSparse * fh3V0K0sEtaPtMassMCRec[fgkiNBinsCent]
eta-pt spectrum of all generated K0s in event
TH2D * fh2QAV0EtaRows[fgkiNQAIndeces]
pseudorapidity
TH1D * fh1EventCent2
number of events for different centralities
TH2D * fh2QAV0PhiPhiK0s[fgkiNQAIndeces]
daughters pseudorapidity vs pseudorapidity
Definition: External.C:196
static const Int_t fgkiNBinsPtJetInit
TH1D * fh1DCAInALambda[fgkiNBinsCent]
TH2D * fh2V0LambdaMCResolMPt[fgkiNBinsCent]
eta_daughter-pt_daughter-pt_V0 reconstructed
TH2D * fh2QAV0EtaPtLambdaPeak[fgkiNQAIndeces]
TH2D * fh2V0K0sEtaPtMCGen[fgkiNBinsCent]
pt spectrum of false reconstructed K0s in event
TH2D * fh2EventCentTracks
number of events for different centralities
THnSparse * fh3V0ALambdaInJetPtMassMCRec[fgkiNBinsCent]