]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/StrangenessInJets/AliAnalysisTaskV0sInJets.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
create static objects without new + add syswatch log to the standalone filtering...
[u/mrichter/AliRoot.git] / PWGJE / StrangenessInJets / AliAnalysisTaskV0sInJets.cxx
CommitLineData
257f0eee 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
7729cc34 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
46ClassImp(AliAnalysisTaskV0sInJets)
47
48// upper edges of centrality bins
49const 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
55const Double_t AliAnalysisTaskV0sInJets::fgkdBinsPtV0[2] = {0, 30};
56const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtV0 = sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtV0)/sizeof((AliAnalysisTaskV0sInJets::fgkdBinsPtV0)[0])-1;
57const 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
59const Double_t AliAnalysisTaskV0sInJets::fgkdBinsPtJet[2] = {0, 100};
60const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJet = sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet)/sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet[0])-1;
61const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[AliAnalysisTaskV0sInJets::fgkiNBinsPtJet]-(AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[0])/5.); // bin width 5 GeV/c
62// axis: K0S invariant mass
63const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsMassK0s = 200;
64const Double_t AliAnalysisTaskV0sInJets::fgkfMassK0sMin = 0.4;
65const Double_t AliAnalysisTaskV0sInJets::fgkfMassK0sMax = 0.6;
66// axis: Lambda invariant mass
67const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsMassLambda = 200;
68const Double_t AliAnalysisTaskV0sInJets::fgkfMassLambdaMin = 1.05;
69const Double_t AliAnalysisTaskV0sInJets::fgkfMassLambdaMax = 1.25;
70
71
72// Default constructor
73AliAnalysisTaskV0sInJets::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 ffCutPtJetMin(0),
93 ffCutPtTrackMin(5),
94 ffRadiusJet(0.4),
95 fbJetSelection(0),
96 fbMCAnalysis(0),
97// fbTreeOutput(0),
98 fRandom(0),
99
100 ffCutVertexZ(10),
101 ffCutVertexR2(1),
102 ffCutCentLow(0),
103 ffCutCentHigh(80),
104
105/*
106 fBranchV0Rec(0),
107 fBranchV0Gen(0),
108 fBranchJet(0),
109 fEventInfo(0),
110*/
111 ffCentrality(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
319AliAnalysisTaskV0sInJets::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 ffCutPtJetMin(0),
339 ffCutPtTrackMin(5),
340 ffRadiusJet(0.4),
341 fbJetSelection(0),
342 fbMCAnalysis(0),
343// fbTreeOutput(0),
344 fRandom(0),
345
346 ffCutVertexZ(10),
347 ffCutVertexR2(1),
348 ffCutCentLow(0),
349 ffCutCentHigh(80),
350/*
351 fBranchV0Rec(0),
352 fBranchV0Gen(0),
353 fBranchJet(0),
354 fEventInfo(0),
355*/
356 ffCentrality(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
573AliAnalysisTaskV0sInJets::~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
597void 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 fPtV0Min = fgkdBinsPtV0[0];
705 Double_t fPtV0Max = fgkdBinsPtV0[fgkiNBinsPtV0];
706 Int_t iNJetPtBins = fgkiNBinsPtJetInit;
707 Double_t fJetPtMin = fgkdBinsPtJet[0];
708 Double_t fJetPtMax = fgkdBinsPtJet[fgkiNBinsPtJet];
709
710 fh2CCK0s = new TH2D("fh2CCK0s","K0s candidates in Lambda peak",fgkiNBinsMassK0s,fgkfMassK0sMin,fgkfMassK0sMax,iNBinsPtV0,fPtV0Min,fPtV0Max);
711 fh2CCLambda = new TH2D("fh2CCLambda","Lambda candidates in K0s peak",fgkiNBinsMassLambda,fgkfMassLambdaMin,fgkfMassLambdaMax,iNBinsPtV0,fPtV0Min,fPtV0Max);
712 Int_t binsCorrel[3] = {fgkiNBinsMassK0s, fgkiNBinsMassLambda, iNBinsPtV0};
713 Double_t xminCorrel[3] = {fgkfMassK0sMin, fgkfMassLambdaMin, fPtV0Min};
714 Double_t xmaxCorrel[3] = {fgkfMassK0sMax, fgkfMassLambdaMax, fPtV0Max};
715// Int_t binsCorrel[3] = {200, 200, iNBinsPtV0};
716// Double_t xminCorrel[3] = {0, 0, fPtV0Min};
717// Double_t xmaxCorrel[3] = {2, 2, fPtV0Max};
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 fStepEtaV0 = 0.025;
728 Double_t fRangeEtaV0Max = 0.8;
729 const Int_t iNBinsEtaV0 = 2*Int_t(fRangeEtaV0Max/fStepEtaV0);
730 // inclusive
731 const Int_t iNDimIncl = 3;
732 Int_t binsKIncl[iNDimIncl] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
733 Double_t xminKIncl[iNDimIncl] = {fgkfMassK0sMin, fPtV0Min, -fRangeEtaV0Max};
734 Double_t xmaxKIncl[iNDimIncl] = {fgkfMassK0sMax, fPtV0Max, fRangeEtaV0Max};
735 Int_t binsLIncl[iNDimIncl] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
736 Double_t xminLIncl[iNDimIncl] = {fgkfMassLambdaMin, fPtV0Min, -fRangeEtaV0Max};
737 Double_t xmaxLIncl[iNDimIncl] = {fgkfMassLambdaMax, fPtV0Max, fRangeEtaV0Max};
738 // binning in jets
739 const Int_t iNDimInJC = 4;
740 Int_t binsKInJC[iNDimInJC] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
741 Double_t xminKInJC[iNDimInJC] = {fgkfMassK0sMin, fPtV0Min, -fRangeEtaV0Max, fJetPtMin};
742 Double_t xmaxKInJC[iNDimInJC] = {fgkfMassK0sMax, fPtV0Max, fRangeEtaV0Max, fJetPtMax};
743 Int_t binsLInJC[iNDimInJC] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
744 Double_t xminLInJC[iNDimInJC] = {fgkfMassLambdaMin, fPtV0Min, -fRangeEtaV0Max, fJetPtMin};
745 Double_t xmaxLInJC[iNDimInJC] = {fgkfMassLambdaMax, fPtV0Max, fRangeEtaV0Max, fJetPtMax};
746
747 // binning eff inclusive vs eta-pT
748 Double_t fStepDeltaEta = 0.1;
749 Double_t fRangeDeltaEtaMax = 0.5;
750 const Int_t iNBinsDeltaEta = 2*Int_t(fRangeDeltaEtaMax/fStepDeltaEta);
751 Int_t binsEtaK[3] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
752 Double_t xminEtaK[3] = {fgkfMassK0sMin, fPtV0Min, -fRangeEtaV0Max};
753 Double_t xmaxEtaK[3] = {fgkfMassK0sMax, fPtV0Max, fRangeEtaV0Max};
754 Int_t binsEtaL[3] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
755 Double_t xminEtaL[3] = {fgkfMassLambdaMin, fPtV0Min, -fRangeEtaV0Max};
756 Double_t xmaxEtaL[3] = {fgkfMassLambdaMax, fPtV0Max, fRangeEtaV0Max};
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] = {fgkfMassK0sMin, fPtV0Min, -fRangeEtaV0Max, fJetPtMin, -fRangeDeltaEtaMax};
761 Double_t xmaxEtaKInRec[5] = {fgkfMassK0sMax, fPtV0Max, fRangeEtaV0Max, fJetPtMax, fRangeDeltaEtaMax};
762 Int_t binsEtaLInRec[5] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
763 Double_t xminEtaLInRec[5] = {fgkfMassLambdaMin, fPtV0Min, -fRangeEtaV0Max, fJetPtMin, -fRangeDeltaEtaMax};
764 Double_t xmaxEtaLInRec[5] = {fgkfMassLambdaMax, fPtV0Max, fRangeEtaV0Max, fJetPtMax, fRangeDeltaEtaMax};
765 // generated
766 Int_t binsEtaInGen[4] = {iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
767 Double_t xminEtaInGen[4] = {fPtV0Min, -fRangeEtaV0Max, fJetPtMin, -fRangeDeltaEtaMax};
768 Double_t xmaxEtaInGen[4] = {fPtV0Max, fRangeEtaV0Max, fJetPtMax, fRangeDeltaEtaMax};
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, fPtV0Min, -fRangeEtaV0Max, fPtV0Min, fJetPtMin};
773 Double_t xmaxEtaDaughter[iNDimEtaD] = {2, 1, fPtV0Max, fRangeEtaV0Max, fPtV0Max, fJetPtMax};
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,fgkfMassK0sMin,fgkfMassK0sMax);
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,fgkfMassLambdaMin,fgkfMassLambdaMax);
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,fgkfMassLambdaMin,fgkfMassLambdaMax);
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,fJetPtMin,fJetPtMax,100,0,ffRadiusJet+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,fJetPtMin,fJetPtMax,100,0,ffRadiusJet+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,fJetPtMin,fJetPtMax,100,0,ffRadiusJet+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,fJetPtMin,fJetPtMax);
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,fJetPtMin,fJetPtMax);
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*ffCutVertexZ,1.5*ffCutVertexZ);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,fgkiNBinsMassK0s,fgkfMassK0sMin,fgkfMassK0sMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsEtaV0,-fRangeEtaV0Max,fRangeEtaV0Max);
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,fPtV0Min,fPtV0Max,iNJetPtBins,fJetPtMin,fJetPtMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsPtV0,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,fgkiNBinsMassLambda,fgkfMassLambdaMin,fgkfMassLambdaMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsEtaV0,-fRangeEtaV0Max,fRangeEtaV0Max);
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,fPtV0Min,fPtV0Max,iNJetPtBins,fJetPtMin,fJetPtMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsPtV0,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,fgkiNBinsMassLambda,fgkfMassLambdaMin,fgkfMassLambdaMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsEtaV0,-fRangeEtaV0Max,fRangeEtaV0Max);
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,fPtV0Min,fPtV0Max,iNJetPtBins,fJetPtMin,fJetPtMax);
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max,iNBinsPtV0,fPtV0Min,fPtV0Max);
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] = {fPtV0Min, dPtXiMin, fJetPtMin};
961 Double_t xmaxFD[iNDimFD] = {fPtV0Max, dPtXiMax, fJetPtMax};
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,fPtV0Min,fPtV0Max);
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,fPtV0Min,fPtV0Max);
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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassK0sMin,fgkfMassK0sMax,100,0,2);
1058 fh2CutCos[i] = new TH2D(Form("fh2CutCos_%d",i),"Cuts: CPA;#it{m}_{inv} (GeV/#it{c}^{2});CPA",fgkiNBinsMassK0s,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassK0sMin,fgkfMassK0sMax,1500,0,150);
1060 fh2CutEtaK0s[i] = new TH2D(Form("fh2CutEtaK0s_%d",i),"Cuts: K0s: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta",fgkiNBinsMassK0s,fgkfMassK0sMin,fgkfMassK0sMax,200,-2,2);
1061 fh2CutEtaLambda[i] = new TH2D(Form("fh2CutEtaLambda_%d",i),"Cuts: Lambda: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta",fgkiNBinsMassLambda,fgkfMassLambdaMin,fgkfMassLambdaMax,200,-2,2);
1062 fh2CutRapK0s[i] = new TH2D(Form("fh2CutRapK0s_%d",i),"Cuts: K0s: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y",fgkiNBinsMassK0s,fgkfMassK0sMin,fgkfMassK0sMax,200,-2,2);
1063 fh2CutRapLambda[i] = new TH2D(Form("fh2CutRapLambda_%d",i),"Cuts: Lambda: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y",fgkiNBinsMassLambda,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassK0sMin,fgkfMassK0sMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassLambdaMin,fgkfMassLambdaMax,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,fgkfMassK0sMin,fgkfMassK0sMax);
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,fgkfMassLambdaMin,fgkfMassLambdaMax);
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,fgkfMassLambdaMin,fgkfMassLambdaMax);
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
1196void 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 printf("TaskV0sInJets: No output AOD found\n");
1223 return;
1224 }
1225 if (!fAODIn)
1226 {
1227 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 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 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 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,ffCutVertexZ,ffCutVertexR2,ffCutCentLow,ffCutCentHigh,1,0.1)) // cut on |delta z| in 2011 data between SPD vertex and nominal primary vertex
1278// if (!IsSelectedForJets(fAODIn,ffCutVertexZ,ffCutVertexR2,ffCutCentLow,ffCutCentHigh)) // no need for cutting in 2010 data
1279 {
1280 if(fDebug>5) printf("TaskV0sInJets: Event rejected\n");
1281 return;
1282 }
1283
1284// ffCentrality = fAODIn->GetHeader()->GetCentrality(); // event centrality
1285 ffCentrality = fAODIn->GetHeader()->GetCentralityP()->GetCentralityPercentile("V0M"); // event centrality
1286 Int_t iCentIndex = GetCentralityBinIndex(ffCentrality); // 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 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(ffCentrality);
1305 fh2EventCentTracks->Fill(ffCentrality,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// AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE); // enable AOD output
1321// AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE); // enable AOD output
1322
1323 AliAODv0* v0 = 0; // pointer to V0 candidates
1324// AliV0Object* objectV0 = 0;
1325 TVector3 vecV0Momentum; // 3D vector of V0 momentum
1326 Double_t fMassV0K0s = 0; // invariant mass of the K0s candidate
1327 Double_t fMassV0Lambda = 0; // invariant mass of the Lambda candidate
1328 Double_t fMassV0ALambda = 0; // invariant mass of the Lambda candidate
1329 Int_t iNV0CandTot = 0; // counter of all V0 candidates at the beginning
1330 Int_t iNV0CandK0s = 0; // counter of K0s candidates at the end
1331 Int_t iNV0CandLambda = 0; // counter of Lambda candidates at the end
1332 Int_t iNV0CandALambda = 0; // counter of Lambda candidates at the end
1333
1334 Bool_t bUseOldCuts = 0; // old reconstruction cuts
1335 Bool_t bUseAliceCuts = 0; // cuts used by Alice Zimmermann
1336 Bool_t bUseIouriCuts = 0; // cuts used by Iouri
1337 Bool_t bPrintCuts = 1; // print out which cuts are applied
1338 Bool_t bPrintJetSelection = 0; // print out which jets are selected
1339
1340 // Values of V0 reconstruction cuts:
1341 // Daughter tracks
1342 Int_t iRefit = AliAODTrack::kTPCrefit; // TPC refit for daughter tracks
1343 Double_t fDCAToPrimVtxMin = fdCutDCAToPrimVtxMin; // 0.1; // [cm] min DCA of daughters to the prim vtx
1344 Double_t fDCADaughtersMax = fdCutDCADaughtersMax; // 1.; // [sigma of TPC tracking] max DCA between daughters
1345 Double_t fEtaDaughterMax = 0.8; // max |pseudorapidity| of daughter tracks
1346 Double_t fNSigmadEdxMax = fdCutNSigmadEdxMax;// 3.; // [sigma dE/dx] max difference between measured and expected signal of dE/dx in the TPC
1347 Double_t fPtProtonPIDMax = 1.; // [GeV/c] maxium pT of proton for applying PID cut
1348 // V0 candidate
1349 Bool_t bOnFly = 0; // on-the-fly (yes) or offline (no) reconstructed
1350 Double_t fCPAMin = fdCutCPAMin;// 0.998; // min cosine of the pointing angle
1351 Double_t fRadiusDecayMin = 5.; // [cm] min radial distance of the decay vertex
1352 Double_t fRadiusDecayMax = 100.; // [cm] max radial distance of the decay vertex
1353 Double_t fEtaMax = 0.7; // max |pseudorapidity| of V0
1354 Double_t fNTauMax = fdCutNTauMax; // 5.0; // [tau] max proper lifetime in multiples of the mean lifetime
1355
1356 // Old cuts Start
1357 Double_t fNCrossedRowsTPCMin = 70.; // min number of crossed TPC rows (turned off)
1358// Double_t fCrossedRowsOverFindMin = 0.8; // min ratio crossed rows / findable clusters (turned off)
1359// Double_t fCrossedRowsOverFindMax = 1e3; // max ratio crossed rows / findable clusters (turned off)
1360 Double_t fPtDaughterMin = 0.150; // [GeV/c] min transverse momentum of daughter tracks (turned off)
1361 Double_t fRapMax = 0.75; // max |rapidity| of V0 (turned off)
1362 // Old cuts End
1363
1364 // Other cuts
1365 Double_t fNSigmaMassMax = 3.; // [sigma m] max difference between candidate mass and real particle mass (used only for mass peak method of signal extraction)
1366 Double_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
1367
1368 // Selection of active cuts
1369 Bool_t bCutEtaDaughter = 1; // daughter pseudorapidity
1370 Bool_t bCutRapV0 = 0; // V0 rapidity
1371 Bool_t bCutEtaV0 = 1; // V0 pseudorapidity
1372 Bool_t bCutTau = 1; // V0 lifetime
1373 Bool_t bCutPid = 1; // PID (TPC dE/dx)
1374 Bool_t bCutArmPod = 1; // Armenteros-Podolanski for K0S
1375// Bool_t bCutCross = 0; // cross contamination
1376
1377 if (bUseOldCuts)
1378 {
1379 bCutRapV0 = 1;
1380 fEtaMax = 0.75;
1381 fNTauMax = 3.0;
1382 }
1383 else if (bUseAliceCuts)
1384 {
1385// bOnFly = 1;
1386 fEtaMax = 0.75;
1387 fNTauMax = 5.0;
1388 }
1389 else if (bUseIouriCuts)
1390 {
1391 bCutRapV0 = 1;
1392 bCutEtaV0 = 0;
1393 fNTauMax = 3.0;
1394 fRapMax = 0.5;
1395 }
1396
1397 Double_t fCTauK0s = 2.6844; // [cm] c tau of K0S
1398 Double_t fCTauLambda = 7.89; // [cm] c tau of Lambda
1399
1400 // Load PDG values of particle masses
1401 Double_t fMassK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1402 Double_t fMassLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1403
1404 // PDG codes of used particles
1405 Int_t iPdgCodePion = 211;
1406 Int_t iPdgCodeProton = 2212;
1407 Int_t iPdgCodeK0s = 310;
1408 Int_t iPdgCodeLambda = 3122;
1409
1410 // Jet selection: ffCutPtJetMin, ffCutPtTrackMin
1411 Double_t fJetEtaWindow = fEtaMax-ffRadiusJet; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1412 Double_t fCutJetAreaMin = 0.6*TMath::Pi()*ffRadiusJet*ffRadiusJet; // minimum jet area
1413 Double_t dRadiusExcludeCone = 2*ffRadiusJet; // radius of cones around jets excluded for V0 outside jets
1414 Bool_t bLeadingJetOnly = 0;
1415
1416 if (bUseAliceCuts)
1417 {
1418 ffCutPtJetMin = 5;
1419 ffCutPtTrackMin = 5;
1420 fCutJetAreaMin = 0;
1421 bLeadingJetOnly = 0;
1422 }
1423
1424// Int_t iNJetAll = 0; // number of reconstructed jets in fBranchJet
1425// iNJetAll = fBranchJet->GetEntriesFast(); // number of reconstructed jets
1426 TClonesArray* jetArray = 0; // object where the input jets are stored
1427 Int_t iNJet = 0; // number of reconstructed jets in the input
1428 TClonesArray* jetArraySel = new TClonesArray("AliAODJet",0); // object where the selected jets are copied
1429 Int_t iNJetSel = 0; // number of selected reconstructed jets
1430// iNJetSel = jetArraySel->GetEntriesFast(); // number of selected reconstructed jets
1431 TClonesArray* jetArrayPerp = new TClonesArray("AliAODJet",0); // object where the perp. cones are stored
1432 Int_t iNJetPerp = 0; // number of perpendicular cones
1433
1434 AliAODJet* jet = 0; // pointer to a jet
1435// AliJetObject* objectJet = 0;
1436 AliAODJet* jetPerp = 0; // pointer to a perp. cone
1437 AliAODJet* jetRnd = 0; // pointer to a rand. cone
1438 TVector3 vecJetMomentum; // 3D vector of jet momentum
1439// TVector3 vecPerpConeMomentum; // 3D vector of perpendicular cone momentum
1440// TVector3 vecRndConeMomentum; // 3D vector of random cone momentum
1441 Bool_t bJetEventGood = kTRUE; // indicator of good jet events
1442
1443// printf("iNJetAll, iNJetSel: %d %d\n",iNJetAll,iNJetSel);
1444
1445 if (fbJetSelection) // analysis of V0s in jets is switched on
1446 {
1447 jetArray = (TClonesArray*)(fAODOut->FindListObject(fsJetBranchName.Data())); // find object with jets in the output AOD
1448 if (!jetArray)
1449 {
1450 printf("TaskV0sInJets: No array of name: %s\n",fsJetBranchName.Data());
1451 bJetEventGood = kFALSE;
1452 }
1453 if (bJetEventGood)
1454 iNJet = jetArray->GetEntriesFast();
1455 if (bJetEventGood && !iNJet) // check whether there are some jets
1456 {
1457 printf("TaskV0sInJets: No jets in array\n");
1458 bJetEventGood = kFALSE;
1459 }
1460 }
1461 else // no in-jet analysis
1462 bJetEventGood = kFALSE;
1463
1464 // select good jets and copy them to another array
1465 if (bJetEventGood)
1466 {
1467 if (bLeadingJetOnly)
1468 iNJet = 1; // only leading jets
1469 printf("TaskV0sInJets: Jet selection for %d jets\n",iNJet);
1470 for (Int_t iJet = 0; iJet<iNJet; iJet++)
1471 {
1472 AliAODJet* jetSel = (AliAODJet*)jetArray->At(iJet); // load a jet in the list
1473 if (!jetSel)
1474 {
1475 printf("TaskV0sInJets: Cannot load jet %d\n",iJet);
1476 continue;
1477 }
1478 if (bPrintJetSelection)
1479 printf("jet: i = %d, pT = %f, eta = %f, phi = %f, pt lead tr = %f ",iJet,jetSel->Pt(),jetSel->Eta(),jetSel->Phi(),jetSel->GetPtLeading());
1480// printf("TaskV0sInJets: Checking pt > %.2f for jet %d with pt %.2f\n",ffCutPtJetMin,iJet,jetSel->Pt());
1481 if (jetSel->Pt() < ffCutPtJetMin) // selection of high-pt jets
1482 {
1483 if (bPrintJetSelection)
1484 printf("rejected (pt)\n");
1485 continue;
1486 }
1487// printf("TaskV0sInJets: Checking |eta| < %.2f for jet %d with |eta| %.2f\n",fJetEtaWindow,iJet,TMath::Abs(jetSel->Eta()));
1488 if (TMath::Abs(jetSel->Eta()) > fJetEtaWindow) // selection of jets in the chosen pseudorapidity range
1489 {
1490 if (bPrintJetSelection)
1491 printf("rejected (eta)\n");
1492 continue;
1493 }
1494 if (!bUseOldCuts)
1495 {
1496 if (jetSel->EffectiveAreaCharged() < fCutJetAreaMin)
1497 continue;
1498 }
1499 Int_t iNTracksInJet = 0;
1500 Double_t fPtLeadTrack = 0; // pt of the leading track
1501// Int_t iLeadTrack = 0;
1502 iNTracksInJet = jetSel->GetRefTracks()->GetEntriesFast(); // number od tracks that constitute the jet
1503// printf("TaskV0sInJets: Searching for leading track from %d tracks in jet %d\n",iNTracksInJet,iJet);
1504 if (ffCutPtTrackMin > 0) // a positive min leading track pt is set
1505 {
1506 for (Int_t j = 0; j < iNTracksInJet; j++) // find the track with the highest pt
1507 {
1508 AliAODTrack* track = (AliAODTrack*)jetSel->GetTrack(j); // is this the leading track?
1509 if (!track)
1510 continue;
1511// printf("TaskV0sInJets: %d: %.2f\n",j,track->Pt());
1512 if (track->Pt() > fPtLeadTrack)
1513 {
1514 fPtLeadTrack = track->Pt();
1515// iLeadTrack = j;
1516 }
1517 }
1518// printf("Leading track pT: my: %f, ali: %f\n",fPtLeadTrack,jetSel->GetPtLeading());
1519// printf("TaskV0sInJets: Checking leading track pt > %.2f for pt %.2f of track %d in jet %d\n",ffCutPtTrackMin,fPtLeadTrack,iLeadTrack,iJet);
1520 if (fPtLeadTrack < ffCutPtTrackMin) // selection of high-pt jet-track events
1521 {
1522 if (bPrintJetSelection)
1523 printf("rejected (track pt)\n");
1524 continue;
1525 }
1526 }
1527 if (bPrintJetSelection)
1528 printf("accepted\n");
1529 if(fDebug>5) printf("TaskV0sInJets: Jet %d with pt %.2f passed selection\n",iJet,jetSel->Pt());
1530/*
1531 if (fbTreeOutput)
1532 {
1533// new ((*fBranchJet)[iNJetAll++]) AliAODJet(*((AliAODJet*)jetSel));
1534 objectJet = new ((*fBranchJet)[iNJetAll++]) AliJetObject(jetSel); // copy selected jet to the array
1535// objectJet->SetPtLeadingTrack(fPtLeadTrack);
1536 objectJet->SetRadius(ffRadiusJet);
1537 objectJet = 0;
1538 }
1539*/
1540 TLorentzVector vecPerpPlus(*(jetSel->MomentumVector()));
1541 vecPerpPlus.RotateZ(TMath::Pi()/2.); // rotate vector by 90 deg around z
1542 TLorentzVector vecPerpMinus(*(jetSel->MomentumVector()));
1543 vecPerpMinus.RotateZ(-TMath::Pi()/2.); // rotate vector by -90 deg around z
1544// AliAODJet jetTmp = AliAODJet(vecPerp);
1545 if(fDebug>5) printf("TaskV0sInJets: Adding perp. cones number %d, %d\n",iNJetPerp,iNJetPerp+1);
1546// 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());
1547 new ((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpPlus); // write perp. cone to the array
1548 new ((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpMinus); // write perp. cone to the array
1549 if(fDebug>5) printf("TaskV0sInJets: Adding jet number %d\n",iNJetSel);
1550// printf("TaskV0sInJets: Adding jet number %d: pT %f, phi %f, eta %f\n",iNJetSel,jetSel->Pt(),jetSel->Phi(),jetSel->Eta());
1551 new ((*jetArraySel)[iNJetSel++]) AliAODJet(*((AliAODJet*)jetSel)); // copy selected jet to the array
1552 }
1553 if(fDebug>5) printf("TaskV0sInJets: Added jets: %d\n",iNJetSel);
1554 iNJetSel = jetArraySel->GetEntriesFast();
1555 printf("TaskV0sInJets: Selected jets in array: %d\n",iNJetSel);
1556 fh1NJetPerEvent[iCentIndex]->Fill(iNJetSel);
1557 // fill jet spectra
1558 for (Int_t iJet = 0; iJet<iNJetSel; iJet++)
1559 {
1560 jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
1561 fh1PtJet[iCentIndex]->Fill(jet->Pt()); // pt spectrum of selected jets
1562 fh1EtaJet[iCentIndex]->Fill(jet->Eta()); // eta spectrum of selected jets
1563 fh2EtaPtJet[iCentIndex]->Fill(jet->Eta(),jet->Pt()); // eta-pT spectrum of selected jets
1564 fh1PhiJet[iCentIndex]->Fill(jet->Phi()); // phi spectrum of selected jets
1565 Double_t dAreaExcluded = TMath::Pi()*dRadiusExcludeCone*dRadiusExcludeCone; // area of the cone
1566 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fEtaMax-jet->Eta()); // positive eta overhang
1567 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fEtaMax+jet->Eta()); // negative eta overhang
1568 fh1AreaExcluded->Fill(iCentIndex,dAreaExcluded);
1569 }
1570 jet = 0;
1571 }
1572
1573 if (bJetEventGood) // there should be some reconstructed jets
1574 {
1575 fh1EventCounterCut->Fill(4); // events with jet(s)
1576 fh1EventCounterCutCent[iCentIndex]->Fill(4); // events with jet(s)
1577 if (iNJetSel)
1578 {
1579 fh1EventCounterCut->Fill(5); // events with selected jets
1580 fh1EventCounterCutCent[iCentIndex]->Fill(5);
1581 }
1582 }
1583
1584 if (iNJetSel)
1585 {
1586 jetRnd = GetRandomCone(jetArraySel,fJetEtaWindow,2*ffRadiusJet);
1587 if (jetRnd)
1588 {
1589 fh1NRndConeCent->Fill(iCentIndex);
1590 fh2EtaPhiRndCone[iCentIndex]->Fill(jetRnd->Eta(),jetRnd->Phi());
1591 }
1592 }
1593
1594 // Loading primary vertex info
1595 AliAODVertex* primVtx = fAODIn->GetPrimaryVertex(); // get the primary vertex
1596 Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
1597 primVtx->GetXYZ(dPrimVtxPos);
1598 fh1VtxZ[iCentIndex]->Fill(dPrimVtxPos[2]);
1599 fh2VtxXY[iCentIndex]->Fill(dPrimVtxPos[0],dPrimVtxPos[1]);
1600
1601 /*===== Start of loop over V0 candidates =====*/
1602 printf("TaskV0sInJets: Start of V0 loop\n");
1603 for (Int_t iV0 = 0; iV0 < iNV0s; iV0++)
1604 {
1605 v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1606 if (!v0)
1607 continue;
1608
1609 iNV0CandTot++;
1610
1611 // Initialization of status indicators
1612 Bool_t bIsCandidateK0s = kTRUE; // candidate for K0s
1613 Bool_t bIsCandidateLambda = kTRUE; // candidate for Lambda
1614 Bool_t bIsCandidateALambda = kTRUE; // candidate for Lambda
1615 Bool_t bIsInPeakK0s = kFALSE; // candidate within the K0s mass peak
1616 Bool_t bIsInPeakLambda = kFALSE; // candidate within the Lambda mass peak
1617 Bool_t bIsInConeJet = kFALSE; // candidate within the jet cones
1618 Bool_t bIsInConePerp = kFALSE; // candidate within the perpendicular cone
1619 Bool_t bIsInConeRnd = kFALSE; // candidate within the random cone
1620 Bool_t bIsOutsideCones = kFALSE; // candidate outside excluded cones
1621
1622 // Invariant mass calculation
1623 fMassV0K0s = v0->MassK0Short();
1624 fMassV0Lambda = v0->MassLambda();
1625 fMassV0ALambda = v0->MassAntiLambda();
1626
1627 Int_t iCutIndex = 0; // indicator of current selection step
1628 // 0
1629 // All V0 candidates
1630 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1631 iCutIndex++;
1632
1633 // Skip candidates outside the histogram range
1634 if ( (fMassV0K0s < fgkfMassK0sMin) || (fMassV0K0s >= fgkfMassK0sMax) )
1635 bIsCandidateK0s = kFALSE;
1636 if ( (fMassV0Lambda < fgkfMassLambdaMin) || (fMassV0Lambda >= fgkfMassLambdaMax) )
1637 bIsCandidateLambda = kFALSE;
1638 if ( (fMassV0ALambda < fgkfMassLambdaMin) || (fMassV0ALambda >= fgkfMassLambdaMax) )
1639 bIsCandidateALambda = kFALSE;
1640 if (!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1641 continue;
1642
1643 Double_t fPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
1644 vecV0Momentum = TVector3(v0->Px(),v0->Py(),v0->Pz()); // set the vector of V0 momentum
1645
1646 // Sigma of the mass peak window
1647 Double_t fMassPeakWindowK0s = fNSigmaMassMax*MassPeakSigmaOld(fPtV0,0);
1648 Double_t fMassPeakWindowLambda = fNSigmaMassMax*MassPeakSigmaOld(fPtV0,1);
1649// Double_t fMassPeakWindowK0s = fNSigmaMassMax*MassPeakSigma(iCentIndex,fPtV0,0);
1650// Double_t fMassPeakWindowLambda = fNSigmaMassMax*MassPeakSigma(iCentIndex,fPtV0,1);
1651
1652 // Invariant mass peak selection
1653 if (TMath::Abs(fMassV0K0s-fMassK0s) < fMassPeakWindowK0s)
1654 bIsInPeakK0s = kTRUE;
1655 if (TMath::Abs(fMassV0Lambda-fMassLambda) < fMassPeakWindowLambda)
1656 bIsInPeakLambda = kTRUE;
1657
1658 // Retrieving all relevant properties of the V0 candidate
1659 Bool_t bOnFlyStatus = v0->GetOnFlyStatus(); // online (on fly) reconstructed vs offline reconstructed
1660 const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1661 const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
1662 Double_t fPtDaughterPos = trackPos->Pt(); // transverse momentum of a daughter track
1663 Double_t fPtDaughterNeg = trackNeg->Pt();
1664 Double_t fNRowsPos = trackPos->GetTPCClusterInfo(2,1); // crossed TPC pad rows of a daughter track
1665 Double_t fNRowsNeg = trackNeg->GetTPCClusterInfo(2,1);
1666 Double_t fDCAToPrimVtxPos = TMath::Abs(v0->DcaPosToPrimVertex()); // dca of a daughter to the primary vertex
1667 Double_t fDCAToPrimVtxNeg = TMath::Abs(v0->DcaNegToPrimVertex());
1668 Double_t fDCADaughters = v0->DcaV0Daughters(); // dca between daughters
1669 Double_t fCPA = v0->CosPointingAngle(primVtx); // cosine of the pointing angle
1670 Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
1671// Double_t dSecVtxPos[3] = {v0->DecayVertexV0X(),v0->DecayVertexV0Y(),v0->DecayVertexV0Z()}; // V0 vertex position
1672 v0->GetSecondaryVtx(dSecVtxPos);
1673 Double_t fRadiusDecay = TMath::Sqrt(dSecVtxPos[0]*dSecVtxPos[0] + dSecVtxPos[1]*dSecVtxPos[1]); // distance of the V0 vertex from the z-axis
1674 Double_t fEtaDaughterNeg = trackNeg->Eta(); // = v0->EtaProng(1), pseudorapidity of a daughter track
1675 Double_t fEtaDaughterPos = trackPos->Eta(); // = v0->EtaProng(0)
1676 Double_t fRapK0s = v0->RapK0Short(); // rapidity calculated for K0s assumption
1677 Double_t fRapLambda = v0->RapLambda(); // rapidity calculated for Lambda assumption
1678 Double_t fEta = v0->Eta(); // V0 pseudorapidity
1679// Double_t fPhi = v0->Phi(); // V0 pseudorapidity
1680 Double_t dDecayPath[3];
1681 for (Int_t iPos = 0; iPos < 3; iPos++)
1682 dDecayPath[iPos] = dSecVtxPos[iPos]-dPrimVtxPos[iPos]; // vector of the V0 path
1683 Double_t fDecLen = TMath::Sqrt(dDecayPath[0]*dDecayPath[0]+dDecayPath[1]*dDecayPath[1]+dDecayPath[2]*dDecayPath[2]); // path length L
1684 Double_t fDecLen2D = TMath::Sqrt(dDecayPath[0]*dDecayPath[0]+dDecayPath[1]*dDecayPath[1]); // transverse path length R
1685 Double_t fLOverP = fDecLen/v0->P(); // L/p
1686 Double_t fROverPt = fDecLen2D/fPtV0; // R/pT
1687 Double_t fMLOverPK0s = fMassK0s*fLOverP; // m*L/p = c*(proper lifetime)
1688// Double_t fMLOverPLambda = fMassLambda*fLOverP; // m*L/p
1689 Double_t fMROverPtK0s = fMassK0s*fROverPt; // m*R/pT
1690 Double_t fMROverPtLambda = fMassLambda*fROverPt; // m*R/pT
1691 Double_t fNSigmaPosPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos,AliPID::kPion)); // difference between measured and expected signal of the dE/dx in the TPC
1692 Double_t fNSigmaPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos,AliPID::kProton));
1693 Double_t fNSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg,AliPID::kPion));
1694 Double_t fNSigmaNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg,AliPID::kProton));
1695 Double_t fAlpha = v0->AlphaV0(); // Armenteros-Podolanski alpha
1696 Double_t fPtArm = v0->PtArmV0(); // Armenteros-Podolanski pT
1697 AliAODVertex* prodVtxDaughterPos = (AliAODVertex*)(trackPos->GetProdVertex()); // production vertex of the positive daughter track
1698 Char_t cTypeVtxProdPos = prodVtxDaughterPos->GetType(); // type of the production vertex
1699 AliAODVertex* prodVtxDaughterNeg = (AliAODVertex*)(trackNeg->GetProdVertex()); // production vertex of the negative daughter track
1700 Char_t cTypeVtxProdNeg = prodVtxDaughterNeg->GetType(); // type of the production vertex
1701
1702 fh2Tau3DVs2D[0]->Fill(fPtV0,fLOverP/fROverPt);
1703
1704 // QA histograms before cuts
1705 FillQAHistogramV0(primVtx,v0,0,bIsCandidateK0s,bIsCandidateLambda,bIsInPeakK0s,bIsInPeakLambda);
1706 // Cut vs mass histograms before cuts
1707 if (bIsCandidateK0s)
1708 {
1709 fh2CutTPCRowsK0s[0]->Fill(fMassV0K0s,fNRowsPos);
1710 fh2CutTPCRowsK0s[0]->Fill(fMassV0K0s,fNRowsNeg);
1711 fh2CutPtPosK0s[0]->Fill(fMassV0K0s,fPtDaughterPos);
1712 fh2CutPtNegK0s[0]->Fill(fMassV0K0s,fPtDaughterNeg);
1713 fh2CutDCAVtx[0]->Fill(fMassV0K0s,fDCAToPrimVtxPos);
1714 fh2CutDCAVtx[0]->Fill(fMassV0K0s,fDCAToPrimVtxNeg);
1715 fh2CutDCAV0[0]->Fill(fMassV0K0s,fDCADaughters);
1716 fh2CutCos[0]->Fill(fMassV0K0s,fCPA);
1717 fh2CutR[0]->Fill(fMassV0K0s,fRadiusDecay);
1718 fh2CutEtaK0s[0]->Fill(fMassV0K0s,fEtaDaughterPos);
1719 fh2CutEtaK0s[0]->Fill(fMassV0K0s,fEtaDaughterNeg);
1720 fh2CutRapK0s[0]->Fill(fMassV0K0s,fRapK0s);
1721 fh2CutCTauK0s[0]->Fill(fMassV0K0s,fMROverPtK0s/fCTauK0s);
1722 fh2CutPIDPosK0s[0]->Fill(fMassV0K0s,fNSigmaPosPion);
1723 fh2CutPIDNegK0s[0]->Fill(fMassV0K0s,fNSigmaNegPion);
1724 }
1725 if (bIsCandidateLambda)
1726 {
1727 fh2CutTPCRowsLambda[0]->Fill(fMassV0Lambda,fNRowsPos);
1728 fh2CutTPCRowsLambda[0]->Fill(fMassV0Lambda,fNRowsNeg);
1729 fh2CutPtPosLambda[0]->Fill(fMassV0Lambda,fPtDaughterPos);
1730 fh2CutPtNegLambda[0]->Fill(fMassV0Lambda,fPtDaughterNeg);
1731 fh2CutEtaLambda[0]->Fill(fMassV0Lambda,fEtaDaughterPos);
1732 fh2CutEtaLambda[0]->Fill(fMassV0Lambda,fEtaDaughterNeg);
1733 fh2CutRapLambda[0]->Fill(fMassV0Lambda,fRapLambda);
1734 fh2CutCTauLambda[0]->Fill(fMassV0Lambda,fMROverPtLambda/fCTauLambda);
1735 fh2CutPIDPosLambda[0]->Fill(fMassV0Lambda,fNSigmaPosProton);
1736 fh2CutPIDNegLambda[0]->Fill(fMassV0Lambda,fNSigmaNegPion);
1737 }
1738
1739 /*===== Start of reconstruction cutting =====*/
1740
1741 // 1
1742 // All V0 candidates
1743 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1744 iCutIndex++;
1745
1746 /* Start of global cuts */
1747 // 2
1748 // Reconstruction method
1749 if (bPrintCuts) printf("Rec: Applying cut: Reconstruction method: on-the-fly? %s\n",(bOnFly ? "yes" : "no"));
1750 if (bOnFlyStatus!=bOnFly)
1751 continue;
1752 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1753 iCutIndex++;
1754
1755 // 3
1756 // Tracks TPC OK
1757 if (bPrintCuts) printf("Rec: Applying cut: Correct charge of daughters\n");
1758 if ( !trackNeg || !trackPos )
1759 continue;
1760 if (trackNeg->Charge() == trackPos->Charge()) // daughters have different charge?
1761 continue;
1762 if (trackNeg->Charge() != -1) // daughters have expected charge?
1763 continue;
1764 if (trackPos->Charge() != 1) // daughters have expected charge?
1765 continue;
1766
1767 if (bPrintCuts) printf("Rec: Applying cut: TPC refit: %d\n",iRefit);
1768 if (!trackNeg->IsOn(iRefit)) // TPC refit is ON?
1769 continue;
1770 if (bPrintCuts) printf("Rec: Applying cut: Type of production vertex of daughter: Not %d\n",AliAODVertex::kKink);
1771 if(cTypeVtxProdNeg == AliAODVertex::kKink) // kink daughter rejection
1772 continue;
1773 // Old cuts Start
1774 if (bUseOldCuts)
1775 {
1776 if (bPrintCuts) printf("Rec: Applying cut: Number of TPC rows: > %f\n",fNCrossedRowsTPCMin);
1777 if (fNRowsNeg < fNCrossedRowsTPCMin) // Crossed TPC padrows
1778 continue;
1779// Int_t findable = trackNeg->GetTPCNclsF(); // Findable clusters
1780// if (findable <= 0)
1781// continue;
1782// if (fNRowsNeg/findable < fCrossedRowsOverFindMin)
1783// continue;
1784// if (fNRowsNeg/findable > fCrossedRowsOverFindMax)
1785// continue;
1786 }
1787 // Old cuts End
1788
1789 if (!trackPos->IsOn(iRefit))
1790 continue;
1791 if(cTypeVtxProdPos == AliAODVertex::kKink) // kink daughter rejection
1792 continue;
1793 // Old cuts Start
1794 if (bUseOldCuts)
1795 {
1796 if (fNRowsPos < fNCrossedRowsTPCMin)
1797 continue;
1798// findable = trackPos->GetTPCNclsF();
1799// if (findable <= 0)
1800// continue;
1801// if (fNRowsPos/findable < fCrossedRowsOverFindMin)
1802// continue;
1803// if (fNRowsPos/findable > fCrossedRowsOverFindMax)
1804// continue;
1805 }
1806 // Old cuts End
1807
1808 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1809 iCutIndex++;
1810
1811 // 4
1812 // Daughters: transverse momentum cut
1813 if (bUseOldCuts)
1814 {
1815 if (bPrintCuts) printf("Rec: Applying cut: Daughter pT: > %f\n",fPtDaughterMin);
1816 if ( ( fPtDaughterNeg < fPtDaughterMin ) || ( fPtDaughterPos < fPtDaughterMin ) )
1817 continue;
1818 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1819 }
1820 iCutIndex++;
1821
1822 // 5
1823 // Daughters: Impact parameter of daughters to prim vtx
1824 if (bPrintCuts) printf("Rec: Applying cut: Daughter DCA to prim vtx: > %f\n",fDCAToPrimVtxMin);
1825 if ( ( fDCAToPrimVtxNeg < fDCAToPrimVtxMin ) || ( fDCAToPrimVtxPos < fDCAToPrimVtxMin ) )
1826 continue;
1827 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1828 iCutIndex++;
1829
1830 // 6
1831 // Daughters: DCA
1832 if (bPrintCuts) printf("Rec: Applying cut: DCA between daughters: < %f\n",fDCADaughtersMax);
1833 if (fDCADaughters > fDCADaughtersMax)
1834 continue;
1835 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1836 iCutIndex++;
1837
1838 // 7
1839 // V0: Cosine of the pointing angle
1840 if (bPrintCuts) printf("Rec: Applying cut: CPA: > %f\n",fCPAMin);
1841 if (fCPA < fCPAMin)
1842 continue;
1843 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1844 iCutIndex++;
1845
1846 // 8
1847 // V0: Fiducial volume
1848 if (bPrintCuts) printf("Rec: Applying cut: Decay radius: > %f, < %f\n",fRadiusDecayMin,fRadiusDecayMax);
1849 if ( (fRadiusDecay < fRadiusDecayMin) || (fRadiusDecay > fRadiusDecayMax) )
1850 continue;
1851 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1852 iCutIndex++;
1853
1854 // 9
1855 // Daughters: pseudorapidity cut
1856 if (bCutEtaDaughter)
1857 {
1858 if (bPrintCuts) printf("Rec: Applying cut: Daughter |eta|: < %f\n",fEtaDaughterMax);
1859 if ( (TMath::Abs(fEtaDaughterNeg) > fEtaDaughterMax) || (TMath::Abs(fEtaDaughterPos) > fEtaDaughterMax) )
1860 continue;
1861 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1862 }
1863 iCutIndex++;
1864 /* End of global cuts */
1865
1866 /* Start of particle-dependent cuts */
1867 // 10
1868 // V0: rapidity cut & pseudorapidity cut
1869 if (bCutRapV0)
1870 {
1871 if (bPrintCuts) printf("Rec: Applying cut: V0 |y|: < %f\n",fRapMax);
1872 if (TMath::Abs(fRapK0s) > fRapMax)
1873 bIsCandidateK0s = kFALSE;
1874 if (TMath::Abs(fRapLambda) > fRapMax)
1875 {
1876 bIsCandidateLambda = kFALSE;
1877 bIsCandidateALambda = kFALSE;
1878 }
1879 }
1880 if (bCutEtaV0)
1881 {
1882 if (bPrintCuts) printf("Rec: Applying cut: V0 |eta|: < %f\n",fEtaMax);
1883 if (TMath::Abs(fEta) > fEtaMax)
1884 {
1885 bIsCandidateK0s = kFALSE;
1886 bIsCandidateLambda = kFALSE;
1887 bIsCandidateALambda = kFALSE;
1888 }
1889 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1890 }
1891 iCutIndex++;
1892
1893 // 11
1894 // Lifetime cut
1895 if (bCutTau)
1896 {
1897 if (bPrintCuts) printf("Rec: Applying cut: Proper lifetime: < %f\n",fNTauMax);
1898 if (fMROverPtK0s > fNTauMax*fCTauK0s)
1899 bIsCandidateK0s = kFALSE;
1900 if (fMROverPtLambda > fNTauMax*fCTauLambda)
1901 {
1902 bIsCandidateLambda = kFALSE;
1903 bIsCandidateALambda = kFALSE;
1904 }
1905 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1906 }
1907 iCutIndex++;
1908
1909 // 12
1910 // Daughter PID
1911 if (bCutPid)
1912 {
1913 if (bUseOldCuts)
1914 {
1915 if (bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (both daughters): < %f\n",fNSigmadEdxMax);
1916 if (fNSigmaPosPion > fNSigmadEdxMax || fNSigmaNegPion > fNSigmadEdxMax) // pi+, pi-
1917 bIsCandidateK0s = kFALSE;
1918 if (fNSigmaPosProton > fNSigmadEdxMax || fNSigmaNegPion > fNSigmadEdxMax) // p+, pi-
1919 bIsCandidateLambda = kFALSE;
1920 if (fNSigmaNegProton > fNSigmadEdxMax || fNSigmaPosPion > fNSigmadEdxMax) // p-, pi+
1921 bIsCandidateALambda = kFALSE;
1922 }
1923 else
1924 {
1925 if (bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (proton below %f GeV/c): < %f\n",fPtProtonPIDMax,fNSigmadEdxMax);
1926 if ( (fPtDaughterPos < fPtProtonPIDMax) && (fNSigmaPosProton > fNSigmadEdxMax) ) // p+
1927 bIsCandidateLambda = kFALSE;
1928 if ( (fPtDaughterNeg < fPtProtonPIDMax) && (fNSigmaNegProton > fNSigmadEdxMax) ) // p-
1929 bIsCandidateALambda = kFALSE;
1930 }
1931 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1932 }
1933 iCutIndex++;
1934
1935 Double_t valueCorrel[3] = {fMassV0K0s,fMassV0Lambda,fPtV0};
1936 if (bIsCandidateK0s && bIsCandidateLambda)
1937 fh3CCMassCorrelBoth->Fill(valueCorrel); // correlation of mass distribution of candidates selected as both K0s and Lambda
1938 if (bIsCandidateK0s && !bIsCandidateLambda)
1939 fh3CCMassCorrelKNotL->Fill(valueCorrel); // correlation of mass distribution of candidates selected as K0s and not Lambda
1940 if (!bIsCandidateK0s && bIsCandidateLambda)
1941 fh3CCMassCorrelLNotK->Fill(valueCorrel); // correlation of mass distribution of candidates selected as not K0s and Lambda
1942
1943 // 13
1944 // Armenteros-Podolanski cut
1945 if (bCutArmPod)
1946 {
1947 if (bPrintCuts) printf("Rec: Applying cut: Armenteros-Podolanski (K0S): pT > %f * |alpha|\n",0.2);
1948 if(fPtArm < TMath::Abs(0.2*fAlpha))
1949 bIsCandidateK0s = kFALSE;
1950 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1951 }
1952 iCutIndex++;
1953
1954 // Cross contamination
1955 if (bIsInPeakK0s)
1956 {
1957 if (bIsCandidateLambda) // Lambda candidates in K0s peak, excluded from Lambda candidates by CC cut
1958 fh2CCLambda->Fill(fMassV0Lambda,fPtV0);
1959 }
1960 if (bIsInPeakLambda)
1961 {
1962 if (bIsCandidateK0s) // K0s candidates in Lambda peak, excluded from K0s candidates by CC cut
1963 fh2CCK0s->Fill(fMassV0K0s,fPtV0);
1964 }
1965// if (bCutCross)
1966// {
1967// if (bIsInPeakK0s)
1968// bIsCandidateLambda = kFALSE;
1969// if (bIsInPeakLambda)
1970// bIsCandidateK0s = kFALSE;
1971// FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1972// }
1973// iCutIndex++;
1974
1975 /* End of particle-dependent cuts */
1976
1977 /*===== End of reconstruction cutting =====*/
1978
1979 if (!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1980 continue;
1981
1982/*
1983 if(fDebug>5) printf("TaskV0sInJets: Adding selected V0 to branch\n");
1984 // Add selected candidates to the output tree branch
1985 if ((bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda) && fbTreeOutput)
1986 {
1987 objectV0 = new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliV0Object(v0,primVtx);
1988// new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliAODv0(*((AliAODv0*)v0));
1989 objectV0->SetIsCandidateK0S(bIsCandidateK0s);
1990 objectV0->SetIsCandidateLambda(bIsCandidateLambda);
1991 objectV0->SetIsCandidateALambda(bIsCandidateALambda);
1992 objectV0->SetNSigmaPosProton(fNSigmaPosProton);
1993 objectV0->SetNSigmaNegProton(fNSigmaNegProton);
1994 }
1995*/
1996
1997 // Selection of V0s in jet cones, perpendicular cones, random cones, outside cones
1998 if (bJetEventGood && iNJetSel && (bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda))
1999 {
2000 // Selection of V0s in jet cones
2001 if(fDebug>5) printf("TaskV0sInJets: Searching for V0 %d %d in %d jet cones\n",bIsCandidateK0s,bIsCandidateLambda,iNJetSel);
2002 for (Int_t iJet = 0; iJet<iNJetSel; iJet++)
2003 {
2004 jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2005 vecJetMomentum = TVector3(jet->Px(),jet->Py(),jet->Pz()); // set the vector of jet momentum
2006 if(fDebug>5) printf("TaskV0sInJets: Checking if V0 %d %d in jet cone %d\n",bIsCandidateK0s,bIsCandidateLambda,iJet);
2007 if (IsParticleInCone(v0,jet,ffRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2008 {
2009 if(fDebug>5) printf("TaskV0sInJets: V0 %d %d found in jet cone %d\n",bIsCandidateK0s,bIsCandidateLambda,iJet);
2010 bIsInConeJet = kTRUE;
2011 break;
2012 }
2013 }
2014 // Selection of V0s in perp. cones
2015 if(fDebug>5) printf("TaskV0sInJets: Searching for V0 %d %d in %d perp. cones\n",bIsCandidateK0s,bIsCandidateLambda,iNJetSel);
2016 for (Int_t iJet = 0; iJet<iNJetPerp; iJet++)
2017 {
2018 jetPerp = (AliAODJet*)jetArrayPerp->At(iJet); // load a jet in the list
2019 if(fDebug>5) printf("TaskV0sInJets: Checking if V0 %d %d in perp. cone %d\n",bIsCandidateK0s,bIsCandidateLambda,iJet);
2020 if (IsParticleInCone(v0,jetPerp,ffRadiusJet)) // V0 in perp. cone
2021 {
2022 if(fDebug>5) printf("TaskV0sInJets: V0 %d %d found in perp. cone %d\n",bIsCandidateK0s,bIsCandidateLambda,iJet);
2023 bIsInConePerp = kTRUE;
2024 break;
2025 }
2026 }
2027 // Selection of V0s in random cones
2028 if (jetRnd)
2029 {
2030 if(fDebug>5) printf("TaskV0sInJets: Searching for V0 %d %d in the rnd. cone\n",bIsCandidateK0s,bIsCandidateLambda);
2031 if (IsParticleInCone(v0,jetRnd,ffRadiusJet)) // V0 in rnd. cone?
2032 {
2033 if(fDebug>5) printf("TaskV0sInJets: V0 %d %d found in the rnd. cone\n",bIsCandidateK0s,bIsCandidateLambda);
2034 bIsInConeRnd = kTRUE;
2035 }
2036 }
2037 // Selection of V0s outside jet cones
2038 if(fDebug>5) printf("TaskV0sInJets: Searching for V0 %d %d outside jet cones\n",bIsCandidateK0s,bIsCandidateLambda);
2039 if (!OverlapWithJets(jetArraySel,v0,dRadiusExcludeCone)) // V0 oustide jet cones
2040 {
2041 if(fDebug>5) printf("TaskV0sInJets: V0 %d %d found outside jet cones\n",bIsCandidateK0s,bIsCandidateLambda);
2042 bIsOutsideCones = kTRUE;
2043 }
2044 }
2045
2046 // QA histograms after cuts
2047 FillQAHistogramV0(primVtx,v0,1,bIsCandidateK0s,bIsCandidateLambda,bIsInPeakK0s,bIsInPeakLambda);
2048 // Cut vs mass histograms after cuts
2049 if (bIsCandidateK0s)
2050 {
2051 fh2CutTPCRowsK0s[1]->Fill(fMassV0K0s,fNRowsPos);
2052 fh2CutTPCRowsK0s[1]->Fill(fMassV0K0s,fNRowsNeg);
2053 fh2CutPtPosK0s[1]->Fill(fMassV0K0s,fPtDaughterPos);
2054 fh2CutPtNegK0s[1]->Fill(fMassV0K0s,fPtDaughterNeg);
2055 fh2CutDCAVtx[1]->Fill(fMassV0K0s,fDCAToPrimVtxPos);
2056 fh2CutDCAVtx[1]->Fill(fMassV0K0s,fDCAToPrimVtxNeg);
2057 fh2CutDCAV0[1]->Fill(fMassV0K0s,fDCADaughters);
2058 fh2CutCos[1]->Fill(fMassV0K0s,fCPA);
2059 fh2CutR[1]->Fill(fMassV0K0s,fRadiusDecay);
2060 fh2CutEtaK0s[1]->Fill(fMassV0K0s,fEtaDaughterPos);
2061 fh2CutEtaK0s[1]->Fill(fMassV0K0s,fEtaDaughterNeg);
2062 fh2CutRapK0s[1]->Fill(fMassV0K0s,fRapK0s);
2063 fh2CutCTauK0s[1]->Fill(fMassV0K0s,fMROverPtK0s/fCTauK0s);
2064 fh2CutPIDPosK0s[1]->Fill(fMassV0K0s,fNSigmaPosPion);
2065 fh2CutPIDNegK0s[1]->Fill(fMassV0K0s,fNSigmaNegPion);
2066 fh1DeltaZK0s[iCentIndex]->Fill(dDecayPath[2]);
2067 }
2068 if (bIsCandidateLambda)
2069 {
2070 fh2CutTPCRowsLambda[1]->Fill(fMassV0Lambda,fNRowsPos);
2071 fh2CutTPCRowsLambda[1]->Fill(fMassV0Lambda,fNRowsNeg);
2072 fh2CutPtPosLambda[1]->Fill(fMassV0Lambda,fPtDaughterPos);
2073 fh2CutPtNegLambda[1]->Fill(fMassV0Lambda,fPtDaughterNeg);
2074 fh2CutEtaLambda[1]->Fill(fMassV0Lambda,fEtaDaughterPos);
2075 fh2CutEtaLambda[1]->Fill(fMassV0Lambda,fEtaDaughterNeg);
2076 fh2CutRapLambda[1]->Fill(fMassV0Lambda,fRapLambda);
2077 fh2CutCTauLambda[1]->Fill(fMassV0Lambda,fMROverPtLambda/fCTauLambda);
2078 fh2CutPIDPosLambda[1]->Fill(fMassV0Lambda,fNSigmaPosProton);
2079 fh2CutPIDNegLambda[1]->Fill(fMassV0Lambda,fNSigmaNegPion);
2080 fh1DeltaZLambda[iCentIndex]->Fill(dDecayPath[2]);
2081 }
2082
2083 /*===== Start of filling V0 spectra =====*/
2084
2085 Double_t fAngle = TMath::Pi(); // angle between V0 momentum and jet momentum
2086 if (bIsInConeJet)
2087 {
2088 fAngle = vecV0Momentum.Angle(vecJetMomentum);
2089 }
2090
2091 // iCutIndex = 14
2092 if (bIsCandidateK0s)
2093 {
2094 // 14 K0s candidates after cuts
2095// printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,fMassV0K0s,fPtV0,fEta,fPhi);
2096 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex, iCentIndex);
2097 Double_t valueKIncl[3] = {fMassV0K0s,fPtV0,fEta};
2098 fhnV0InclusiveK0s[iCentIndex]->Fill(valueKIncl);
2099 fh1V0InvMassK0sCent[iCentIndex]->Fill(fMassV0K0s);
2100
2101 fh1QACTau2D[1]->Fill(fMROverPtK0s/fCTauK0s);
2102 fh1QACTau3D[1]->Fill(fMLOverPK0s/fCTauK0s);
2103 fh2Tau3DVs2D[1]->Fill(fPtV0,fLOverP/fROverPt);
2104
2105 if (iNJetSel)
2106 {
2107 // 15 K0s in jet events
2108 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex+1, iCentIndex);
2109 }
2110 if (bIsInConeJet)
2111 {
2112 // 16 K0s in jets
2113 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex+2, iCentIndex);
2114 Double_t valueKInJC[4] = {fMassV0K0s,fPtV0,fEta,jet->Pt()};
2115 fhnV0InJetK0s[iCentIndex]->Fill(valueKInJC);
2116 fh2V0PtJetAngleK0s[iCentIndex]->Fill(jet->Pt(),fAngle);
2117 }
2118 if (bIsOutsideCones)
2119 {
2120 Double_t valueKOutJC[3] = {fMassV0K0s,fPtV0,fEta};
2121 fhnV0OutJetK0s[iCentIndex]->Fill(valueKOutJC);
2122 }
2123 if (bIsInConePerp)
2124 {
2125 Double_t valueKInPC[4] = {fMassV0K0s,fPtV0,fEta,jetPerp->Pt()};
2126 fhnV0InPerpK0s[iCentIndex]->Fill(valueKInPC);
2127 }
2128 if (bIsInConeRnd)
2129 {
2130 Double_t valueKInRnd[3] = {fMassV0K0s,fPtV0,fEta};
2131 fhnV0InRndK0s[iCentIndex]->Fill(valueKInRnd);
2132 }
2133 if (!iNJetSel)
2134 {
2135 Double_t valueKNoJet[3] = {fMassV0K0s,fPtV0,fEta};
2136 fhnV0NoJetK0s[iCentIndex]->Fill(valueKNoJet);
2137 }
2138 iNV0CandK0s++;
2139 }
2140 if (bIsCandidateLambda)
2141 {
2142 // 14 Lambda candidates after cuts
2143// printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,fMassV0Lambda,fPtV0,fEta,fPhi);
2144 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex, iCentIndex);
2145 Double_t valueLIncl[3] = {fMassV0Lambda,fPtV0,fEta};
2146 fhnV0InclusiveLambda[iCentIndex]->Fill(valueLIncl);
2147 fh1V0InvMassLambdaCent[iCentIndex]->Fill(fMassV0Lambda);
2148 if (iNJetSel)
2149 {
2150 // 15 Lambda in jet events
2151 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex+1, iCentIndex);
2152 }
2153 if (bIsInConeJet)
2154 {
2155 // 16 Lambda in jets
2156 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex+2, iCentIndex);
2157 Double_t valueLInJC[4] = {fMassV0Lambda,fPtV0,fEta,jet->Pt()};
2158 fhnV0InJetLambda[iCentIndex]->Fill(valueLInJC);
2159 fh2V0PtJetAngleLambda[iCentIndex]->Fill(jet->Pt(),fAngle);
2160 }
2161 if (bIsOutsideCones)
2162 {
2163 Double_t valueLOutJet[3] = {fMassV0Lambda,fPtV0,fEta};
2164 fhnV0OutJetLambda[iCentIndex]->Fill(valueLOutJet);
2165 }
2166 if (bIsInConePerp)
2167 {
2168 Double_t valueLInPC[4] = {fMassV0Lambda,fPtV0,fEta,jetPerp->Pt()};
2169 fhnV0InPerpLambda[iCentIndex]->Fill(valueLInPC);
2170 }
2171 if (bIsInConeRnd)
2172 {
2173 Double_t valueLInRnd[3] = {fMassV0Lambda,fPtV0,fEta};
2174 fhnV0InRndLambda[iCentIndex]->Fill(valueLInRnd);
2175 }
2176 if (!iNJetSel)
2177 {
2178 Double_t valueLNoJet[3] = {fMassV0Lambda,fPtV0,fEta};
2179 fhnV0NoJetLambda[iCentIndex]->Fill(valueLNoJet);
2180 }
2181 iNV0CandLambda++;
2182 }
2183 if (bIsCandidateALambda)
2184 {
2185 // 14 ALambda candidates after cuts
2186// printf("AL: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,fMassV0ALambda,fPtV0,fEta,fPhi);
2187 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex, iCentIndex);
2188 Double_t valueALIncl[3] = {fMassV0ALambda,fPtV0,fEta};
2189 fhnV0InclusiveALambda[iCentIndex]->Fill(valueALIncl);
2190 fh1V0InvMassALambdaCent[iCentIndex]->Fill(fMassV0ALambda);
2191 if (iNJetSel)
2192 {
2193 // 15 ALambda in jet events
2194 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex+1, iCentIndex);
2195 }
2196 if (bIsInConeJet)
2197 {
2198 // 16 ALambda in jets
2199 FillCandidates(fMassV0K0s, fMassV0Lambda, fMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex+2, iCentIndex);
2200 Double_t valueLInJC[4] = {fMassV0ALambda,fPtV0,fEta,jet->Pt()};
2201 fhnV0InJetALambda[iCentIndex]->Fill(valueLInJC);
2202 fh2V0PtJetAngleALambda[iCentIndex]->Fill(jet->Pt(),fAngle);
2203 }
2204 if (bIsOutsideCones)
2205 {
2206 Double_t valueALOutJet[3] = {fMassV0ALambda,fPtV0,fEta};
2207 fhnV0OutJetALambda[iCentIndex]->Fill(valueALOutJet);
2208 }
2209 if (bIsInConePerp)
2210 {
2211 Double_t valueLInPC[4] = {fMassV0ALambda,fPtV0,fEta,jetPerp->Pt()};
2212 fhnV0InPerpALambda[iCentIndex]->Fill(valueLInPC);
2213 }
2214 if (bIsInConeRnd)
2215 {
2216 Double_t valueALInRnd[3] = {fMassV0ALambda,fPtV0,fEta};
2217 fhnV0InRndALambda[iCentIndex]->Fill(valueALInRnd);
2218 }
2219 if (!iNJetSel)
2220 {
2221 Double_t valueALNoJet[3] = {fMassV0ALambda,fPtV0,fEta};
2222 fhnV0NoJetALambda[iCentIndex]->Fill(valueALNoJet);
2223 }
2224 iNV0CandALambda++;
2225 }
2226 /*===== End of filling V0 spectra =====*/
2227
2228
2229 /*===== Association of reconstructed V0 candidates with MC particles =====*/
2230 if (fbMCAnalysis)
2231 {
2232 // Associate selected candidates only
2233// if ( !(bIsCandidateK0s && bIsInPeakK0s) && !(bIsCandidateLambda && bIsInPeakLambda) ) // signal candidates
2234 if ( !(bIsCandidateK0s) && !(bIsCandidateLambda) && !(bIsCandidateALambda) ) // chosen candidates with any mass
2235 continue;
2236
2237 // Get MC labels of reconstructed daughter tracks
2238 Int_t iLabelPos = TMath::Abs(trackPos->GetLabel());
2239 Int_t iLabelNeg = TMath::Abs(trackNeg->GetLabel());
2240
2241 // Make sure MC daughters are in the array range
2242 if ( (iLabelNeg<0) || (iLabelNeg>=iNTracksMC) || (iLabelPos<0) || (iLabelPos>=iNTracksMC) )
2243 continue;
2244
2245 // Get MC particles corresponding to reconstructed daughter tracks
2246 AliAODMCParticle* particleMCDaughterNeg = (AliAODMCParticle*)arrayMC->At(iLabelNeg);
2247 AliAODMCParticle* particleMCDaughterPos = (AliAODMCParticle*)arrayMC->At(iLabelPos);
2248 if (!particleMCDaughterNeg || !particleMCDaughterPos)
2249 continue;
2250
2251 // Make sure MC daughter particles are not physical primary
2252 if ( (particleMCDaughterNeg->IsPhysicalPrimary()) || (particleMCDaughterPos->IsPhysicalPrimary()) )
2253 continue;
2254
2255 // Get identities of MC daughter particles
2256 Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
2257 Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
2258
2259 // Get index of the mother particle for each MC daughter particle
2260 Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
2261 Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
2262
2263 if ( (iIndexMotherNeg<0) || (iIndexMotherNeg>=iNTracksMC) || (iIndexMotherPos<0) || (iIndexMotherPos>=iNTracksMC) )
2264 continue;
2265
2266 // Check whether MC daughter particles have the same mother
2267 if (iIndexMotherNeg != iIndexMotherPos)
2268 continue;
2269
2270 // Get the MC mother particle of both MC daughter particles
2271 AliAODMCParticle* particleMCMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherPos);
2272 if (!particleMCMother)
2273 continue;
2274
2275 // Get identity of the MC mother particle
2276 Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
2277
2278 // Skip not interesting particles
2279 if ( (iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda) )
2280 continue;
2281
2282 // Check identity of the MC mother particle and the decay channel
2283 // Is MC mother particle K0S?
2284 Bool_t bV0MCIsK0s = ( (iPdgCodeMother==iPdgCodeK0s) && (iPdgCodeDaughterPos==+iPdgCodePion) && (iPdgCodeDaughterNeg==-iPdgCodePion) );
2285 // Is MC mother particle Lambda?
2286 Bool_t bV0MCIsLambda = ( (iPdgCodeMother==+iPdgCodeLambda) && (iPdgCodeDaughterPos==+iPdgCodeProton) && (iPdgCodeDaughterNeg==-iPdgCodePion) );
2287 // Is MC mother particle anti Lambda?
2288 Bool_t bV0MCIsALambda = ( (iPdgCodeMother==-iPdgCodeLambda) && (iPdgCodeDaughterPos==+iPdgCodePion) && (iPdgCodeDaughterNeg==-iPdgCodeProton) );
2289
2290 Double_t fPtV0MC = particleMCMother->Pt();
2291// Double_t fRapV0MC = particleMCMother->Y();
2292 Double_t fEtaV0MC = particleMCMother->Eta();
2293// Double_t fPhiV0MC = particleMCMother->Phi();
2294
2295 // Is MC mother particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2296// Bool_t bV0MCIsPrimary = particleMCMother->IsPhysicalPrimary();
2297 // Get the MC mother particle of the MC mother particle
2298 Int_t iIndexMotherOfMother = particleMCMother->GetMother();
2299 AliAODMCParticle* particleMCMotherOfMother = 0;
2300 if (iIndexMotherOfMother >= 0)
2301 particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2302 // Get identity of the MC mother particle of the MC mother particle if it exists
2303 Int_t iPdgCodeMotherOfMother = 0;
2304 if (particleMCMotherOfMother)
2305 iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2306 // 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*-)
2307// Bool_t bV0MCComesFromSigma = kFALSE; // Is MC mother particle daughter of a Sigma?
2308// if ( (particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && ( (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) ) )
2309// bV0MCComesFromSigma = kTRUE;
2310 // Should MC mother particle be considered as primary when it is Lambda?
2311// Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2312 // Check if the MC mother particle of the MC mother particle is a Xi (3322 - Xi0, 3312 - Xi-)
2313 Bool_t bV0MCComesFromXi = ( (particleMCMotherOfMother) && ( (iPdgCodeMotherOfMother==3322) || (iPdgCodeMotherOfMother==3312) ) ); // Is MC mother particle daughter of a Xi?
2314 Bool_t bV0MCComesFromAXi = ( (particleMCMotherOfMother) && ( (iPdgCodeMotherOfMother==-3322) || (iPdgCodeMotherOfMother==-3312) ) ); // Is MC mother particle daughter of a anti-Xi?
2315
2316 // Get the distance between production point of the MC mother particle and the primary vertex
2317 Double_t dx = dPrimVtxMCX-particleMCMother->Xv();
2318 Double_t dy = dPrimVtxMCY-particleMCMother->Yv();
2319 Double_t dz = dPrimVtxMCZ-particleMCMother->Zv();
2320 Double_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2321 Bool_t bV0MCIsPrimaryDist = (fDistPrimary < fDistPrimaryMax); // Is close enough to be considered primary-like?
2322
2323 // phi, eta resolution for V0-reconstruction
2324// Double_t fResolutionV0Eta = particleMCMother->Eta()-v0->Eta();
2325// Double_t fResolutionV0Phi = particleMCMother->Phi()-v0->Phi();
2326
2327/*
2328 if (fbTreeOutput)
2329 {
2330 objectV0->SetPtTrue(fPtV0MC);
2331 objectV0->SetEtaTrue(fEtaV0MC);
2332 objectV0->SetPhiTrue(fPhiV0MC);
2333 objectV0->SetPDGCode(iPdgCodeMother);
2334 objectV0->SetPDGCodeMother(iPdgCodeMotherOfMother);
2335 }
2336*/
2337
2338 // K0s
2339// if (bIsCandidateK0s && bIsInPeakK0s) // selected candidates in peak
2340 if (bIsCandidateK0s) // selected candidates with any mass
2341 {
2342// if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2343 if (bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2344 {
2345// if (fbTreeOutput)
2346// objectV0->SetOrigin(1);
2347 fh2V0K0sPtMassMCRec[iCentIndex]->Fill(fPtV0MC,fMassV0K0s);
2348 Double_t valueEtaK[3] = {fMassV0K0s,fPtV0MC,fEtaV0MC};
2349 fh3V0K0sEtaPtMassMCRec[iCentIndex]->Fill(valueEtaK);
2350
2351 Double_t valueEtaDKNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,0};
2352 fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKNeg);
2353 Double_t valueEtaDKPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,0};
2354 fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKPos);
2355
2356 fh2V0K0sMCResolMPt[iCentIndex]->Fill(fMassV0K0s-fMassK0s,fPtV0);
2357 fh2V0K0sMCPtGenPtRec[iCentIndex]->Fill(fPtV0MC,fPtV0);
2358 if (bIsInConeJet) // true V0 associated to a candidate in jet
2359 {
2360 Double_t valueKInJCMC[4] = {fMassV0K0s,fPtV0MC,fEtaV0MC,jet->Pt()};
2361 fh3V0K0sInJetPtMassMCRec[iCentIndex]->Fill(valueKInJCMC);
2362 Double_t valueEtaKIn[5] = {fMassV0K0s,fPtV0MC,fEtaV0MC,jet->Pt(),fEtaV0MC-jet->Eta()};
2363 fh4V0K0sInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaKIn);
2364
2365 Double_t valueEtaDKJCNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2366 fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCNeg);
2367 Double_t valueEtaDKJCPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2368 fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCPos);
2369 }
2370 }
2371 if (bV0MCIsK0s && !bV0MCIsPrimaryDist) // not primary K0s
2372 {
2373// if (fbTreeOutput)
2374// objectV0->SetOrigin(-1);
2375 fh1V0K0sPtMCRecFalse[iCentIndex]->Fill(fPtV0MC);
2376 }
2377 }
2378 // Lambda
2379// if (bIsCandidateLambda && bIsInPeakLambda) // selected candidates in peak
2380 if (bIsCandidateLambda) // selected candidates with any mass
2381 {
2382// if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2383 if (bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2384 {
2385// if (fbTreeOutput)
2386// objectV0->SetOrigin(1);
2387 fh2V0LambdaPtMassMCRec[iCentIndex]->Fill(fPtV0MC,fMassV0Lambda);
2388 Double_t valueEtaL[3] = {fMassV0Lambda,fPtV0MC,fEtaV0MC};
2389 fh3V0LambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaL);
2390
2391 Double_t valueEtaDLNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,0};
2392 fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLNeg);
2393 Double_t valueEtaDLPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,0};
2394 fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLPos);
2395
2396 fh2V0LambdaMCResolMPt[iCentIndex]->Fill(fMassV0Lambda-fMassLambda,fPtV0);
2397 fh2V0LambdaMCPtGenPtRec[iCentIndex]->Fill(fPtV0MC,fPtV0);
2398 if (bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2399 {
2400 Double_t valueLInJCMC[4] = {fMassV0Lambda,fPtV0MC,fEtaV0MC,jet->Pt()};
2401 fh3V0LambdaInJetPtMassMCRec[iCentIndex]->Fill(valueLInJCMC);
2402 Double_t valueEtaLIn[5] = {fMassV0Lambda,fPtV0MC,fEtaV0MC,jet->Pt(),fEtaV0MC-jet->Eta()};
2403 fh4V0LambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaLIn);
2404
2405 Double_t valueEtaDLJCNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2406 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCNeg);
2407 Double_t valueEtaDLJCPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2408 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCPos);
2409 }
2410 }
2411 // Fill the feed-down histograms
2412 if (bV0MCIsLambda && bV0MCComesFromXi)
2413 {
2414// if (fbTreeOutput)
2415// objectV0->SetOrigin(2);
2416 Double_t valueFDLIncl[3] = {fPtV0MC,particleMCMotherOfMother->Pt(),0.};
2417 fhnV0LambdaInclMCFD[iCentIndex]->Fill(valueFDLIncl);
2418 if (bIsInConeRnd)
2419 {
2420 fhnV0LambdaBulkMCFD[iCentIndex]->Fill(valueFDLIncl);
2421 }
2422 if (bIsInConeJet)
2423 {
2424 Double_t valueFDLInJets[3] = {fPtV0MC,particleMCMotherOfMother->Pt(),jet->Pt()};
2425 fhnV0LambdaInJetsMCFD[iCentIndex]->Fill(valueFDLInJets);
2426 }
2427 }
2428 if (bV0MCIsLambda && !bV0MCIsPrimaryDist && !bV0MCComesFromXi) // not primary Lambda
2429 {
2430// if (fbTreeOutput)
2431// objectV0->SetOrigin(-1);
2432 fh1V0LambdaPtMCRecFalse[iCentIndex]->Fill(fPtV0MC);
2433 }
2434 }
2435 // anti-Lambda
2436// if (bIsCandidateALambda && bIsInPeakALambda) // selected candidates in peak
2437 if (bIsCandidateALambda) // selected candidates with any mass
2438 {
2439// if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2440 if (bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2441 {
2442// if (fbTreeOutput)
2443// objectV0->SetOrigin(1);
2444 fh2V0ALambdaPtMassMCRec[iCentIndex]->Fill(fPtV0MC,fMassV0ALambda);
2445 Double_t valueEtaAL[3] = {fMassV0ALambda,fPtV0MC,fEtaV0MC};
2446 fh3V0ALambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaAL);
2447
2448 Double_t valueEtaDALNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,0};
2449 fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALNeg);
2450 Double_t valueEtaDALPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,0};
2451 fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALPos);
2452
2453 fh2V0ALambdaMCResolMPt[iCentIndex]->Fill(fMassV0ALambda-fMassLambda,fPtV0);
2454 fh2V0ALambdaMCPtGenPtRec[iCentIndex]->Fill(fPtV0MC,fPtV0);
2455 if (bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2456 {
2457 Double_t valueALInJCMC[4] = {fMassV0ALambda,fPtV0MC,fEtaV0MC,jet->Pt()};
2458 fh3V0ALambdaInJetPtMassMCRec[iCentIndex]->Fill(valueALInJCMC);
2459 Double_t valueEtaALIn[5] = {fMassV0ALambda,fPtV0MC,fEtaV0MC,jet->Pt(),fEtaV0MC-jet->Eta()};
2460 fh4V0ALambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaALIn);
2461
2462 Double_t valueEtaDALJCNeg[6] = {0,particleMCDaughterNeg->Eta(),particleMCDaughterNeg->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2463 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCNeg);
2464 Double_t valueEtaDALJCPos[6] = {1,particleMCDaughterPos->Eta(),particleMCDaughterPos->Pt(),fEtaV0MC,fPtV0MC,jet->Pt()};
2465 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCPos);
2466 }
2467 }
2468 // Fill the feed-down histograms
2469 if (bV0MCIsALambda && bV0MCComesFromAXi)
2470 {
2471// if (fbTreeOutput)
2472// objectV0->SetOrigin(2);
2473 Double_t valueFDALIncl[3] = {fPtV0MC,particleMCMotherOfMother->Pt(),0.};
2474 fhnV0ALambdaInclMCFD[iCentIndex]->Fill(valueFDALIncl);
2475 if (bIsInConeRnd)
2476 {
2477 fhnV0ALambdaBulkMCFD[iCentIndex]->Fill(valueFDALIncl);
2478 }
2479 if (bIsInConeJet)
2480 {
2481 Double_t valueFDALInJets[3] = {fPtV0MC,particleMCMotherOfMother->Pt(),jet->Pt()};
2482 fhnV0ALambdaInJetsMCFD[iCentIndex]->Fill(valueFDALInJets);
2483 }
2484 }
2485 if (bV0MCIsALambda && !bV0MCIsPrimaryDist && !bV0MCComesFromAXi) // not primary anti Lambda
2486 {
2487// if (fbTreeOutput)
2488// objectV0->SetOrigin(-1);
2489 fh1V0ALambdaPtMCRecFalse[iCentIndex]->Fill(fPtV0MC);
2490 }
2491 }
2492 }
2493 /*===== End Association of reconstructed V0 candidates with MC particles =====*/
2494 }
2495 /*===== End of V0 loop =====*/
2496
2497 fh1V0CandPerEvent->Fill(iNV0CandTot);
2498 fh1V0CandPerEventCentK0s[iCentIndex]->Fill(iNV0CandK0s);
2499 fh1V0CandPerEventCentLambda[iCentIndex]->Fill(iNV0CandLambda);
2500 fh1V0CandPerEventCentALambda[iCentIndex]->Fill(iNV0CandALambda);
2501
2502 printf("TaskV0sInJets: End of V0 loop\n");
2503
2504 // Spectra of generated particles
2505 if (fbMCAnalysis)
2506 {
2507 for (Int_t iPartMC = 0; iPartMC < iNTracksMC; iPartMC++)
2508 {
2509 // Get MC particle
2510 AliAODMCParticle* particleMC = (AliAODMCParticle*)arrayMC->At(iPartMC);
2511 if(!particleMC)
2512 continue;
2513
2514 // Get identity of MC particle
2515 Int_t iPdgCodeParticleMC = particleMC->GetPdgCode();
2516 // Fill Xi spectrum (3322 - Xi0, 3312 - Xi-)
2517// if ( (iPdgCodeParticleMC==3322) || (iPdgCodeParticleMC==3312) )
2518 if ( (iPdgCodeParticleMC==3312) && (TMath::Abs(particleMC->Y())<0.5) )
2519 {
2520// if (fbTreeOutput)
2521// new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2522 fh1V0XiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2523 }
2524 if ( (iPdgCodeParticleMC==-3312) && (TMath::Abs(particleMC->Y())<0.5) )
2525 {
2526// if (fbTreeOutput)
2527// new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2528 fh1V0AXiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2529 }
2530 // Skip not interesting particles
2531 if ( (iPdgCodeParticleMC != iPdgCodeK0s) && (TMath::Abs(iPdgCodeParticleMC) != iPdgCodeLambda) )
2532 continue;
2533
2534 // Check identity of the MC V0 particle
2535 // Is MC V0 particle K0S?
2536 Bool_t bV0MCIsK0s = (iPdgCodeParticleMC==iPdgCodeK0s);
2537 // Is MC V0 particle Lambda?
2538 Bool_t bV0MCIsLambda = (iPdgCodeParticleMC==+iPdgCodeLambda);
2539 // Is MC V0 particle anti-Lambda?
2540 Bool_t bV0MCIsALambda = (iPdgCodeParticleMC==-iPdgCodeLambda);
2541
2542 Double_t fPtMC = particleMC->Pt();
2543 Double_t fRapMC = particleMC->Y();
2544 Double_t fEtaMC = particleMC->Eta();
2545
2546 // V0 rapidity cut
2547 if (bCutRapV0)
2548 {
2549 if (bPrintCuts) printf("Gen: Applying cut: V0 |y|: < %f\n",fRapMax);
2550 if ( (TMath::Abs(fRapMC) > fRapMax) )
2551 continue;
2552 }
2553 // V0 pseudorapidity cut
2554 if (bCutEtaV0)
2555 {
2556 if (bPrintCuts) printf("Gen: Applying cut: V0 |eta|: < %f\n",fEtaMax);
2557 if ( (TMath::Abs(fEtaMC) > fEtaMax) )
2558 continue;
2559 }
2560 /*
2561 // Is MC V0 particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2562 Bool_t bV0MCIsPrimary = particleMC->IsPhysicalPrimary();
2563
2564 // Get the MC mother particle of the MC V0 particle
2565 Int_t iIndexMotherOfMother = particleMC->GetMother();
2566 AliAODMCParticle* particleMCMotherOfMother = 0;
2567 if (iIndexMotherOfMother >= 0)
2568 particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2569 // Get identity of the MC mother particle of the MC V0 particle if it exists
2570 Int_t iPdgCodeMotherOfMother = 0;
2571 if (particleMCMotherOfMother)
2572 iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2573 // Check if the MC mother particle is a physical primary Sigma
2574 Bool_t bV0MCComesFromSigma = kFALSE;
2575 if ((particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) )
2576 bV0MCComesFromSigma = kTRUE;
2577 // Should the MC V0 particle be considered as primary when it is Lambda?
2578 Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2579 */
2580 // Reject non primary particles
2581// if (!bV0MCIsPrimaryLambda)
2582// continue;
2583
2584 // Get the distance between the production point of the MC V0 particle and the primary vertex
2585 Double_t dx = dPrimVtxMCX-particleMC->Xv();
2586 Double_t dy = dPrimVtxMCY-particleMC->Yv();
2587 Double_t dz = dPrimVtxMCZ-particleMC->Zv();
2588 Double_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2589 Bool_t bV0MCIsPrimaryDist = (fDistPrimary < fDistPrimaryMax); // Is close enough to be considered primary-like?
2590
2591 // Check whether the MC V0 particle is in a MC jet
2592 AliAODJet* jetMC = 0;
2593 Bool_t bIsMCV0InJet = kFALSE;
2594 if (iNJetSel)
2595 {
2596 if(fDebug>5) printf("TaskV0sInJets: Searching for gen V0 in %d MC jets\n",iNJetSel);
2597 for (Int_t iJet = 0; iJet<iNJetSel; iJet++)
2598 {
2599 jetMC = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2600 if(fDebug>5) printf("TaskV0sInJets: Checking if gen V0 in MC jet %d\n",iJet);
2601 if (IsParticleInCone(particleMC,jetMC,ffRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2602 {
2603 if(fDebug>5) printf("TaskV0sInJets: gen V0 found in MC jet %d\n",iJet);
2604 bIsMCV0InJet = kTRUE;
2605 break;
2606 }
2607 }
2608 }
2609
2610 // Select only primary-like MC V0 particles
2611 // K0s
2612// if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2613 if (bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2614 {
2615// if (fbTreeOutput)
2616// new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2617 fh1V0K0sPtMCGen[iCentIndex]->Fill(fPtMC);
2618 fh2V0K0sEtaPtMCGen[iCentIndex]->Fill(fPtMC,fEtaMC);
2619 if (bIsMCV0InJet)
2620 {
2621 fh2V0K0sInJetPtMCGen[iCentIndex]->Fill(fPtMC,jetMC->Pt());
2622 Double_t valueEtaKInGen[4] = {fPtMC,fEtaMC,jetMC->Pt(),fEtaMC-jetMC->Eta()};
2623 fh3V0K0sInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaKInGen);
2624 }
2625 }
2626 // Lambda
2627// if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2628 if (bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2629 {
2630// if (fbTreeOutput)
2631// new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2632 fh1V0LambdaPtMCGen[iCentIndex]->Fill(fPtMC);
2633 fh2V0LambdaEtaPtMCGen[iCentIndex]->Fill(fPtMC,fEtaMC);
2634 if (bIsMCV0InJet)
2635 {
2636 fh2V0LambdaInJetPtMCGen[iCentIndex]->Fill(fPtMC,jetMC->Pt());
2637 Double_t valueEtaLInGen[4] = {fPtMC,fEtaMC,jetMC->Pt(),fEtaMC-jetMC->Eta()};
2638 fh3V0LambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaLInGen);
2639 }
2640 }
2641 // anti-Lambda
2642// if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2643 if (bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2644 {
2645// if (fbTreeOutput)
2646// new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2647 fh1V0ALambdaPtMCGen[iCentIndex]->Fill(fPtMC);
2648 fh2V0ALambdaEtaPtMCGen[iCentIndex]->Fill(fPtMC,fEtaMC);
2649 if (bIsMCV0InJet)
2650 {
2651 fh2V0ALambdaInJetPtMCGen[iCentIndex]->Fill(fPtMC,jetMC->Pt());
2652 Double_t valueEtaALInGen[4] = {fPtMC,fEtaMC,jetMC->Pt(),fEtaMC-jetMC->Eta()};
2653 fh3V0ALambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaALInGen);
2654 }
2655 }
2656 }
2657 }
2658
2659// if (fbTreeOutput)
2660// ftreeOut->Fill();
2661
2662 jetArraySel->Delete();
2663 delete jetArraySel;
2664 jetArrayPerp->Delete();
2665 delete jetArrayPerp;
2666 if (jetRnd)
2667 delete jetRnd;
2668 jetRnd = 0;
2669
2670 PostData(1,fOutputListStd);
2671 PostData(2,fOutputListQA);
2672 PostData(3,fOutputListCuts);
2673 PostData(4,fOutputListMC);
2674// if (fbTreeOutput)
2675// PostData(5,ftreeOut);
2676// if(fDebug>5) printf("TaskV0sInJets: UserExec: End\n");
2677}
2678
2679void AliAnalysisTaskV0sInJets::FillQAHistogramV0(AliAODVertex* vtx, const AliAODv0* vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda)
2680{
2681 if (!IsCandK0s && !IsCandLambda)
2682 return;
2683
2684// Double_t fMassK0s = vZero->MassK0Short();
2685// Double_t fMassLambda = vZero->MassLambda();
2686
2687 fh1QAV0Status[iIndexHisto]->Fill(vZero->GetOnFlyStatus());
2688
2689 AliAODTrack* trackNeg=(AliAODTrack*)vZero->GetDaughter(1); // negative track
2690 AliAODTrack* trackPos=(AliAODTrack*)vZero->GetDaughter(0); // positive track
2691
2692 Short_t fTotalCharge = 0;
2693 for (Int_t i = 0; i < 2; i++)
2694 {
2695 AliAODTrack* track = (AliAODTrack*)vZero->GetDaughter(i); // track
2696 // Tracks TPC OK
2697 fh1QAV0TPCRefit[iIndexHisto]->Fill(track->IsOn(AliAODTrack::kTPCrefit));
2698 Double_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
2699 fh1QAV0TPCRows[iIndexHisto]->Fill(nCrossedRowsTPC);
2700 Int_t findable = track->GetTPCNclsF();
2701 fh1QAV0TPCFindable[iIndexHisto]->Fill(findable);
2702 if (findable != 0)
2703 {
2704 fh1QAV0TPCRowsFind[iIndexHisto]->Fill(nCrossedRowsTPC/findable);
2705 }
2706 // Daughters: pseudo-rapidity cut
2707 fh1QAV0Eta[iIndexHisto]->Fill(track->Eta());
2708 if ( (nCrossedRowsTPC > (160./(250.-85.)*(255.*TMath::Abs(tan(track->Theta()))-85.))+20.) && (track->Eta() < 0) && (track->Pt() > 0.15) )
2709// if (IsCandK0s)
2710 {
2711 fh2QAV0EtaRows[iIndexHisto]->Fill(track->Eta(),nCrossedRowsTPC);
2712 fh2QAV0PtRows[iIndexHisto]->Fill(track->Pt(),nCrossedRowsTPC);
2713 fh2QAV0PhiRows[iIndexHisto]->Fill(track->Phi(),nCrossedRowsTPC);
2714 fh2QAV0NClRows[iIndexHisto]->Fill(findable,nCrossedRowsTPC);
2715 fh2QAV0EtaNCl[iIndexHisto]->Fill(track->Eta(),findable);
2716 }
2717
2718 // Daughters: transverse momentum cut
2719 fh1QAV0Pt[iIndexHisto]->Fill(track->Pt());
2720 fTotalCharge+=track->Charge();
2721 }
2722 fh1QAV0Charge[iIndexHisto]->Fill(fTotalCharge);
2723
2724 // Daughters: Impact parameter of daughters to prim vtx
2725 fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaNegToPrimVertex()));
2726 fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaPosToPrimVertex()));
2727// fh2CutDCAVtx[iIndexHisto]->Fill(fMassK0s,TMath::Abs(vZero->DcaNegToPrimVertex()));
2728
2729 // Daughters: DCA
2730 fh1QAV0DCAV0[iIndexHisto]->Fill(vZero->DcaV0Daughters());
2731// fh2CutDCAV0[iIndexHisto]->Fill(fMassK0s,vZero->DcaV0Daughters());
2732
2733 // V0: Cosine of the pointing angle
2734 fh1QAV0Cos[iIndexHisto]->Fill(vZero->CosPointingAngle(vtx));
2735// fh2CutCos[iIndexHisto]->Fill(fMassK0s,vZero->CosPointingAngle(vtx));
2736
2737 // V0: Fiducial volume
2738 Double_t xyz[3];
2739 vZero->GetSecondaryVtx(xyz);
2740 Double_t r2=xyz[0]*xyz[0] + xyz[1]*xyz[1];
2741 fh1QAV0R[iIndexHisto]->Fill(TMath::Sqrt(r2));
2742
2743 Double_t fAlpha = vZero->AlphaV0();
2744 Double_t fPtArm = vZero->PtArmV0();
2745
2746 if (IsCandK0s)
2747 {
2748 if (IsInPeakK0s)
2749 {
2750// fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2751// fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2752 fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(vZero->Eta(),vZero->Pt());
2753 fh2QAV0PtPtK0sPeak[iIndexHisto]->Fill(trackNeg->Pt(),trackPos->Pt());
2754 fh2ArmPodK0s[iIndexHisto]->Fill(fAlpha,fPtArm);
2755 }
2756 fh2QAV0EtaEtaK0s[iIndexHisto]->Fill(trackNeg->Eta(),trackPos->Eta());
2757 fh2QAV0PhiPhiK0s[iIndexHisto]->Fill(trackNeg->Phi(),trackPos->Phi());
2758 fh1QAV0RapK0s[iIndexHisto]->Fill(vZero->RapK0Short());
2759 }
2760
2761 if (IsCandLambda)
2762 {
2763 if (IsInPeakLambda)
2764 {
2765// fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2766// fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2767 fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(vZero->Eta(),vZero->Pt());
2768 fh2QAV0PtPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(),trackPos->Pt());
2769 fh2ArmPodLambda[iIndexHisto]->Fill(fAlpha,fPtArm);
2770 }
2771 fh2QAV0EtaEtaLambda[iIndexHisto]->Fill(trackNeg->Eta(),trackPos->Eta());
2772 fh2QAV0PhiPhiLambda[iIndexHisto]->Fill(trackNeg->Phi(),trackPos->Phi());
2773 fh1QAV0RapLambda[iIndexHisto]->Fill(vZero->RapLambda());
2774 }
2775
2776 fh2ArmPod[iIndexHisto]->Fill(fAlpha,fPtArm);
2777
2778}
2779
2780void 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*/)
2781{
2782 if (isK)
2783 {
2784 fh1V0CounterCentK0s[iCent]->Fill(iCut);
2785 fh1V0InvMassK0sAll[iCut]->Fill(mK);
2786 }
2787 if (isL)
2788 {
2789 fh1V0CounterCentLambda[iCent]->Fill(iCut);
2790 fh1V0InvMassLambdaAll[iCut]->Fill(mL);
2791 }
2792 if (isAL)
2793 {
2794 fh1V0CounterCentALambda[iCent]->Fill(iCut);
2795 fh1V0InvMassALambdaAll[iCut]->Fill(mAL);
2796 }
2797}
2798
2799Bool_t AliAnalysisTaskV0sInJets::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
2800{
2801// decides whether a particle is inside a jet cone
2802 if (!part1 || !part2)
2803 return kFALSE;
2804
2805 TVector3 vecMom2(part2->Px(),part2->Py(),part2->Pz());
2806 TVector3 vecMom1(part1->Px(),part1->Py(),part1->Pz());
2807 Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
2808 if(dR<dRMax) // momentum vectors of part1 and part2 are closer than dRMax
2809 return kTRUE;
2810 return kFALSE;
2811}
2812
2813Bool_t AliAnalysisTaskV0sInJets::OverlapWithJets(const TClonesArray* array, const AliVParticle* part, Double_t dDistance) const
2814{
2815// decides whether a cone overlaps with other jets
2816 if (!part)
2817 {
2818 printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No part\n");
2819 return kFALSE;
2820 }
2821 if (!array)
2822 {
2823 printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No array\n");
2824 return kFALSE;
2825 }
2826 Int_t iNJets = array->GetEntriesFast();
2827 if (iNJets<=0)
2828 {
2829 printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Warning: No jets\n");
2830 return kFALSE;
2831 }
2832 AliVParticle* jet = 0;
2833 for (Int_t iJet=0; iJet<iNJets; iJet++)
2834 {
2835 jet = (AliVParticle*)array->At(iJet);
2836 if (!jet)
2837 {
2838 printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: Failed to load jet %d/%d\n",iJet,iNJets);
2839 continue;
2840 }
2841 if (IsParticleInCone(part,jet,dDistance))
2842 return kTRUE;
2843 }
2844 return kFALSE;
2845}
2846
2847AliAODJet* AliAnalysisTaskV0sInJets::GetRandomCone(const TClonesArray* array, Double_t dEtaMax, Double_t dDistance) const
2848{
2849// generate a random cone which does not overlap with selected jets
2850// printf("Generating random cone...\n");
2851 TLorentzVector vecCone;
2852 AliAODJet* part = 0;
2853 Double_t dEta, dPhi;
2854 Int_t iNTrialsMax = 10;
2855 Bool_t bStatus = kFALSE;
2856 for (Int_t iTry=0; iTry<iNTrialsMax; iTry++)
2857 {
2858// printf("Try %d\n",iTry);
2859 dEta = dEtaMax*(2*fRandom->Rndm()-1.); // random eta in [-dEtaMax,+dEtaMax]
2860 dPhi = TMath::TwoPi()*fRandom->Rndm(); // random phi in [0,2*Pi]
2861 vecCone.SetPtEtaPhiM(1.,dEta,dPhi,0.);
2862 part = new AliAODJet(vecCone);
2863 if (!OverlapWithJets(array,part,dDistance))
2864 {
2865 bStatus = kTRUE;
2866// printf("Success\n");
2867 break;
2868 }
2869 else
2870 delete part;
2871 }
2872 if (!bStatus)
2873 part = 0;
2874 return part;
2875}
2876
2877Double_t AliAnalysisTaskV0sInJets::AreaCircSegment(Double_t dRadius, Double_t dDistance) const
2878{
2879// calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
2880 Double_t dEpsilon = 1e-2;
2881 Double_t dR = dRadius;
2882 Double_t dD = dDistance;
2883 if (TMath::Abs(dR)<dEpsilon)
2884 {
2885 printf("AliAnalysisTaskV0sInJets::AreaCircSegment: Error: Too small radius: %f < %f\n",dR,dEpsilon);
2886 return 0.;
2887 }
2888 if (dD>=dR)
2889 return 0.;
2890 if (dD<=-dR)
2891 return TMath::Pi()*dR*dR;
2892 return dR*dR*TMath::ACos(dD/dR)-dD*TMath::Sqrt(dR*dR-dD*dD);
2893}
2894
2895Bool_t AliAnalysisTaskV0sInJets::IsSelectedForJets(AliAODEvent* fAOD,Double_t fVtxZCut,Double_t fVtxR2Cut,Double_t fCentCutLo,Double_t fCentCutUp,Bool_t bCutDeltaZ,Double_t fDeltaZMax)
2896{
2897// event selection
2898 AliAODVertex* vertex = fAOD->GetPrimaryVertex();
2899 if (!vertex)
2900 return kFALSE;
2901 if (vertex->GetNContributors() < 3)
2902 return kFALSE;
2903 TString vtxTitle(vertex->GetTitle());
2904 if (vtxTitle.Contains("TPCVertex"))
2905 return kFALSE;
2906 Double_t zVertex = vertex->GetZ();
2907 if (TMath::Abs(zVertex) > fVtxZCut)
2908 return kFALSE;
2909 if (bCutDeltaZ)
2910 {
2911 AliAODVertex* vertexSPD = fAOD->GetPrimaryVertexSPD();
2912 if (!vertexSPD)
2913 {
2914// printf("IsSelectedForJets: Error: No SPD vertex\n");
2915 return kFALSE;
2916 }
2917 Double_t zVertexSPD = vertexSPD->GetZ();
2918 if (TMath::Abs(zVertex-zVertexSPD) > fDeltaZMax)
2919 {
2920// printf("IsSelectedForJets: Rejecting event due to delta z = %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
2921 return kFALSE;
2922 }
2923// printf("IsSelectedForJets: Event OK: %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
2924 }
2925 Double_t xVertex = vertex->GetX();
2926 Double_t yVertex = vertex->GetY();
2927 Double_t radiusSq = yVertex*yVertex+xVertex*xVertex;
2928 if (radiusSq > fVtxR2Cut)
2929 return kFALSE;
2930 Double_t centrality;
2931// centrality = fAOD->GetHeader()->GetCentrality();
2932 centrality = fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentile("V0M");
2933 if (centrality < 0)
2934 return kFALSE;
2935 if( (fCentCutUp < 0) || (fCentCutLo < 0) || (fCentCutUp > 100) || (fCentCutLo > 100) || (fCentCutLo > fCentCutUp) )
2936 return kFALSE;
2937 if ( (centrality < fCentCutLo) || (centrality > fCentCutUp) )
2938 return kFALSE;
2939 return kTRUE;
2940}
2941
2942Int_t AliAnalysisTaskV0sInJets::GetCentralityBinIndex(Double_t centrality)
2943{
2944// returns index of the centrality bin corresponding to the provided value of centrality
2945 if (centrality < 0 || centrality > fgkiCentBinRanges[fgkiNBinsCent-1])
2946 return -1;
2947 for (Int_t i = 0; i < fgkiNBinsCent; i++)
2948 {
2949 if (centrality <= fgkiCentBinRanges[i])
2950 return i;
2951 }
2952 return -1;
2953}
2954
2955Int_t AliAnalysisTaskV0sInJets::GetCentralityBinEdge(Int_t index)
2956{
2957// returns the upper edge of the centrality bin corresponding to the provided value of index
2958 if (index < 0 || index >= fgkiNBinsCent)
2959 return -1;
2960 return fgkiCentBinRanges[index];
2961}
2962
2963TString AliAnalysisTaskV0sInJets::GetCentBinLabel(Int_t index)
2964{
2965// get string with centrality range for given bin
2966 TString lowerEdge = ( (index == 0) ? "0" : Form("%d",GetCentralityBinEdge(index-1)));
2967 TString upperEdge = Form("%d",GetCentralityBinEdge(index));
2968 return Form("%s-%s %%",lowerEdge.Data(),upperEdge.Data());
2969}
2970
2971Double_t AliAnalysisTaskV0sInJets::MassPeakSigmaOld(Double_t pt, Int_t particle)
2972{
2973// estimation of the sigma of the invariant-mass peak as a function of pT and particle type
2974 switch (particle)
2975 {
2976 case 0: // K0S
2977 return 0.0044 + 0.0004*(pt - 1.);
2978 break;
2979 case 1: // Lambda
2980 return 0.0023 + 0.00034*(pt - 1.);
2981 break;
2982 default:
2983 return 0;
2984 break;
2985 }
2986}