]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/STRANGENESS/LambdaK0/AliAnalysisTaskExtractPerformanceV0.cxx
Added y shift to the histograms. Needed for the acquisiton of efficiencies.
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0 / AliAnalysisTaskExtractPerformanceV0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
17 //
18 // Modified version of AliAnalysisTaskCheckCascade.cxx.
19 // This is a 'hybrid' output version, in that it uses a classic TTree
20 // ROOT object to store the candidates, plus a couple of histograms filled on
21 // a per-event basis for storing variables too numerous to put in a tree. 
22 //
23 // --- Adapted to look for lambdas as well, using code from 
24 //        AliAnalysisTaskCheckPerformanceStrange.cxx
25 //
26 //  --- Algorithm Description 
27 //   1. Loop over primaries in stack to acquire generated charged Xi
28 //   2. Loop over stack to find V0s, fill TH3Fs "PrimRawPt"s for Efficiency
29 //   3. Perform Physics Selection
30 //   4. Perform Primary Vertex |z|<10cm selection
31 //   5. Perform Primary Vertex NoTPCOnly vertexing selection (>0 contrib.)
32 //   6. Perform Pileup Rejection
33 //   7. Analysis Loops: 
34 //    7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
35 //    7b. Fill TTree object with V0 information, candidates
36 //
37 //  Please Report Any Bugs! 
38 //
39 //   --- David Dobrigkeit Chinellato
40 //        (david.chinellato@gmail.com)
41 //
42 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
43
44 class TTree;
45 class TParticle;
46 class TVector3;
47
48 //class AliMCEventHandler;
49 //class AliMCEvent;
50 //class AliStack;
51
52 class AliESDVertex;
53 class AliAODVertex;
54 class AliESDv0;
55 class AliAODv0;
56
57 #include <Riostream.h>
58 #include "TList.h"
59 #include "TH1.h"
60 #include "TH2.h"
61 #include "TH3.h"
62 #include "TFile.h"
63 #include "THnSparse.h"
64 #include "TVector3.h"
65 #include "TCanvas.h"
66 #include "TMath.h"
67 #include "TLegend.h"
68 #include "AliLog.h"
69
70 #include "AliESDEvent.h"
71 #include "AliAODEvent.h"
72 #include "AliV0vertexer.h"
73 #include "AliCascadeVertexer.h"
74 #include "AliESDpid.h"
75 #include "AliESDtrack.h"
76 #include "AliESDtrackCuts.h"
77 #include "AliInputEventHandler.h"
78 #include "AliAnalysisManager.h"
79 #include "AliMCEventHandler.h"
80 #include "AliMCEvent.h"
81 #include "AliStack.h"
82
83 #include "AliCFContainer.h"
84 #include "AliMultiplicity.h"
85 #include "AliAODMCParticle.h"
86 #include "AliESDcascade.h"
87 #include "AliAODcascade.h"
88 #include "AliESDUtils.h"
89 #include "AliGenEventHeader.h"
90
91 #include "AliAnalysisTaskExtractPerformanceV0.h"
92
93 using std::cout;
94 using std::endl;
95
96 ClassImp(AliAnalysisTaskExtractPerformanceV0)
97
98 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0() 
99   : AliAnalysisTaskSE(), fListHistV0(0), fTree(0), fPIDResponse(0), fESDtrackCuts(0),
100    fkIsNuclear   ( kFALSE ), 
101    fkSwitchINT7  ( kFALSE ),
102    fkUseOnTheFly ( kFALSE ),
103    fkTakeAllTracks ( kFALSE ),
104    fpArapidityShift ( 0.465 ),
105 //------------------------------------------------
106 // Tree Variables 
107
108   fTreeVariablePrimaryStatus(0),
109   fTreeVariablePrimaryStatusMother(0),
110   fTreeVariableChi2V0(0),
111         fTreeVariableDcaV0Daughters(0),
112         fTreeVariableDcaV0ToPrimVertex(0),
113         fTreeVariableDcaPosToPrimVertex(0),
114         fTreeVariableDcaNegToPrimVertex(0),
115         fTreeVariableV0CosineOfPointingAngle(0),
116         fTreeVariableV0Radius(0),
117         fTreeVariablePt(0),
118         fTreeVariablePtMC(0),
119         fTreeVariableRapK0Short(0),
120         fTreeVariableRapLambda(0),
121         fTreeVariableRapMC(0),
122         fTreeVariableInvMassK0s(0),
123         fTreeVariableInvMassLambda(0),
124         fTreeVariableInvMassAntiLambda(0),
125         fTreeVariableAlphaV0(0),
126         fTreeVariablePtArmV0(0),
127         fTreeVariableNegTotMomentum(0),
128         fTreeVariablePosTotMomentum(0),
129         fTreeVariableNegTransvMomentum(0),
130         fTreeVariablePosTransvMomentum(0),
131         fTreeVariableNegTransvMomentumMC(0),
132         fTreeVariablePosTransvMomentumMC(0),
133    
134         fTreeVariableNSigmasPosProton(0),
135         fTreeVariableNSigmasPosPion(0),
136         fTreeVariableNSigmasNegProton(0),
137         fTreeVariableNSigmasNegPion(0),
138
139         fTreeVariablePtMother(0),
140         fTreeVariableV0CreationRadius(0),
141   fTreeVariablePID(0),
142   fTreeVariablePIDPositive(0),
143   fTreeVariablePIDNegative(0),
144   fTreeVariablePIDMother(0),
145   fTreeVariableIndexStatus(0),
146   fTreeVariableIndexStatusMother(0),
147
148   fTreeVariableRunNumber(0),
149   fTreeVariableEventNumber(0),
150
151         fTreeVariableDistOverTotMom(0),
152
153         fTreeVariablePosEta(0),
154         fTreeVariableNegEta(0),
155
156         fTreeVariableVertexZ(0),
157
158   fTreeVariableLeastNbrCrossedRows(0),
159   fTreeVariableLeastRatioCrossedRowsOverFindable(0),
160   fTreeVariableMultiplicity(0),
161   fTreeVariableMultiplicityMC(0),
162
163   fTreeVariableV0x(0),
164   fTreeVariableV0y(0),
165   fTreeVariableV0z(0),
166
167   fTreeVariableV0Px(0),
168   fTreeVariableV0Py(0),
169   fTreeVariableV0Pz(0),
170
171   fTreeVariableMCV0x(0),
172   fTreeVariableMCV0y(0),
173   fTreeVariableMCV0z(0),
174
175   fTreeVariableMCV0Px(0),
176   fTreeVariableMCV0Py(0),
177   fTreeVariableMCV0Pz(0),
178
179   fTreeVariablePVx(0),
180   fTreeVariablePVy(0),
181   fTreeVariablePVz(0),
182
183   fTreeVariableMCPVx(0),
184   fTreeVariableMCPVy(0),
185   fTreeVariableMCPVz(0),
186
187   fTreeVariableIsNonInjected(0),
188
189 //------------------------------------------------
190 // HISTOGRAMS
191 // --- Filled on an Event-by-event basis
192 //------------------------------------------------
193    fHistV0MultiplicityBeforeTrigSel(0),
194    fHistV0MultiplicityForTrigEvt(0), 
195    fHistV0MultiplicityForSelEvt(0),
196    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
197    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
198    fHistMultiplicityBeforeTrigSel(0),
199    fHistMultiplicityForTrigEvt(0),
200    fHistMultiplicity(0),
201    fHistMultiplicityNoTPCOnly(0),
202    fHistMultiplicityNoTPCOnlyNoPileup(0),
203
204         f2dHistMultiplicityVsTrueBeforeTrigSel(0),
205         f2dHistMultiplicityVsTrueForTrigEvt(0),
206         f2dHistMultiplicityVsTrue(0),
207         f2dHistMultiplicityVsTrueNoTPCOnly(0),
208         f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup(0),
209
210   //Raw Data for Vertex Z position estimator change
211         f2dHistMultiplicityVsVertexZBeforeTrigSel(0),
212         f2dHistMultiplicityVsVertexZForTrigEvt(0),
213         f2dHistMultiplicityVsVertexZ(0),
214         f2dHistMultiplicityVsVertexZNoTPCOnly(0),
215         f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup(0),
216
217   fHistGenVertexZBeforeTrigSel(0),     
218   fHistGenVertexZForTrigEvt(0),
219   fHistGenVertexZ(0),
220   fHistGenVertexZNoTPCOnly(0),
221   fHistGenVertexZNoTPCOnlyNoPileup(0),
222
223 //------------------------------------------------
224 // PARTICLE HISTOGRAMS
225 // --- Filled on a Particle-by-Particle basis
226 //------------------------------------------------
227    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
228    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
229    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
230    f3dHistPrimAnalysisPtVsYCMSVsMultLambda(0),
231    f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda(0),
232    f3dHistPrimAnalysisPtVsYCMSVsMultK0Short(0),
233    f3dHistPrimRawPtVsYVsMultLambda(0),
234    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
235    f3dHistPrimRawPtVsYVsMultK0Short(0),
236    f3dHistPrimRawPtVsYCMSVsMultLambda(0),
237    f3dHistPrimRawPtVsYCMSVsMultAntiLambda(0),
238    f3dHistPrimRawPtVsYCMSVsMultK0Short(0),
239    f3dHistPrimRawPtVsYVsMultNonInjLambda(0),
240    f3dHistPrimRawPtVsYVsMultNonInjAntiLambda(0),
241    f3dHistPrimRawPtVsYVsMultNonInjK0Short(0),
242    f3dHistPrimRawPtVsYVsMultMCLambda(0),
243    f3dHistPrimRawPtVsYVsMultMCAntiLambda(0),
244    f3dHistPrimRawPtVsYVsMultMCK0Short(0),
245    f3dHistPrimRawPtVsYVsVertexZLambda(0),
246    f3dHistPrimRawPtVsYVsVertexZAntiLambda(0),
247    f3dHistPrimRawPtVsYVsVertexZK0Short(0),
248    f3dHistPrimCloseToPVPtVsYVsMultLambda(0),
249    f3dHistPrimCloseToPVPtVsYVsMultAntiLambda(0),
250    f3dHistPrimCloseToPVPtVsYVsMultK0Short(0),
251    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
252    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
253    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
254    f3dHistGenPtVsYVsMultXiMinus(0),
255    f3dHistGenPtVsYVsMultXiPlus(0),
256    f3dHistGenPtVsYVsMultOmegaMinus(0),
257    f3dHistGenPtVsYVsMultOmegaPlus(0),
258    f3dHistGenSelectedPtVsYVsMultXiMinus(0),
259    f3dHistGenSelectedPtVsYVsMultXiPlus(0),
260    f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
261    f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
262    fHistPVx(0),
263    fHistPVy(0),
264    fHistPVz(0),
265    fHistPVxAnalysis(0),
266    fHistPVyAnalysis(0),
267    fHistPVzAnalysis(0),
268    fHistPVxAnalysisHasHighPtLambda(0),
269    fHistPVyAnalysisHasHighPtLambda(0),
270    fHistPVzAnalysisHasHighPtLambda(0),
271    fHistSwappedV0Counter(0)
272 {
273   // Dummy Constructor
274 }
275
276 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0(const char *name) 
277   : AliAnalysisTaskSE(name), fListHistV0(0), fTree(0), fPIDResponse(0), fESDtrackCuts(0),
278    fkIsNuclear   ( kFALSE ), 
279    fkSwitchINT7  ( kFALSE ),
280    fkUseOnTheFly ( kFALSE ),
281    fkTakeAllTracks ( kFALSE ),
282    fpArapidityShift ( 0.465 ),
283 //------------------------------------------------
284 // Tree Variables 
285
286   fTreeVariablePrimaryStatus(0),
287   fTreeVariablePrimaryStatusMother(0),
288   fTreeVariableChi2V0(0),
289         fTreeVariableDcaV0Daughters(0),
290         fTreeVariableDcaV0ToPrimVertex(0),
291         fTreeVariableDcaPosToPrimVertex(0),
292         fTreeVariableDcaNegToPrimVertex(0),
293         fTreeVariableV0CosineOfPointingAngle(0),
294         fTreeVariableV0Radius(0),
295         fTreeVariablePt(0),
296         fTreeVariablePtMC(0),
297         fTreeVariableRapK0Short(0),
298         fTreeVariableRapLambda(0),
299         fTreeVariableRapMC(0),
300         fTreeVariableInvMassK0s(0),
301         fTreeVariableInvMassLambda(0),
302         fTreeVariableInvMassAntiLambda(0),
303         fTreeVariableAlphaV0(0),
304         fTreeVariablePtArmV0(0),
305         fTreeVariableNegTotMomentum(0),
306         fTreeVariablePosTotMomentum(0),
307         fTreeVariableNegTransvMomentum(0),
308         fTreeVariablePosTransvMomentum(0),
309         fTreeVariableNegTransvMomentumMC(0),
310         fTreeVariablePosTransvMomentumMC(0),
311    
312         fTreeVariableNSigmasPosProton(0),
313         fTreeVariableNSigmasPosPion(0),
314         fTreeVariableNSigmasNegProton(0),
315         fTreeVariableNSigmasNegPion(0),
316
317         fTreeVariablePtMother(0),
318         fTreeVariableV0CreationRadius(0),
319   fTreeVariablePID(0),
320   fTreeVariablePIDPositive(0),
321   fTreeVariablePIDNegative(0),
322   fTreeVariablePIDMother(0),
323   fTreeVariableIndexStatus(0),
324   fTreeVariableIndexStatusMother(0),
325
326   fTreeVariableRunNumber(0),
327   fTreeVariableEventNumber(0),
328
329         fTreeVariableDistOverTotMom(0),
330
331         fTreeVariablePosEta(0),
332         fTreeVariableNegEta(0),
333
334         fTreeVariableVertexZ(0),
335
336   fTreeVariableLeastNbrCrossedRows(0),
337   fTreeVariableLeastRatioCrossedRowsOverFindable(0),
338   fTreeVariableMultiplicity(0),
339   fTreeVariableMultiplicityMC(0),
340
341   fTreeVariableV0x(0),
342   fTreeVariableV0y(0),
343   fTreeVariableV0z(0),
344
345   fTreeVariableV0Px(0),
346   fTreeVariableV0Py(0),
347   fTreeVariableV0Pz(0),
348
349   fTreeVariableMCV0x(0),
350   fTreeVariableMCV0y(0),
351   fTreeVariableMCV0z(0),
352
353   fTreeVariableMCV0Px(0),
354   fTreeVariableMCV0Py(0),
355   fTreeVariableMCV0Pz(0),
356
357   fTreeVariablePVx(0),
358   fTreeVariablePVy(0),
359   fTreeVariablePVz(0),
360
361   fTreeVariableMCPVx(0),
362   fTreeVariableMCPVy(0),
363   fTreeVariableMCPVz(0),
364
365   fTreeVariableIsNonInjected(0),
366
367
368 //------------------------------------------------
369 // HISTOGRAMS
370 // --- Filled on an Event-by-event basis
371 //------------------------------------------------
372    fHistV0MultiplicityBeforeTrigSel(0),
373    fHistV0MultiplicityForTrigEvt(0), 
374    fHistV0MultiplicityForSelEvt(0),
375    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
376    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
377    fHistMultiplicityBeforeTrigSel(0),
378    fHistMultiplicityForTrigEvt(0),
379    fHistMultiplicity(0),
380    fHistMultiplicityNoTPCOnly(0),
381    fHistMultiplicityNoTPCOnlyNoPileup(0),
382
383         f2dHistMultiplicityVsTrueBeforeTrigSel(0),
384         f2dHistMultiplicityVsTrueForTrigEvt(0),
385         f2dHistMultiplicityVsTrue(0),
386         f2dHistMultiplicityVsTrueNoTPCOnly(0),
387         f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup(0),
388
389   //Raw Data for Vertex Z position estimator change
390         f2dHistMultiplicityVsVertexZBeforeTrigSel(0),
391         f2dHistMultiplicityVsVertexZForTrigEvt(0),
392         f2dHistMultiplicityVsVertexZ(0),
393         f2dHistMultiplicityVsVertexZNoTPCOnly(0),
394         f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup(0),
395
396   fHistGenVertexZBeforeTrigSel(0),     
397   fHistGenVertexZForTrigEvt(0),
398   fHistGenVertexZ(0),
399   fHistGenVertexZNoTPCOnly(0),
400   fHistGenVertexZNoTPCOnlyNoPileup(0),
401
402 //------------------------------------------------
403 // PARTICLE HISTOGRAMS
404 // --- Filled on a Particle-by-Particle basis
405 //------------------------------------------------
406    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
407    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
408    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
409    f3dHistPrimRawPtVsYVsMultLambda(0),
410    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
411    f3dHistPrimRawPtVsYVsMultK0Short(0),
412    f3dHistPrimRawPtVsYVsMultNonInjLambda(0),
413    f3dHistPrimRawPtVsYVsMultNonInjAntiLambda(0),
414    f3dHistPrimRawPtVsYVsMultNonInjK0Short(0),
415    f3dHistPrimRawPtVsYVsMultMCLambda(0),
416    f3dHistPrimRawPtVsYVsMultMCAntiLambda(0),
417    f3dHistPrimRawPtVsYVsMultMCK0Short(0),
418    f3dHistPrimRawPtVsYVsVertexZLambda(0),
419    f3dHistPrimRawPtVsYVsVertexZAntiLambda(0),
420    f3dHistPrimRawPtVsYVsVertexZK0Short(0),
421    f3dHistPrimCloseToPVPtVsYVsMultLambda(0),
422    f3dHistPrimCloseToPVPtVsYVsMultAntiLambda(0),
423    f3dHistPrimCloseToPVPtVsYVsMultK0Short(0),
424    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
425    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
426    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
427    f3dHistGenPtVsYVsMultXiMinus(0),
428    f3dHistGenPtVsYVsMultXiPlus(0),
429    f3dHistGenPtVsYVsMultOmegaMinus(0),
430    f3dHistGenPtVsYVsMultOmegaPlus(0),
431    f3dHistGenSelectedPtVsYVsMultXiMinus(0),
432    f3dHistGenSelectedPtVsYVsMultXiPlus(0),
433    f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
434    f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
435    fHistPVx(0),
436    fHistPVy(0),
437    fHistPVz(0),
438    fHistPVxAnalysis(0),
439    fHistPVyAnalysis(0),
440    fHistPVzAnalysis(0),
441    fHistPVxAnalysisHasHighPtLambda(0),
442    fHistPVyAnalysisHasHighPtLambda(0),
443    fHistPVzAnalysisHasHighPtLambda(0),
444    fHistSwappedV0Counter(0)
445 {
446    // Constructor
447    // Output slot #0 writes into a TList container (Cascade)
448    DefineOutput(1, TList::Class());
449    DefineOutput(2, TTree::Class());
450 }
451
452
453 AliAnalysisTaskExtractPerformanceV0::~AliAnalysisTaskExtractPerformanceV0()
454 {
455 //------------------------------------------------
456 // DESTRUCTOR
457 //------------------------------------------------
458
459    if (fListHistV0){
460       delete fListHistV0;
461       fListHistV0 = 0x0;
462    }
463    if (fTree){
464       delete fTree;
465       fTree = 0x0;
466    }
467     //cleanup esd track cuts object too...
468    if (fESDtrackCuts){
469     delete fESDtrackCuts;
470     fESDtrackCuts = 0x0; 
471   }
472 }
473
474 //________________________________________________________________________
475 void AliAnalysisTaskExtractPerformanceV0::UserCreateOutputObjects()
476 {
477
478    OpenFile(2); 
479    // Called once
480
481 //------------------------------------------------
482
483    fTree = new TTree("fTree","V0Candidates");
484
485 //------------------------------------------------
486 // fTree Branch definitions - V0 Tree
487 //------------------------------------------------
488
489 //-----------BASIC-INFO---------------------------
490 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatus",&fTreeVariablePrimaryStatus,"fTreeVariablePrimaryStatus/I");        
491 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatusMother",&fTreeVariablePrimaryStatusMother,"fTreeVariablePrimaryStatusMother/I");      
492 /* 2*/   fTree->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"Chi2V0/F");
493 /* 3*/   fTree->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
494 /* 4*/   fTree->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
495 /* 5*/   fTree->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
496 /* 6*/   fTree->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
497 /* 7*/   fTree->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
498 /* 7*/   fTree->Branch("fTreeVariablePtMC",&fTreeVariablePtMC,"fTreeVariablePtMC/F");
499 /* 8*/   fTree->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
500 /* 9*/   fTree->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
501 /*10*/   fTree->Branch("fTreeVariableRapMC",&fTreeVariableRapMC,"fTreeVariableRapMC/F");
502 /*11*/   fTree->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
503 /*12*/   fTree->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
504 /*13*/   fTree->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
505 /*14*/   fTree->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
506 /*15*/   fTree->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
507 /*16*/   fTree->Branch("fTreeVariableNegTransvMomentum",&fTreeVariableNegTransvMomentum,"fTreeVariableNegTransvMomentum/F");
508 /*17*/   fTree->Branch("fTreeVariablePosTransvMomentum",&fTreeVariablePosTransvMomentum,"fTreeVariablePosTransvMomentum/F");
509 /*18*/   fTree->Branch("fTreeVariableNegTransvMomentumMC",&fTreeVariableNegTransvMomentumMC,"fTreeVariableNegTransvMomentumMC/F");
510 /*19*/   fTree->Branch("fTreeVariablePosTransvMomentumMC",&fTreeVariablePosTransvMomentumMC,"fTreeVariablePosTransvMomentumMC/F");
511 /*20*/   fTree->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
512 /*21*/   fTree->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
513 /*22*/   fTree->Branch("fTreeVariablePID",&fTreeVariablePID,"fTreeVariablePID/I");
514 /*23*/   fTree->Branch("fTreeVariablePIDPositive",&fTreeVariablePIDPositive,"fTreeVariablePIDPositive/I");
515 /*24*/   fTree->Branch("fTreeVariablePIDNegative",&fTreeVariablePIDNegative,"fTreeVariablePIDNegative/I");
516 /*25*/   fTree->Branch("fTreeVariablePIDMother",&fTreeVariablePIDMother,"fTreeVariablePIDMother/I");
517 /*26*/   fTree->Branch("fTreeVariablePtXiMother",&fTreeVariablePtMother,"fTreeVariablePtMother/F");
518 /*27*/   fTree->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
519 //-----------MULTIPLICITY-INFO--------------------
520 /*28*/   fTree->Branch("fTreeVariableMultiplicity",&fTreeVariableMultiplicity,"fTreeVariableMultiplicity/I");
521 /*28*/   fTree->Branch("fTreeVariableMultiplicityMC",&fTreeVariableMultiplicityMC,"fTreeVariableMultiplicityMC/I");
522 //------------------------------------------------
523 /*29*/   fTree->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
524 /*30*/   fTree->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
525 /*31*/   fTree->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
526 /*32*/   fTree->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
527 /*33*/   fTree->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
528 //------------------------------------------------
529 /*34*/   fTree->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
530 /*35*/   fTree->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
531 /*36*/   fTree->Branch("fTreeVariableV0CreationRadius",&fTreeVariableV0CreationRadius,"fTreeVariableV0CreationRadius/F");
532 /*37*/   fTree->Branch("fTreeVariableIndexStatus",&fTreeVariableIndexStatus,"fTreeVariableIndexStatus/I");
533 /*38*/   fTree->Branch("fTreeVariableIndexStatusMother",&fTreeVariableIndexStatusMother,"fTreeVariableIndexStatusMother/I");
534
535 /*39*/   fTree->Branch("fTreeVariableRunNumber",&fTreeVariableRunNumber,"fTreeVariableRunNumber/I");
536 /*40*/   fTree->Branch("fTreeVariableEventNumber",&fTreeVariableEventNumber,"fTreeVariableEventNumber/l");
537
538 /*34*/   fTree->Branch("fTreeVariableVertexZ",&fTreeVariableVertexZ,"fTreeVariableVertexZ/F");
539
540 //-----------FOR CTAU DEBUGGING: Full Phase Space + Decay Position Info 
541         fTree->Branch("fTreeVariableV0x",&fTreeVariableV0x,"fTreeVariableV0x/F");
542         fTree->Branch("fTreeVariableV0y",&fTreeVariableV0y,"fTreeVariableV0y/F");
543         fTree->Branch("fTreeVariableV0z",&fTreeVariableV0z,"fTreeVariableV0z/F");
544
545         fTree->Branch("fTreeVariableV0Px",&fTreeVariableV0Px,"fTreeVariableV0Px/F");
546         fTree->Branch("fTreeVariableV0Py",&fTreeVariableV0Py,"fTreeVariableV0Py/F");
547         fTree->Branch("fTreeVariableV0Pz",&fTreeVariableV0Pz,"fTreeVariableV0Pz/F");
548
549 //-----------FOR CTAU DEBUGGING: Full Phase Space + Decay Position Info, perfect info from MC
550         fTree->Branch("fTreeVariableMCV0x",&fTreeVariableMCV0x,"fTreeVariableMCV0x/F");
551         fTree->Branch("fTreeVariableMCV0y",&fTreeVariableMCV0y,"fTreeVariableMCV0y/F");
552         fTree->Branch("fTreeVariableMCV0z",&fTreeVariableMCV0z,"fTreeVariableMCV0z/F");
553
554         fTree->Branch("fTreeVariableMCV0Px",&fTreeVariableMCV0Px,"fTreeVariableMCV0Px/F");
555         fTree->Branch("fTreeVariableMCV0Py",&fTreeVariableMCV0Py,"fTreeVariableMCV0Py/F");
556         fTree->Branch("fTreeVariableMCV0Pz",&fTreeVariableMCV0Pz,"fTreeVariableMCV0Pz/F");
557
558 //-----------FOR CTAU DEBUGGING: Primary vertex info 
559         fTree->Branch("fTreeVariablePVx",&fTreeVariablePVx,"fTreeVariablePVx/F");
560         fTree->Branch("fTreeVariablePVy",&fTreeVariablePVy,"fTreeVariablePVy/F");
561         fTree->Branch("fTreeVariablePVz",&fTreeVariablePVz,"fTreeVariablePVz/F");
562
563         fTree->Branch("fTreeVariableMCPVx",&fTreeVariableMCPVx,"fTreeVariableMCPVx/F");
564         fTree->Branch("fTreeVariableMCPVy",&fTreeVariableMCPVy,"fTreeVariableMCPVy/F");
565         fTree->Branch("fTreeVariableMCPVz",&fTreeVariableMCPVz,"fTreeVariableMCPVz/F");
566
567         fTree->Branch("fTreeVariableIsNonInjected",&fTreeVariableIsNonInjected,"fTreeVariableIsNonInjected/O"); //O for bOOlean... 
568 //------------------------------------------------
569 // Particle Identification Setup
570 //------------------------------------------------
571
572    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
573    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
574    fPIDResponse = inputHandler->GetPIDResponse();
575
576   // Multiplicity 
577
578     if(! fESDtrackCuts ){
579           fESDtrackCuts = new AliESDtrackCuts();
580     }
581
582 //------------------------------------------------
583 // V0 Multiplicity Histograms
584 //------------------------------------------------
585
586    // Create histograms
587    OpenFile(1);
588    fListHistV0 = new TList();
589    fListHistV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
590
591
592    if(! fHistV0MultiplicityBeforeTrigSel) {
593       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", 
594          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", 
595          25, 0, 25);
596       fListHistV0->Add(fHistV0MultiplicityBeforeTrigSel);
597    }
598            
599    if(! fHistV0MultiplicityForTrigEvt) {
600       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", 
601          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", 
602          25, 0, 25);
603       fListHistV0->Add(fHistV0MultiplicityForTrigEvt);
604    }
605
606    if(! fHistV0MultiplicityForSelEvt) {
607       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", 
608          "V0s per event;Nbr of V0s/Evt;Events", 
609          25, 0, 25);
610       fListHistV0->Add(fHistV0MultiplicityForSelEvt);
611    }
612
613    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
614       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", 
615          "V0s per event;Nbr of V0s/Evt;Events", 
616          25, 0, 25);
617       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
618    }
619    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
620       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", 
621          "V0s per event;Nbr of V0s/Evt;Events", 
622          25, 0, 25);
623       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
624    }
625
626 //------------------------------------------------
627 // Track Multiplicity Histograms
628 //------------------------------------------------
629
630    if(! fHistMultiplicityBeforeTrigSel) {
631       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", 
632          "Tracks per event;Nbr of Tracks;Events", 
633          200, 0, 200);          
634       fListHistV0->Add(fHistMultiplicityBeforeTrigSel);
635    }
636    if(! fHistMultiplicityForTrigEvt) {
637       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", 
638          "Tracks per event;Nbr of Tracks;Events", 
639          200, 0, 200);          
640       fListHistV0->Add(fHistMultiplicityForTrigEvt);
641    }
642    if(! fHistMultiplicity) {
643       fHistMultiplicity = new TH1F("fHistMultiplicity", 
644          "Tracks per event;Nbr of Tracks;Events", 
645          200, 0, 200);          
646       fListHistV0->Add(fHistMultiplicity);
647    }
648    if(! fHistMultiplicityNoTPCOnly) {
649       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", 
650          "Tracks per event;Nbr of Tracks;Events", 
651          200, 0, 200);          
652       fListHistV0->Add(fHistMultiplicityNoTPCOnly);
653    }
654    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
655       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", 
656          "Tracks per event;Nbr of Tracks;Events", 
657          200, 0, 200);          
658       fListHistV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
659    }
660
661   //Raw Data for J/Psi paper Technique
662         //TH2F    *f2dHistMultiplicityVsTrueBeforeTrigSel;              //! multiplicity distribution    
663         //TH2F    *f2dHistMultiplicityVsTrueForTrigEvt;                         //! multiplicity distribution
664         //TH2F    *f2dHistMultiplicityVsTrue;                                                   //! multiplicity distribution
665         //TH2F    *f2dHistMultiplicityVsTrueNoTPCOnly;                          //! multiplicity distribution
666         //TH2F    *f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup;                  //! multiplicity distribution
667
668    if(! f2dHistMultiplicityVsTrueBeforeTrigSel) {
669       f2dHistMultiplicityVsTrueBeforeTrigSel = new TH2F("f2dHistMultiplicityVsTrueBeforeTrigSel", 
670          "Tracks per event", 200, 0, 200, 200, 0, 200);                 
671       fListHistV0->Add(f2dHistMultiplicityVsTrueBeforeTrigSel);
672    }
673    if(! f2dHistMultiplicityVsTrueForTrigEvt) {
674       f2dHistMultiplicityVsTrueForTrigEvt = new TH2F("f2dHistMultiplicityVsTrueForTrigEvt", 
675          "Tracks per event", 200, 0, 200, 200, 0, 200);                 
676       fListHistV0->Add(f2dHistMultiplicityVsTrueForTrigEvt);
677    }
678    if(! f2dHistMultiplicityVsTrue) {
679       f2dHistMultiplicityVsTrue = new TH2F("f2dHistMultiplicityVsTrue", 
680          "Tracks per event", 200, 0, 200, 200, 0, 200);                 
681       fListHistV0->Add(f2dHistMultiplicityVsTrue);
682    }
683    if(! f2dHistMultiplicityVsTrueNoTPCOnly) {
684       f2dHistMultiplicityVsTrueNoTPCOnly = new TH2F("f2dHistMultiplicityVsTrueNoTPCOnly", 
685          "Tracks per event", 200, 0, 200, 200, 0, 200);                 
686       fListHistV0->Add(f2dHistMultiplicityVsTrueNoTPCOnly);
687    }
688    if(! f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup) {
689       f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup = new TH2F("f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup", 
690          "Tracks per event", 200, 0, 200, 200, 0, 200);                 
691       fListHistV0->Add(f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup);
692    }
693
694
695   //Raw Data for Vertex Z position estimator change
696         //TH2F    *f2dHistMultiplicityVsVertexZBeforeTrigSel;           //! multiplicity distribution    
697         //TH2F    *f2dHistMultiplicityVsVertexZForTrigEvt;                      //! multiplicity distribution
698         //TH2F    *f2dHistMultiplicityVsVertexZ;                                                //! multiplicity distribution
699         //TH2F    *f2dHistMultiplicityVsVertexZNoTPCOnly;                               //! multiplicity distribution
700         //TH2F    *f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup;                       //! multiplicity distribution
701
702    if(! f2dHistMultiplicityVsVertexZBeforeTrigSel) {
703       f2dHistMultiplicityVsVertexZBeforeTrigSel = new TH2F("f2dHistMultiplicityVsVertexZBeforeTrigSel", 
704          "Tracks per event", 200, 0, 200,400, -20, 20);                 
705       fListHistV0->Add(f2dHistMultiplicityVsVertexZBeforeTrigSel);
706    }
707    if(! f2dHistMultiplicityVsVertexZForTrigEvt) {
708       f2dHistMultiplicityVsVertexZForTrigEvt = new TH2F("f2dHistMultiplicityVsVertexZForTrigEvt", 
709          "Tracks per event", 200, 0, 200, 400, -20, 20);                
710       fListHistV0->Add(f2dHistMultiplicityVsVertexZForTrigEvt);
711    }
712    if(! f2dHistMultiplicityVsVertexZ) {
713       f2dHistMultiplicityVsVertexZ = new TH2F("f2dHistMultiplicityVsVertexZ", 
714          "Tracks per event", 200, 0, 200, 400, -20, 20);                
715       fListHistV0->Add(f2dHistMultiplicityVsVertexZ);
716    }
717    if(! f2dHistMultiplicityVsVertexZNoTPCOnly) {
718       f2dHistMultiplicityVsVertexZNoTPCOnly = new TH2F("f2dHistMultiplicityVsVertexZNoTPCOnly", 
719          "Tracks per event", 200, 0, 200, 400, -20, 20);                
720       fListHistV0->Add(f2dHistMultiplicityVsVertexZNoTPCOnly);
721    }
722    if(! f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup) {
723       f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup = new TH2F("f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup", 
724          "Tracks per event", 200, 0, 200, 400, -20, 20);                
725       fListHistV0->Add(f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup);
726    }
727
728 //Generated PVz test
729
730 //   TH1F      *fHistGenVertexZBeforeTrigSel;     //! multiplicity distribution      
731 //   TH1F      *fHistGenVertexZForTrigEvt;        //! multiplicity distribution
732 //   TH1F      *fHistGenVertexZ;                  //! multiplicity distribution
733 //   TH1F      *fHistGenVertexZNoTPCOnly;         //! multiplicity distribution
734 //   TH1F      *fHistGenVertexZNoTPCOnlyNoPileup; //! multiplicity distribution
735
736    if(! fHistGenVertexZBeforeTrigSel) {
737          fHistGenVertexZBeforeTrigSel = new TH1F("fHistGenVertexZBeforeTrigSel", 
738             "PV z position;Nbr of Evts;z", 
739             400, -20, 20);       
740       fListHistV0->Add(fHistGenVertexZBeforeTrigSel);
741    }
742    if(! fHistGenVertexZForTrigEvt) {
743          fHistGenVertexZForTrigEvt = new TH1F("fHistGenVertexZForTrigEvt", 
744             "PV z position;Nbr of Evts;z", 
745             400, -20, 20);       
746       fListHistV0->Add(fHistGenVertexZForTrigEvt);
747    }
748    if(! fHistGenVertexZ) {
749          fHistGenVertexZ = new TH1F("fHistGenVertexZ", 
750             "PV z position;Nbr of Evts;z", 
751             400, -20, 20);       
752       fListHistV0->Add(fHistGenVertexZ);
753    }
754    if(! fHistGenVertexZNoTPCOnly) {
755          fHistGenVertexZNoTPCOnly = new TH1F("fHistGenVertexZNoTPCOnly", 
756             "PV z position;Nbr of Evts;z", 
757             400, -20, 20);       
758       fListHistV0->Add(fHistGenVertexZNoTPCOnly);
759    }
760    if(! fHistGenVertexZNoTPCOnlyNoPileup) {
761          fHistGenVertexZNoTPCOnlyNoPileup = new TH1F("fHistGenVertexZNoTPCOnlyNoPileup", 
762             "PV z position;Nbr of Evts;z", 
763             400, -20, 20);       
764       fListHistV0->Add(fHistGenVertexZNoTPCOnlyNoPileup);
765    }
766
767
768 //------------------------------------------------
769 // Generated Particle Histograms
770 //------------------------------------------------
771
772    Int_t lCustomNBins = 200; 
773    Double_t lCustomPtUpperLimit = 20; 
774    Int_t lCustomNBinsMultiplicity = 100;
775
776 //----------------------------------
777 // Raw Generated (Pre-physics-selection)
778 //----------------------------------
779
780 //--- 3D Histo (Pt, Y, Multiplicity)  
781
782    if(! f3dHistPrimRawPtVsYVsMultLambda) {
783       f3dHistPrimRawPtVsYVsMultLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
784       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultLambda);
785    }
786    if(! f3dHistPrimRawPtVsYVsMultAntiLambda) {
787       f3dHistPrimRawPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
788       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultAntiLambda);
789    }
790    if(! f3dHistPrimRawPtVsYVsMultK0Short) {
791       f3dHistPrimRawPtVsYVsMultK0Short = new TH3F( "f3dHistPrimRawPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
792       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultK0Short);
793    }
794
795    if(! f3dHistPrimRawPtVsYCMSVsMultLambda) {
796       f3dHistPrimRawPtVsYCMSVsMultLambda = new TH3F( "f3dHistPrimRawPtVsYCMSVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
797       fListHistV0->Add(f3dHistPrimRawPtVsYCMSVsMultLambda);
798    }
799    if(! f3dHistPrimRawPtVsYCMSVsMultAntiLambda) {
800       f3dHistPrimRawPtVsYCMSVsMultAntiLambda = new TH3F( "f3dHistPrimRawPtVsYCMSVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
801       fListHistV0->Add(f3dHistPrimRawPtVsYCMSVsMultAntiLambda);
802    }
803    if(! f3dHistPrimRawPtVsYCMSVsMultK0Short) {
804       f3dHistPrimRawPtVsYCMSVsMultK0Short = new TH3F( "f3dHistPrimRawPtVsYCMSVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
805       fListHistV0->Add(f3dHistPrimRawPtVsYCMSVsMultK0Short);
806    }
807
808 //---> Non-injected particles
809
810    if(! f3dHistPrimRawPtVsYVsMultNonInjLambda) {
811       f3dHistPrimRawPtVsYVsMultNonInjLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultNonInjLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
812       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultNonInjLambda);
813    }
814    if(! f3dHistPrimRawPtVsYVsMultNonInjAntiLambda) {
815       f3dHistPrimRawPtVsYVsMultNonInjAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultNonInjAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
816       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultNonInjAntiLambda);
817    }
818    if(! f3dHistPrimRawPtVsYVsMultNonInjK0Short) {
819       f3dHistPrimRawPtVsYVsMultNonInjK0Short = new TH3F( "f3dHistPrimRawPtVsYVsMultNonInjK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
820       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultNonInjK0Short);
821    }
822
823 //--- 3D Histo (Pt, Y, MultiplicityMC)  
824
825    if(! f3dHistPrimRawPtVsYVsMultMCLambda) {
826       f3dHistPrimRawPtVsYVsMultMCLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultMCLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
827       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCLambda);
828    }
829    if(! f3dHistPrimRawPtVsYVsMultMCAntiLambda) {
830       f3dHistPrimRawPtVsYVsMultMCAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultMCAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
831       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCAntiLambda);
832    }
833    if(! f3dHistPrimRawPtVsYVsMultMCK0Short) {
834       f3dHistPrimRawPtVsYVsMultMCK0Short = new TH3F( "f3dHistPrimRawPtVsYVsMultMCK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
835       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCK0Short);
836    }
837
838 //--- 3D Histo (Pt, Y, VertexZ)  
839
840    if(! f3dHistPrimRawPtVsYVsVertexZLambda) {
841       f3dHistPrimRawPtVsYVsVertexZLambda = new TH3F( "f3dHistPrimRawPtVsYVsVertexZLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs VertexZiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
842       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZLambda);
843    }
844    if(! f3dHistPrimRawPtVsYVsVertexZAntiLambda) {
845       f3dHistPrimRawPtVsYVsVertexZAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsVertexZAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs VertexZiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
846       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZAntiLambda);
847    }
848    if(! f3dHistPrimRawPtVsYVsVertexZK0Short) {
849       f3dHistPrimRawPtVsYVsVertexZK0Short = new TH3F( "f3dHistPrimRawPtVsYVsVertexZK0Short", "Pt_{K0S} Vs Y_{K0S} Vs VertexZiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
850       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZK0Short);
851    }
852
853 //--- 3D Histo (Pt, Y, Multiplicity), close to PV criterion
854
855    if(! f3dHistPrimCloseToPVPtVsYVsMultLambda) {
856       f3dHistPrimCloseToPVPtVsYVsMultLambda = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
857       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultLambda);
858    }
859    if(! f3dHistPrimCloseToPVPtVsYVsMultAntiLambda) {
860       f3dHistPrimCloseToPVPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
861       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultAntiLambda);
862    }
863    if(! f3dHistPrimCloseToPVPtVsYVsMultK0Short) {
864       f3dHistPrimCloseToPVPtVsYVsMultK0Short = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
865       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultK0Short);
866    }
867
868
869 //--- 3D Histo (Pt, Y, Proper Decay Length)
870
871    if(! f3dHistPrimRawPtVsYVsDecayLengthLambda) {
872       f3dHistPrimRawPtVsYVsDecayLengthLambda = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs DecayLength; Pt_{lambda} (GeV/c); Y_{#Lambda} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
873       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthLambda);
874    }
875    if(! f3dHistPrimRawPtVsYVsDecayLengthAntiLambda) {
876       f3dHistPrimRawPtVsYVsDecayLengthAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs DecayLength; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
877       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthAntiLambda);
878    }
879    if(! f3dHistPrimRawPtVsYVsDecayLengthK0Short) {
880       f3dHistPrimRawPtVsYVsDecayLengthK0Short = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthK0Short", "Pt_{K0S} Vs Y_{K0S} Vs DecayLength; Pt_{K0S} (GeV/c); Y_{K0S} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
881       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthK0Short);
882    }
883
884 //--------------------------------------------------------------------------------------
885 //--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, all generated
886
887    if(! f3dHistGenPtVsYVsMultXiMinus) {
888       f3dHistGenPtVsYVsMultXiMinus = new TH3F( "f3dHistGenPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
889       fListHistV0->Add(f3dHistGenPtVsYVsMultXiMinus);
890    }
891    if(! f3dHistGenPtVsYVsMultXiPlus) {
892       f3dHistGenPtVsYVsMultXiPlus = new TH3F( "f3dHistGenPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
893       fListHistV0->Add(f3dHistGenPtVsYVsMultXiPlus);
894    }
895 //--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
896
897    if(! f3dHistGenPtVsYVsMultOmegaMinus) {
898       f3dHistGenPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
899       fListHistV0->Add(f3dHistGenPtVsYVsMultOmegaMinus);
900    }
901    if(! f3dHistGenPtVsYVsMultOmegaPlus) {
902       f3dHistGenPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
903       fListHistV0->Add(f3dHistGenPtVsYVsMultOmegaPlus);
904    }
905
906 //--------------------------------------------------------------------------------------
907 //--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, at selected analysis evts
908
909    if(! f3dHistGenSelectedPtVsYVsMultXiMinus) {
910       f3dHistGenSelectedPtVsYVsMultXiMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
911       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultXiMinus);
912    }
913    if(! f3dHistGenSelectedPtVsYVsMultXiPlus) {
914       f3dHistGenSelectedPtVsYVsMultXiPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
915       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultXiPlus);
916    }
917 //--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
918
919    if(! f3dHistGenSelectedPtVsYVsMultOmegaMinus) {
920       f3dHistGenSelectedPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
921       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultOmegaMinus);
922    }
923    if(! f3dHistGenSelectedPtVsYVsMultOmegaPlus) {
924       f3dHistGenSelectedPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
925       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultOmegaPlus);
926    }
927
928
929 //----------------------------------
930 // Histos at analysis level 
931 //----------------------------------
932
933    if(! f3dHistPrimAnalysisPtVsYVsMultLambda) {
934       f3dHistPrimAnalysisPtVsYVsMultLambda = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
935       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultLambda);
936    }
937    if(! f3dHistPrimAnalysisPtVsYVsMultAntiLambda) {
938       f3dHistPrimAnalysisPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
939       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultAntiLambda);
940    }
941    if(! f3dHistPrimAnalysisPtVsYVsMultK0Short) {
942       f3dHistPrimAnalysisPtVsYVsMultK0Short = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
943       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultK0Short);
944    }
945
946    if(! f3dHistPrimAnalysisPtVsYCMSVsMultLambda) {
947       f3dHistPrimAnalysisPtVsYCMSVsMultLambda = new TH3F( "f3dHistPrimAnalysisPtVsYCMSVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
948       fListHistV0->Add(f3dHistPrimAnalysisPtVsYCMSVsMultLambda);
949    }
950    if(! f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda) {
951       f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda = new TH3F( "f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
952       fListHistV0->Add(f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda);
953    }
954    if(! f3dHistPrimAnalysisPtVsYCMSVsMultK0Short) {
955       f3dHistPrimAnalysisPtVsYCMSVsMultK0Short = new TH3F( "f3dHistPrimAnalysisPtVsYCMSVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
956       fListHistV0->Add(f3dHistPrimAnalysisPtVsYCMSVsMultK0Short);
957    }
958
959 //----------------------------------
960 // Primary Vertex Position Histos
961 //----------------------------------
962
963    if(! fHistPVx) {
964          fHistPVx = new TH1F("fHistPVx", 
965             "PV x position;Nbr of Evts;x", 
966             2000, -0.5, 0.5);       
967       fListHistV0->Add(fHistPVx);
968    }
969    if(! fHistPVy) {
970          fHistPVy = new TH1F("fHistPVy", 
971             "PV y position;Nbr of Evts;y", 
972             2000, -0.5, 0.5);       
973       fListHistV0->Add(fHistPVy);
974    }
975    if(! fHistPVz) {
976          fHistPVz = new TH1F("fHistPVz", 
977             "PV z position;Nbr of Evts;z", 
978             400, -20, 20);       
979       fListHistV0->Add(fHistPVz);
980    }
981
982    if(! fHistPVxAnalysis) {
983          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis", 
984             "PV x position;Nbr of Evts;x", 
985             2000, -0.5, 0.5);       
986       fListHistV0->Add(fHistPVxAnalysis);
987    }
988    if(! fHistPVyAnalysis) {
989          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis", 
990             "PV y position;Nbr of Evts;y", 
991             2000, -0.5, 0.5);       
992       fListHistV0->Add(fHistPVyAnalysis);
993    }
994    if(! fHistPVzAnalysis) {
995          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis", 
996             "PV z position;Nbr of Evts;z", 
997             400, -20, 20);       
998       fListHistV0->Add(fHistPVzAnalysis);
999    }
1000
1001    if(! fHistPVxAnalysisHasHighPtLambda) {
1002          fHistPVxAnalysisHasHighPtLambda = new TH1F("fHistPVxAnalysisHasHighPtLambda", 
1003             "PV x position;Nbr of Evts;x", 
1004             2000, -0.5, 0.5);       
1005       fListHistV0->Add(fHistPVxAnalysisHasHighPtLambda);
1006    }
1007    if(! fHistPVyAnalysisHasHighPtLambda) {
1008          fHistPVyAnalysisHasHighPtLambda = new TH1F("fHistPVyAnalysisHasHighPtLambda", 
1009             "PV y position;Nbr of Evts;y", 
1010             2000, -0.5, 0.5);       
1011       fListHistV0->Add(fHistPVyAnalysisHasHighPtLambda);
1012    }
1013    if(! fHistPVzAnalysisHasHighPtLambda) {
1014          fHistPVzAnalysisHasHighPtLambda = new TH1F("fHistPVzAnalysisHasHighPtLambda", 
1015             "PV z position;Nbr of Evts;z", 
1016             400, -20, 20);       
1017       fListHistV0->Add(fHistPVzAnalysisHasHighPtLambda);
1018    }
1019    if(! fHistSwappedV0Counter) {
1020       fHistSwappedV0Counter = new TH1F("fHistSwappedV0Counter", 
1021          "Swap or not histo;Swapped (1) or not (0); count", 
1022          2, 0, 2);              
1023       fListHistV0->Add(fHistSwappedV0Counter);
1024    }
1025
1026    //List of Histograms: Normal
1027    PostData(1, fListHistV0);
1028
1029    //TTree Object: Saved to base directory. Should cache to disk while saving. 
1030    //(Important to avoid excessive memory usage, particularly when merging)
1031    PostData(2, fTree);
1032
1033 }// end UserCreateOutputObjects
1034
1035
1036 //________________________________________________________________________
1037 void AliAnalysisTaskExtractPerformanceV0::UserExec(Option_t *) 
1038 {
1039   // Main loop
1040   // Called for each event
1041
1042    AliESDEvent *lESDevent = 0x0;
1043    AliMCEvent  *lMCevent  = 0x0; 
1044    AliStack    *lMCstack  = 0x0; 
1045
1046    Int_t    lNumberOfV0s                      = -1;
1047    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
1048    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
1049    Double_t lMagneticField                 = -10.;
1050    
1051   // Connect to the InputEvent   
1052   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
1053
1054    // Appropriate for ESD analysis! 
1055       
1056    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1057    if (!lESDevent) {
1058       AliWarning("ERROR: lESDevent not available \n");
1059       return;
1060    }
1061         
1062    fTreeVariableRunNumber = lESDevent->GetRunNumber();
1063    fTreeVariableEventNumber =  
1064     ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
1065       ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
1066         ((ULong64_t)lESDevent->GetBunchCrossNumber() )  );
1067
1068    lMCevent = MCEvent();
1069    if (!lMCevent) {
1070       Printf("ERROR: Could not retrieve MC event \n");
1071       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;   
1072       return;
1073    }
1074
1075    lMCstack = lMCevent->Stack();
1076    if (!lMCstack) {
1077       Printf("ERROR: Could not retrieve MC stack \n");
1078       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1079       return;
1080    }
1081    TArrayF mcPrimaryVtx;
1082    AliGenEventHeader* mcHeader=lMCevent->GenEventHeader();
1083    if(!mcHeader) return;
1084    mcHeader->PrimaryVertex(mcPrimaryVtx);
1085         
1086 //------------------------------------------------
1087 // Multiplicity Information Acquistion
1088 //------------------------------------------------
1089
1090    //REVISED multiplicity estimator after 'multiplicity day' (2011)
1091    Int_t lMultiplicity = -100; 
1092
1093    //testing purposes
1094    if(fkIsNuclear == kFALSE) lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1095
1096    //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
1097    //---> Warning: Experimental
1098    if(fkIsNuclear == kTRUE){ 
1099       AliCentrality* centrality;
1100       centrality = lESDevent->GetCentrality();
1101       lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
1102       if (centrality->GetQuality()>1) {
1103         PostData(1, fListHistV0);
1104         PostData(2, fTree);
1105         return;
1106       }
1107    }
1108   
1109    //Set variable for filling tree afterwards!
1110    //---> pp case......: GetReferenceMultiplicity
1111    //---> Pb-Pb case...: Centrality by V0M
1112    fTreeVariableMultiplicity = lMultiplicity;
1113
1114    fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
1115    fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
1116         
1117 //------------------------------------------------
1118 // MC Information Acquistion
1119 //------------------------------------------------
1120
1121    Int_t iNumberOfPrimaries = -1;
1122    iNumberOfPrimaries = lMCstack->GetNprimary();
1123    if(iNumberOfPrimaries < 1) return; 
1124    Bool_t lHasHighPtLambda = kFALSE;
1125
1126 //------------------------------------------------
1127 // Variable Definition
1128 //------------------------------------------------
1129
1130    Int_t lNbMCPrimary        = 0;
1131
1132    Int_t lPdgcodeCurrentPart = 0;
1133    Double_t lRapCurrentPart  = 0;
1134    Double_t lPtCurrentPart   = 0;
1135   
1136    //Int_t lComeFromSigma      = 0;
1137
1138    // current mc particle 's mother
1139    //Int_t iCurrentMother  = 0;
1140    lNbMCPrimary = lMCstack->GetNprimary();
1141
1142 //------------------------------------------------
1143 // Pre-Physics Selection
1144 //------------------------------------------------
1145
1146 //----- Loop on primary Xi, Omega --------------------------------------------------------------
1147    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) 
1148    {// This is the begining of the loop on primaries
1149       
1150       TParticle* lCurrentParticlePrimary = 0x0; 
1151       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
1152       if(!lCurrentParticlePrimary){
1153          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1154          continue;
1155       }
1156       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) { 
1157          Double_t lRapXiMCPrimary = -100;
1158          if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) { 
1159            if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
1160              lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
1161            }
1162          }
1163
1164          //=================================================================================
1165          // Xi Histograms
1166          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){ 
1167             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1168             f3dHistGenPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1169          }
1170          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){ 
1171             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1172             f3dHistGenPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1173          }
1174          // Omega Histograms
1175          if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){ 
1176             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1177             f3dHistGenPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1178          }
1179          if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){ 
1180             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1181             f3dHistGenPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1182          }
1183       } 
1184    }
1185 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
1186
1187 //--------- GENERATED NUMBER OF CHARGED PARTICLES
1188 // ---> Set Variables to Zero again
1189 // ---> Variable Definition
1190
1191   Long_t lNumberOfCharged = 0; 
1192
1193 //----- Loop on Stack ----------------------------------------------------------------
1194    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++) 
1195    {// This is the begining of the loop on tracks
1196       TParticle* particleOne = lMCstack->Particle(iCurrentLabelStack);
1197       if(!particleOne) continue;
1198       if(!particleOne->GetPDG()) continue;
1199       Double_t lThisCharge = particleOne->GetPDG()->Charge()/3.;
1200       if(TMath::Abs(lThisCharge)<0.001) continue;
1201       if(! (lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) ) continue;
1202       
1203       Double_t gpt = particleOne -> Pt();
1204       Double_t geta = particleOne -> Eta(); 
1205
1206       if( TMath::Abs(geta) < 0.5) lNumberOfCharged++; 
1207    }//End of loop on tracks
1208 //----- End Loop on Stack ------------------------------------------------------------
1209
1210    //Double_t lpArapidityShift = 0.465;
1211    Bool_t lStackNatural = kTRUE;
1212 //----- Loop on Lambda, K0Short ----------------------------------------------------------------
1213    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++) 
1214    {// This is the begining of the loop on tracks
1215       
1216       TParticle* lCurrentParticleForLambdaCheck = 0x0; 
1217       lCurrentParticleForLambdaCheck = lMCstack->Particle( iCurrentLabelStack );
1218       if(!lCurrentParticleForLambdaCheck){
1219          Printf("V0s loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1220          continue;
1221       }
1222
1223       //=================================================================================
1224       //Single-Strange checks
1225       // Keep only K0s, Lambda and AntiLambda:
1226       lPdgcodeCurrentPart = lCurrentParticleForLambdaCheck->GetPdgCode();             
1227
1228       if ( (lCurrentParticleForLambdaCheck->GetPdgCode() == 310   ) ||
1229            (lCurrentParticleForLambdaCheck->GetPdgCode() == 3122  ) ||
1230            (lCurrentParticleForLambdaCheck->GetPdgCode() == -3122 ) )
1231            {
1232          lRapCurrentPart   = MyRapidity(lCurrentParticleForLambdaCheck->Energy(),lCurrentParticleForLambdaCheck->Pz());
1233          lPtCurrentPart    = lCurrentParticleForLambdaCheck->Pt();
1234
1235           //Use Close to PV for filling CloseToPV histograms!
1236          Double_t dx, dy, dz; 
1237
1238          dx = ( (mcPrimaryVtx.At(0)) - (lCurrentParticleForLambdaCheck->Vx()) ); 
1239          dy = ( (mcPrimaryVtx.At(1)) - (lCurrentParticleForLambdaCheck->Vy()) );
1240          dz = ( (mcPrimaryVtx.At(2)) - (lCurrentParticleForLambdaCheck->Vz()) );
1241          Double_t lDistToPV = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
1242          if( lDistToPV <= 0.001){ 
1243            if( lPdgcodeCurrentPart == 3122 ){
1244               f3dHistPrimCloseToPVPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1245            }
1246            if( lPdgcodeCurrentPart == -3122 ){
1247               f3dHistPrimCloseToPVPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1248            }
1249            if( lPdgcodeCurrentPart == 310 ){
1250               f3dHistPrimCloseToPVPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1251            }
1252          }
1253
1254          //Use Physical Primaries only for filling PrimRaw Histograms!
1255          if ( lMCstack->IsPhysicalPrimary(iCurrentLabelStack)!=kTRUE ) continue;
1256
1257           lStackNatural = lMCevent->IsFromBGEvent(iCurrentLabelStack); //Is it? 
1258           if (!lStackNatural){
1259             if (!(lCurrentParticleForLambdaCheck->GetFirstMother()<0)) 
1260               {lStackNatural = kTRUE;} // because there are primaries (ALICE definition) not produced in the collision
1261           }
1262
1263          if( lPdgcodeCurrentPart == 3122 ){
1264             f3dHistPrimRawPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1265             f3dHistPrimRawPtVsYCMSVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1266             if(lStackNatural){f3dHistPrimRawPtVsYVsMultNonInjLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);}
1267             f3dHistPrimRawPtVsYVsMultMCLambda->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1268             f3dHistPrimRawPtVsYVsVertexZLambda->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1269             if( TMath::Abs( lCurrentParticleForLambdaCheck->Eta() )<1.2 && lPtCurrentPart>2 ){
1270                lHasHighPtLambda = kTRUE; //Keep track of events with Lambda within |eta|<1.2 and pt>2
1271             }
1272          }
1273          if( lPdgcodeCurrentPart == -3122 ){
1274             f3dHistPrimRawPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1275             f3dHistPrimRawPtVsYCMSVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1276             if(lStackNatural){f3dHistPrimRawPtVsYVsMultNonInjAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);}
1277             f3dHistPrimRawPtVsYVsMultMCAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1278             f3dHistPrimRawPtVsYVsVertexZAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1279          }
1280          if( lPdgcodeCurrentPart == 310 ){
1281             f3dHistPrimRawPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1282             f3dHistPrimRawPtVsYCMSVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1283             if(lStackNatural){f3dHistPrimRawPtVsYVsMultNonInjK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);}
1284             f3dHistPrimRawPtVsYVsMultMCK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1285             f3dHistPrimRawPtVsYVsVertexZK0Short->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1286          }
1287          //Decay Length Acquisition=====================================================
1288          Double_t decaylength = -1; 
1289          Double_t lV0Mass = -1; 
1290           
1291          if( !(lCurrentParticleForLambdaCheck->GetDaughter(0) < 0) ) {
1292             TParticle* lDght0ofV0 = lMCstack->Particle(  lCurrentParticleForLambdaCheck->GetDaughter(0) ); //get first daughter
1293             if(lDght0ofV0){ // skip if not defined. 
1294                decaylength = TMath::Sqrt(
1295                                         TMath::Power( lCurrentParticleForLambdaCheck->Vx() - lDght0ofV0->Vx() , 2) + 
1296                                         TMath::Power( lCurrentParticleForLambdaCheck->Vy() - lDght0ofV0->Vy() , 2) + 
1297                                         TMath::Power( lCurrentParticleForLambdaCheck->Vz() - lDght0ofV0->Vz() , 2)
1298                );
1299                //Need to correct for relativitity! Involves multiplying by mass and dividing by momentum. 
1300                if(TMath::Abs( lPdgcodeCurrentPart ) == 3122 ) { lV0Mass = 1.115683; }
1301                if(TMath::Abs( lPdgcodeCurrentPart ) == 310 ) { lV0Mass = 0.497614; }
1302                if( lCurrentParticleForLambdaCheck->P() + 1e-10 != 0 ) decaylength = ( lV0Mass * decaylength ) / ( lCurrentParticleForLambdaCheck->P() + 1e-10 );
1303                if( lCurrentParticleForLambdaCheck->P() + 1e-10 == 0 ) decaylength = 1e+5;
1304             }
1305          }
1306          if( lPdgcodeCurrentPart == 3122) f3dHistPrimRawPtVsYVsDecayLengthLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
1307          if( lPdgcodeCurrentPart == -3122) f3dHistPrimRawPtVsYVsDecayLengthAntiLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
1308          if( lPdgcodeCurrentPart == 310) f3dHistPrimRawPtVsYVsDecayLengthK0Short ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
1309       }
1310    }//End of loop on tracks
1311 //----- End Loop on Lambda, K0Short ------------------------------------------------------------
1312
1313
1314    f2dHistMultiplicityVsTrueBeforeTrigSel->Fill ( lMultiplicity , lNumberOfCharged );
1315
1316     fTreeVariableMultiplicityMC = lNumberOfCharged;
1317
1318    fHistGenVertexZBeforeTrigSel->Fill( (mcPrimaryVtx.At(2)) );
1319
1320    lPdgcodeCurrentPart = 0;
1321    lRapCurrentPart  = 0;
1322    lPtCurrentPart   = 0;
1323
1324 //------------------------------------------------
1325 // Physics Selection
1326 //------------------------------------------------
1327
1328    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1329    Bool_t isSelected = 0;
1330    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1331
1332    //pA triggering: CINT7
1333    if( fkSwitchINT7 ) isSelected = (maskIsSelected & AliVEvent::kINT7) == AliVEvent::kINT7;
1334
1335    //Standard Min-Bias Selection
1336    if ( ! isSelected ) { 
1337       PostData(1, fListHistV0);
1338       PostData(2, fTree);
1339       return;
1340    }
1341
1342    f2dHistMultiplicityVsTrueForTrigEvt->Fill ( lMultiplicity , lNumberOfCharged );
1343    fHistGenVertexZForTrigEvt->Fill( mcPrimaryVtx.At(2) );
1344 //------------------------------------------------
1345 // After Trigger Selection
1346 //------------------------------------------------
1347
1348    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
1349   
1350    //Set variable for filling tree afterwards!
1351    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
1352    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
1353
1354 //------------------------------------------------
1355 // Getting: Primary Vertex + MagField Info
1356 //------------------------------------------------
1357
1358    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
1359    // get the vtx stored in ESD found with tracks
1360    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
1361         
1362    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();      
1363    // get the best primary vertex available for the event
1364    // As done in AliCascadeVertexer, we keep the one which is the best one available.
1365    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1366    // This one will be used for next calculations (DCA essentially)
1367    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
1368
1369    Double_t lPrimaryVtxPosition[3];
1370    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
1371    lPrimaryVtxPosition[0] = primaryVtx->GetX();
1372    lPrimaryVtxPosition[1] = primaryVtx->GetY();
1373    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
1374    fHistPVx->Fill( lPrimaryVtxPosition[0] );
1375    fHistPVy->Fill( lPrimaryVtxPosition[1] );
1376    fHistPVz->Fill( lPrimaryVtxPosition[2] );
1377
1378    f2dHistMultiplicityVsVertexZForTrigEvt->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1379
1380 //------------------------------------------------
1381 // Primary Vertex Z position: SKIP
1382 //------------------------------------------------
1383
1384    if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) { 
1385       AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !"); 
1386       PostData(1, fListHistV0);
1387       PostData(2, fTree);
1388       return; 
1389    }
1390
1391    f2dHistMultiplicityVsVertexZ->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1392
1393    lMagneticField = lESDevent->GetMagneticField( );
1394    fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
1395    fHistMultiplicity->Fill(lMultiplicity);
1396    f2dHistMultiplicityVsTrue->Fill ( lMultiplicity , lNumberOfCharged );
1397    fHistGenVertexZ->Fill( (mcPrimaryVtx.At(2)) );
1398 //------------------------------------------------
1399 // SKIP: Events with well-established PVtx
1400 //------------------------------------------------
1401         
1402    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
1403    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
1404    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
1405       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1406       PostData(1, fListHistV0);
1407       PostData(2, fTree);
1408       return;
1409    }
1410
1411    f2dHistMultiplicityVsVertexZNoTPCOnly->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1412    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
1413    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
1414    f2dHistMultiplicityVsTrueNoTPCOnly->Fill ( lMultiplicity , lNumberOfCharged );
1415    fHistGenVertexZNoTPCOnly->Fill( (mcPrimaryVtx.At(2)) );
1416 //------------------------------------------------
1417 // Pileup Rejection Studies
1418 //------------------------------------------------
1419
1420    // FIXME : quality selection regarding pile-up rejection 
1421    if(lESDevent->IsPileupFromSPD() && !fkIsNuclear ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1422       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
1423       PostData(1, fListHistV0);
1424       PostData(2, fTree);
1425       return;
1426    }
1427    f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1428    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
1429    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
1430    f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup->Fill ( lMultiplicity , lNumberOfCharged );
1431    fHistGenVertexZNoTPCOnlyNoPileup->Fill( (mcPrimaryVtx.At(2)) );
1432    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
1433    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){ 
1434       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
1435       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
1436       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
1437       if ( lHasHighPtLambda == kTRUE ){ 
1438          fHistPVxAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[0] );
1439          fHistPVyAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[1] );
1440          fHistPVzAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[2] );
1441       }
1442    }
1443
1444   fTreeVariableVertexZ = lPrimaryVtxPosition[2];
1445
1446   fTreeVariablePVx = lPrimaryVtxPosition[0];
1447   fTreeVariablePVy = lPrimaryVtxPosition[1];
1448   fTreeVariablePVz = lPrimaryVtxPosition[2];
1449
1450   fTreeVariableMCPVx = (mcPrimaryVtx.At(0));
1451   fTreeVariableMCPVy = (mcPrimaryVtx.At(1));
1452   fTreeVariableMCPVz = (mcPrimaryVtx.At(2));
1453
1454 //------------------------------------------------
1455 // stack loop starts here
1456 //------------------------------------------------
1457
1458 //---> Loop over ALL PARTICLES
1459  
1460    for (Int_t iMc = 0; iMc < (lMCstack->GetNtrack()); iMc++) {  
1461       TParticle *p0 = lMCstack->Particle(iMc); 
1462       if (!p0) {
1463          //Printf("ERROR: particle with label %d not found in lMCstack (mc loop)", iMc);
1464          continue;
1465       }
1466       lPdgcodeCurrentPart = p0->GetPdgCode();
1467
1468       // Keep only K0s, Lambda and AntiLambda:
1469       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
1470         
1471       lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
1472       lPtCurrentPart    = p0->Pt();
1473
1474         //Use Physical Primaries only for filling PrimRaw Histograms!
1475       if ( lMCstack->IsPhysicalPrimary(iMc)!=kTRUE ) continue;
1476
1477       if( lPdgcodeCurrentPart == 3122 ){
1478          f3dHistPrimAnalysisPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1479          f3dHistPrimAnalysisPtVsYCMSVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1480       }
1481       if( lPdgcodeCurrentPart == -3122 ){
1482          f3dHistPrimAnalysisPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1483          f3dHistPrimAnalysisPtVsYCMSVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1484       }
1485       if( lPdgcodeCurrentPart == 310 ){
1486          f3dHistPrimAnalysisPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1487          f3dHistPrimAnalysisPtVsYCMSVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart+fpArapidityShift, lMultiplicity);
1488       }
1489    }
1490
1491 //----- Loop on primary Xi, Omega --------------------------------------------------------------
1492    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) 
1493    {// This is the begining of the loop on primaries
1494       
1495       TParticle* lCurrentParticlePrimary = 0x0; 
1496       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
1497       if(!lCurrentParticlePrimary){
1498          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1499          continue;
1500       }
1501       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) { 
1502          Double_t lRapXiMCPrimary = -100;
1503          if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) { 
1504            if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
1505              lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
1506            }
1507          }
1508
1509          //=================================================================================
1510          // Xi Histograms
1511          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){ 
1512             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1513             f3dHistGenSelectedPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1514          }
1515          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){ 
1516             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1517             f3dHistGenSelectedPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1518          }
1519          // Omega Histograms
1520          if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){ 
1521             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1522             f3dHistGenSelectedPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1523          }
1524          if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){ 
1525             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1526             f3dHistGenSelectedPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1527          }
1528       } 
1529    }
1530 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
1531
1532 //------------------------------------------------
1533 // MAIN LAMBDA LOOP STARTS HERE
1534 //------------------------------------------------
1535
1536    //Variable definition
1537    Int_t    lOnFlyStatus = 0;
1538    Double_t lChi2V0 = 0;
1539    Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
1540    Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
1541    Double_t lV0CosineOfPointingAngle = 0;
1542    Double_t lV0Radius = 0, lPt = 0;
1543    Double_t lRapK0Short = 0, lRapLambda = 0;
1544    Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
1545    Double_t lAlphaV0 = 0, lPtArmV0 = 0;
1546    Double_t fMinV0Pt = 0; 
1547    Double_t fMaxV0Pt = 100; 
1548
1549    Int_t nv0s = 0;
1550    nv0s = lESDevent->GetNumberOfV0s();
1551    
1552    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) 
1553         {// This is the begining of the V0 loop
1554       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
1555       if (!v0) continue;
1556
1557       //---> Fix On-the-Fly candidates
1558       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
1559         fHistSwappedV0Counter -> Fill( 1 );
1560       }else{
1561         fHistSwappedV0Counter -> Fill( 0 ); 
1562       }
1563       if ( fkUseOnTheFly ) CheckChargeV0(v0); 
1564
1565
1566       Double_t tV0mom[3];
1567       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] ); 
1568       Double_t lV0TotalMomentum = TMath::Sqrt(
1569          tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
1570
1571       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]); 
1572       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
1573       lPt = v0->Pt();
1574       lRapK0Short = v0->RapK0Short();
1575       lRapLambda  = v0->RapLambda();
1576
1577       //Set Variables for later filling
1578       fTreeVariableV0x = tDecayVertexV0[0];
1579       fTreeVariableV0y = tDecayVertexV0[1];
1580       fTreeVariableV0z = tDecayVertexV0[2];
1581
1582       //Set Variables for later filling
1583       fTreeVariableV0Px = tV0mom[0];
1584       fTreeVariableV0Py = tV0mom[1];
1585       fTreeVariableV0Pz = tV0mom[2];
1586
1587       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
1588
1589       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
1590       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
1591
1592       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1593       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1594
1595       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
1596       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
1597       if (!pTrack || !nTrack) {
1598          Printf("ERROR: Could not retreive one of the daughter track");
1599          continue;
1600       }
1601
1602       fTreeVariableNegEta = nTrack->Eta();
1603       fTreeVariablePosEta = pTrack->Eta();
1604
1605       // Filter like-sign V0 (next: add counter and distribution)
1606       if ( pTrack->GetSign() == nTrack->GetSign()){
1607          continue;
1608       } 
1609
1610       //________________________________________________________________________
1611       // Track quality cuts 
1612       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
1613       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
1614       fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
1615       if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
1616          fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
1617
1618       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
1619       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;      
1620       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
1621
1622       if ( ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) )&&(fkTakeAllTracks==kFALSE) ) continue;
1623         
1624       //GetKinkIndex condition
1625       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
1626
1627       //Findable clusters > 0 condition
1628       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
1629
1630       //Compute ratio Crossed Rows / Findable clusters
1631       //Note: above test avoids division by zero! 
1632       Float_t lPosTrackCrossedRowsOverFindable = -1;
1633       Float_t lNegTrackCrossedRowsOverFindable = -1;
1634       if ( ((double)(pTrack->GetTPCNclsF()) ) != 0 ) lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); 
1635       if ( ((double)(nTrack->GetTPCNclsF()) ) != 0 ) lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); 
1636
1637       fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
1638       if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
1639          fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
1640
1641       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
1642       if ( (fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8)&&(fkTakeAllTracks==kFALSE) ) continue;
1643
1644       //End track Quality Cuts
1645       //________________________________________________________________________
1646
1647       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(lPrimaryVtxPosition[0],
1648                                                         lPrimaryVtxPosition[1],
1649                                                         lMagneticField) );
1650
1651       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(lPrimaryVtxPosition[0],
1652                                                         lPrimaryVtxPosition[1],
1653                                                         lMagneticField) );
1654
1655       lOnFlyStatus = v0->GetOnFlyStatus();
1656       lChi2V0 = v0->GetChi2V0();
1657       lDcaV0Daughters = v0->GetDcaV0Daughters();
1658       lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1659       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1660       fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
1661
1662       // Getting invariant mass infos directly from ESD
1663       v0->ChangeMassHypothesis(310);
1664       lInvMassK0s = v0->GetEffMass();
1665       v0->ChangeMassHypothesis(3122);
1666       lInvMassLambda = v0->GetEffMass();
1667       v0->ChangeMassHypothesis(-3122);
1668       lInvMassAntiLambda = v0->GetEffMass();
1669       lAlphaV0 = v0->AlphaV0();
1670       lPtArmV0 = v0->PtArmV0();
1671
1672       //fTreeVariableOnFlyStatus = lOnFlyStatus;
1673       //fHistV0OnFlyStatus->Fill(lOnFlyStatus);
1674
1675 //===============================================
1676 // Monte Carlo Association starts here
1677 //===============================================
1678
1679       //---> Set Everything to "I don't know" before starting
1680
1681       fTreeVariablePIDPositive = 0;
1682       fTreeVariablePIDNegative = 0;
1683
1684       fTreeVariableIndexStatus = 0;
1685       fTreeVariableIndexStatusMother = 0;
1686
1687       fTreeVariablePtMother = -1;
1688       fTreeVariablePtMC = -1;
1689       fTreeVariableRapMC = -100;
1690
1691       fTreeVariablePID = -1; 
1692       fTreeVariablePIDMother = -1;
1693
1694       fTreeVariablePrimaryStatus = 0; 
1695       fTreeVariablePrimaryStatusMother = 0; 
1696       fTreeVariableV0CreationRadius = -1;
1697
1698       Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrack->GetLabel() );  
1699       Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrack->GetLabel() );
1700                 
1701       TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1702       TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1703             
1704       fTreeVariablePosTransvMomentumMC = mcPosV0Dghter->Pt();
1705       fTreeVariableNegTransvMomentumMC = mcNegV0Dghter->Pt();
1706
1707       Int_t lPIDPositive = mcPosV0Dghter -> GetPdgCode();
1708       Int_t lPIDNegative = mcNegV0Dghter -> GetPdgCode();
1709
1710       fTreeVariablePIDPositive = lPIDPositive;
1711       fTreeVariablePIDNegative = lPIDNegative;
1712
1713       Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
1714       Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
1715
1716       if( lblMotherPosV0Dghter == lblMotherNegV0Dghter && lblMotherPosV0Dghter > -1 ){
1717          //either label is fine, they're equal at this stage
1718          TParticle* pThisV0 = lMCstack->Particle( lblMotherPosV0Dghter ); 
1719          //Set tree variables
1720          fTreeVariablePID   = pThisV0->GetPdgCode(); //PDG Code
1721          fTreeVariablePtMC  = pThisV0->Pt(); //Perfect Pt
1722
1723           fTreeVariableIsNonInjected = lMCevent->IsFromBGEvent(lblMotherPosV0Dghter); //Is it? 
1724           if (!fTreeVariableIsNonInjected){
1725             if (!(pThisV0->GetFirstMother()<0)) 
1726               {fTreeVariableIsNonInjected = kTRUE;} // because there are primaries (ALICE definition) not produced in the collision
1727           }
1728
1729          //Set Variables for later filling
1730          //Be careful: Vx, Vy, Vz: Creation vertex. So decay position is the 
1731          //Creation vertex of any one of the daughters!
1732          fTreeVariableMCV0x = mcPosV0Dghter->Vx();
1733          fTreeVariableMCV0y = mcPosV0Dghter->Vy();
1734          fTreeVariableMCV0z = mcPosV0Dghter->Vz();
1735
1736          //Set Variables for later filling
1737          fTreeVariableMCV0Px = pThisV0->Px();
1738          fTreeVariableMCV0Py = pThisV0->Py();
1739          fTreeVariableMCV0Pz = pThisV0->Pz();
1740
1741          //Only Interested if it's a Lambda, AntiLambda or K0s 
1742          //Avoid the Junction Bug! PYTHIA has particles with Px=Py=Pz=E=0 occasionally, 
1743          //having particle code 88 (unrecognized by PDG), for documentation purposes.
1744          //Even ROOT's TParticle::Y() is not prepared to deal with that exception!
1745          //Note that TParticle::Pt() is immune (that would just return 0)...
1746          //Though granted that that should be extremely rare in this precise condition...
1747          if( TMath::Abs(fTreeVariablePID) == 3122 || fTreeVariablePID==310 ){
1748             fTreeVariableRapMC = pThisV0->Y(); //Perfect Y
1749          }
1750          fTreeVariableV0CreationRadius = TMath::Sqrt(
1751           TMath::Power(  ( (mcPrimaryVtx.At(0)) - (pThisV0->Vx()) ) , 2) + 
1752           TMath::Power(  ( (mcPrimaryVtx.At(1)) - (pThisV0->Vy()) ) , 2) + 
1753           TMath::Power(  ( (mcPrimaryVtx.At(2)) - (pThisV0->Vz()) ) , 2) 
1754          );
1755          if( lblMotherPosV0Dghter  < lNbMCPrimary ) fTreeVariableIndexStatus = 1; //looks primary
1756          if( lblMotherPosV0Dghter >= lNbMCPrimary ) fTreeVariableIndexStatus = 2; //looks secondary
1757          if( lMCstack->IsPhysicalPrimary       (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 1; //Is Primary!
1758          if( lMCstack->IsSecondaryFromWeakDecay(lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 2; //Weak Decay!
1759          if( lMCstack->IsSecondaryFromMaterial (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 3; //Material Int!
1760          
1761          //Now we try to acquire the V0 parent particle, if possible
1762          Int_t lblThisV0Parent = pThisV0->GetFirstMother();
1763          if ( lblThisV0Parent > -1 ){ //if it has a parent, get it and store specs
1764             TParticle* pThisV0Parent = lMCstack->Particle( lblThisV0Parent );
1765             fTreeVariablePIDMother   = pThisV0Parent->GetPdgCode(); //V0 Mother PDG
1766             fTreeVariablePtMother    = pThisV0Parent->Pt();         //V0 Mother Pt
1767             //Primary Status for the V0 Mother particle 
1768             if( lblThisV0Parent  < lNbMCPrimary ) fTreeVariableIndexStatusMother = 1; //looks primary
1769             if( lblThisV0Parent >= lNbMCPrimary ) fTreeVariableIndexStatusMother = 2; //looks secondary
1770             if( lMCstack->IsPhysicalPrimary       (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 1; //Is Primary!
1771             if( lMCstack->IsSecondaryFromWeakDecay(lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 2; //Weak Decay!
1772             if( lMCstack->IsSecondaryFromMaterial (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 3; //Material Int!
1773          }
1774       }
1775
1776       fTreeVariablePt = v0->Pt();
1777       fTreeVariableChi2V0 = lChi2V0; 
1778       fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
1779       fTreeVariableDcaV0Daughters = lDcaV0Daughters;
1780       fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle; 
1781       fTreeVariableV0Radius = lV0Radius;
1782       fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
1783       fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
1784       fTreeVariableInvMassK0s = lInvMassK0s;
1785       fTreeVariableInvMassLambda = lInvMassLambda;
1786       fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
1787       fTreeVariableRapK0Short = lRapK0Short;
1788
1789       fTreeVariableRapLambda = lRapLambda;
1790       fTreeVariableAlphaV0 = lAlphaV0;
1791       fTreeVariablePtArmV0 = lPtArmV0;
1792
1793       //Official means of acquiring N-sigmas 
1794       fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
1795       fTreeVariableNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
1796       fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
1797       fTreeVariableNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
1798
1799 //tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]
1800       Double_t lDistanceTravelled = TMath::Sqrt(
1801                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
1802                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
1803                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
1804                                         );
1805       fTreeVariableDistOverTotMom = 1e+5;
1806       if( lV0TotalMomentum + 1e-10 != 0 ) fTreeVariableDistOverTotMom = lDistanceTravelled / (lV0TotalMomentum + 1e-10); //avoid division by zero, to be sure
1807
1808       Double_t lMomentumPosTemp[3];
1809       pTrack->GetPxPyPz(lMomentumPosTemp);
1810       Double_t lPtPosTemporary = sqrt(pow(lMomentumPosTemp[0],2) + pow(lMomentumPosTemp[1],2));
1811
1812       Double_t lMomentumNegTemp[3];
1813       nTrack->GetPxPyPz(lMomentumNegTemp);
1814       Double_t lPtNegTemporary = sqrt(pow(lMomentumNegTemp[0],2) + pow(lMomentumNegTemp[1],2));
1815
1816       fTreeVariablePosTransvMomentum = lPtPosTemporary;
1817       fTreeVariableNegTransvMomentum = lPtNegTemporary;
1818
1819
1820 //------------------------------------------------
1821 // Fill Tree! 
1822 //------------------------------------------------
1823
1824       // The conditionals are meant to decrease excessive
1825       // memory usage! 
1826
1827       //Modified version: Keep only OnFlyStatus == 0
1828       //Keep only if included in a parametric InvMass Region 20 sigmas away from peak
1829
1830       //First Selection: Reject OnFly
1831       if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){
1832          //Second Selection: rough 20-sigma band, parametric. 
1833          //K0Short: Enough to parametrize peak broadening with linear function.    
1834          Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt; 
1835          Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
1836          //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
1837          //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
1838          Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt); 
1839          Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
1840          //Do Selection      
1841          if( (fTreeVariableInvMassLambda     < lUpperLimitLambda  && fTreeVariableInvMassLambda     > lLowerLimitLambda     ) || 
1842              (fTreeVariableInvMassAntiLambda < lUpperLimitLambda  && fTreeVariableInvMassAntiLambda > lLowerLimitLambda     ) || 
1843              (fTreeVariableInvMassK0s        < lUpperLimitK0Short && fTreeVariableInvMassK0s        > lLowerLimitK0Short    ) ){
1844              //Pre-selection in case this is AA...
1845              if( fkIsNuclear == kFALSE ) fTree->Fill();
1846              if( fkIsNuclear == kTRUE){ 
1847              //If this is a nuclear collision___________________
1848              // ... pre-filter with daughter eta selection only (not TPC)
1849                if ( TMath::Abs(fTreeVariableNegEta)<0.8 && TMath::Abs(fTreeVariablePosEta)<0.8 ) fTree->Fill();
1850              }//end nuclear_____________________________________
1851          }
1852       }
1853
1854 //------------------------------------------------
1855 // Fill tree over.
1856 //------------------------------------------------
1857
1858
1859    }// This is the end of the V0 loop
1860
1861 //------------------------------------------------
1862
1863    // Post output data.
1864    PostData(1, fListHistV0);
1865    PostData(2, fTree);
1866 }
1867
1868 //________________________________________________________________________
1869 void AliAnalysisTaskExtractPerformanceV0::Terminate(Option_t *)
1870 {
1871    // Draw result to the screen
1872    // Called once at the end of the query
1873
1874    TList *cRetrievedList = 0x0;
1875    cRetrievedList = (TList*)GetOutputData(1);
1876    if(!cRetrievedList){
1877       Printf("ERROR - AliAnalysisTaskExtractV0 : ouput data container list not available\n");
1878       return;
1879    }    
1880         
1881    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1882    if (!fHistV0MultiplicityForTrigEvt) {
1883       Printf("ERROR - AliAnalysisTaskExtractV0 : fHistV0MultiplicityForTrigEvt not available");
1884       return;
1885    }
1886   
1887    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractV0","V0 Multiplicity",10,10,510,510);
1888    canCheck->cd(1)->SetLogy();
1889
1890    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1891    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1892 }
1893
1894 //----------------------------------------------------------------------------
1895
1896 Double_t AliAnalysisTaskExtractPerformanceV0::MyRapidity(Double_t rE, Double_t rPz) const
1897 {
1898    // Local calculation for rapidity
1899    Double_t ReturnValue = -100;
1900    if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){ 
1901       ReturnValue =  0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1902    }
1903    return ReturnValue;
1904
1905
1906 //________________________________________________________________________
1907 void AliAnalysisTaskExtractPerformanceV0::CheckChargeV0(AliESDv0 *v0)
1908 {
1909    // This function checks charge of negative and positive daughter tracks. 
1910    // If incorrectly defined (onfly vertexer), swaps out. 
1911    if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
1912       //V0 daughter track swapping is required! Note: everything is swapped here... P->N, N->P
1913       Long_t lCorrectNidx = v0->GetPindex();
1914       Long_t lCorrectPidx = v0->GetNindex();
1915       Double32_t        lCorrectNmom[3];
1916       Double32_t        lCorrectPmom[3];
1917       v0->GetPPxPyPz( lCorrectNmom[0], lCorrectNmom[1], lCorrectNmom[2] );
1918       v0->GetNPxPyPz( lCorrectPmom[0], lCorrectPmom[1], lCorrectPmom[2] );
1919
1920       AliExternalTrackParam     lCorrectParamN(
1921         v0->GetParamP()->GetX() , 
1922         v0->GetParamP()->GetAlpha() , 
1923         v0->GetParamP()->GetParameter() , 
1924         v0->GetParamP()->GetCovariance() 
1925       );
1926       AliExternalTrackParam     lCorrectParamP(
1927         v0->GetParamN()->GetX() , 
1928         v0->GetParamN()->GetAlpha() , 
1929         v0->GetParamN()->GetParameter() , 
1930         v0->GetParamN()->GetCovariance() 
1931       );
1932       lCorrectParamN.SetMostProbablePt( v0->GetParamP()->GetMostProbablePt() );
1933       lCorrectParamP.SetMostProbablePt( v0->GetParamN()->GetMostProbablePt() );
1934
1935       //Get Variables___________________________________________________
1936       Double_t lDcaV0Daughters = v0 -> GetDcaV0Daughters();
1937       Double_t lCosPALocal     = v0 -> GetV0CosineOfPointingAngle(); 
1938       Bool_t lOnFlyStatusLocal = v0 -> GetOnFlyStatus();
1939
1940       //Create Replacement Object_______________________________________
1941       AliESDv0 *v0correct = new AliESDv0(lCorrectParamN,lCorrectNidx,lCorrectParamP,lCorrectPidx);
1942       v0correct->SetDcaV0Daughters          ( lDcaV0Daughters   );
1943       v0correct->SetV0CosineOfPointingAngle ( lCosPALocal       );
1944       v0correct->ChangeMassHypothesis       ( kK0Short          );
1945       v0correct->SetOnFlyStatus             ( lOnFlyStatusLocal );
1946
1947       //Reverse Cluster info..._________________________________________
1948       v0correct->SetClusters( v0->GetClusters( 1 ), v0->GetClusters ( 0 ) );
1949
1950       *v0 = *v0correct;
1951       //Proper cleanup..._______________________________________________
1952       v0correct->Delete();
1953       v0correct = 0x0;
1954
1955       //Just another cross-check and output_____________________________
1956       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ) {
1957         AliWarning("Found Swapped Charges, tried to correct but something FAILED!");
1958       }else{
1959         //AliWarning("Found Swapped Charges and fixed.");
1960       }
1961       //________________________________________________________________
1962    }else{
1963       //Don't touch it! ---
1964       //Printf("Ah, nice. Charges are already ordered...");
1965    }
1966    return;
1967
1968