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