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