1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // task for analysis of V0s (K0S, (anti-)Lambda) in charged jets
17 // Author: Vit Kucera (vit.kucera@cern.ch)
23 #include "THnSparse.h"
26 #include "AliAnalysisTask.h"
27 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliAODEvent.h"
31 #include "AliAODTrack.h"
32 #include <TDatabasePDG.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"
44 #include "AliAnalysisTaskV0sInJets.h"
46 ClassImp(AliAnalysisTaskV0sInJets)
48 // upper edges of centrality bins
49 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10, 30, 50, 80}; // Alice Zimmermann
50 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10, 20, 40, 60, 80}; // Vit Kucera, initial binning
51 //const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {5, 10, 20, 40, 60, 80}; // Iouri Belikov, LF analysis
52 const Int_t AliAnalysisTaskV0sInJets::fgkiCentBinRanges[AliAnalysisTaskV0sInJets::fgkiNBinsCent] = {10}; // only central
55 const Double_t AliAnalysisTaskV0sInJets::fgkdBinsPtV0[2] = {0, 12};
56 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtV0 = sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtV0) / sizeof((AliAnalysisTaskV0sInJets::fgkdBinsPtV0)[0]) - 1;
57 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtV0Init = int(((AliAnalysisTaskV0sInJets::fgkdBinsPtV0)[AliAnalysisTaskV0sInJets::fgkiNBinsPtV0] - (AliAnalysisTaskV0sInJets::fgkdBinsPtV0)[0]) / 0.1); // bin width 0.1 GeV/c
59 const Double_t AliAnalysisTaskV0sInJets::fgkdBinsPtJet[2] = {0, 100};
60 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJet = sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet) / sizeof(AliAnalysisTaskV0sInJets::fgkdBinsPtJet[0]) - 1;
61 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsPtJetInit = int(((AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[AliAnalysisTaskV0sInJets::fgkiNBinsPtJet] - (AliAnalysisTaskV0sInJets::fgkdBinsPtJet)[0]) / 5.); // bin width 5 GeV/c
62 // axis: K0S invariant mass
63 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsMassK0s = 300;
64 const Double_t AliAnalysisTaskV0sInJets::fgkdMassK0sMin = 0.35;
65 const Double_t AliAnalysisTaskV0sInJets::fgkdMassK0sMax = 0.65;
66 // axis: Lambda invariant mass
67 const Int_t AliAnalysisTaskV0sInJets::fgkiNBinsMassLambda = 200;
68 const Double_t AliAnalysisTaskV0sInJets::fgkdMassLambdaMin = 1.05;
69 const Double_t AliAnalysisTaskV0sInJets::fgkdMassLambdaMax = 1.25;
72 // Default constructor
73 AliAnalysisTaskV0sInJets::AliAnalysisTaskV0sInJets():
86 fdCutDCAToPrimVtxMin(0.1),
87 fdCutDCADaughtersMax(1.),
88 fdCutNSigmadEdxMax(3),
115 fh1EventCounterCut(0),
118 fh1EventCent2Jets(0),
119 fh1EventCent2NoJets(0),
120 fh2EventCentTracks(0),
121 fh1V0CandPerEvent(0),
128 fh3CCMassCorrelBoth(0),
129 fh3CCMassCorrelKNotL(0),
130 fh3CCMassCorrelLNotK(0)
132 for(Int_t i = 0; i < fgkiNQAIndeces; i++)
134 fh1QAV0Status[i] = 0;
135 fh1QAV0TPCRefit[i] = 0;
136 fh1QAV0TPCRows[i] = 0;
137 fh1QAV0TPCFindable[i] = 0;
138 fh1QAV0TPCRowsFind[i] = 0;
140 fh2QAV0EtaRows[i] = 0;
141 fh2QAV0PtRows[i] = 0;
142 fh2QAV0PhiRows[i] = 0;
143 fh2QAV0NClRows[i] = 0;
144 fh2QAV0EtaNCl[i] = 0;
146 fh2QAV0EtaPtK0sPeak[i] = 0;
147 fh2QAV0EtaEtaK0s[i] = 0;
148 fh2QAV0PhiPhiK0s[i] = 0;
149 fh1QAV0RapK0s[i] = 0;
150 fh2QAV0PtPtK0sPeak[i] = 0;
153 fh2QAV0EtaPtLambdaPeak[i] = 0;
154 fh2QAV0EtaEtaLambda[i] = 0;
155 fh2QAV0PhiPhiLambda[i] = 0;
156 fh1QAV0RapLambda[i] = 0;
157 fh2QAV0PtPtLambdaPeak[i] = 0;
158 fh2ArmPodLambda[i] = 0;
160 fh2QAV0EtaPtALambdaPeak[i] = 0;
161 fh2QAV0EtaEtaALambda[i] = 0;
162 fh2QAV0PhiPhiALambda[i] = 0;
163 fh1QAV0RapALambda[i] = 0;
164 fh2QAV0PtPtALambdaPeak[i] = 0;
165 fh2ArmPodALambda[i] = 0;
168 fh1QAV0Charge[i] = 0;
169 fh1QAV0DCAVtx[i] = 0;
179 fh2CutTPCRowsK0s[i] = 0;
180 fh2CutTPCRowsLambda[i] = 0;
181 fh2CutPtPosK0s[i] = 0;
182 fh2CutPtNegK0s[i] = 0;
183 fh2CutPtPosLambda[i] = 0;
184 fh2CutPtNegLambda[i] = 0;
190 fh2CutEtaLambda[i] = 0;
192 fh2CutRapLambda[i] = 0;
193 fh2CutCTauK0s[i] = 0;
194 fh2CutCTauLambda[i] = 0;
195 fh2CutPIDPosK0s[i] = 0;
196 fh2CutPIDNegK0s[i] = 0;
197 fh2CutPIDPosLambda[i] = 0;
198 fh2CutPIDNegLambda[i] = 0;
203 for(Int_t i = 0; i < fgkiNCategV0; i++)
205 fh1V0InvMassK0sAll[i] = 0;
206 fh1V0InvMassLambdaAll[i] = 0;
207 fh1V0InvMassALambdaAll[i] = 0;
209 for(Int_t i = 0; i < fgkiNBinsCent; i++)
211 fh1EventCounterCutCent[i] = 0;
212 fh1V0CounterCentK0s[i] = 0;
213 fh1V0CounterCentLambda[i] = 0;
214 fh1V0CounterCentALambda[i] = 0;
215 fh1V0CandPerEventCentK0s[i] = 0;
216 fh1V0CandPerEventCentLambda[i] = 0;
217 fh1V0CandPerEventCentALambda[i] = 0;
218 fh1V0InvMassK0sCent[i] = 0;
219 fh1V0InvMassLambdaCent[i] = 0;
220 fh1V0InvMassALambdaCent[i] = 0;
221 fh1V0K0sPtMCGen[i] = 0;
222 fh2V0K0sPtMassMCRec[i] = 0;
223 fh1V0K0sPtMCRecFalse[i] = 0;
224 fh2V0K0sEtaPtMCGen[i] = 0;
225 fh3V0K0sEtaPtMassMCRec[i] = 0;
226 fh2V0K0sInJetPtMCGen[i] = 0;
227 fh3V0K0sInJetPtMassMCRec[i] = 0;
228 fh3V0K0sInJetEtaPtMCGen[i] = 0;
229 fh4V0K0sInJetEtaPtMassMCRec[i] = 0;
230 fh2V0K0sMCResolMPt[i] = 0;
231 fh2V0K0sMCPtGenPtRec[i] = 0;
232 fh1V0LambdaPtMCGen[i] = 0;
233 fh2V0LambdaPtMassMCRec[i] = 0;
234 fh1V0LambdaPtMCRecFalse[i] = 0;
235 fh2V0LambdaEtaPtMCGen[i] = 0;
236 fh3V0LambdaEtaPtMassMCRec[i] = 0;
237 fh2V0LambdaInJetPtMCGen[i] = 0;
238 fh3V0LambdaInJetPtMassMCRec[i] = 0;
239 fh3V0LambdaInJetEtaPtMCGen[i] = 0;
240 fh4V0LambdaInJetEtaPtMassMCRec[i] = 0;
241 fh2V0LambdaMCResolMPt[i] = 0;
242 fh2V0LambdaMCPtGenPtRec[i] = 0;
243 fhnV0LambdaInclMCFD[i] = 0;
244 fhnV0LambdaInJetsMCFD[i] = 0;
245 fhnV0LambdaBulkMCFD[i] = 0;
246 fh1V0XiPtMCGen[i] = 0;
247 fh1V0ALambdaPt[i] = 0;
248 fh1V0ALambdaPtMCGen[i] = 0;
249 fh1V0ALambdaPtMCRec[i] = 0;
250 fh2V0ALambdaPtMassMCRec[i] = 0;
251 fh1V0ALambdaPtMCRecFalse[i] = 0;
252 fh2V0ALambdaEtaPtMCGen[i] = 0;
253 fh3V0ALambdaEtaPtMassMCRec[i] = 0;
254 fh2V0ALambdaInJetPtMCGen[i] = 0;
255 fh2V0ALambdaInJetPtMCRec[i] = 0;
256 fh3V0ALambdaInJetPtMassMCRec[i] = 0;
257 fh3V0ALambdaInJetEtaPtMCGen[i] = 0;
258 fh4V0ALambdaInJetEtaPtMassMCRec[i] = 0;
259 fh2V0ALambdaMCResolMPt[i] = 0;
260 fh2V0ALambdaMCPtGenPtRec[i] = 0;
261 fhnV0ALambdaInclMCFD[i] = 0;
262 fhnV0ALambdaInJetsMCFD[i] = 0;
263 fhnV0ALambdaBulkMCFD[i] = 0;
264 fh1V0AXiPtMCGen[i] = 0;
267 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
268 fhnV0K0sInclDaughterEtaPtPtMCRec[i] = 0;
269 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
270 fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = 0;
271 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
272 fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = 0;
273 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
274 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
275 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
276 fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = 0;
277 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
278 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
281 fhnV0InclusiveK0s[i] = 0;
282 fhnV0InclusiveLambda[i] = 0;
283 fhnV0InclusiveALambda[i] = 0;
285 fhnV0InJetK0s[i] = 0;
286 fhnV0InPerpK0s[i] = 0;
287 fhnV0InRndK0s[i] = 0;
288 fhnV0InMedK0s[i] = 0;
289 fhnV0OutJetK0s[i] = 0;
290 fhnV0NoJetK0s[i] = 0;
291 fhnV0InJetLambda[i] = 0;
292 fhnV0InPerpLambda[i] = 0;
293 fhnV0InRndLambda[i] = 0;
294 fhnV0InMedLambda[i] = 0;
295 fhnV0OutJetLambda[i] = 0;
296 fhnV0NoJetLambda[i] = 0;
297 fhnV0InJetALambda[i] = 0;
298 fhnV0InPerpALambda[i] = 0;
299 fhnV0InRndALambda[i] = 0;
300 fhnV0InMedALambda[i] = 0;
301 fhnV0OutJetALambda[i] = 0;
302 fhnV0NoJetALambda[i] = 0;
304 fh2V0PtJetAngleK0s[i] = 0;
305 fh2V0PtJetAngleLambda[i] = 0;
306 fh2V0PtJetAngleALambda[i] = 0;
308 fh1DCAInLambda[i] = 0;
309 fh1DCAInALambda[i] = 0;
311 fh1DCAOutLambda[i] = 0;
312 fh1DCAOutALambda[i] = 0;
313 // fh1DeltaZK0s[i] = 0;
314 // fh1DeltaZLambda[i] = 0;
315 // fh1DeltaZALambda[i] = 0;
321 fh1NJetPerEvent[i] = 0;
322 fh2EtaPhiRndCone[i] = 0;
323 fh2EtaPhiMedCone[i] = 0;
331 AliAnalysisTaskV0sInJets::AliAnalysisTaskV0sInJets(const char* name):
332 AliAnalysisTaskSE(name),
344 fdCutDCAToPrimVtxMin(0.1),
345 fdCutDCADaughtersMax(1.),
346 fdCutNSigmadEdxMax(3),
351 fsJetBgBranchName(0),
372 fh1EventCounterCut(0),
375 fh1EventCent2Jets(0),
376 fh1EventCent2NoJets(0),
377 fh2EventCentTracks(0),
378 fh1V0CandPerEvent(0),
385 fh3CCMassCorrelBoth(0),
386 fh3CCMassCorrelKNotL(0),
387 fh3CCMassCorrelLNotK(0)
389 for(Int_t i = 0; i < fgkiNQAIndeces; i++)
391 fh1QAV0Status[i] = 0;
392 fh1QAV0TPCRefit[i] = 0;
393 fh1QAV0TPCRows[i] = 0;
394 fh1QAV0TPCFindable[i] = 0;
395 fh1QAV0TPCRowsFind[i] = 0;
397 fh2QAV0EtaRows[i] = 0;
398 fh2QAV0PtRows[i] = 0;
399 fh2QAV0PhiRows[i] = 0;
400 fh2QAV0NClRows[i] = 0;
401 fh2QAV0EtaNCl[i] = 0;
403 fh2QAV0EtaPtK0sPeak[i] = 0;
404 fh2QAV0EtaEtaK0s[i] = 0;
405 fh2QAV0PhiPhiK0s[i] = 0;
406 fh1QAV0RapK0s[i] = 0;
407 fh2QAV0PtPtK0sPeak[i] = 0;
410 fh2QAV0EtaPtLambdaPeak[i] = 0;
411 fh2QAV0EtaEtaLambda[i] = 0;
412 fh2QAV0PhiPhiLambda[i] = 0;
413 fh1QAV0RapLambda[i] = 0;
414 fh2QAV0PtPtLambdaPeak[i] = 0;
415 fh2ArmPodLambda[i] = 0;
417 fh2QAV0EtaPtALambdaPeak[i] = 0;
418 fh2QAV0EtaEtaALambda[i] = 0;
419 fh2QAV0PhiPhiALambda[i] = 0;
420 fh1QAV0RapALambda[i] = 0;
421 fh2QAV0PtPtALambdaPeak[i] = 0;
422 fh2ArmPodALambda[i] = 0;
425 fh1QAV0Charge[i] = 0;
426 fh1QAV0DCAVtx[i] = 0;
436 fh2CutTPCRowsK0s[i] = 0;
437 fh2CutTPCRowsLambda[i] = 0;
438 fh2CutPtPosK0s[i] = 0;
439 fh2CutPtNegK0s[i] = 0;
440 fh2CutPtPosLambda[i] = 0;
441 fh2CutPtNegLambda[i] = 0;
447 fh2CutEtaLambda[i] = 0;
449 fh2CutRapLambda[i] = 0;
450 fh2CutCTauK0s[i] = 0;
451 fh2CutCTauLambda[i] = 0;
452 fh2CutPIDPosK0s[i] = 0;
453 fh2CutPIDNegK0s[i] = 0;
454 fh2CutPIDPosLambda[i] = 0;
455 fh2CutPIDNegLambda[i] = 0;
460 for(Int_t i = 0; i < fgkiNCategV0; i++)
462 fh1V0InvMassK0sAll[i] = 0;
463 fh1V0InvMassLambdaAll[i] = 0;
464 fh1V0InvMassALambdaAll[i] = 0;
466 for(Int_t i = 0; i < fgkiNBinsCent; i++)
468 fh1EventCounterCutCent[i] = 0;
469 fh1V0CounterCentK0s[i] = 0;
470 fh1V0CounterCentLambda[i] = 0;
471 fh1V0CounterCentALambda[i] = 0;
472 fh1V0CandPerEventCentK0s[i] = 0;
473 fh1V0CandPerEventCentLambda[i] = 0;
474 fh1V0CandPerEventCentALambda[i] = 0;
475 fh1V0InvMassK0sCent[i] = 0;
476 fh1V0InvMassLambdaCent[i] = 0;
477 fh1V0InvMassALambdaCent[i] = 0;
478 fh1V0K0sPtMCGen[i] = 0;
479 fh2V0K0sPtMassMCRec[i] = 0;
480 fh1V0K0sPtMCRecFalse[i] = 0;
481 fh2V0K0sEtaPtMCGen[i] = 0;
482 fh3V0K0sEtaPtMassMCRec[i] = 0;
483 fh2V0K0sInJetPtMCGen[i] = 0;
484 fh3V0K0sInJetPtMassMCRec[i] = 0;
485 fh3V0K0sInJetEtaPtMCGen[i] = 0;
486 fh4V0K0sInJetEtaPtMassMCRec[i] = 0;
487 fh2V0K0sMCResolMPt[i] = 0;
488 fh2V0K0sMCPtGenPtRec[i] = 0;
489 fh1V0LambdaPtMCGen[i] = 0;
490 fh2V0LambdaPtMassMCRec[i] = 0;
491 fh1V0LambdaPtMCRecFalse[i] = 0;
492 fh2V0LambdaEtaPtMCGen[i] = 0;
493 fh3V0LambdaEtaPtMassMCRec[i] = 0;
494 fh2V0LambdaInJetPtMCGen[i] = 0;
495 fh3V0LambdaInJetPtMassMCRec[i] = 0;
496 fh3V0LambdaInJetEtaPtMCGen[i] = 0;
497 fh4V0LambdaInJetEtaPtMassMCRec[i] = 0;
498 fh2V0LambdaMCResolMPt[i] = 0;
499 fh2V0LambdaMCPtGenPtRec[i] = 0;
500 fhnV0LambdaInclMCFD[i] = 0;
501 fhnV0LambdaInJetsMCFD[i] = 0;
502 fhnV0LambdaBulkMCFD[i] = 0;
503 fh1V0XiPtMCGen[i] = 0;
504 fh1V0ALambdaPt[i] = 0;
505 fh1V0ALambdaPtMCGen[i] = 0;
506 fh1V0ALambdaPtMCRec[i] = 0;
507 fh2V0ALambdaPtMassMCRec[i] = 0;
508 fh1V0ALambdaPtMCRecFalse[i] = 0;
509 fh2V0ALambdaEtaPtMCGen[i] = 0;
510 fh3V0ALambdaEtaPtMassMCRec[i] = 0;
511 fh2V0ALambdaInJetPtMCGen[i] = 0;
512 fh2V0ALambdaInJetPtMCRec[i] = 0;
513 fh3V0ALambdaInJetPtMassMCRec[i] = 0;
514 fh3V0ALambdaInJetEtaPtMCGen[i] = 0;
515 fh4V0ALambdaInJetEtaPtMassMCRec[i] = 0;
516 fh2V0ALambdaMCResolMPt[i] = 0;
517 fh2V0ALambdaMCPtGenPtRec[i] = 0;
518 fhnV0ALambdaInclMCFD[i] = 0;
519 fhnV0ALambdaInJetsMCFD[i] = 0;
520 fhnV0ALambdaBulkMCFD[i] = 0;
521 fh1V0AXiPtMCGen[i] = 0;
524 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = 0;
525 fhnV0K0sInclDaughterEtaPtPtMCRec[i] = 0;
526 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = 0;
527 fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = 0;
528 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = 0;
529 fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = 0;
530 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
531 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
532 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = 0;
533 fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = 0;
534 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = 0;
535 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = 0;
538 fhnV0InclusiveK0s[i] = 0;
539 fhnV0InclusiveLambda[i] = 0;
540 fhnV0InclusiveALambda[i] = 0;
542 fhnV0InJetK0s[i] = 0;
543 fhnV0InPerpK0s[i] = 0;
544 fhnV0InRndK0s[i] = 0;
545 fhnV0InMedK0s[i] = 0;
546 fhnV0OutJetK0s[i] = 0;
547 fhnV0NoJetK0s[i] = 0;
548 fhnV0InJetLambda[i] = 0;
549 fhnV0InPerpLambda[i] = 0;
550 fhnV0InRndLambda[i] = 0;
551 fhnV0InMedLambda[i] = 0;
552 fhnV0OutJetLambda[i] = 0;
553 fhnV0NoJetLambda[i] = 0;
554 fhnV0InJetALambda[i] = 0;
555 fhnV0InPerpALambda[i] = 0;
556 fhnV0InRndALambda[i] = 0;
557 fhnV0InMedALambda[i] = 0;
558 fhnV0OutJetALambda[i] = 0;
559 fhnV0NoJetALambda[i] = 0;
561 fh2V0PtJetAngleK0s[i] = 0;
562 fh2V0PtJetAngleLambda[i] = 0;
563 fh2V0PtJetAngleALambda[i] = 0;
565 fh1DCAInLambda[i] = 0;
566 fh1DCAInALambda[i] = 0;
568 fh1DCAOutLambda[i] = 0;
569 fh1DCAOutALambda[i] = 0;
570 // fh1DeltaZK0s[i] = 0;
571 // fh1DeltaZLambda[i] = 0;
572 // fh1DeltaZALambda[i] = 0;
578 fh1NJetPerEvent[i] = 0;
579 fh2EtaPhiRndCone[i] = 0;
580 fh2EtaPhiMedCone[i] = 0;
585 // Define input and output slots here
586 // Input slot #0 works with a TChain
587 DefineInput(0, TChain::Class());
588 // Output slot #0 id reserved by the base class for AOD
589 // Output slot #1 writes into a TList container
590 DefineOutput(1, TList::Class());
591 DefineOutput(2, TList::Class());
592 DefineOutput(3, TList::Class());
593 DefineOutput(4, TList::Class());
594 DefineOutput(5, TTree::Class());
597 AliAnalysisTaskV0sInJets::~AliAnalysisTaskV0sInJets()
601 fBranchV0Rec->Delete();
605 fBranchV0Gen->Delete();
609 fBranchJet->Delete();
613 fEventInfo->Delete();
621 void AliAnalysisTaskV0sInJets::UserCreateOutputObjects()
626 fRandom = new TRandom3(0);
629 if (!fBranchV0Rec && fbTreeOutput)
631 // fBranchV0Rec = new TClonesArray("AliAODv0",0);
632 fBranchV0Rec = new TClonesArray("AliV0Object",0);
633 fBranchV0Rec->SetName("branch_V0Rec");
635 if (!fBranchV0Gen && fbTreeOutput)
637 fBranchV0Gen = new TClonesArray("AliAODMCParticle",0);
638 fBranchV0Gen->SetName("branch_V0Gen");
640 if (!fBranchJet && fbTreeOutput)
642 // fBranchJet = new TClonesArray("AliAODJet",0);
643 fBranchJet = new TClonesArray("AliJetObject",0);
644 fBranchJet->SetName("branch_Jet");
646 if (!fEventInfo && fbTreeOutput)
648 fEventInfo = new AliEventInfoObject();
649 fEventInfo->SetName("eventInfo");
651 Int_t dSizeBuffer = 32000; // default 32000
654 ftreeOut = new TTree("treeV0","Tree V0");
655 ftreeOut->Branch("branch_V0Rec",&fBranchV0Rec,dSizeBuffer,2);
656 ftreeOut->Branch("branch_V0Gen",&fBranchV0Gen,dSizeBuffer,2);
657 ftreeOut->Branch("branch_Jet",&fBranchJet,dSizeBuffer,2);
658 ftreeOut->Branch("eventInfo",&fEventInfo,dSizeBuffer,2);
662 fOutputListStd = new TList();
663 fOutputListStd->SetOwner();
664 fOutputListQA = new TList();
665 fOutputListQA->SetOwner();
666 fOutputListCuts = new TList();
667 fOutputListCuts->SetOwner();
668 fOutputListMC = new TList();
669 fOutputListMC->SetOwner();
672 const Int_t iNCategEvent = 6;
673 TString categEvent[iNCategEvent] = {"coll. candid.", "AOD OK", "vtx & cent", "with V0", "with jets", "jet selection"};
674 // labels for stages of V0 selection
675 TString categV0[fgkiNCategV0] = {"all"/*0*/, "mass range"/*1*/, "rec. method"/*2*/, "tracks TPC"/*3*/, "track pt"/*4*/, "DCA prim v"/*5*/, "DCA daughters"/*6*/, "CPA"/*7*/, "volume"/*8*/, "track #it{#eta}"/*9*/, "V0 #it{y} & #it{#eta}"/*10*/, "lifetime"/*11*/, "PID"/*12*/, "Arm.-Pod."/*13*/, "inclusive"/*14*/, "in jet event"/*15*/, "in jet"/*16*/};
677 fh1EventCounterCut = new TH1D("fh1EventCounterCut", "Number of events after filtering;selection filter;counts", iNCategEvent, 0, iNCategEvent);
678 for(Int_t i = 0; i < iNCategEvent; i++)
679 fh1EventCounterCut->GetXaxis()->SetBinLabel(i + 1, categEvent[i].Data());
680 fh1EventCent2 = new TH1D("fh1EventCent2", "Number of events vs centrality;centrality;counts", 100, 0, 100);
681 fh1EventCent2Jets = new TH1D("fh1EventCent2Jets", "Number of sel.-jet events vs centrality;centrality;counts", 100, 0, 100);
682 fh1EventCent2NoJets = new TH1D("fh1EventCent2NoJets", "Number of no-jet events vs centrality;centrality;counts", 100, 0, 100);
683 fh2EventCentTracks = new TH2D("fh2EventCentTracks", "Number of tracks vs centrality;centrality;tracks;counts", 100, 0, 100, 150, 0, 15e3);
684 fh1EventCent = new TH1D("fh1EventCent", "Number of events in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
685 for(Int_t i = 0; i < fgkiNBinsCent; i++)
686 fh1EventCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
687 fh1NRndConeCent = new TH1D("fh1NRndConeCent", "Number of rnd. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
688 for(Int_t i = 0; i < fgkiNBinsCent; i++)
689 fh1NRndConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
690 fh1NMedConeCent = new TH1D("fh1NMedConeCent", "Number of med.-cl. cones in centrality bins;centrality;counts", fgkiNBinsCent, 0, fgkiNBinsCent);
691 for(Int_t i = 0; i < fgkiNBinsCent; i++)
692 fh1NMedConeCent->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
693 fh1AreaExcluded = new TH1D("fh1AreaExcluded", "Area of excluded cones in centrality bins;centrality;area", fgkiNBinsCent, 0, fgkiNBinsCent);
694 for(Int_t i = 0; i < fgkiNBinsCent; i++)
695 fh1AreaExcluded->GetXaxis()->SetBinLabel(i + 1, GetCentBinLabel(i).Data());
696 fOutputListStd->Add(fh1EventCounterCut);
697 fOutputListStd->Add(fh1EventCent);
698 fOutputListStd->Add(fh1EventCent2);
699 fOutputListStd->Add(fh1EventCent2Jets);
700 fOutputListStd->Add(fh1EventCent2NoJets);
701 fOutputListStd->Add(fh1NRndConeCent);
702 fOutputListStd->Add(fh1NMedConeCent);
703 fOutputListStd->Add(fh1AreaExcluded);
704 fOutputListStd->Add(fh2EventCentTracks);
706 fh1V0CandPerEvent = new TH1D("fh1V0CandPerEvent", "Number of all V0 candidates per event;candidates;events", 1000, 0, 1000);
707 fOutputListStd->Add(fh1V0CandPerEvent);
709 for(Int_t i = 0; i < fgkiNBinsCent; i++)
711 fh1EventCounterCutCent[i] = new TH1D(Form("fh1EventCounterCutCent_%d", i), Form("Number of events after filtering, cent %s;selection filter;counts", GetCentBinLabel(i).Data()), iNCategEvent, 0, iNCategEvent);
712 for(Int_t j = 0; j < iNCategEvent; j++)
713 fh1EventCounterCutCent[i]->GetXaxis()->SetBinLabel(j + 1, categEvent[j].Data());
714 fh1V0CandPerEventCentK0s[i] = new TH1D(Form("fh1V0CandPerEventCentK0s_%d", i), Form("Number of selected K0s candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
715 fh1V0CandPerEventCentLambda[i] = new TH1D(Form("fh1V0CandPerEventCentLambda_%d", i), Form("Number of selected Lambda candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
716 fh1V0CandPerEventCentALambda[i] = new TH1D(Form("fh1V0CandPerEventCentALambda_%d", i), Form("Number of selected ALambda candidates per event, cent %s;candidates;events", GetCentBinLabel(i).Data()), 100, 0, 100);
717 fh1V0CounterCentK0s[i] = new TH1D(Form("fh1V0CounterCentK0s_%d", i), Form("Number of K0s candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
718 fh1V0CounterCentLambda[i] = new TH1D(Form("fh1V0CounterCentLambda_%d", i), Form("Number of Lambda candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
719 fh1V0CounterCentALambda[i] = new TH1D(Form("fh1V0CounterCentALambda_%d", i), Form("Number of ALambda candidates after cuts, cent %s;cut;counts", GetCentBinLabel(i).Data()), fgkiNCategV0, 0, fgkiNCategV0);
720 for(Int_t j = 0; j < fgkiNCategV0; j++)
722 fh1V0CounterCentK0s[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
723 fh1V0CounterCentLambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
724 fh1V0CounterCentALambda[i]->GetXaxis()->SetBinLabel(j + 1, categV0[j].Data());
726 fOutputListStd->Add(fh1EventCounterCutCent[i]);
727 fOutputListStd->Add(fh1V0CandPerEventCentK0s[i]);
728 fOutputListStd->Add(fh1V0CandPerEventCentLambda[i]);
729 fOutputListStd->Add(fh1V0CandPerEventCentALambda[i]);
730 fOutputListStd->Add(fh1V0CounterCentK0s[i]);
731 fOutputListStd->Add(fh1V0CounterCentLambda[i]);
732 fOutputListStd->Add(fh1V0CounterCentALambda[i]);
734 // pt binning for V0 and jets
735 Int_t iNBinsPtV0 = fgkiNBinsPtV0Init;
736 Double_t dPtV0Min = fgkdBinsPtV0[0];
737 Double_t dPtV0Max = fgkdBinsPtV0[fgkiNBinsPtV0];
738 Int_t iNJetPtBins = fgkiNBinsPtJetInit;
739 Double_t dJetPtMin = fgkdBinsPtJet[0];
740 Double_t dJetPtMax = fgkdBinsPtJet[fgkiNBinsPtJet];
742 fh2CCK0s = new TH2D("fh2CCK0s", "K0s candidates in Lambda peak", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
743 fh2CCLambda = new TH2D("fh2CCLambda", "Lambda candidates in K0s peak", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, iNBinsPtV0, dPtV0Min, dPtV0Max);
744 Int_t binsCorrel[3] = {fgkiNBinsMassK0s, fgkiNBinsMassLambda, iNBinsPtV0};
745 Double_t xminCorrel[3] = {fgkdMassK0sMin, fgkdMassLambdaMin, dPtV0Min};
746 Double_t xmaxCorrel[3] = {fgkdMassK0sMax, fgkdMassLambdaMax, dPtV0Max};
747 // Int_t binsCorrel[3] = {200, 200, iNBinsPtV0};
748 // Double_t xminCorrel[3] = {0, 0, dPtV0Min};
749 // Double_t xmaxCorrel[3] = {2, 2, dPtV0Max};
750 fh3CCMassCorrelBoth = new THnSparseD("fh3CCMassCorrelBoth", "Mass correlation: K0S && Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
751 fh3CCMassCorrelKNotL = new THnSparseD("fh3CCMassCorrelKNotL", "Mass correlation: K0S, not Lambda;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
752 fh3CCMassCorrelLNotK = new THnSparseD("fh3CCMassCorrelLNotK", "Mass correlation: Lambda, not K0S;m K0S;m Lambda;pT", 3, binsCorrel, xminCorrel, xmaxCorrel);
753 fOutputListQA->Add(fh2CCK0s);
754 fOutputListQA->Add(fh2CCLambda);
755 fOutputListQA->Add(fh3CCMassCorrelBoth);
756 fOutputListQA->Add(fh3CCMassCorrelKNotL);
757 fOutputListQA->Add(fh3CCMassCorrelLNotK);
759 Double_t dStepEtaV0 = 0.025;
760 Double_t dRangeEtaV0Max = 0.8;
761 const Int_t iNBinsEtaV0 = 2 * Int_t(dRangeEtaV0Max / dStepEtaV0);
763 const Int_t iNDimIncl = 3;
764 Int_t binsKIncl[iNDimIncl] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
765 Double_t xminKIncl[iNDimIncl] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
766 Double_t xmaxKIncl[iNDimIncl] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
767 Int_t binsLIncl[iNDimIncl] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
768 Double_t xminLIncl[iNDimIncl] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
769 Double_t xmaxLIncl[iNDimIncl] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
771 const Int_t iNDimInJC = 4;
772 Int_t binsKInJC[iNDimInJC] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
773 Double_t xminKInJC[iNDimInJC] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
774 Double_t xmaxKInJC[iNDimInJC] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
775 Int_t binsLInJC[iNDimInJC] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins};
776 Double_t xminLInJC[iNDimInJC] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin};
777 Double_t xmaxLInJC[iNDimInJC] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax};
779 // binning eff inclusive vs eta-pT
780 Double_t dStepDeltaEta = 0.1;
781 Double_t dRangeDeltaEtaMax = 0.5;
782 const Int_t iNBinsDeltaEta = 2 * Int_t(dRangeDeltaEtaMax / dStepDeltaEta);
783 Int_t binsEtaK[3] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0};
784 Double_t xminEtaK[3] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max};
785 Double_t xmaxEtaK[3] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max};
786 Int_t binsEtaL[3] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0};
787 Double_t xminEtaL[3] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max};
788 Double_t xmaxEtaL[3] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max};
789 // binning eff in jets vs eta-pT
791 Int_t binsEtaKInRec[5] = {fgkiNBinsMassK0s, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
792 Double_t xminEtaKInRec[5] = {fgkdMassK0sMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
793 Double_t xmaxEtaKInRec[5] = {fgkdMassK0sMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
794 Int_t binsEtaLInRec[5] = {fgkiNBinsMassLambda, iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
795 Double_t xminEtaLInRec[5] = {fgkdMassLambdaMin, dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
796 Double_t xmaxEtaLInRec[5] = {fgkdMassLambdaMax, dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
798 Int_t binsEtaInGen[4] = {iNBinsPtV0, iNBinsEtaV0, iNJetPtBins, iNBinsDeltaEta};
799 Double_t xminEtaInGen[4] = {dPtV0Min, -dRangeEtaV0Max, dJetPtMin, -dRangeDeltaEtaMax};
800 Double_t xmaxEtaInGen[4] = {dPtV0Max, dRangeEtaV0Max, dJetPtMax, dRangeDeltaEtaMax};
801 // daughter eta: charge-etaD-ptD-etaV0-ptV0-ptJet
802 const Int_t iNDimEtaD = 6;
803 Int_t binsEtaDaughter[iNDimEtaD] = {2, 20, iNBinsPtV0, iNBinsEtaV0, iNBinsPtV0, iNJetPtBins};
804 Double_t xminEtaDaughter[iNDimEtaD] = {0, -1, dPtV0Min, -dRangeEtaV0Max, dPtV0Min, dJetPtMin};
805 Double_t xmaxEtaDaughter[iNDimEtaD] = {2, 1, dPtV0Max, dRangeEtaV0Max, dPtV0Max, dJetPtMax};
807 for(Int_t i = 0; i < fgkiNBinsCent; i++)
809 fh1V0InvMassK0sCent[i] = new TH1D(Form("fh1V0InvMassK0sCent_%d", i), Form("K0s: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
810 fh1V0InvMassLambdaCent[i] = new TH1D(Form("fh1V0InvMassLambdaCent_%d", i), Form("Lambda: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
811 fh1V0InvMassALambdaCent[i] = new TH1D(Form("fh1V0InvMassALambdaCent_%d", i), Form("ALambda: V0 invariant mass, cent %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", GetCentBinLabel(i).Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
812 fOutputListStd->Add(fh1V0InvMassK0sCent[i]);
813 fOutputListStd->Add(fh1V0InvMassLambdaCent[i]);
814 fOutputListStd->Add(fh1V0InvMassALambdaCent[i]);
816 fhnV0InclusiveK0s[i] = new THnSparseD(Form("fhnV0InclusiveK0s_C%d", i), "K0s: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
817 fhnV0InclusiveLambda[i] = new THnSparseD(Form("fhnV0InclusiveLambda_C%d", i), "Lambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
818 fhnV0InclusiveALambda[i] = new THnSparseD(Form("fhnV0InclusiveALambda_C%d", i), "ALambda: V0 invariant mass vs pt;#it{m}_{inv} (GeV/#it{c}^{2});pt (GeV/#it{c});counts", iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
819 fOutputListStd->Add(fhnV0InclusiveK0s[i]);
820 fOutputListStd->Add(fhnV0InclusiveLambda[i]);
821 fOutputListStd->Add(fhnV0InclusiveALambda[i]);
823 fhnV0InJetK0s[i] = new THnSparseD(Form("fhnV0InJetK0s_%d", i), Form("K0s: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
824 fOutputListStd->Add(fhnV0InJetK0s[i]);
825 fhnV0InPerpK0s[i] = new THnSparseD(Form("fhnV0InPerpK0s_%d", i), Form("K0s: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
826 fOutputListStd->Add(fhnV0InPerpK0s[i]);
827 fhnV0InRndK0s[i] = new THnSparseD(Form("fhnV0InRndK0s_%d", i), Form("K0s: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
828 fOutputListStd->Add(fhnV0InRndK0s[i]);
829 fhnV0InMedK0s[i] = new THnSparseD(Form("fhnV0InMedK0s_%d", i), Form("K0s: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
830 fOutputListStd->Add(fhnV0InMedK0s[i]);
831 fhnV0OutJetK0s[i] = new THnSparseD(Form("fhnV0OutJetK0s_%d", i), Form("K0s: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
832 fOutputListStd->Add(fhnV0OutJetK0s[i]);
833 fhnV0NoJetK0s[i] = new THnSparseD(Form("fhnV0NoJetK0s_%d", i), Form("K0s: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsKIncl, xminKIncl, xmaxKIncl);
834 fOutputListStd->Add(fhnV0NoJetK0s[i]);
835 fhnV0InJetLambda[i] = new THnSparseD(Form("fhnV0InJetLambda_%d", i), Form("Lambda: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
836 fOutputListStd->Add(fhnV0InJetLambda[i]);
837 fhnV0InPerpLambda[i] = new THnSparseD(Form("fhnV0InPerpLambda_%d", i), Form("Lambda: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
838 fOutputListStd->Add(fhnV0InPerpLambda[i]);
839 fhnV0InRndLambda[i] = new THnSparseD(Form("fhnV0InRndLambda_%d", i), Form("Lambda: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
840 fOutputListStd->Add(fhnV0InRndLambda[i]);
841 fhnV0InMedLambda[i] = new THnSparseD(Form("fhnV0InMedLambda_%d", i), Form("Lambda: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
842 fOutputListStd->Add(fhnV0InMedLambda[i]);
843 fhnV0OutJetLambda[i] = new THnSparseD(Form("fhnV0OutJetLambda_%d", i), Form("Lambda: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
844 fOutputListStd->Add(fhnV0OutJetLambda[i]);
845 fhnV0NoJetLambda[i] = new THnSparseD(Form("fhnV0NoJetLambda_%d", i), Form("Lambda: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
846 fOutputListStd->Add(fhnV0NoJetLambda[i]);
847 fhnV0InJetALambda[i] = new THnSparseD(Form("fhnV0InJetALambda_%d", i), Form("ALambda: Mass vs Pt in jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
848 fOutputListStd->Add(fhnV0InJetALambda[i]);
849 fhnV0InPerpALambda[i] = new THnSparseD(Form("fhnV0InPerpALambda_%d", i), Form("ALambda: Mass vs Pt in perp. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
850 fOutputListStd->Add(fhnV0InPerpALambda[i]);
851 fhnV0InRndALambda[i] = new THnSparseD(Form("fhnV0InRndALambda_%d", i), Form("ALambda: Mass vs Pt in rnd. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
852 fOutputListStd->Add(fhnV0InRndALambda[i]);
853 fhnV0InMedALambda[i] = new THnSparseD(Form("fhnV0InMedALambda_%d", i), Form("ALambda: Mass vs Pt in med.-cl. cones, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
854 fOutputListStd->Add(fhnV0InMedALambda[i]);
855 fhnV0OutJetALambda[i] = new THnSparseD(Form("fhnV0OutJetALambda_%d", i), Form("ALambda: Pt outside jets, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
856 fOutputListStd->Add(fhnV0OutJetALambda[i]);
857 fhnV0NoJetALambda[i] = new THnSparseD(Form("fhnV0NoJetALambda_%d", i), Form("ALambda: Pt in jet-less events, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimIncl, binsLIncl, xminLIncl, xmaxLIncl);
858 fOutputListStd->Add(fhnV0NoJetALambda[i]);
860 fh2V0PtJetAngleK0s[i] = new TH2D(Form("fh2V0PtJetAngleK0s_%d", i), Form("K0s: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
861 fOutputListStd->Add(fh2V0PtJetAngleK0s[i]);
862 fh2V0PtJetAngleLambda[i] = new TH2D(Form("fh2V0PtJetAngleLambda_%d", i), Form("Lambda: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
863 fOutputListStd->Add(fh2V0PtJetAngleLambda[i]);
864 fh2V0PtJetAngleALambda[i] = new TH2D(Form("fh2V0PtJetAngleALambda_%d", i), Form("ALambda: #it{p}_{T}^{jet} vs angle V0-jet, cent: %s;#it{p}_{T}^{jet};#it{#alpha}", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax, 100, 0, fdRadiusJet + 0.1);
865 fOutputListStd->Add(fh2V0PtJetAngleALambda[i]);
867 fh1DCAInK0s[i] = new TH1D(Form("fh1DCAInK0s_%d", i), Form("K0s in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
868 fOutputListQA->Add(fh1DCAInK0s[i]);
869 fh1DCAInLambda[i] = new TH1D(Form("fh1DCAInLambda_%d", i), Form("Lambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
870 fOutputListQA->Add(fh1DCAInLambda[i]);
871 fh1DCAInALambda[i] = new TH1D(Form("fh1DCAInALambda_%d", i), Form("ALambda in jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
872 fOutputListQA->Add(fh1DCAInALambda[i]);
874 fh1DCAOutK0s[i] = new TH1D(Form("fh1DCAOutK0s_%d", i), Form("K0s outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
875 fOutputListQA->Add(fh1DCAOutK0s[i]);
876 fh1DCAOutLambda[i] = new TH1D(Form("fh1DCAOutLambda_%d", i), Form("Lambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
877 fOutputListQA->Add(fh1DCAOutLambda[i]);
878 fh1DCAOutALambda[i] = new TH1D(Form("fh1DCAOutALambda_%d", i), Form("ALambda outside jets: DCA daughters, cent %s;DCA (#sigma)", GetCentBinLabel(i).Data()), 50, 0, 1);
879 fOutputListQA->Add(fh1DCAOutALambda[i]);
882 fh1DeltaZK0s[i] = new TH1D(Form("fh1DeltaZK0s_%d", i), Form("K0s: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
883 fOutputListQA->Add(fh1DeltaZK0s[i]);
884 fh1DeltaZLambda[i] = new TH1D(Form("fh1DeltaZLambda_%d", i), Form("Lambda: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
885 fOutputListQA->Add(fh1DeltaZLambda[i]);
886 fh1DeltaZALambda[i] = new TH1D(Form("fh1DeltaZALambda_%d", i), Form("ALambda: #Delta#it{z} vertices, cent %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 50, -10, 10);
887 fOutputListQA->Add(fh1DeltaZALambda[i]);
891 fh1PtJet[i] = new TH1D(Form("fh1PtJet_%d", i), Form("Jet pt spectrum, cent: %s;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), iNJetPtBins, dJetPtMin, dJetPtMax);
892 fOutputListStd->Add(fh1PtJet[i]);
893 fh1EtaJet[i] = new TH1D(Form("fh1EtaJet_%d", i), Form("Jet eta spectrum, cent: %s;#it{#eta} jet", GetCentBinLabel(i).Data()), 80, -1., 1.);
894 fOutputListStd->Add(fh1EtaJet[i]);
895 fh2EtaPtJet[i] = new TH2D(Form("fh2EtaPtJet_%d", i), Form("Jet eta vs pT spectrum, cent: %s;#it{#eta} jet;#it{p}_{T} jet (GeV/#it{c})", GetCentBinLabel(i).Data()), 80, -1., 1., iNJetPtBins, dJetPtMin, dJetPtMax);
896 fOutputListStd->Add(fh2EtaPtJet[i]);
897 fh2EtaPhiRndCone[i] = new TH2D(Form("fh2EtaPhiRndCone_%d", i), Form("Rnd. cones: eta vs phi, cent: %s;#it{#eta} cone;#it{#phi} cone", GetCentBinLabel(i).Data()), 80, -1., 1., 100, 0., TMath::TwoPi());
898 fOutputListStd->Add(fh2EtaPhiRndCone[i]);
899 fh2EtaPhiMedCone[i] = new TH2D(Form("fh2EtaPhiMedCone_%d", i), Form("Med.-cl. cones: eta vs phi, cent: %s;#it{#eta} cone;#it{#phi} cone", GetCentBinLabel(i).Data()), 80, -1., 1., 100, 0., TMath::TwoPi());
900 fOutputListStd->Add(fh2EtaPhiMedCone[i]);
901 fh1PhiJet[i] = new TH1D(Form("fh1PhiJet_%d", i), Form("Jet phi spectrum, cent: %s;#it{#phi} jet", GetCentBinLabel(i).Data()), 100, 0., TMath::TwoPi());
902 fOutputListStd->Add(fh1PhiJet[i]);
903 fh1NJetPerEvent[i] = new TH1D(Form("fh1NJetPerEvent_%d", i), Form("Number of selected jets per event, cent: %s;# jets;# events", GetCentBinLabel(i).Data()), 100, 0., 100.);
904 fOutputListStd->Add(fh1NJetPerEvent[i]);
906 fh1VtxZ[i] = new TH1D(Form("fh1VtxZ_%d", i), Form("#it{z} coordinate of the primary vertex, cent: %s;#it{z} (cm)", GetCentBinLabel(i).Data()), 150, -1.5 * fdCutVertexZ, 1.5 * fdCutVertexZ);
907 fOutputListQA->Add(fh1VtxZ[i]);
908 fh2VtxXY[i] = new TH2D(Form("fh2VtxXY_%d", i), Form("#it{xy} coordinate of the primary vertex, cent: %s;#it{x} (cm);#it{y} (cm)", GetCentBinLabel(i).Data()), 200, -0.2, 0.2, 500, -0.5, 0.5);
909 fOutputListQA->Add(fh2VtxXY[i]);
910 // fOutputListStd->Add([i]);
914 fh1V0K0sPtMCGen[i] = new TH1D(Form("fh1V0K0sPtMCGen_%d", i), Form("MC K0s generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
915 fOutputListMC->Add(fh1V0K0sPtMCGen[i]);
916 fh2V0K0sPtMassMCRec[i] = new TH2D(Form("fh2V0K0sPtMassMCRec_%d", i), Form("MC K0s associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
917 fOutputListMC->Add(fh2V0K0sPtMassMCRec[i]);
918 fh1V0K0sPtMCRecFalse[i] = new TH1D(Form("fh1V0K0sPtMCRecFalse_%d", i), Form("MC K0s false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
919 fOutputListMC->Add(fh1V0K0sPtMCRecFalse[i]);
921 fh2V0K0sEtaPtMCGen[i] = new TH2D(Form("fh2V0K0sEtaPtMCGen_%d", i), Form("MC K0s generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
922 fOutputListMC->Add(fh2V0K0sEtaPtMCGen[i]);
923 fh3V0K0sEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0K0sEtaPtMassMCRec_%d", i), Form("MC K0s associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaK, xminEtaK, xmaxEtaK);
924 fOutputListMC->Add(fh3V0K0sEtaPtMassMCRec[i]);
926 fh2V0K0sInJetPtMCGen[i] = new TH2D(Form("fh2V0K0sInJetPtMCGen_%d", i), Form("MC K0s in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
927 fOutputListMC->Add(fh2V0K0sInJetPtMCGen[i]);
928 fh3V0K0sInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0K0sInJetPtMassMCRec_%d", i), Form("MC K0s in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
929 fOutputListMC->Add(fh3V0K0sInJetPtMassMCRec[i]);
931 fh3V0K0sInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0K0sInJetEtaPtMCGen_%d", i), Form("MC K0s generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
932 fOutputListMC->Add(fh3V0K0sInJetEtaPtMCGen[i]);
933 fh4V0K0sInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0K0sInJetEtaPtMassMCRec_%d", i), Form("MC K0s associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaKInRec, xminEtaKInRec, xmaxEtaKInRec);
934 fOutputListMC->Add(fh4V0K0sInJetEtaPtMassMCRec[i]);
936 fh2V0K0sMCResolMPt[i] = new TH2D(Form("fh2V0K0sMCResolMPt_%d", i), Form("MC K0s associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
937 fOutputListMC->Add(fh2V0K0sMCResolMPt[i]);
938 fh2V0K0sMCPtGenPtRec[i] = new TH2D(Form("fh2V0K0sMCPtGenPtRec_%d", i), Form("MC K0s associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
939 fOutputListMC->Add(fh2V0K0sMCPtGenPtRec[i]);
942 fh1V0LambdaPtMCGen[i] = new TH1D(Form("fh1V0LambdaPtMCGen_%d", i), Form("MC Lambda generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
943 fOutputListMC->Add(fh1V0LambdaPtMCGen[i]);
944 fh2V0LambdaPtMassMCRec[i] = new TH2D(Form("fh2V0LambdaPtMassMCRec_%d", i), Form("MC Lambda associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
945 fOutputListMC->Add(fh2V0LambdaPtMassMCRec[i]);
946 fh1V0LambdaPtMCRecFalse[i] = new TH1D(Form("fh1V0LambdaPtMCRecFalse_%d", i), Form("MC Lambda false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
947 fOutputListMC->Add(fh1V0LambdaPtMCRecFalse[i]);
949 fh2V0LambdaEtaPtMCGen[i] = new TH2D(Form("fh2V0LambdaEtaPtMCGen_%d", i), Form("MC Lambda generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
950 fOutputListMC->Add(fh2V0LambdaEtaPtMCGen[i]);
951 fh3V0LambdaEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0LambdaEtaPtMassMCRec_%d", i), Form("MC Lambda associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaL, xminEtaL, xmaxEtaL);
952 fOutputListMC->Add(fh3V0LambdaEtaPtMassMCRec[i]);
954 fh2V0LambdaInJetPtMCGen[i] = new TH2D(Form("fh2V0LambdaInJetPtMCGen_%d", i), Form("MC Lambda in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
955 fOutputListMC->Add(fh2V0LambdaInJetPtMCGen[i]);
956 fh3V0LambdaInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0LambdaInJetPtMassMCRec_%d", i), Form("MC Lambda in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
957 fOutputListMC->Add(fh3V0LambdaInJetPtMassMCRec[i]);
959 fh3V0LambdaInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0LambdaInJetEtaPtMCGen_%d", i), Form("MC Lambda generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
960 fOutputListMC->Add(fh3V0LambdaInJetEtaPtMCGen[i]);
961 fh4V0LambdaInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0LambdaInJetEtaPtMassMCRec_%d", i), Form("MC Lambda associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaLInRec, xminEtaLInRec, xmaxEtaLInRec);
962 fOutputListMC->Add(fh4V0LambdaInJetEtaPtMassMCRec[i]);
964 fh2V0LambdaMCResolMPt[i] = new TH2D(Form("fh2V0LambdaMCResolMPt_%d", i), Form("MC Lambda associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
965 fOutputListMC->Add(fh2V0LambdaMCResolMPt[i]);
966 fh2V0LambdaMCPtGenPtRec[i] = new TH2D(Form("fh2V0LambdaMCPtGenPtRec_%d", i), Form("MC Lambda associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
967 fOutputListMC->Add(fh2V0LambdaMCPtGenPtRec[i]);
970 fh1V0ALambdaPtMCGen[i] = new TH1D(Form("fh1V0ALambdaPtMCGen_%d", i), Form("MC ALambda generated: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
971 fOutputListMC->Add(fh1V0ALambdaPtMCGen[i]);
972 fh2V0ALambdaPtMassMCRec[i] = new TH2D(Form("fh2V0ALambdaPtMassMCRec_%d", i), Form("MC ALambda associated: pt-m spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{m}_{inv} (GeV/#it{c}^{2})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
973 fOutputListMC->Add(fh2V0ALambdaPtMassMCRec[i]);
974 fh1V0ALambdaPtMCRecFalse[i] = new TH1D(Form("fh1V0ALambdaPtMCRecFalse_%d", i), Form("MC ALambda false: pt spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max);
975 fOutputListMC->Add(fh1V0ALambdaPtMCRecFalse[i]);
977 fh2V0ALambdaEtaPtMCGen[i] = new TH2D(Form("fh2V0ALambdaEtaPtMCGen_%d", i), Form("MC ALambda generated: pt-eta spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsEtaV0, -dRangeEtaV0Max, dRangeEtaV0Max);
978 fOutputListMC->Add(fh2V0ALambdaEtaPtMCGen[i]);
979 fh3V0ALambdaEtaPtMassMCRec[i] = new THnSparseD(Form("fh3V0ALambdaEtaPtMassMCRec_%d", i), Form("MC ALambda associated: m-pt-eta spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta", GetCentBinLabel(i).Data()), 3, binsEtaL, xminEtaL, xmaxEtaL);
980 fOutputListMC->Add(fh3V0ALambdaEtaPtMassMCRec[i]);
982 fh2V0ALambdaInJetPtMCGen[i] = new TH2D(Form("fh2V0ALambdaInJetPtMCGen_%d", i), Form("MC ALambda in jet generated: pt-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNJetPtBins, dJetPtMin, dJetPtMax);
983 fOutputListMC->Add(fh2V0ALambdaInJetPtMCGen[i]);
984 fh3V0ALambdaInJetPtMassMCRec[i] = new THnSparseD(Form("fh3V0ALambdaInJetPtMassMCRec_%d", i), Form("MC ALambda in jet associated: m-pt-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
985 fOutputListMC->Add(fh3V0ALambdaInJetPtMassMCRec[i]);
987 fh3V0ALambdaInJetEtaPtMCGen[i] = new THnSparseD(Form("fh3V0ALambdaInJetEtaPtMCGen_%d", i), Form("MC ALambda generated: pt-eta-ptJet spectrum, cent: %s;MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 4, binsEtaInGen, xminEtaInGen, xmaxEtaInGen);
988 fOutputListMC->Add(fh3V0ALambdaInJetEtaPtMCGen[i]);
989 fh4V0ALambdaInJetEtaPtMassMCRec[i] = new THnSparseD(Form("fh4V0ALambdaInJetEtaPtMassMCRec_%d", i), Form("MC ALambda associated: m-pt-eta-ptJet spectrum, cent: %s;#it{m}_{inv} (GeV/#it{c}^{2});MC #it{p}_{T} (GeV/#it{c});#eta;#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), 5, binsEtaLInRec, xminEtaLInRec, xmaxEtaLInRec);
990 fOutputListMC->Add(fh4V0ALambdaInJetEtaPtMassMCRec[i]);
992 fh2V0ALambdaMCResolMPt[i] = new TH2D(Form("fh2V0ALambdaMCResolMPt_%d", i), Form("MC ALambda associated: #Delta#it{m} vs pt, cent %s;#Delta#it{m} = #it{m}_{inv} - #it{m}_{true} (GeV/#it{c}^{2});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), 100, -0.02, 0.02, iNBinsPtV0, dPtV0Min, dPtV0Max);
993 fOutputListMC->Add(fh2V0ALambdaMCResolMPt[i]);
994 fh2V0ALambdaMCPtGenPtRec[i] = new TH2D(Form("fh2V0ALambdaMCPtGenPtRec_%d", i), Form("MC ALambda associated: pt gen vs pt rec, cent %s;#it{p}_{T}^{gen} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtV0, dPtV0Min, dPtV0Max, iNBinsPtV0, dPtV0Min, dPtV0Max);
995 fOutputListMC->Add(fh2V0ALambdaMCPtGenPtRec[i]);
997 Int_t iNBinsPtXi = 80;
998 Double_t dPtXiMin = 0;
999 Double_t dPtXiMax = 8;
1000 const Int_t iNDimFD = 3;
1001 Int_t binsFD[iNDimFD] = {iNBinsPtV0, iNBinsPtXi, iNJetPtBins};
1002 Double_t xminFD[iNDimFD] = {dPtV0Min, dPtXiMin, dJetPtMin};
1003 Double_t xmaxFD[iNDimFD] = {dPtV0Max, dPtXiMax, dJetPtMax};
1004 fhnV0LambdaInclMCFD[i] = new THnSparseD(Form("fhnV0LambdaInclMCFD_%d", i), Form("MC Lambda associated, inclusive, from Xi: pt-pt, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1005 fOutputListMC->Add(fhnV0LambdaInclMCFD[i]);
1006 fhnV0LambdaInJetsMCFD[i] = new THnSparseD(Form("fhnV0LambdaInJetsMCFD_%d", i), Form("MC Lambda associated, in JC, from Xi: pt-pt-ptJet, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1007 fOutputListMC->Add(fhnV0LambdaInJetsMCFD[i]);
1008 fhnV0LambdaBulkMCFD[i] = new THnSparseD(Form("fhnV0LambdaBulkMCFD_%d", i), Form("MC Lambda associated, in no jet events, from Xi: pt-pt, cent %s;#it{p}_{T}^{#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1009 fOutputListMC->Add(fhnV0LambdaBulkMCFD[i]);
1010 fh1V0XiPtMCGen[i] = new TH1D(Form("fh1V0XiPtMCGen_%d", i), Form("MC Xi^{-} generated: Pt spectrum, cent %s;#it{p}_{T}^{#Xi^{-},gen.} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtXi, dPtXiMin, dPtXiMax);
1011 fOutputListMC->Add(fh1V0XiPtMCGen[i]);
1012 fhnV0ALambdaInclMCFD[i] = new THnSparseD(Form("fhnV0ALambdaInclMCFD_%d", i), Form("MC ALambda associated, from AXi: pt-pt, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1013 fOutputListMC->Add(fhnV0ALambdaInclMCFD[i]);
1014 fhnV0ALambdaInJetsMCFD[i] = new THnSparseD(Form("fhnV0ALambdaInJetsMCFD_%d", i), Form("MC ALambda associated, in JC, from AXi: pt-pt-ptJet, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1015 fOutputListMC->Add(fhnV0ALambdaInJetsMCFD[i]);
1016 fhnV0ALambdaBulkMCFD[i] = new THnSparseD(Form("fhnV0ALambdaBulkMCFD_%d", i), Form("MC ALambda associated, in no jet events, from AXi: pt-pt-ptJet, cent %s;#it{p}_{T}^{A#Lambda,gen.} (GeV/#it{c});#it{p}_{T}^{A#Xi,gen.} (GeV/#it{c});#it{p}_{T}^{jet} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNDimFD, binsFD, xminFD, xmaxFD);
1017 fOutputListMC->Add(fhnV0ALambdaBulkMCFD[i]);
1018 fh1V0AXiPtMCGen[i] = new TH1D(Form("fh1V0AXiPtMCGen_%d", i), Form("MC AXi^{-} generated: Pt spectrum, cent %s;#it{p}_{T}^{A#Xi^{-},gen.} (GeV/#it{c})", GetCentBinLabel(i).Data()), iNBinsPtXi, dPtXiMin, dPtXiMax);
1019 fOutputListMC->Add(fh1V0AXiPtMCGen[i]);
1022 // fhnV0K0sInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1023 fhnV0K0sInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInclDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1024 // fhnV0K0sInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0K0sInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC K0S, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1025 fhnV0K0sInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0K0sInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC K0S, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1026 // fhnV0LambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1027 fhnV0LambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1028 // fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0LambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC Lambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1029 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0LambdaInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC Lambda, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1030 // fhnV0ALambdaInclDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, inclusive, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1031 fhnV0ALambdaInclDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInclDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, inclusive, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1032 // fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i] = new THnSparseD(Form("fhnV0ALambdaInJetsDaughterEtaPtPtMCGen_%d",i),Form("MC ALambda, in JC, gen., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet",GetCentBinLabel(i).Data()),iNDimEtaD,binsEtaDaughter,xminEtaDaughter,xmaxEtaDaughter);
1033 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i] = new THnSparseD(Form("fhnV0ALambdaInJetsDaughterEtaPtPtMCRec_%d", i), Form("MC ALambda, in JC, assoc., daughters: charge-etaD-ptD-etaV0-ptV0-ptJet, cent: %s;charge;eta daughter;pT daughter;eta V0;pT V0;pT jet", GetCentBinLabel(i).Data()), iNDimEtaD, binsEtaDaughter, xminEtaDaughter, xmaxEtaDaughter);
1035 // fOutputListMC->Add(fhnV0K0sInclDaughterEtaPtPtMCGen[i]);
1036 fOutputListMC->Add(fhnV0K0sInclDaughterEtaPtPtMCRec[i]);
1037 // fOutputListMC->Add(fhnV0K0sInJetsDaughterEtaPtPtMCGen[i]);
1038 fOutputListMC->Add(fhnV0K0sInJetsDaughterEtaPtPtMCRec[i]);
1039 // fOutputListMC->Add(fhnV0LambdaInclDaughterEtaPtPtMCGen[i]);
1040 fOutputListMC->Add(fhnV0LambdaInclDaughterEtaPtPtMCRec[i]);
1041 // fOutputListMC->Add(fhnV0LambdaInJetsDaughterEtaPtPtMCGen[i]);
1042 fOutputListMC->Add(fhnV0LambdaInJetsDaughterEtaPtPtMCRec[i]);
1043 // fOutputListMC->Add(fhnV0ALambdaInclDaughterEtaPtPtMCGen[i]);
1044 fOutputListMC->Add(fhnV0ALambdaInclDaughterEtaPtPtMCRec[i]);
1045 // fOutputListMC->Add(fhnV0ALambdaInJetsDaughterEtaPtPtMCGen[i]);
1046 fOutputListMC->Add(fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[i]);
1051 for(Int_t i = 0; i < fgkiNQAIndeces; i++)
1053 // [i] = new TH1D(Form("%d",i),";;Counts",,,);
1054 fh1QAV0Status[i] = new TH1D(Form("fh1QAV0Status_%d", i), "QA: V0 status", 2, 0, 2);
1055 fh1QAV0TPCRefit[i] = new TH1D(Form("fh1QAV0TPCRefit_%d", i), "QA: TPC refit", 2, 0, 2);
1056 fh1QAV0TPCRows[i] = new TH1D(Form("fh1QAV0TPCRows_%d", i), "QA: TPC Rows", 160, 0, 160);
1057 fh1QAV0TPCFindable[i] = new TH1D(Form("fh1QAV0TPCFindable_%d", i), "QA: TPC Findable", 160, 0, 160);
1058 fh1QAV0TPCRowsFind[i] = new TH1D(Form("fh1QAV0TPCRowsFind_%d", i), "QA: TPC Rows/Findable", 100, 0, 2);
1059 fh1QAV0Eta[i] = new TH1D(Form("fh1QAV0Eta_%d", i), "QA: Daughter Eta", 200, -2, 2);
1060 fh2QAV0EtaRows[i] = new TH2D(Form("fh2QAV0EtaRows_%d", i), "QA: Daughter Eta vs TPC rows;#eta;TPC rows", 200, -2, 2, 160, 0, 160);
1061 fh2QAV0PtRows[i] = new TH2D(Form("fh2QAV0PtRows_%d", i), "QA: Daughter Pt vs TPC rows;pt;TPC rows", 100, 0, 10, 160, 0, 160);
1062 fh2QAV0PhiRows[i] = new TH2D(Form("fh2QAV0PhiRows_%d", i), "QA: Daughter Phi vs TPC rows;#phi;TPC rows", 100, 0, TMath::TwoPi(), 160, 0, 160);
1063 fh2QAV0NClRows[i] = new TH2D(Form("fh2QAV0NClRows_%d", i), "QA: Daughter NCl vs TPC rows;findable clusters;TPC rows", 100, 0, 160, 160, 0, 160);
1064 fh2QAV0EtaNCl[i] = new TH2D(Form("fh2QAV0EtaNCl_%d", i), "QA: Daughter Eta vs NCl;#eta;findable clusters", 200, -2, 2, 160, 0, 160);
1066 fh2QAV0EtaPtK0sPeak[i] = new TH2D(Form("fh2QAV0EtaPtK0sPeak_%d", i), "QA: K0s: Daughter Eta vs V0 pt, peak;track eta;V0 pt", 200, -2, 2, iNBinsPtV0, dPtV0Min, dPtV0Max);
1067 fh2QAV0EtaEtaK0s[i] = new TH2D(Form("fh2QAV0EtaEtaK0s_%d", i), "QA: K0s: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1068 fh2QAV0PhiPhiK0s[i] = new TH2D(Form("fh2QAV0PhiPhiK0s_%d", i), "QA: K0s: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1069 fh1QAV0RapK0s[i] = new TH1D(Form("fh1QAV0RapK0s_%d", i), "QA: K0s: V0 Rapidity", 200, -2, 2);
1070 fh2QAV0PtPtK0sPeak[i] = new TH2D(Form("fh2QAV0PtPtK0sPeak_%d", i), "QA: K0s: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1072 fh2QAV0EtaPtLambdaPeak[i] = new TH2D(Form("fh2QAV0EtaPtLambdaPeak_%d", i), "QA: Lambda: Daughter Eta vs V0 pt, peak;track eta;V0 pt", 200, -2, 2, iNBinsPtV0, dPtV0Min, dPtV0Max);
1073 fh2QAV0EtaEtaLambda[i] = new TH2D(Form("fh2QAV0EtaEtaLambda_%d", i), "QA: Lambda: Eta vs Eta Daughter", 200, -2, 2, 200, -2, 2);
1074 fh2QAV0PhiPhiLambda[i] = new TH2D(Form("fh2QAV0PhiPhiLambda_%d", i), "QA: Lambda: Phi vs Phi Daughter", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
1075 fh1QAV0RapLambda[i] = new TH1D(Form("fh1QAV0RapLambda_%d", i), "QA: Lambda: V0 Rapidity", 200, -2, 2);
1076 fh2QAV0PtPtLambdaPeak[i] = new TH2D(Form("fh2QAV0PtPtLambdaPeak_%d", i), "QA: Lambda: Daughter Pt vs Pt;neg pt;pos pt", 100, 0, 5, 100, 0, 5);
1078 fh1QAV0Pt[i] = new TH1D(Form("fh1QAV0Pt_%d", i), "QA: Daughter Pt", 100, 0, 5);
1079 fh1QAV0Charge[i] = new TH1D(Form("fh1QAV0Charge_%d", i), "QA: V0 Charge", 3, -1, 2);
1080 fh1QAV0DCAVtx[i] = new TH1D(Form("fh1QAV0DCAVtx_%d", i), "QA: DCA daughters to primary vertex", 100, 0, 10);
1081 fh1QAV0DCAV0[i] = new TH1D(Form("fh1QAV0DCAV0_%d", i), "QA: DCA daughters", 100, 0, 2);
1082 fh1QAV0Cos[i] = new TH1D(Form("fh1QAV0Cos_%d", i), "QA: CPA", 10000, 0.9, 1);
1083 fh1QAV0R[i] = new TH1D(Form("fh1QAV0R_%d", i), "QA: R", 1500, 0, 150);
1084 fh1QACTau2D[i] = new TH1D(Form("fh1QACTau2D_%d", i), "QA: K0s: c#tau 2D;mR/pt#tau", 100, 0, 10);
1085 fh1QACTau3D[i] = new TH1D(Form("fh1QACTau3D_%d", i), "QA: K0s: c#tau 3D;mL/p#tau", 100, 0, 10);
1087 fh2ArmPod[i] = new TH2D(Form("fh2ArmPod_%d", i), "Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1088 fh2ArmPodK0s[i] = new TH2D(Form("fh2ArmPodK0s_%d", i), "K0s: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1089 fh2ArmPodLambda[i] = new TH2D(Form("fh2ArmPodLambda_%d", i), "Lambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1090 fh2ArmPodALambda[i] = new TH2D(Form("fh2ArmPodALambda_%d", i), "ALambda: Armenteros-Podolanski;#alpha;#it{p}_{T}^{Arm}", 100, -1., 1., 50, 0., 0.25);
1092 fOutputListQA->Add(fh1QAV0Status[i]);
1093 fOutputListQA->Add(fh1QAV0TPCRefit[i]);
1094 fOutputListQA->Add(fh1QAV0TPCRows[i]);
1095 fOutputListQA->Add(fh1QAV0TPCFindable[i]);
1096 fOutputListQA->Add(fh1QAV0TPCRowsFind[i]);
1097 fOutputListQA->Add(fh1QAV0Eta[i]);
1098 fOutputListQA->Add(fh2QAV0EtaRows[i]);
1099 fOutputListQA->Add(fh2QAV0PtRows[i]);
1100 fOutputListQA->Add(fh2QAV0PhiRows[i]);
1101 fOutputListQA->Add(fh2QAV0NClRows[i]);
1102 fOutputListQA->Add(fh2QAV0EtaNCl[i]);
1104 fOutputListQA->Add(fh2QAV0EtaPtK0sPeak[i]);
1105 fOutputListQA->Add(fh2QAV0EtaEtaK0s[i]);
1106 fOutputListQA->Add(fh2QAV0PhiPhiK0s[i]);
1107 fOutputListQA->Add(fh1QAV0RapK0s[i]);
1108 fOutputListQA->Add(fh2QAV0PtPtK0sPeak[i]);
1110 fOutputListQA->Add(fh2QAV0EtaPtLambdaPeak[i]);
1111 fOutputListQA->Add(fh2QAV0EtaEtaLambda[i]);
1112 fOutputListQA->Add(fh2QAV0PhiPhiLambda[i]);
1113 fOutputListQA->Add(fh1QAV0RapLambda[i]);
1114 fOutputListQA->Add(fh2QAV0PtPtLambdaPeak[i]);
1116 fOutputListQA->Add(fh1QAV0Pt[i]);
1117 fOutputListQA->Add(fh1QAV0Charge[i]);
1118 fOutputListQA->Add(fh1QAV0DCAVtx[i]);
1119 fOutputListQA->Add(fh1QAV0DCAV0[i]);
1120 fOutputListQA->Add(fh1QAV0Cos[i]);
1121 fOutputListQA->Add(fh1QAV0R[i]);
1122 fOutputListQA->Add(fh1QACTau2D[i]);
1123 fOutputListQA->Add(fh1QACTau3D[i]);
1125 fOutputListQA->Add(fh2ArmPod[i]);
1126 fOutputListQA->Add(fh2ArmPodK0s[i]);
1127 fOutputListQA->Add(fh2ArmPodLambda[i]);
1128 fOutputListQA->Add(fh2ArmPodALambda[i]);
1131 fh2CutTPCRowsK0s[i] = new TH2D(Form("fh2CutTPCRowsK0s_%d", i), "Cuts: K0s: TPC Rows vs mass;#it{m}_{inv} (GeV/#it{c}^{2});TPC rows", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 160, 0, 160);
1132 fh2CutTPCRowsLambda[i] = new TH2D(Form("fh2CutTPCRowsLambda_%d", i), "Cuts: Lambda: TPC Rows vs mass;#it{m}_{inv} (GeV/#it{c}^{2});TPC rows", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 160, 0, 160);
1133 fh2CutPtPosK0s[i] = new TH2D(Form("fh2CutPtPosK0s_%d", i), "Cuts: K0s: Pt pos;#it{m}_{inv} (GeV/#it{c}^{2});pt pos", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 5);
1134 fh2CutPtNegK0s[i] = new TH2D(Form("fh2CutPtNegK0s_%d", i), "Cuts: K0s: Pt neg;#it{m}_{inv} (GeV/#it{c}^{2});pt neg", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 5);
1135 fh2CutPtPosLambda[i] = new TH2D(Form("fh2CutPtPosLambda_%d", i), "Cuts: Lambda: Pt pos;#it{m}_{inv} (GeV/#it{c}^{2});pt pos", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 5);
1136 fh2CutPtNegLambda[i] = new TH2D(Form("fh2CutPtNegLambda_%d", i), "Cuts: Lambda: Pt neg;#it{m}_{inv} (GeV/#it{c}^{2});pt neg", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 5);
1137 fh2CutDCAVtx[i] = new TH2D(Form("fh2CutDCAVtx_%d", i), "Cuts: DCA daughters to prim. vtx.;#it{m}_{inv} (GeV/#it{c}^{2});DCA daughter to prim. vtx. (cm)", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1138 fh2CutDCAV0[i] = new TH2D(Form("fh2CutDCAV0_%d", i), "Cuts: DCA daughters;#it{m}_{inv} (GeV/#it{c}^{2});DCA daughters / #sigma_{TPC}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 2);
1139 fh2CutCos[i] = new TH2D(Form("fh2CutCos_%d", i), "Cuts: CPA;#it{m}_{inv} (GeV/#it{c}^{2});CPA", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 10000, 0.9, 1);
1140 fh2CutR[i] = new TH2D(Form("fh2CutR_%d", i), "Cuts: R;#it{m}_{inv} (GeV/#it{c}^{2});R (cm)", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 1500, 0, 150);
1141 fh2CutEtaK0s[i] = new TH2D(Form("fh2CutEtaK0s_%d", i), "Cuts: K0s: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 200, -2, 2);
1142 fh2CutEtaLambda[i] = new TH2D(Form("fh2CutEtaLambda_%d", i), "Cuts: Lambda: Eta;#it{m}_{inv} (GeV/#it{c}^{2});#eta", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 200, -2, 2);
1143 fh2CutRapK0s[i] = new TH2D(Form("fh2CutRapK0s_%d", i), "Cuts: K0s: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 200, -2, 2);
1144 fh2CutRapLambda[i] = new TH2D(Form("fh2CutRapLambda_%d", i), "Cuts: Lambda: Rapidity;#it{m}_{inv} (GeV/#it{c}^{2});y", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 200, -2, 2);
1145 fh2CutCTauK0s[i] = new TH2D(Form("fh2CutCTauK0s_%d", i), "Cuts: K0s: #it{c#tau};#it{m}_{inv} (GeV/#it{c}^{2});#it{mL/p#tau}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1146 fh2CutCTauLambda[i] = new TH2D(Form("fh2CutCTauLambda_%d", i), "Cuts: Lambda: #it{c#tau};#it{m}_{inv} (GeV/#it{c}^{2});#it{mL/p#tau}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1147 fh2CutPIDPosK0s[i] = new TH2D(Form("fh2CutPIDPosK0s_%d", i), "Cuts: K0s: PID pos;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1148 fh2CutPIDNegK0s[i] = new TH2D(Form("fh2CutPIDNegK0s_%d", i), "Cuts: K0s: PID neg;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax, 100, 0, 10);
1149 fh2CutPIDPosLambda[i] = new TH2D(Form("fh2CutPIDPosLambda_%d", i), "Cuts: Lambda: PID pos;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1150 fh2CutPIDNegLambda[i] = new TH2D(Form("fh2CutPIDNegLambda_%d", i), "Cuts: Lambda: PID neg;#it{m}_{inv} (GeV/#it{c}^{2});##sigma_{d#it{E}/d#it{x}}", fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax, 100, 0, 10);
1151 fh2Tau3DVs2D[i] = new TH2D(Form("fh2Tau3DVs2D_%d", i), "Decay 3D vs 2D;pt;3D/2D", 100, 0, 10, 200, 0.5, 1.5);
1153 fOutputListCuts->Add(fh2CutTPCRowsK0s[i]);
1154 fOutputListCuts->Add(fh2CutTPCRowsLambda[i]);
1155 fOutputListCuts->Add(fh2CutPtPosK0s[i]);
1156 fOutputListCuts->Add(fh2CutPtNegK0s[i]);
1157 fOutputListCuts->Add(fh2CutPtPosLambda[i]);
1158 fOutputListCuts->Add(fh2CutPtNegLambda[i]);
1159 fOutputListCuts->Add(fh2CutDCAVtx[i]);
1160 fOutputListCuts->Add(fh2CutDCAV0[i]);
1161 fOutputListCuts->Add(fh2CutCos[i]);
1162 fOutputListCuts->Add(fh2CutR[i]);
1163 fOutputListCuts->Add(fh2CutEtaK0s[i]);
1164 fOutputListCuts->Add(fh2CutEtaLambda[i]);
1165 fOutputListCuts->Add(fh2CutRapK0s[i]);
1166 fOutputListCuts->Add(fh2CutRapLambda[i]);
1167 fOutputListCuts->Add(fh2CutCTauK0s[i]);
1168 fOutputListCuts->Add(fh2CutCTauLambda[i]);
1169 fOutputListCuts->Add(fh2CutPIDPosK0s[i]);
1170 fOutputListCuts->Add(fh2CutPIDNegK0s[i]);
1171 fOutputListCuts->Add(fh2CutPIDPosLambda[i]);
1172 fOutputListCuts->Add(fh2CutPIDNegLambda[i]);
1173 fOutputListCuts->Add(fh2Tau3DVs2D[i]);
1177 for(Int_t i = 0; i < fgkiNCategV0; i++)
1179 fh1V0InvMassK0sAll[i] = new TH1D(Form("fh1V0InvMassK0sAll_%d", i), Form("K0s: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassK0s, fgkdMassK0sMin, fgkdMassK0sMax);
1180 fh1V0InvMassLambdaAll[i] = new TH1D(Form("fh1V0InvMassLambdaAll_%d", i), Form("Lambda: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
1181 fh1V0InvMassALambdaAll[i] = new TH1D(Form("fh1V0InvMassALambdaAll_%d", i), Form("ALambda: V0 invariant mass, %s;#it{m}_{inv} (GeV/#it{c}^{2});counts", categV0[i].Data()), fgkiNBinsMassLambda, fgkdMassLambdaMin, fgkdMassLambdaMax);
1182 fOutputListStd->Add(fh1V0InvMassK0sAll[i]);
1183 fOutputListStd->Add(fh1V0InvMassLambdaAll[i]);
1184 fOutputListStd->Add(fh1V0InvMassALambdaAll[i]);
1187 for(Int_t i = 0; i < fOutputListStd->GetEntries(); ++i)
1189 TH1* h1 = dynamic_cast<TH1*>(fOutputListStd->At(i));
1195 THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListStd->At(i));
1198 for(Int_t i = 0; i < fOutputListQA->GetEntries(); ++i)
1200 TH1* h1 = dynamic_cast<TH1*>(fOutputListQA->At(i));
1206 THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListQA->At(i));
1209 for(Int_t i = 0; i < fOutputListCuts->GetEntries(); ++i)
1211 TH1* h1 = dynamic_cast<TH1*>(fOutputListCuts->At(i));
1217 THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListCuts->At(i));
1220 for(Int_t i = 0; i < fOutputListMC->GetEntries(); ++i)
1222 TH1* h1 = dynamic_cast<TH1*>(fOutputListMC->At(i));
1228 THnSparse* hn = dynamic_cast<THnSparse*>(fOutputListMC->At(i));
1232 PostData(1, fOutputListStd);
1233 PostData(2, fOutputListQA);
1234 PostData(3, fOutputListCuts);
1235 PostData(4, fOutputListMC);
1236 // if (fbTreeOutput)
1237 // PostData(5,ftreeOut);
1240 void AliAnalysisTaskV0sInJets::UserExec(Option_t*)
1242 // Main loop, called for each event
1243 if(fDebug > 5) printf("TaskV0sInJets: UserExec: Start\n");
1245 // reset branches for each event
1247 fBranchV0Rec->Clear();
1249 fBranchV0Gen->Clear();
1251 fBranchJet->Clear();
1253 fEventInfo->Reset();
1258 if(fDebug > 2) printf("TaskV0sInJets: AOD analysis\n");
1259 fh1EventCounterCut->Fill(0); // all available selected events (collision candidates)
1261 if(fDebug > 5) printf("TaskV0sInJets: UserExec: Loading AOD\n");
1262 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); // input AOD
1263 fAODOut = AODEvent(); // output AOD
1266 if(fDebug > 0) printf("TaskV0sInJets: No output AOD found\n");
1271 if(fDebug > 0) printf("TaskV0sInJets: No input AOD found\n");
1274 if(fDebug > 5) printf("TaskV0sInJets: UserExec: Loading AOD OK\n");
1276 TClonesArray* arrayMC = 0; // array particles in the MC event
1277 AliAODMCHeader* headerMC = 0; // MC header
1278 Int_t iNTracksMC = 0; // number of MC tracks
1279 Double_t dPrimVtxMCX = 0., dPrimVtxMCY = 0., dPrimVtxMCZ = 0.; // position of the MC primary vertex
1284 arrayMC = (TClonesArray*)fAODIn->FindListObject(AliAODMCParticle::StdBranchName());
1287 if(fDebug > 0) printf("TaskV0sInJets: No MC array found\n");
1290 if(fDebug > 5) printf("TaskV0sInJets: MC array found\n");
1291 iNTracksMC = arrayMC->GetEntriesFast();
1292 if(fDebug > 5) printf("TaskV0sInJets: There are %d MC tracks in this event\n", iNTracksMC);
1295 headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
1298 if(fDebug > 0) printf("TaskV0sInJets: No MC header found\n");
1301 // get position of the MC primary vertex
1302 dPrimVtxMCX = headerMC->GetVtxX();
1303 dPrimVtxMCY = headerMC->GetVtxY();
1304 dPrimVtxMCZ = headerMC->GetVtxZ();
1307 // PID Response Task object
1308 AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1309 AliInputEventHandler* inputHandler = (AliInputEventHandler*)mgr->GetInputEventHandler();
1310 AliPIDResponse* fPIDResponse = inputHandler->GetPIDResponse();
1313 if(fDebug > 0) printf("TaskV0sInJets: No PID response object found\n");
1318 fh1EventCounterCut->Fill(1);
1321 if(!IsSelectedForJets(fAODIn, fdCutVertexZ, fdCutVertexR2, fdCutCentLow, fdCutCentHigh, 1, 0.1)) // cut on |delta z| in 2011 data between SPD vertex and nominal primary vertex
1322 // if (!IsSelectedForJets(fAODIn,fdCutVertexZ,fdCutVertexR2,fdCutCentLow,fdCutCentHigh)) // no need for cutting in 2010 data
1324 if(fDebug > 5) printf("TaskV0sInJets: Event rejected\n");
1328 // fdCentrality = ((AliVAODHeader*)fAODIn->GetHeader())->GetCentrality(); // event centrality
1329 fdCentrality = ((AliVAODHeader*)fAODIn->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M"); // event centrality
1332 Int_t iCentIndex = GetCentralityBinIndex(fdCentrality); // get index of centrality bin
1335 if(fDebug > 5) printf("TaskV0sInJets: Event is out of histogram range\n");
1338 fh1EventCounterCut->Fill(2); // selected events (vertex, centrality)
1339 fh1EventCounterCutCent[iCentIndex]->Fill(2);
1341 UInt_t iNTracks = fAODIn->GetNumberOfTracks(); // get number of tracks in event
1342 if(fDebug > 5) printf("TaskV0sInJets: There are %d tracks in this event\n", iNTracks);
1346 Int_t iNV0s = fAODIn->GetNumberOfV0s(); // get the number of V0 candidates
1349 if(fDebug > 2) printf("TaskV0sInJets: No V0s found in event\n");
1353 //===== Event is OK for the analysis =====
1354 fh1EventCent->Fill(iCentIndex);
1355 fh1EventCent2->Fill(fdCentrality);
1356 fh2EventCentTracks->Fill(fdCentrality, iNTracks);
1358 // if (fbTreeOutput)
1359 // fEventInfo->SetAODEvent(fAODIn);
1360 // printf("V0sInJets: EventInfo: Centrality: %f\n",fEventInfo->GetCentrality());
1364 fh1EventCounterCut->Fill(3); // events with V0s
1365 fh1EventCounterCutCent[iCentIndex]->Fill(3);
1368 // Int_t iNV0SelV0Rec = 0;
1369 // Int_t iNV0SelV0Gen = 0;
1371 AliAODv0* v0 = 0; // pointer to V0 candidates
1372 // AliV0Object* objectV0 = 0;
1373 TVector3 vecV0Momentum; // 3D vector of V0 momentum
1374 Double_t dMassV0K0s = 0; // invariant mass of the K0s candidate
1375 Double_t dMassV0Lambda = 0; // invariant mass of the Lambda candidate
1376 Double_t dMassV0ALambda = 0; // invariant mass of the Lambda candidate
1377 Int_t iNV0CandTot = 0; // counter of all V0 candidates at the beginning
1378 Int_t iNV0CandK0s = 0; // counter of K0s candidates at the end
1379 Int_t iNV0CandLambda = 0; // counter of Lambda candidates at the end
1380 Int_t iNV0CandALambda = 0; // counter of Lambda candidates at the end
1382 Bool_t bUseOldCuts = 0; // old reconstruction cuts
1383 Bool_t bUseAliceCuts = 0; // cuts used by Alice Zimmermann
1384 Bool_t bUseIouriCuts = 0; // cuts used by Iouri
1385 Bool_t bPrintCuts = 0; // print out which cuts are applied
1386 Bool_t bPrintJetSelection = 0; // print out which jets are selected
1388 // Values of V0 reconstruction cuts:
1390 Int_t iRefit = AliAODTrack::kTPCrefit; // TPC refit for daughter tracks
1391 Double_t dDCAToPrimVtxMin = fdCutDCAToPrimVtxMin; // 0.1; // [cm] min DCA of daughters to the prim vtx
1392 Double_t dDCADaughtersMax = fdCutDCADaughtersMax; // 1.; // [sigma of TPC tracking] max DCA between daughters
1393 Double_t dEtaDaughterMax = 0.8; // max |pseudorapidity| of daughter tracks
1394 Double_t dNSigmadEdxMax = fdCutNSigmadEdxMax;// 3.; // [sigma dE/dx] max difference between measured and expected signal of dE/dx in the TPC
1395 Double_t dPtProtonPIDMax = 1.; // [GeV/c] maxium pT of proton for applying PID cut
1397 Bool_t bOnFly = 0; // on-the-fly (yes) or offline (no) reconstructed
1398 Double_t dCPAMin = fdCutCPAMin;// 0.998; // min cosine of the pointing angle
1399 Double_t dRadiusDecayMin = 5.; // [cm] min radial distance of the decay vertex
1400 Double_t dRadiusDecayMax = 100.; // [cm] max radial distance of the decay vertex
1401 Double_t dEtaMax = 0.7; // max |pseudorapidity| of V0
1402 Double_t dNTauMax = fdCutNTauMax; // 5.0; // [tau] max proper lifetime in multiples of the mean lifetime
1405 Double_t dNCrossedRowsTPCMin = 70.; // min number of crossed TPC rows (turned off)
1406 // Double_t dCrossedRowsOverFindMin = 0.8; // min ratio crossed rows / findable clusters (turned off)
1407 // Double_t dCrossedRowsOverFindMax = 1e3; // max ratio crossed rows / findable clusters (turned off)
1408 Double_t dPtDaughterMin = 0.150; // [GeV/c] min transverse momentum of daughter tracks (turned off)
1409 Double_t dRapMax = 0.75; // max |rapidity| of V0 (turned off)
1413 Double_t dNSigmaMassMax = 3.; // [sigma m] max difference between candidate mass and real particle mass (used only for mass peak method of signal extraction)
1414 Double_t dDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
1416 // Selection of active cuts
1417 Bool_t bCutEtaDaughter = 1; // daughter pseudorapidity
1418 Bool_t bCutRapV0 = 0; // V0 rapidity
1419 Bool_t bCutEtaV0 = 1; // V0 pseudorapidity
1420 Bool_t bCutTau = 1; // V0 lifetime
1421 Bool_t bCutPid = 1; // PID (TPC dE/dx)
1422 Bool_t bCutArmPod = 1; // Armenteros-Podolanski for K0S
1423 // Bool_t bCutCross = 0; // cross contamination
1431 else if(bUseAliceCuts)
1437 else if(bUseIouriCuts)
1445 Double_t dCTauK0s = 2.6844; // [cm] c tau of K0S
1446 Double_t dCTauLambda = 7.89; // [cm] c tau of Lambda
1448 // Load PDG values of particle masses
1449 Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1450 Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1452 // PDG codes of used particles
1453 Int_t iPdgCodePion = 211;
1454 Int_t iPdgCodeProton = 2212;
1455 Int_t iPdgCodeK0s = 310;
1456 Int_t iPdgCodeLambda = 3122;
1458 // Jet selection: fdCutPtJetMin, fdCutPtTrackMin
1459 Double_t dJetEtaWindow = dEtaMax - fdRadiusJet; // max jet |pseudorapidity|, to make sure that V0s can appear in the entire jet area
1460 Double_t dCutJetAreaMin = 0.6 * TMath::Pi() * fdRadiusJet * fdRadiusJet; // minimum jet area
1461 Double_t dRadiusExcludeCone = 2 * fdRadiusJet; // radius of cones around jets excluded for V0 outside jets
1462 Bool_t bLeadingJetOnly = 0;
1467 fdCutPtTrackMin = 5;
1469 bLeadingJetOnly = 0;
1472 // Int_t iNJetAll = 0; // number of reconstructed jets in fBranchJet
1473 // iNJetAll = fBranchJet->GetEntriesFast(); // number of reconstructed jets
1474 TClonesArray* jetArray = 0; // object where the input jets are stored
1475 TClonesArray* jetArrayBg = 0; // object where the kt clusters are stored
1476 Int_t iNJet = 0; // number of reconstructed jets in the input
1477 TClonesArray* jetArraySel = new TClonesArray("AliAODJet", 0); // object where the selected jets are copied
1478 Int_t iNJetSel = 0; // number of selected reconstructed jets
1479 // iNJetSel = jetArraySel->GetEntriesFast(); // number of selected reconstructed jets
1480 TClonesArray* jetArrayPerp = new TClonesArray("AliAODJet", 0); // object where the perp. cones are stored
1481 Int_t iNJetPerp = 0; // number of perpendicular cones
1483 AliAODJet* jet = 0; // pointer to a jet
1484 // AliJetObject* objectJet = 0;
1485 AliAODJet* jetPerp = 0; // pointer to a perp. cone
1486 AliAODJet* jetRnd = 0; // pointer to a rand. cone
1487 AliAODJet* jetMed = 0; // pointer to a median cluster
1488 TVector3 vecJetMomentum; // 3D vector of jet momentum
1489 // TVector3 vecPerpConeMomentum; // 3D vector of perpendicular cone momentum
1490 // TVector3 vecRndConeMomentum; // 3D vector of random cone momentum
1491 Bool_t bJetEventGood = kTRUE; // indicator of good jet events
1493 // printf("iNJetAll, iNJetSel: %d %d\n",iNJetAll,iNJetSel);
1495 if(fbJetSelection) // analysis of V0s in jets is switched on
1497 jetArray = (TClonesArray*)(fAODOut->FindListObject(fsJetBranchName.Data())); // find object with jets in the output AOD
1500 if(fDebug > 0) printf("TaskV0sInJets: No array of name: %s\n", fsJetBranchName.Data());
1501 bJetEventGood = kFALSE;
1504 iNJet = jetArray->GetEntriesFast();
1505 if(bJetEventGood && !iNJet) // check whether there are some jets
1507 if(fDebug > 2) printf("TaskV0sInJets: No jets in array\n");
1508 bJetEventGood = kFALSE;
1512 // printf("TaskV0sInJets: Loading bg array of name: %s\n",fsJetBgBranchName.Data());
1513 jetArrayBg = (TClonesArray*)(fAODOut->FindListObject(fsJetBgBranchName.Data())); // find object with jets in the output AOD
1516 if(fDebug > 0) printf("TaskV0sInJets: No bg array of name: %s\n", fsJetBgBranchName.Data());
1517 // bJetEventGood = kFALSE;
1521 else // no in-jet analysis
1522 bJetEventGood = kFALSE;
1524 // select good jets and copy them to another array
1528 iNJet = 1; // only leading jets
1529 if(fDebug > 5) printf("TaskV0sInJets: Jet selection for %d jets\n", iNJet);
1530 for(Int_t iJet = 0; iJet < iNJet; iJet++)
1532 AliAODJet* jetSel = (AliAODJet*)jetArray->At(iJet); // load a jet in the list
1535 if(fDebug > 0) printf("TaskV0sInJets: Cannot load jet %d\n", iJet);
1538 if(bPrintJetSelection)
1539 if(fDebug > 7) printf("jet: i = %d, pT = %f, eta = %f, phi = %f, pt lead tr = %f ", iJet, jetSel->Pt(), jetSel->Eta(), jetSel->Phi(), jetSel->GetPtLeading());
1540 // printf("TaskV0sInJets: Checking pt > %.2f for jet %d with pt %.2f\n",fdCutPtJetMin,iJet,jetSel->Pt());
1541 if(jetSel->Pt() < fdCutPtJetMin) // selection of high-pt jets
1543 if(bPrintJetSelection)
1544 if(fDebug > 7) printf("rejected (pt)\n");
1547 // printf("TaskV0sInJets: Checking |eta| < %.2f for jet %d with |eta| %.2f\n",dJetEtaWindow,iJet,TMath::Abs(jetSel->Eta()));
1548 if(TMath::Abs(jetSel->Eta()) > dJetEtaWindow) // selection of jets in the chosen pseudorapidity range
1550 if(bPrintJetSelection)
1551 if(fDebug > 7) printf("rejected (eta)\n");
1556 if(jetSel->EffectiveAreaCharged() < dCutJetAreaMin)
1559 Int_t iNTracksInJet = 0;
1560 Double_t dPtLeadTrack = 0; // pt of the leading track
1561 // Int_t iLeadTrack = 0;
1562 iNTracksInJet = jetSel->GetRefTracks()->GetEntriesFast(); // number od tracks that constitute the jet
1563 // printf("TaskV0sInJets: Searching for leading track from %d tracks in jet %d\n",iNTracksInJet,iJet);
1564 if(fdCutPtTrackMin > 0) // a positive min leading track pt is set
1566 for(Int_t j = 0; j < iNTracksInJet; j++) // find the track with the highest pt
1568 AliAODTrack* track = (AliAODTrack*)jetSel->GetTrack(j); // is this the leading track?
1571 // printf("TaskV0sInJets: %d: %.2f\n",j,track->Pt());
1572 if(track->Pt() > dPtLeadTrack)
1574 dPtLeadTrack = track->Pt();
1578 // printf("Leading track pT: my: %f, ali: %f\n",dPtLeadTrack,jetSel->GetPtLeading());
1579 // printf("TaskV0sInJets: Checking leading track pt > %.2f for pt %.2f of track %d in jet %d\n",fdCutPtTrackMin,dPtLeadTrack,iLeadTrack,iJet);
1580 if(dPtLeadTrack < fdCutPtTrackMin) // selection of high-pt jet-track events
1582 if(bPrintJetSelection)
1583 if(fDebug > 7) printf("rejected (track pt)\n");
1587 if(bPrintJetSelection)
1588 if(fDebug > 7) printf("accepted\n");
1589 if(fDebug > 5) printf("TaskV0sInJets: Jet %d with pt %.2f passed selection\n", iJet, jetSel->Pt());
1593 // new ((*fBranchJet)[iNJetAll++]) AliAODJet(*((AliAODJet*)jetSel));
1594 objectJet = new ((*fBranchJet)[iNJetAll++]) AliJetObject(jetSel); // copy selected jet to the array
1595 // objectJet->SetPtLeadingTrack(dPtLeadTrack);
1596 objectJet->SetRadius(fdRadiusJet);
1600 TLorentzVector vecPerpPlus(*(jetSel->MomentumVector()));
1601 vecPerpPlus.RotateZ(TMath::Pi() / 2.); // rotate vector by 90 deg around z
1602 TLorentzVector vecPerpMinus(*(jetSel->MomentumVector()));
1603 vecPerpMinus.RotateZ(-TMath::Pi() / 2.); // rotate vector by -90 deg around z
1604 // AliAODJet jetTmp = AliAODJet(vecPerp);
1605 if(fDebug > 5) printf("TaskV0sInJets: Adding perp. cones number %d, %d\n", iNJetPerp, iNJetPerp + 1);
1606 // printf("TaskV0sInJets: Adding perp. cone number %d: pT %f, phi %f, eta %f, pT %f, phi %f, eta %f\n",iNJetSel,vecPerp.Pt(),vecPerp.Phi(),vecPerp.Eta(),jetTmp.Pt(),jetTmp.Phi(),jetTmp.Eta());
1607 new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpPlus); // write perp. cone to the array
1608 new((*jetArrayPerp)[iNJetPerp++]) AliAODJet(vecPerpMinus); // write perp. cone to the array
1609 if(fDebug > 5) printf("TaskV0sInJets: Adding jet number %d\n", iNJetSel);
1610 // printf("TaskV0sInJets: Adding jet number %d: pT %f, phi %f, eta %f\n",iNJetSel,jetSel->Pt(),jetSel->Phi(),jetSel->Eta());
1611 new((*jetArraySel)[iNJetSel++]) AliAODJet(*((AliAODJet*)jetSel)); // copy selected jet to the array
1613 if(fDebug > 5) printf("TaskV0sInJets: Added jets: %d\n", iNJetSel);
1614 iNJetSel = jetArraySel->GetEntriesFast();
1615 if(fDebug > 2) printf("TaskV0sInJets: Selected jets in array: %d\n", iNJetSel);
1616 fh1NJetPerEvent[iCentIndex]->Fill(iNJetSel);
1618 for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
1620 jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
1621 fh1PtJet[iCentIndex]->Fill(jet->Pt()); // pt spectrum of selected jets
1622 fh1EtaJet[iCentIndex]->Fill(jet->Eta()); // eta spectrum of selected jets
1623 fh2EtaPtJet[iCentIndex]->Fill(jet->Eta(), jet->Pt()); // eta-pT spectrum of selected jets
1624 fh1PhiJet[iCentIndex]->Fill(jet->Phi()); // phi spectrum of selected jets
1625 Double_t dAreaExcluded = TMath::Pi() * dRadiusExcludeCone * dRadiusExcludeCone; // area of the cone
1626 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax - jet->Eta()); // positive eta overhang
1627 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone, dEtaMax + jet->Eta()); // negative eta overhang
1628 fh1AreaExcluded->Fill(iCentIndex, dAreaExcluded);
1633 if(bJetEventGood) // there should be some reconstructed jets
1635 fh1EventCounterCut->Fill(4); // events with jet(s)
1636 fh1EventCounterCutCent[iCentIndex]->Fill(4); // events with jet(s)
1639 fh1EventCounterCut->Fill(5); // events with selected jets
1640 fh1EventCounterCutCent[iCentIndex]->Fill(5);
1644 fh1EventCent2Jets->Fill(fdCentrality);
1646 fh1EventCent2NoJets->Fill(fdCentrality);
1650 jetRnd = GetRandomCone(jetArraySel, dJetEtaWindow, 2 * fdRadiusJet);
1653 fh1NRndConeCent->Fill(iCentIndex);
1654 fh2EtaPhiRndCone[iCentIndex]->Fill(jetRnd->Eta(), jetRnd->Phi());
1656 jetMed = GetMedianCluster(jetArrayBg, dJetEtaWindow);
1659 fh1NMedConeCent->Fill(iCentIndex);
1660 fh2EtaPhiMedCone[iCentIndex]->Fill(jetMed->Eta(), jetMed->Phi());
1664 // Loading primary vertex info
1665 AliAODVertex* primVtx = fAODIn->GetPrimaryVertex(); // get the primary vertex
1666 Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
1667 primVtx->GetXYZ(dPrimVtxPos);
1668 fh1VtxZ[iCentIndex]->Fill(dPrimVtxPos[2]);
1669 fh2VtxXY[iCentIndex]->Fill(dPrimVtxPos[0], dPrimVtxPos[1]);
1671 //===== Start of loop over V0 candidates =====
1672 if(fDebug > 2) printf("TaskV0sInJets: Start of V0 loop\n");
1673 for(Int_t iV0 = 0; iV0 < iNV0s; iV0++)
1675 v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1681 // Initialization of status indicators
1682 Bool_t bIsCandidateK0s = kTRUE; // candidate for K0s
1683 Bool_t bIsCandidateLambda = kTRUE; // candidate for Lambda
1684 Bool_t bIsCandidateALambda = kTRUE; // candidate for Lambda
1685 Bool_t bIsInPeakK0s = kFALSE; // candidate within the K0s mass peak
1686 Bool_t bIsInPeakLambda = kFALSE; // candidate within the Lambda mass peak
1687 Bool_t bIsInConeJet = kFALSE; // candidate within the jet cones
1688 Bool_t bIsInConePerp = kFALSE; // candidate within the perpendicular cone
1689 Bool_t bIsInConeRnd = kFALSE; // candidate within the random cone
1690 Bool_t bIsInConeMed = kFALSE; // candidate within the median-cluster cone
1691 Bool_t bIsOutsideCones = kFALSE; // candidate outside excluded cones
1693 // Invariant mass calculation
1694 dMassV0K0s = v0->MassK0Short();
1695 dMassV0Lambda = v0->MassLambda();
1696 dMassV0ALambda = v0->MassAntiLambda();
1698 Int_t iCutIndex = 0; // indicator of current selection step
1700 // All V0 candidates
1701 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1704 // Skip candidates outside the histogram range
1705 if((dMassV0K0s < fgkdMassK0sMin) || (dMassV0K0s >= fgkdMassK0sMax))
1706 bIsCandidateK0s = kFALSE;
1707 if((dMassV0Lambda < fgkdMassLambdaMin) || (dMassV0Lambda >= fgkdMassLambdaMax))
1708 bIsCandidateLambda = kFALSE;
1709 if((dMassV0ALambda < fgkdMassLambdaMin) || (dMassV0ALambda >= fgkdMassLambdaMax))
1710 bIsCandidateALambda = kFALSE;
1711 if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
1714 Double_t dPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
1715 vecV0Momentum = TVector3(v0->Px(), v0->Py(), v0->Pz()); // set the vector of V0 momentum
1717 // Sigma of the mass peak window
1718 Double_t dMassPeakWindowK0s = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 0);
1719 Double_t dMassPeakWindowLambda = dNSigmaMassMax * MassPeakSigmaOld(dPtV0, 1);
1720 // Double_t dMassPeakWindowK0s = dNSigmaMassMax*MassPeakSigma(iCentIndex,dPtV0,0);
1721 // Double_t dMassPeakWindowLambda = dNSigmaMassMax*MassPeakSigma(iCentIndex,dPtV0,1);
1723 // Invariant mass peak selection
1724 if(TMath::Abs(dMassV0K0s - dMassPDGK0s) < dMassPeakWindowK0s)
1725 bIsInPeakK0s = kTRUE;
1726 if(TMath::Abs(dMassV0Lambda - dMassPDGLambda) < dMassPeakWindowLambda)
1727 bIsInPeakLambda = kTRUE;
1729 // Retrieving all relevant properties of the V0 candidate
1730 Bool_t bOnFlyStatus = v0->GetOnFlyStatus(); // online (on fly) reconstructed vs offline reconstructed
1731 const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1732 const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
1733 Double_t dPtDaughterPos = trackPos->Pt(); // transverse momentum of a daughter track
1734 Double_t dPtDaughterNeg = trackNeg->Pt();
1735 Double_t dNRowsPos = trackPos->GetTPCClusterInfo(2, 1); // crossed TPC pad rows of a daughter track
1736 Double_t dNRowsNeg = trackNeg->GetTPCClusterInfo(2, 1);
1737 Double_t dDCAToPrimVtxPos = TMath::Abs(v0->DcaPosToPrimVertex()); // dca of a daughter to the primary vertex
1738 Double_t dDCAToPrimVtxNeg = TMath::Abs(v0->DcaNegToPrimVertex());
1739 Double_t dDCADaughters = v0->DcaV0Daughters(); // dca between daughters
1740 Double_t dCPA = v0->CosPointingAngle(primVtx); // cosine of the pointing angle
1741 Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
1742 // Double_t dSecVtxPos[3] = {v0->DecayVertexV0X(),v0->DecayVertexV0Y(),v0->DecayVertexV0Z()}; // V0 vertex position
1743 v0->GetSecondaryVtx(dSecVtxPos);
1744 Double_t dRadiusDecay = TMath::Sqrt(dSecVtxPos[0] * dSecVtxPos[0] + dSecVtxPos[1] * dSecVtxPos[1]); // distance of the V0 vertex from the z-axis
1745 Double_t dEtaDaughterNeg = trackNeg->Eta(); // = v0->EtaProng(1), pseudorapidity of a daughter track
1746 Double_t dEtaDaughterPos = trackPos->Eta(); // = v0->EtaProng(0)
1747 Double_t dRapK0s = v0->RapK0Short(); // rapidity calculated for K0s assumption
1748 Double_t dRapLambda = v0->RapLambda(); // rapidity calculated for Lambda assumption
1749 Double_t dEtaV0 = v0->Eta(); // V0 pseudorapidity
1750 // Double_t dPhiV0 = v0->Phi(); // V0 pseudorapidity
1751 Double_t dDecayPath[3];
1752 for(Int_t iPos = 0; iPos < 3; iPos++)
1753 dDecayPath[iPos] = dSecVtxPos[iPos] - dPrimVtxPos[iPos]; // vector of the V0 path
1754 Double_t dDecLen = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]); // path length L
1755 Double_t dDecLen2D = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]); // transverse path length R
1756 Double_t dLOverP = dDecLen / v0->P(); // L/p
1757 Double_t dROverPt = dDecLen2D / dPtV0; // R/pT
1758 Double_t dMLOverPK0s = dMassPDGK0s * dLOverP; // m*L/p = c*(proper lifetime)
1759 // Double_t dMLOverPLambda = dMassPDGLambda*dLOverP; // m*L/p
1760 Double_t dMROverPtK0s = dMassPDGK0s * dROverPt; // m*R/pT
1761 Double_t dMROverPtLambda = dMassPDGLambda * dROverPt; // m*R/pT
1762 Double_t dNSigmaPosPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kPion)); // difference between measured and expected signal of the dE/dx in the TPC
1763 Double_t dNSigmaPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackPos, AliPID::kProton));
1764 Double_t dNSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kPion));
1765 Double_t dNSigmaNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(trackNeg, AliPID::kProton));
1766 Double_t dAlpha = v0->AlphaV0(); // Armenteros-Podolanski alpha
1767 Double_t dPtArm = v0->PtArmV0(); // Armenteros-Podolanski pT
1768 AliAODVertex* prodVtxDaughterPos = (AliAODVertex*)(trackPos->GetProdVertex()); // production vertex of the positive daughter track
1769 Char_t cTypeVtxProdPos = prodVtxDaughterPos->GetType(); // type of the production vertex
1770 AliAODVertex* prodVtxDaughterNeg = (AliAODVertex*)(trackNeg->GetProdVertex()); // production vertex of the negative daughter track
1771 Char_t cTypeVtxProdNeg = prodVtxDaughterNeg->GetType(); // type of the production vertex
1773 // fh2Tau3DVs2D[0]->Fill(dPtV0, dLOverP / dROverPt);
1775 // QA histograms before cuts
1776 FillQAHistogramV0(primVtx, v0, 0, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
1777 // Cut vs mass histograms before cuts
1781 fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsPos);
1782 fh2CutTPCRowsK0s[0]->Fill(dMassV0K0s, dNRowsNeg);
1783 fh2CutPtPosK0s[0]->Fill(dMassV0K0s, dPtDaughterPos);
1784 fh2CutPtNegK0s[0]->Fill(dMassV0K0s, dPtDaughterNeg);
1785 fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
1786 fh2CutDCAVtx[0]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
1787 fh2CutDCAV0[0]->Fill(dMassV0K0s, dDCADaughters);
1788 fh2CutCos[0]->Fill(dMassV0K0s, dCPA);
1789 fh2CutR[0]->Fill(dMassV0K0s, dRadiusDecay);
1790 fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterPos);
1791 fh2CutEtaK0s[0]->Fill(dMassV0K0s, dEtaDaughterNeg);
1792 fh2CutRapK0s[0]->Fill(dMassV0K0s, dRapK0s);
1793 fh2CutCTauK0s[0]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
1794 fh2CutPIDPosK0s[0]->Fill(dMassV0K0s, dNSigmaPosPion);
1795 fh2CutPIDNegK0s[0]->Fill(dMassV0K0s, dNSigmaNegPion);
1797 if(bIsCandidateLambda)
1799 fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsPos);
1800 fh2CutTPCRowsLambda[0]->Fill(dMassV0Lambda, dNRowsNeg);
1801 fh2CutPtPosLambda[0]->Fill(dMassV0Lambda, dPtDaughterPos);
1802 fh2CutPtNegLambda[0]->Fill(dMassV0Lambda, dPtDaughterNeg);
1803 fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterPos);
1804 fh2CutEtaLambda[0]->Fill(dMassV0Lambda, dEtaDaughterNeg);
1805 fh2CutRapLambda[0]->Fill(dMassV0Lambda, dRapLambda);
1806 fh2CutCTauLambda[0]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
1807 fh2CutPIDPosLambda[0]->Fill(dMassV0Lambda, dNSigmaPosProton);
1808 fh2CutPIDNegLambda[0]->Fill(dMassV0Lambda, dNSigmaNegPion);
1812 //===== Start of reconstruction cutting =====
1815 // All V0 candidates
1816 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1819 // Start of global cuts
1821 // Reconstruction method
1822 if(bPrintCuts) printf("Rec: Applying cut: Reconstruction method: on-the-fly? %s\n", (bOnFly ? "yes" : "no"));
1823 if(bOnFlyStatus != bOnFly)
1825 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1830 if(bPrintCuts) printf("Rec: Applying cut: Correct charge of daughters\n");
1831 if(!trackNeg || !trackPos)
1833 if(trackNeg->Charge() == trackPos->Charge()) // daughters have different charge?
1835 if(trackNeg->Charge() != -1) // daughters have expected charge?
1837 if(trackPos->Charge() != 1) // daughters have expected charge?
1840 if(bPrintCuts) printf("Rec: Applying cut: TPC refit: %d\n", iRefit);
1841 if(!trackNeg->IsOn(iRefit)) // TPC refit is ON?
1843 if(bPrintCuts) printf("Rec: Applying cut: Type of production vertex of daughter: Not %d\n", AliAODVertex::kKink);
1844 if(cTypeVtxProdNeg == AliAODVertex::kKink) // kink daughter rejection
1849 if(bPrintCuts) printf("Rec: Applying cut: Number of TPC rows: > %f\n", dNCrossedRowsTPCMin);
1850 if(dNRowsNeg < dNCrossedRowsTPCMin) // Crossed TPC padrows
1852 // Int_t findable = trackNeg->GetTPCNclsF(); // Findable clusters
1853 // if (findable <= 0)
1855 // if (dNRowsNeg/findable < dCrossedRowsOverFindMin)
1857 // if (dNRowsNeg/findable > dCrossedRowsOverFindMax)
1862 if(!trackPos->IsOn(iRefit))
1864 if(cTypeVtxProdPos == AliAODVertex::kKink) // kink daughter rejection
1869 if(dNRowsPos < dNCrossedRowsTPCMin)
1871 // findable = trackPos->GetTPCNclsF();
1872 // if (findable <= 0)
1874 // if (dNRowsPos/findable < dCrossedRowsOverFindMin)
1876 // if (dNRowsPos/findable > dCrossedRowsOverFindMax)
1881 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1885 // Daughters: transverse momentum cut
1888 if(bPrintCuts) printf("Rec: Applying cut: Daughter pT: > %f\n", dPtDaughterMin);
1889 if((dPtDaughterNeg < dPtDaughterMin) || (dPtDaughterPos < dPtDaughterMin))
1891 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1896 // Daughters: Impact parameter of daughters to prim vtx
1897 if(bPrintCuts) printf("Rec: Applying cut: Daughter DCA to prim vtx: > %f\n", dDCAToPrimVtxMin);
1898 if((dDCAToPrimVtxNeg < dDCAToPrimVtxMin) || (dDCAToPrimVtxPos < dDCAToPrimVtxMin))
1900 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1905 if(bPrintCuts) printf("Rec: Applying cut: DCA between daughters: < %f\n", dDCADaughtersMax);
1906 if(dDCADaughters > dDCADaughtersMax)
1908 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1912 // V0: Cosine of the pointing angle
1913 if(bPrintCuts) printf("Rec: Applying cut: CPA: > %f\n", dCPAMin);
1916 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1920 // V0: Fiducial volume
1921 if(bPrintCuts) printf("Rec: Applying cut: Decay radius: > %f, < %f\n", dRadiusDecayMin, dRadiusDecayMax);
1922 if((dRadiusDecay < dRadiusDecayMin) || (dRadiusDecay > dRadiusDecayMax))
1924 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1928 // Daughters: pseudorapidity cut
1931 if(bPrintCuts) printf("Rec: Applying cut: Daughter |eta|: < %f\n", dEtaDaughterMax);
1932 if((TMath::Abs(dEtaDaughterNeg) > dEtaDaughterMax) || (TMath::Abs(dEtaDaughterPos) > dEtaDaughterMax))
1934 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1937 // End of global cuts
1939 // Start of particle-dependent cuts
1941 // V0: rapidity cut & pseudorapidity cut
1944 if(bPrintCuts) printf("Rec: Applying cut: V0 |y|: < %f\n", dRapMax);
1945 if(TMath::Abs(dRapK0s) > dRapMax)
1946 bIsCandidateK0s = kFALSE;
1947 if(TMath::Abs(dRapLambda) > dRapMax)
1949 bIsCandidateLambda = kFALSE;
1950 bIsCandidateALambda = kFALSE;
1955 if(bPrintCuts) printf("Rec: Applying cut: V0 |eta|: < %f\n", dEtaMax);
1956 if(TMath::Abs(dEtaV0) > dEtaMax)
1958 bIsCandidateK0s = kFALSE;
1959 bIsCandidateLambda = kFALSE;
1960 bIsCandidateALambda = kFALSE;
1962 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1970 if(bPrintCuts) printf("Rec: Applying cut: Proper lifetime: < %f\n", dNTauMax);
1971 if(dMROverPtK0s > dNTauMax * dCTauK0s)
1972 bIsCandidateK0s = kFALSE;
1973 if(dMROverPtLambda > dNTauMax * dCTauLambda)
1975 bIsCandidateLambda = kFALSE;
1976 bIsCandidateALambda = kFALSE;
1978 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
1988 if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (both daughters): < %f\n", dNSigmadEdxMax);
1989 if(dNSigmaPosPion > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax) // pi+, pi-
1990 bIsCandidateK0s = kFALSE;
1991 if(dNSigmaPosProton > dNSigmadEdxMax || dNSigmaNegPion > dNSigmadEdxMax) // p+, pi-
1992 bIsCandidateLambda = kFALSE;
1993 if(dNSigmaNegProton > dNSigmadEdxMax || dNSigmaPosPion > dNSigmadEdxMax) // p-, pi+
1994 bIsCandidateALambda = kFALSE;
1998 if(bPrintCuts) printf("Rec: Applying cut: Delta dE/dx (proton below %f GeV/c): < %f\n", dPtProtonPIDMax, dNSigmadEdxMax);
1999 if((dPtDaughterPos < dPtProtonPIDMax) && (dNSigmaPosProton > dNSigmadEdxMax)) // p+
2000 bIsCandidateLambda = kFALSE;
2001 if((dPtDaughterNeg < dPtProtonPIDMax) && (dNSigmaNegProton > dNSigmadEdxMax)) // p-
2002 bIsCandidateALambda = kFALSE;
2004 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2008 Double_t valueCorrel[3] = {dMassV0K0s, dMassV0Lambda, dPtV0};
2009 if(bIsCandidateK0s && bIsCandidateLambda)
2010 fh3CCMassCorrelBoth->Fill(valueCorrel); // correlation of mass distribution of candidates selected as both K0s and Lambda
2011 if(bIsCandidateK0s && !bIsCandidateLambda)
2012 fh3CCMassCorrelKNotL->Fill(valueCorrel); // correlation of mass distribution of candidates selected as K0s and not Lambda
2013 if(!bIsCandidateK0s && bIsCandidateLambda)
2014 fh3CCMassCorrelLNotK->Fill(valueCorrel); // correlation of mass distribution of candidates selected as not K0s and Lambda
2017 // Armenteros-Podolanski cut
2020 if(bPrintCuts) printf("Rec: Applying cut: Armenteros-Podolanski (K0S): pT > %f * |alpha|\n", 0.2);
2021 if(dPtArm < TMath::Abs(0.2 * dAlpha))
2022 bIsCandidateK0s = kFALSE;
2023 // if(dPtArm < 0.025)
2025 // bIsCandidateLambda = kFALSE;
2026 // bIsCandidateALambda = kFALSE;
2028 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2032 // Cross contamination
2035 if(bIsCandidateLambda) // Lambda candidates in K0s peak, excluded from Lambda candidates by CC cut
2036 fh2CCLambda->Fill(dMassV0Lambda, dPtV0);
2040 if(bIsCandidateK0s) // K0s candidates in Lambda peak, excluded from K0s candidates by CC cut
2041 fh2CCK0s->Fill(dMassV0K0s, dPtV0);
2045 // if (bIsInPeakK0s)
2046 // bIsCandidateLambda = kFALSE;
2047 // if (bIsInPeakLambda)
2048 // bIsCandidateK0s = kFALSE;
2049 // FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex, iCentIndex);
2053 // End of particle-dependent cuts
2055 //===== End of reconstruction cutting =====
2057 if(!bIsCandidateK0s && !bIsCandidateLambda && !bIsCandidateALambda)
2061 if(fDebug>5) printf("TaskV0sInJets: Adding selected V0 to branch\n");
2062 // Add selected candidates to the output tree branch
2063 if ((bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda) && fbTreeOutput)
2065 objectV0 = new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliV0Object(v0,primVtx);
2066 // new ((*fBranchV0Rec)[iNV0SelV0Rec++]) AliAODv0(*((AliAODv0*)v0));
2067 objectV0->SetIsCandidateK0S(bIsCandidateK0s);
2068 objectV0->SetIsCandidateLambda(bIsCandidateLambda);
2069 objectV0->SetIsCandidateALambda(bIsCandidateALambda);
2070 objectV0->SetNSigmaPosProton(dNSigmaPosProton);
2071 objectV0->SetNSigmaNegProton(dNSigmaNegProton);
2075 // Selection of V0s in jet cones, perpendicular cones, random cones, outside cones
2076 if(bJetEventGood && iNJetSel && (bIsCandidateK0s || bIsCandidateLambda || bIsCandidateALambda))
2078 // Selection of V0s in jet cones
2079 if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in %d jet cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetSel);
2080 for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2082 jet = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2083 vecJetMomentum = TVector3(jet->Px(), jet->Py(), jet->Pz()); // set the vector of jet momentum
2084 if(fDebug > 5) printf("TaskV0sInJets: Checking if V0 %d %d in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2085 if(IsParticleInCone(v0, jet, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2087 if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in jet cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2088 bIsInConeJet = kTRUE;
2092 // Selection of V0s in perp. cones
2093 if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in %d perp. cones\n", bIsCandidateK0s, bIsCandidateLambda, iNJetSel);
2094 for(Int_t iJet = 0; iJet < iNJetPerp; iJet++)
2096 jetPerp = (AliAODJet*)jetArrayPerp->At(iJet); // load a jet in the list
2097 if(fDebug > 5) printf("TaskV0sInJets: Checking if V0 %d %d in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2098 if(IsParticleInCone(v0, jetPerp, fdRadiusJet)) // V0 in perp. cone
2100 if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in perp. cone %d\n", bIsCandidateK0s, bIsCandidateLambda, iJet);
2101 bIsInConePerp = kTRUE;
2105 // Selection of V0s in random cones
2108 if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2109 if(IsParticleInCone(v0, jetRnd, fdRadiusJet)) // V0 in rnd. cone?
2111 if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in the rnd. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2112 bIsInConeRnd = kTRUE;
2115 // Selection of V0s in median-cluster cones
2118 if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2119 if(IsParticleInCone(v0, jetMed, fdRadiusJet)) // V0 in med. cone?
2121 if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found in the med. cone\n", bIsCandidateK0s, bIsCandidateLambda);
2122 bIsInConeMed = kTRUE;
2125 // Selection of V0s outside jet cones
2126 if(fDebug > 5) printf("TaskV0sInJets: Searching for V0 %d %d outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2127 if(!OverlapWithJets(jetArraySel, v0, dRadiusExcludeCone)) // V0 oustide jet cones
2129 if(fDebug > 5) printf("TaskV0sInJets: V0 %d %d found outside jet cones\n", bIsCandidateK0s, bIsCandidateLambda);
2130 bIsOutsideCones = kTRUE;
2134 // QA histograms after cuts
2135 FillQAHistogramV0(primVtx, v0, 1, bIsCandidateK0s, bIsCandidateLambda, bIsInPeakK0s, bIsInPeakLambda);
2136 // Cut vs mass histograms after cuts
2140 fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsPos);
2141 fh2CutTPCRowsK0s[1]->Fill(dMassV0K0s, dNRowsNeg);
2142 fh2CutPtPosK0s[1]->Fill(dMassV0K0s, dPtDaughterPos);
2143 fh2CutPtNegK0s[1]->Fill(dMassV0K0s, dPtDaughterNeg);
2144 fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxPos);
2145 fh2CutDCAVtx[1]->Fill(dMassV0K0s, dDCAToPrimVtxNeg);
2146 fh2CutDCAV0[1]->Fill(dMassV0K0s, dDCADaughters);
2147 fh2CutCos[1]->Fill(dMassV0K0s, dCPA);
2148 fh2CutR[1]->Fill(dMassV0K0s, dRadiusDecay);
2149 fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterPos);
2150 fh2CutEtaK0s[1]->Fill(dMassV0K0s, dEtaDaughterNeg);
2151 fh2CutRapK0s[1]->Fill(dMassV0K0s, dRapK0s);
2152 fh2CutCTauK0s[1]->Fill(dMassV0K0s, dMROverPtK0s / dCTauK0s);
2153 fh2CutPIDPosK0s[1]->Fill(dMassV0K0s, dNSigmaPosPion);
2154 fh2CutPIDNegK0s[1]->Fill(dMassV0K0s, dNSigmaNegPion);
2155 fh1DeltaZK0s[iCentIndex]->Fill(dDecayPath[2]);
2157 if(bIsCandidateLambda)
2159 fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsPos);
2160 fh2CutTPCRowsLambda[1]->Fill(dMassV0Lambda, dNRowsNeg);
2161 fh2CutPtPosLambda[1]->Fill(dMassV0Lambda, dPtDaughterPos);
2162 fh2CutPtNegLambda[1]->Fill(dMassV0Lambda, dPtDaughterNeg);
2163 fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterPos);
2164 fh2CutEtaLambda[1]->Fill(dMassV0Lambda, dEtaDaughterNeg);
2165 fh2CutRapLambda[1]->Fill(dMassV0Lambda, dRapLambda);
2166 fh2CutCTauLambda[1]->Fill(dMassV0Lambda, dMROverPtLambda / dCTauLambda);
2167 fh2CutPIDPosLambda[1]->Fill(dMassV0Lambda, dNSigmaPosProton);
2168 fh2CutPIDNegLambda[1]->Fill(dMassV0Lambda, dNSigmaNegPion);
2169 fh1DeltaZLambda[iCentIndex]->Fill(dDecayPath[2]);
2173 //===== Start of filling V0 spectra =====
2175 Double_t dAngle = TMath::Pi(); // angle between V0 momentum and jet momentum
2178 dAngle = vecV0Momentum.Angle(vecJetMomentum);
2184 // 14 K0s candidates after cuts
2185 // printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0K0s,dPtV0,dEtaV0,dPhiV0);
2186 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex, iCentIndex);
2187 Double_t valueKIncl[3] = {dMassV0K0s, dPtV0, dEtaV0};
2188 fhnV0InclusiveK0s[iCentIndex]->Fill(valueKIncl);
2189 fh1V0InvMassK0sCent[iCentIndex]->Fill(dMassV0K0s);
2191 fh1QACTau2D[1]->Fill(dMROverPtK0s / dCTauK0s);
2192 fh1QACTau3D[1]->Fill(dMLOverPK0s / dCTauK0s);
2193 // fh2Tau3DVs2D[1]->Fill(dPtV0, dLOverP / dROverPt);
2197 // 15 K0s in jet events
2198 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 1, iCentIndex);
2203 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, bIsCandidateK0s, kFALSE, kFALSE, iCutIndex + 2, iCentIndex);
2204 Double_t valueKInJC[4] = {dMassV0K0s, dPtV0, dEtaV0, jet->Pt()};
2205 fhnV0InJetK0s[iCentIndex]->Fill(valueKInJC);
2206 fh2V0PtJetAngleK0s[iCentIndex]->Fill(jet->Pt(), dAngle);
2210 Double_t valueKOutJC[3] = {dMassV0K0s, dPtV0, dEtaV0};
2211 fhnV0OutJetK0s[iCentIndex]->Fill(valueKOutJC);
2215 Double_t valueKInPC[4] = {dMassV0K0s, dPtV0, dEtaV0, jetPerp->Pt()};
2216 fhnV0InPerpK0s[iCentIndex]->Fill(valueKInPC);
2220 Double_t valueKInRnd[3] = {dMassV0K0s, dPtV0, dEtaV0};
2221 fhnV0InRndK0s[iCentIndex]->Fill(valueKInRnd);
2225 Double_t valueKInMed[3] = {dMassV0K0s, dPtV0, dEtaV0};
2226 fhnV0InMedK0s[iCentIndex]->Fill(valueKInMed);
2230 Double_t valueKNoJet[3] = {dMassV0K0s, dPtV0, dEtaV0};
2231 fhnV0NoJetK0s[iCentIndex]->Fill(valueKNoJet);
2235 if(bIsCandidateLambda)
2237 // 14 Lambda candidates after cuts
2238 // printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0Lambda,dPtV0,dEtaV0,dPhiV0);
2239 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex, iCentIndex);
2240 Double_t valueLIncl[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2241 fhnV0InclusiveLambda[iCentIndex]->Fill(valueLIncl);
2242 fh1V0InvMassLambdaCent[iCentIndex]->Fill(dMassV0Lambda);
2245 // 15 Lambda in jet events
2246 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 1, iCentIndex);
2250 // 16 Lambda in jets
2251 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, bIsCandidateLambda, kFALSE, iCutIndex + 2, iCentIndex);
2252 Double_t valueLInJC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jet->Pt()};
2253 fhnV0InJetLambda[iCentIndex]->Fill(valueLInJC);
2254 fh2V0PtJetAngleLambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2258 Double_t valueLOutJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2259 fhnV0OutJetLambda[iCentIndex]->Fill(valueLOutJet);
2263 Double_t valueLInPC[4] = {dMassV0Lambda, dPtV0, dEtaV0, jetPerp->Pt()};
2264 fhnV0InPerpLambda[iCentIndex]->Fill(valueLInPC);
2268 Double_t valueLInRnd[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2269 fhnV0InRndLambda[iCentIndex]->Fill(valueLInRnd);
2273 Double_t valueLInMed[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2274 fhnV0InMedLambda[iCentIndex]->Fill(valueLInMed);
2278 Double_t valueLNoJet[3] = {dMassV0Lambda, dPtV0, dEtaV0};
2279 fhnV0NoJetLambda[iCentIndex]->Fill(valueLNoJet);
2283 if(bIsCandidateALambda)
2285 // 14 ALambda candidates after cuts
2286 // printf("AL: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",iV0,dMassV0ALambda,dPtV0,dEtaV0,dPhiV0);
2287 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex, iCentIndex);
2288 Double_t valueALIncl[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2289 fhnV0InclusiveALambda[iCentIndex]->Fill(valueALIncl);
2290 fh1V0InvMassALambdaCent[iCentIndex]->Fill(dMassV0ALambda);
2293 // 15 ALambda in jet events
2294 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 1, iCentIndex);
2298 // 16 ALambda in jets
2299 FillCandidates(dMassV0K0s, dMassV0Lambda, dMassV0ALambda, kFALSE, kFALSE, bIsCandidateALambda, iCutIndex + 2, iCentIndex);
2300 Double_t valueLInJC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jet->Pt()};
2301 fhnV0InJetALambda[iCentIndex]->Fill(valueLInJC);
2302 fh2V0PtJetAngleALambda[iCentIndex]->Fill(jet->Pt(), dAngle);
2306 Double_t valueALOutJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2307 fhnV0OutJetALambda[iCentIndex]->Fill(valueALOutJet);
2311 Double_t valueLInPC[4] = {dMassV0ALambda, dPtV0, dEtaV0, jetPerp->Pt()};
2312 fhnV0InPerpALambda[iCentIndex]->Fill(valueLInPC);
2316 Double_t valueALInRnd[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2317 fhnV0InRndALambda[iCentIndex]->Fill(valueALInRnd);
2321 Double_t valueALInMed[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2322 fhnV0InMedALambda[iCentIndex]->Fill(valueALInMed);
2326 Double_t valueALNoJet[3] = {dMassV0ALambda, dPtV0, dEtaV0};
2327 fhnV0NoJetALambda[iCentIndex]->Fill(valueALNoJet);
2331 //===== End of filling V0 spectra =====
2334 //===== Association of reconstructed V0 candidates with MC particles =====
2337 // Associate selected candidates only
2338 // if ( !(bIsCandidateK0s && bIsInPeakK0s) && !(bIsCandidateLambda && bIsInPeakLambda) ) // signal candidates
2339 if(!(bIsCandidateK0s) && !(bIsCandidateLambda) && !(bIsCandidateALambda)) // chosen candidates with any mass
2342 // Get MC labels of reconstructed daughter tracks
2343 Int_t iLabelPos = TMath::Abs(trackPos->GetLabel());
2344 Int_t iLabelNeg = TMath::Abs(trackNeg->GetLabel());
2346 // Make sure MC daughters are in the array range
2347 if((iLabelNeg < 0) || (iLabelNeg >= iNTracksMC) || (iLabelPos < 0) || (iLabelPos >= iNTracksMC))
2350 // Get MC particles corresponding to reconstructed daughter tracks
2351 AliAODMCParticle* particleMCDaughterNeg = (AliAODMCParticle*)arrayMC->At(iLabelNeg);
2352 AliAODMCParticle* particleMCDaughterPos = (AliAODMCParticle*)arrayMC->At(iLabelPos);
2353 if(!particleMCDaughterNeg || !particleMCDaughterPos)
2356 // Make sure MC daughter particles are not physical primary
2357 if((particleMCDaughterNeg->IsPhysicalPrimary()) || (particleMCDaughterPos->IsPhysicalPrimary()))
2360 // Get identities of MC daughter particles
2361 Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
2362 Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
2364 // Get index of the mother particle for each MC daughter particle
2365 Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
2366 Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
2368 if((iIndexMotherNeg < 0) || (iIndexMotherNeg >= iNTracksMC) || (iIndexMotherPos < 0) || (iIndexMotherPos >= iNTracksMC))
2371 // Check whether MC daughter particles have the same mother
2372 if(iIndexMotherNeg != iIndexMotherPos)
2375 // Get the MC mother particle of both MC daughter particles
2376 AliAODMCParticle* particleMCMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherPos);
2377 if(!particleMCMother)
2380 // Get identity of the MC mother particle
2381 Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
2383 // Skip not interesting particles
2384 if((iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda))
2387 // Check identity of the MC mother particle and the decay channel
2388 // Is MC mother particle K0S?
2389 Bool_t bV0MCIsK0s = ((iPdgCodeMother == iPdgCodeK0s) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2390 // Is MC mother particle Lambda?
2391 Bool_t bV0MCIsLambda = ((iPdgCodeMother == +iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodeProton) && (iPdgCodeDaughterNeg == -iPdgCodePion));
2392 // Is MC mother particle anti-Lambda?
2393 Bool_t bV0MCIsALambda = ((iPdgCodeMother == -iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodeProton));
2395 Double_t dPtV0Gen = particleMCMother->Pt();
2396 // Double_t dRapV0MC = particleMCMother->Y();
2397 Double_t dEtaV0Gen = particleMCMother->Eta();
2398 // Double_t dPhiV0Gen = particleMCMother->Phi();
2400 // Is MC mother particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2401 // Bool_t bV0MCIsPrimary = particleMCMother->IsPhysicalPrimary();
2402 // Get the MC mother particle of the MC mother particle
2403 Int_t iIndexMotherOfMother = particleMCMother->GetMother();
2404 AliAODMCParticle* particleMCMotherOfMother = 0;
2405 if(iIndexMotherOfMother >= 0)
2406 particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2407 // Get identity of the MC mother particle of the MC mother particle if it exists
2408 Int_t iPdgCodeMotherOfMother = 0;
2409 if(particleMCMotherOfMother)
2410 iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2411 // Check if the MC mother particle of the MC mother particle is a physical primary Sigma (3212 - Sigma0, 3224 - Sigma*+, 3214 - Sigma*0, 3114 - Sigma*-)
2412 // Bool_t bV0MCComesFromSigma = kFALSE; // Is MC mother particle daughter of a Sigma?
2413 // if ( (particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && ( (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) ) )
2414 // bV0MCComesFromSigma = kTRUE;
2415 // Should MC mother particle be considered as primary when it is Lambda?
2416 // Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2417 // Check if the MC mother particle of the MC mother particle is a Xi (3322 - Xi0, 3312 - Xi-)
2418 Bool_t bV0MCComesFromXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == 3322) || (iPdgCodeMotherOfMother == 3312))); // Is MC mother particle daughter of a Xi?
2419 Bool_t bV0MCComesFromAXi = ((particleMCMotherOfMother) && ((iPdgCodeMotherOfMother == -3322) || (iPdgCodeMotherOfMother == -3312))); // Is MC mother particle daughter of a anti-Xi?
2421 // Get the distance between production point of the MC mother particle and the primary vertex
2422 Double_t dx = dPrimVtxMCX - particleMCMother->Xv();
2423 Double_t dy = dPrimVtxMCY - particleMCMother->Yv();
2424 Double_t dz = dPrimVtxMCZ - particleMCMother->Zv();
2425 Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2426 Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2428 // phi, eta resolution for V0-reconstruction
2429 // Double_t dResolutionV0Eta = particleMCMother->Eta()-v0->Eta();
2430 // Double_t dResolutionV0Phi = particleMCMother->Phi()-v0->Phi();
2435 objectV0->SetPtTrue(dPtV0Gen);
2436 objectV0->SetEtaTrue(dEtaV0Gen);
2437 objectV0->SetPhiTrue(dPhiV0Gen);
2438 objectV0->SetPDGCode(iPdgCodeMother);
2439 objectV0->SetPDGCodeMother(iPdgCodeMotherOfMother);
2444 // if (bIsCandidateK0s && bIsInPeakK0s) // selected candidates in peak
2445 if(bIsCandidateK0s) // selected candidates with any mass
2447 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2448 if(bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2450 // if (fbTreeOutput)
2451 // objectV0->SetOrigin(1);
2452 fh2V0K0sPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0K0s);
2453 Double_t valueEtaK[3] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen};
2454 fh3V0K0sEtaPtMassMCRec[iCentIndex]->Fill(valueEtaK);
2456 Double_t valueEtaDKNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2457 fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKNeg);
2458 Double_t valueEtaDKPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2459 fhnV0K0sInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKPos);
2461 fh2V0K0sMCResolMPt[iCentIndex]->Fill(dMassV0K0s - dMassPDGK0s, dPtV0);
2462 fh2V0K0sMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2463 if(bIsInConeJet) // true V0 associated to a candidate in jet
2465 Double_t valueKInJCMC[4] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2466 fh3V0K0sInJetPtMassMCRec[iCentIndex]->Fill(valueKInJCMC);
2467 Double_t valueEtaKIn[5] = {dMassV0K0s, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2468 fh4V0K0sInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaKIn);
2470 Double_t valueEtaDKJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2471 fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCNeg);
2472 Double_t valueEtaDKJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2473 fhnV0K0sInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDKJCPos);
2476 if(bV0MCIsK0s && !bV0MCIsPrimaryDist) // not primary K0s
2478 // if (fbTreeOutput)
2479 // objectV0->SetOrigin(-1);
2480 fh1V0K0sPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2484 // if (bIsCandidateLambda && bIsInPeakLambda) // selected candidates in peak
2485 if(bIsCandidateLambda) // selected candidates with any mass
2487 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2488 if(bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2490 // if (fbTreeOutput)
2491 // objectV0->SetOrigin(1);
2492 fh2V0LambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0Lambda);
2493 Double_t valueEtaL[3] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen};
2494 fh3V0LambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaL);
2496 Double_t valueEtaDLNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2497 fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLNeg);
2498 Double_t valueEtaDLPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2499 fhnV0LambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLPos);
2501 fh2V0LambdaMCResolMPt[iCentIndex]->Fill(dMassV0Lambda - dMassPDGLambda, dPtV0);
2502 fh2V0LambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2503 if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2505 Double_t valueLInJCMC[4] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2506 fh3V0LambdaInJetPtMassMCRec[iCentIndex]->Fill(valueLInJCMC);
2507 Double_t valueEtaLIn[5] = {dMassV0Lambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2508 fh4V0LambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaLIn);
2510 Double_t valueEtaDLJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2511 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCNeg);
2512 Double_t valueEtaDLJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2513 fhnV0LambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDLJCPos);
2516 // Fill the feed-down histograms
2517 if(bV0MCIsLambda && bV0MCComesFromXi)
2519 // if (fbTreeOutput)
2520 // objectV0->SetOrigin(2);
2521 Double_t valueFDLIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2522 fhnV0LambdaInclMCFD[iCentIndex]->Fill(valueFDLIncl);
2525 fhnV0LambdaBulkMCFD[iCentIndex]->Fill(valueFDLIncl);
2529 Double_t valueFDLInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2530 fhnV0LambdaInJetsMCFD[iCentIndex]->Fill(valueFDLInJets);
2533 if(bV0MCIsLambda && !bV0MCIsPrimaryDist && !bV0MCComesFromXi) // not primary Lambda
2535 // if (fbTreeOutput)
2536 // objectV0->SetOrigin(-1);
2537 fh1V0LambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2541 // if (bIsCandidateALambda && bIsInPeakALambda) // selected candidates in peak
2542 if(bIsCandidateALambda) // selected candidates with any mass
2544 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2545 if(bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2547 // if (fbTreeOutput)
2548 // objectV0->SetOrigin(1);
2549 fh2V0ALambdaPtMassMCRec[iCentIndex]->Fill(dPtV0Gen, dMassV0ALambda);
2550 Double_t valueEtaAL[3] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen};
2551 fh3V0ALambdaEtaPtMassMCRec[iCentIndex]->Fill(valueEtaAL);
2553 Double_t valueEtaDALNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2554 fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALNeg);
2555 Double_t valueEtaDALPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, 0};
2556 fhnV0ALambdaInclDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALPos);
2558 fh2V0ALambdaMCResolMPt[iCentIndex]->Fill(dMassV0ALambda - dMassPDGLambda, dPtV0);
2559 fh2V0ALambdaMCPtGenPtRec[iCentIndex]->Fill(dPtV0Gen, dPtV0);
2560 if(bIsInConeJet) // true V0 associated to a reconstructed candidate in jet
2562 Double_t valueALInJCMC[4] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt()};
2563 fh3V0ALambdaInJetPtMassMCRec[iCentIndex]->Fill(valueALInJCMC);
2564 Double_t valueEtaALIn[5] = {dMassV0ALambda, dPtV0Gen, dEtaV0Gen, jet->Pt(), dEtaV0Gen - jet->Eta()};
2565 fh4V0ALambdaInJetEtaPtMassMCRec[iCentIndex]->Fill(valueEtaALIn);
2567 Double_t valueEtaDALJCNeg[6] = {0, particleMCDaughterNeg->Eta(), particleMCDaughterNeg->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2568 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCNeg);
2569 Double_t valueEtaDALJCPos[6] = {1, particleMCDaughterPos->Eta(), particleMCDaughterPos->Pt(), dEtaV0Gen, dPtV0Gen, jet->Pt()};
2570 fhnV0ALambdaInJetsDaughterEtaPtPtMCRec[iCentIndex]->Fill(valueEtaDALJCPos);
2573 // Fill the feed-down histograms
2574 if(bV0MCIsALambda && bV0MCComesFromAXi)
2576 // if (fbTreeOutput)
2577 // objectV0->SetOrigin(2);
2578 Double_t valueFDALIncl[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), 0.};
2579 fhnV0ALambdaInclMCFD[iCentIndex]->Fill(valueFDALIncl);
2582 fhnV0ALambdaBulkMCFD[iCentIndex]->Fill(valueFDALIncl);
2586 Double_t valueFDALInJets[3] = {dPtV0Gen, particleMCMotherOfMother->Pt(), jet->Pt()};
2587 fhnV0ALambdaInJetsMCFD[iCentIndex]->Fill(valueFDALInJets);
2590 if(bV0MCIsALambda && !bV0MCIsPrimaryDist && !bV0MCComesFromAXi) // not primary anti-Lambda
2592 // if (fbTreeOutput)
2593 // objectV0->SetOrigin(-1);
2594 fh1V0ALambdaPtMCRecFalse[iCentIndex]->Fill(dPtV0Gen);
2598 //===== End Association of reconstructed V0 candidates with MC particles =====
2600 //===== End of V0 loop =====
2602 fh1V0CandPerEvent->Fill(iNV0CandTot);
2603 fh1V0CandPerEventCentK0s[iCentIndex]->Fill(iNV0CandK0s);
2604 fh1V0CandPerEventCentLambda[iCentIndex]->Fill(iNV0CandLambda);
2605 fh1V0CandPerEventCentALambda[iCentIndex]->Fill(iNV0CandALambda);
2607 if(fDebug > 2) printf("TaskV0sInJets: End of V0 loop\n");
2609 // Spectra of generated particles
2612 for(Int_t iPartMC = 0; iPartMC < iNTracksMC; iPartMC++)
2615 AliAODMCParticle* particleMC = (AliAODMCParticle*)arrayMC->At(iPartMC);
2619 // Get identity of MC particle
2620 Int_t iPdgCodeParticleMC = particleMC->GetPdgCode();
2621 // Fill Xi spectrum (3322 - Xi0, 3312 - Xi-)
2622 // if ( (iPdgCodeParticleMC==3322) || (iPdgCodeParticleMC==3312) )
2623 if((iPdgCodeParticleMC == 3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2625 // if (fbTreeOutput)
2626 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2627 fh1V0XiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2629 if((iPdgCodeParticleMC == -3312) && (TMath::Abs(particleMC->Y()) < 0.5))
2631 // if (fbTreeOutput)
2632 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2633 fh1V0AXiPtMCGen[iCentIndex]->Fill(particleMC->Pt());
2635 // Skip not interesting particles
2636 if((iPdgCodeParticleMC != iPdgCodeK0s) && (TMath::Abs(iPdgCodeParticleMC) != iPdgCodeLambda))
2639 // Check identity of the MC V0 particle
2640 // Is MC V0 particle K0S?
2641 Bool_t bV0MCIsK0s = (iPdgCodeParticleMC == iPdgCodeK0s);
2642 // Is MC V0 particle Lambda?
2643 Bool_t bV0MCIsLambda = (iPdgCodeParticleMC == +iPdgCodeLambda);
2644 // Is MC V0 particle anti-Lambda?
2645 Bool_t bV0MCIsALambda = (iPdgCodeParticleMC == -iPdgCodeLambda);
2647 Double_t dPtV0Gen = particleMC->Pt();
2648 Double_t dRapV0Gen = particleMC->Y();
2649 Double_t dEtaV0Gen = particleMC->Eta();
2654 if(bPrintCuts) printf("Gen: Applying cut: V0 |y|: < %f\n", dRapMax);
2655 if((TMath::Abs(dRapV0Gen) > dRapMax))
2658 // V0 pseudorapidity cut
2661 if(bPrintCuts) printf("Gen: Applying cut: V0 |eta|: < %f\n", dEtaMax);
2662 if((TMath::Abs(dEtaV0Gen) > dEtaMax))
2666 // Is MC V0 particle physical primary? Attention!! Definition of IsPhysicalPrimary may change!!
2667 Bool_t bV0MCIsPrimary = particleMC->IsPhysicalPrimary();
2669 // Get the MC mother particle of the MC V0 particle
2670 Int_t iIndexMotherOfMother = particleMC->GetMother();
2671 AliAODMCParticle* particleMCMotherOfMother = 0;
2672 if (iIndexMotherOfMother >= 0)
2673 particleMCMotherOfMother = (AliAODMCParticle*)arrayMC->At(iIndexMotherOfMother);
2674 // Get identity of the MC mother particle of the MC V0 particle if it exists
2675 Int_t iPdgCodeMotherOfMother = 0;
2676 if (particleMCMotherOfMother)
2677 iPdgCodeMotherOfMother = particleMCMotherOfMother->GetPdgCode();
2678 // Check if the MC mother particle is a physical primary Sigma
2679 Bool_t bV0MCComesFromSigma = kFALSE;
2680 if ((particleMCMotherOfMother && particleMCMotherOfMother->IsPhysicalPrimary()) && (TMath::Abs(iPdgCodeMotherOfMother)==3212) || (TMath::Abs(iPdgCodeMotherOfMother)==3224) || (TMath::Abs(iPdgCodeMotherOfMother)==3214) || (TMath::Abs(iPdgCodeMotherOfMother)==3114) )
2681 bV0MCComesFromSigma = kTRUE;
2682 // Should the MC V0 particle be considered as primary when it is Lambda?
2683 Bool_t bV0MCIsPrimaryLambda = (bV0MCIsPrimary || bV0MCComesFromSigma);
2685 // Reject non primary particles
2686 // if (!bV0MCIsPrimaryLambda)
2689 // Get the distance between the production point of the MC V0 particle and the primary vertex
2690 Double_t dx = dPrimVtxMCX - particleMC->Xv();
2691 Double_t dy = dPrimVtxMCY - particleMC->Yv();
2692 Double_t dz = dPrimVtxMCZ - particleMC->Zv();
2693 Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
2694 Bool_t bV0MCIsPrimaryDist = (dDistPrimary < dDistPrimaryMax); // Is close enough to be considered primary-like?
2696 // Check whether the MC V0 particle is in a MC jet
2697 AliAODJet* jetMC = 0;
2698 Bool_t bIsMCV0InJet = kFALSE;
2701 if(fDebug > 5) printf("TaskV0sInJets: Searching for gen V0 in %d MC jets\n", iNJetSel);
2702 for(Int_t iJet = 0; iJet < iNJetSel; iJet++)
2704 jetMC = (AliAODJet*)jetArraySel->At(iJet); // load a jet in the list
2705 if(fDebug > 5) printf("TaskV0sInJets: Checking if gen V0 in MC jet %d\n", iJet);
2706 if(IsParticleInCone(particleMC, jetMC, fdRadiusJet)) // If good jet in event, find out whether V0 is in that jet
2708 if(fDebug > 5) printf("TaskV0sInJets: gen V0 found in MC jet %d\n", iJet);
2709 bIsMCV0InJet = kTRUE;
2715 // Select only primary-like MC V0 particles
2717 // if (bV0MCIsK0s && bV0MCIsPrimary) // well reconstructed candidates
2718 if(bV0MCIsK0s && bV0MCIsPrimaryDist) // well reconstructed candidates
2720 // if (fbTreeOutput)
2721 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2722 fh1V0K0sPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2723 fh2V0K0sEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2726 fh2V0K0sInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2727 Double_t valueEtaKInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2728 fh3V0K0sInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaKInGen);
2732 // if (bV0MCIsLambda && bV0MCIsPrimaryLambda) // well reconstructed candidates
2733 if(bV0MCIsLambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2735 // if (fbTreeOutput)
2736 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2737 fh1V0LambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2738 fh2V0LambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2741 fh2V0LambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2742 Double_t valueEtaLInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2743 fh3V0LambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaLInGen);
2747 // if (bV0MCIsALambda && bV0MCIsPrimaryALambda) // well reconstructed candidates
2748 if(bV0MCIsALambda && bV0MCIsPrimaryDist) // well reconstructed candidates
2750 // if (fbTreeOutput)
2751 // new ((*fBranchV0Gen)[iNV0SelV0Gen++]) AliAODMCParticle(*((AliAODMCParticle*)particleMC));
2752 fh1V0ALambdaPtMCGen[iCentIndex]->Fill(dPtV0Gen);
2753 fh2V0ALambdaEtaPtMCGen[iCentIndex]->Fill(dPtV0Gen, dEtaV0Gen);
2756 fh2V0ALambdaInJetPtMCGen[iCentIndex]->Fill(dPtV0Gen, jetMC->Pt());
2757 Double_t valueEtaALInGen[4] = {dPtV0Gen, dEtaV0Gen, jetMC->Pt(), dEtaV0Gen - jetMC->Eta()};
2758 fh3V0ALambdaInJetEtaPtMCGen[iCentIndex]->Fill(valueEtaALInGen);
2764 // if (fbTreeOutput)
2765 // ftreeOut->Fill();
2767 jetArraySel->Delete();
2769 jetArrayPerp->Delete();
2770 delete jetArrayPerp;
2775 PostData(1, fOutputListStd);
2776 PostData(2, fOutputListQA);
2777 PostData(3, fOutputListCuts);
2778 PostData(4, fOutputListMC);
2779 // if (fbTreeOutput)
2780 // PostData(5,ftreeOut);
2781 // if(fDebug>5) printf("TaskV0sInJets: UserExec: End\n");
2784 void AliAnalysisTaskV0sInJets::FillQAHistogramV0(AliAODVertex* vtx, const AliAODv0* vZero, Int_t iIndexHisto, Bool_t IsCandK0s, Bool_t IsCandLambda, Bool_t IsInPeakK0s, Bool_t IsInPeakLambda)
2786 if(!IsCandK0s && !IsCandLambda)
2789 // Double_t dMassK0s = vZero->MassK0Short();
2790 // Double_t dMassLambda = vZero->MassLambda();
2792 fh1QAV0Status[iIndexHisto]->Fill(vZero->GetOnFlyStatus());
2794 AliAODTrack* trackNeg = (AliAODTrack*)vZero->GetDaughter(1); // negative track
2795 AliAODTrack* trackPos = (AliAODTrack*)vZero->GetDaughter(0); // positive track
2797 Short_t fTotalCharge = 0;
2798 for(Int_t i = 0; i < 2; i++)
2800 AliAODTrack* track = (AliAODTrack*)vZero->GetDaughter(i); // track
2802 fh1QAV0TPCRefit[iIndexHisto]->Fill(track->IsOn(AliAODTrack::kTPCrefit));
2803 Double_t nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
2804 fh1QAV0TPCRows[iIndexHisto]->Fill(nCrossedRowsTPC);
2805 Int_t findable = track->GetTPCNclsF();
2806 fh1QAV0TPCFindable[iIndexHisto]->Fill(findable);
2809 fh1QAV0TPCRowsFind[iIndexHisto]->Fill(nCrossedRowsTPC / findable);
2811 // Daughters: pseudo-rapidity cut
2812 fh1QAV0Eta[iIndexHisto]->Fill(track->Eta());
2813 if((nCrossedRowsTPC > (160. / (250. - 85.) * (255.*TMath::Abs(tan(track->Theta())) - 85.)) + 20.) && (track->Eta() < 0) && (track->Pt() > 0.15))
2816 fh2QAV0EtaRows[iIndexHisto]->Fill(track->Eta(), nCrossedRowsTPC);
2817 fh2QAV0PtRows[iIndexHisto]->Fill(track->Pt(), nCrossedRowsTPC);
2818 fh2QAV0PhiRows[iIndexHisto]->Fill(track->Phi(), nCrossedRowsTPC);
2819 fh2QAV0NClRows[iIndexHisto]->Fill(findable, nCrossedRowsTPC);
2820 fh2QAV0EtaNCl[iIndexHisto]->Fill(track->Eta(), findable);
2823 // Daughters: transverse momentum cut
2824 fh1QAV0Pt[iIndexHisto]->Fill(track->Pt());
2825 fTotalCharge += track->Charge();
2827 fh1QAV0Charge[iIndexHisto]->Fill(fTotalCharge);
2829 // Daughters: Impact parameter of daughters to prim vtx
2830 fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaNegToPrimVertex()));
2831 fh1QAV0DCAVtx[iIndexHisto]->Fill(TMath::Abs(vZero->DcaPosToPrimVertex()));
2832 // fh2CutDCAVtx[iIndexHisto]->Fill(dMassK0s,TMath::Abs(vZero->DcaNegToPrimVertex()));
2835 fh1QAV0DCAV0[iIndexHisto]->Fill(vZero->DcaV0Daughters());
2836 // fh2CutDCAV0[iIndexHisto]->Fill(dMassK0s,vZero->DcaV0Daughters());
2838 // V0: Cosine of the pointing angle
2839 fh1QAV0Cos[iIndexHisto]->Fill(vZero->CosPointingAngle(vtx));
2840 // fh2CutCos[iIndexHisto]->Fill(dMassK0s,vZero->CosPointingAngle(vtx));
2842 // V0: Fiducial volume
2844 vZero->GetSecondaryVtx(xyz);
2845 Double_t r2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
2846 fh1QAV0R[iIndexHisto]->Fill(TMath::Sqrt(r2));
2848 Double_t dAlpha = vZero->AlphaV0();
2849 Double_t dPtArm = vZero->PtArmV0();
2855 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2856 // fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2857 fh2QAV0EtaPtK0sPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2858 fh2QAV0PtPtK0sPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2859 fh2ArmPodK0s[iIndexHisto]->Fill(dAlpha, dPtArm);
2861 fh2QAV0EtaEtaK0s[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2862 fh2QAV0PhiPhiK0s[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2863 fh1QAV0RapK0s[iIndexHisto]->Fill(vZero->RapK0Short());
2870 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Eta(),vZero->Pt());
2871 // fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(trackPos->Eta(),vZero->Pt());
2872 fh2QAV0EtaPtLambdaPeak[iIndexHisto]->Fill(vZero->Eta(), vZero->Pt());
2873 fh2QAV0PtPtLambdaPeak[iIndexHisto]->Fill(trackNeg->Pt(), trackPos->Pt());
2874 fh2ArmPodLambda[iIndexHisto]->Fill(dAlpha, dPtArm);
2876 fh2QAV0EtaEtaLambda[iIndexHisto]->Fill(trackNeg->Eta(), trackPos->Eta());
2877 fh2QAV0PhiPhiLambda[iIndexHisto]->Fill(trackNeg->Phi(), trackPos->Phi());
2878 fh1QAV0RapLambda[iIndexHisto]->Fill(vZero->RapLambda());
2881 fh2ArmPod[iIndexHisto]->Fill(dAlpha, dPtArm);
2885 void AliAnalysisTaskV0sInJets::FillCandidates(Double_t mK, Double_t mL, Double_t mAL, Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut/*cut index*/, Int_t iCent/*cent index*/)
2889 fh1V0CounterCentK0s[iCent]->Fill(iCut);
2890 fh1V0InvMassK0sAll[iCut]->Fill(mK);
2894 fh1V0CounterCentLambda[iCent]->Fill(iCut);
2895 fh1V0InvMassLambdaAll[iCut]->Fill(mL);
2899 fh1V0CounterCentALambda[iCent]->Fill(iCut);
2900 fh1V0InvMassALambdaAll[iCut]->Fill(mAL);
2904 Bool_t AliAnalysisTaskV0sInJets::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
2906 // decides whether a particle is inside a jet cone
2907 if(!part1 || !part2)
2910 TVector3 vecMom2(part2->Px(), part2->Py(), part2->Pz());
2911 TVector3 vecMom1(part1->Px(), part1->Py(), part1->Pz());
2912 Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
2913 if(dR < dRMax) // momentum vectors of part1 and part2 are closer than dRMax
2918 Bool_t AliAnalysisTaskV0sInJets::OverlapWithJets(const TClonesArray* array, const AliVParticle* part, Double_t dDistance) const
2920 // decides whether a cone overlaps with other jets
2923 if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No part\n");
2928 if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: No array\n");
2931 Int_t iNJets = array->GetEntriesFast();
2934 if(fDebug > 2) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Warning: No jets\n");
2937 AliVParticle* jet = 0;
2938 for(Int_t iJet = 0; iJet < iNJets; iJet++)
2940 jet = (AliVParticle*)array->At(iJet);
2943 if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::OverlapWithJets: Error: Failed to load jet %d/%d\n", iJet, iNJets);
2946 if(IsParticleInCone(part, jet, dDistance))
2952 AliAODJet* AliAnalysisTaskV0sInJets::GetRandomCone(const TClonesArray* array, Double_t dEtaConeMax, Double_t dDistance) const
2954 // generate a random cone which does not overlap with selected jets
2955 // printf("Generating random cone...\n");
2956 TLorentzVector vecCone;
2957 AliAODJet* part = 0;
2958 Double_t dEta, dPhi;
2959 Int_t iNTrialsMax = 10;
2960 Bool_t bStatus = kFALSE;
2961 for(Int_t iTry = 0; iTry < iNTrialsMax; iTry++)
2963 // printf("Try %d\n",iTry);
2964 dEta = dEtaConeMax * (2 * fRandom->Rndm() - 1.); // random eta in [-dEtaConeMax,+dEtaConeMax]
2965 dPhi = TMath::TwoPi() * fRandom->Rndm(); // random phi in [0,2*Pi]
2966 vecCone.SetPtEtaPhiM(1., dEta, dPhi, 0.);
2967 part = new AliAODJet(vecCone);
2968 if(!OverlapWithJets(array, part, dDistance))
2971 // printf("Success\n");
2982 AliAODJet* AliAnalysisTaskV0sInJets::GetMedianCluster(const TClonesArray* array, Double_t dEtaConeMax) const
2984 // sort kt clusters by pT/area and return the middle one, based on code in AliAnalysisTaskJetChem
2987 if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Error: No array\n");
2990 Int_t iNClTot = array->GetEntriesFast(); // number of all clusters in the array
2991 Int_t iNCl = 0; // number of accepted clusters
2993 // get list of densities
2994 std::vector<std::vector<Double_t> > vecListClusters; // vector that contains pairs [ index, density ]
2995 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Loop over %d clusters.\n", iNClTot);
2996 for(Int_t ij = 0; ij < iNClTot; ij++)
2998 AliAODJet* clusterBg = (AliAODJet*)(array->At(ij));
3001 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: cluster %d/%d not ok\n", ij, iNClTot);
3004 if(TMath::Abs(clusterBg->Eta()) > 0.9 - fdRadiusJetBg)
3006 // fh2EtaPhiMedCone[0]->Fill(clusterBg->Eta(), clusterBg->Phi());
3007 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Cluster %d/%d used as accepted cluster %d.\n", ij, iNClTot, int(vecListClusters.size()));
3008 Double_t dPtBg = clusterBg->Pt();
3009 Double_t dAreaBg = clusterBg->EffectiveAreaCharged();
3010 Double_t dDensityBg = 0;
3012 dDensityBg = dPtBg / dAreaBg;
3013 std::vector<Double_t> vecCluster;
3014 vecCluster.push_back(ij);
3015 vecCluster.push_back(dDensityBg);
3016 vecListClusters.push_back(vecCluster);
3018 iNCl = vecListClusters.size();
3019 if(iNCl < 3) // need at least 3 clusters (skipping 2 highest)
3021 // if(fDebug > 2) printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Warning: Too little clusters\n");
3025 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Original lists:\n");
3026 // for(Int_t i = 0; i < iNCl; i++)
3027 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3029 // sort list of indeces by density in descending order
3030 std::sort(vecListClusters.begin(), vecListClusters.end(), CompareClusters);
3032 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Sorted lists:\n");
3033 // for(Int_t i = 0; i < iNCl; i++)
3034 // printf("%g %g\n", (vecListClusters[i])[0], (vecListClusters[i])[1]);
3036 // get median cluster with median density
3037 AliAODJet* clusterMed = 0;
3038 Int_t iIndex = 0; // index of the median cluster in the sorted list
3039 Int_t iIndexMed = 0; // index of the median cluster in the original array
3040 if(TMath::Odd(iNCl)) // odd number of clusters
3042 iIndex = (Int_t)(0.5 * (iNCl + 1)); // = (n - skip + 1)/2 + 1, skip = 2
3043 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Odd, median index = %d/%d\n", iIndex, iNCl);
3045 else // even number: picking randomly one of the two closest to median
3047 Int_t iIndex1 = (Int_t)(0.5 * iNCl); // = (n - skip)/2 + 1, skip = 2
3048 Int_t iIndex2 = (Int_t)(0.5 * iNCl + 1); // = (n - skip)/2 + 1 + 1, skip = 2
3049 iIndex = ((fRandom->Rndm() > 0.5) ? iIndex1 : iIndex2);
3050 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: Even, median index = %d or %d -> %d/%d\n", iIndex1, iIndex2, iIndex, iNCl);
3052 iIndexMed = Int_t((vecListClusters[iIndex])[0]);
3054 // printf("AliAnalysisTaskV0sInJets::GetMedianCluster: getting median cluster %d/%d ok, rho = %g\n", iIndexMed, iNClTot, (vecListClusters[iIndex])[1]);
3055 clusterMed = (AliAODJet*)(array->At(iIndexMed));
3057 if(TMath::Abs(clusterMed->Eta()) > dEtaConeMax)
3063 Double_t AliAnalysisTaskV0sInJets::AreaCircSegment(Double_t dRadius, Double_t dDistance) const
3065 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
3066 Double_t dEpsilon = 1e-2;
3067 Double_t dR = dRadius;
3068 Double_t dD = dDistance;
3069 if(TMath::Abs(dR) < dEpsilon)
3071 if(fDebug > 0) printf("AliAnalysisTaskV0sInJets::AreaCircSegment: Error: Too small radius: %f < %f\n", dR, dEpsilon);
3077 return TMath::Pi() * dR * dR;
3078 return dR * dR * TMath::ACos(dD / dR) - dD * TMath::Sqrt(dR * dR - dD * dD);
3081 Bool_t AliAnalysisTaskV0sInJets::IsSelectedForJets(AliAODEvent* fAOD, Double_t dVtxZCut, Double_t dVtxR2Cut, Double_t dCentCutLo, Double_t dCentCutUp, Bool_t bCutDeltaZ, Double_t dDeltaZMax)
3084 AliAODVertex* vertex = fAOD->GetPrimaryVertex();
3087 Int_t iNContribMin = 3;
3090 if(vertex->GetNContributors() < iNContribMin)
3092 TString vtxTitle(vertex->GetTitle());
3093 if(vtxTitle.Contains("TPCVertex"))
3095 Double_t zVertex = vertex->GetZ();
3096 if(TMath::Abs(zVertex) > dVtxZCut)
3100 AliAODVertex* vertexSPD = fAOD->GetPrimaryVertexSPD();
3103 // printf("IsSelectedForJets: Error: No SPD vertex\n");
3106 Double_t zVertexSPD = vertexSPD->GetZ();
3107 if(TMath::Abs(zVertex - zVertexSPD) > dDeltaZMax)
3109 // printf("IsSelectedForJets: Rejecting event due to delta z = %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
3112 // printf("IsSelectedForJets: Event OK: %f - %f = %f\n",zVertex,zVertexSPD,zVertex-zVertexSPD);
3114 Double_t xVertex = vertex->GetX();
3115 Double_t yVertex = vertex->GetY();
3116 Double_t radiusSq = yVertex * yVertex + xVertex * xVertex;
3117 if(radiusSq > dVtxR2Cut)
3119 Double_t centrality;
3120 // centrality = fAOD->GetHeader()->GetCentrality();
3121 centrality = ((AliVAODHeader*)fAOD->GetHeader())->GetCentralityP()->GetCentralityPercentile("V0M");
3126 if((dCentCutUp < 0) || (dCentCutLo < 0) || (dCentCutUp > 100) || (dCentCutLo > 100) || (dCentCutLo > dCentCutUp))
3128 if((centrality < dCentCutLo) || (centrality > dCentCutUp))
3133 if(centrality != -1)
3139 Int_t AliAnalysisTaskV0sInJets::GetCentralityBinIndex(Double_t centrality)
3141 // returns index of the centrality bin corresponding to the provided value of centrality
3142 if(centrality < 0 || centrality > fgkiCentBinRanges[fgkiNBinsCent - 1])
3144 for(Int_t i = 0; i < fgkiNBinsCent; i++)
3146 if(centrality <= fgkiCentBinRanges[i])
3152 Int_t AliAnalysisTaskV0sInJets::GetCentralityBinEdge(Int_t index)
3154 // returns the upper edge of the centrality bin corresponding to the provided value of index
3155 if(index < 0 || index >= fgkiNBinsCent)
3157 return fgkiCentBinRanges[index];
3160 TString AliAnalysisTaskV0sInJets::GetCentBinLabel(Int_t index)
3162 // get string with centrality range for given bin
3163 TString lowerEdge = ((index == 0) ? "0" : Form("%d", GetCentralityBinEdge(index - 1)));
3164 TString upperEdge = Form("%d", GetCentralityBinEdge(index));
3165 return Form("%s-%s %%", lowerEdge.Data(), upperEdge.Data());
3168 Double_t AliAnalysisTaskV0sInJets::MassPeakSigmaOld(Double_t pt, Int_t particle)
3170 // estimation of the sigma of the invariant-mass peak as a function of pT and particle type
3174 return 0.0044 + 0.0004 * (pt - 1.);
3177 return 0.0023 + 0.00034 * (pt - 1.);
3185 bool AliAnalysisTaskV0sInJets::CompareClusters(const std::vector<Double_t> cluster1, const std::vector<Double_t> cluster2)
3187 return (cluster1[1] > cluster2[1]);