Compatibility with the Root trunk
[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
90 #include "AliAnalysisTaskExtractPerformanceV0.h"
91
92 using std::cout;
93 using std::endl;
94
95 ClassImp(AliAnalysisTaskExtractPerformanceV0)
96
97 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0() 
98   : AliAnalysisTaskSE(), fListHistV0(0), fTree(0), fPIDResponse(0),
99    fkIsNuclear   ( kFALSE ), 
100    fkLowEnergyPP ( kFALSE ),
101    fkUseOnTheFly ( kFALSE ),
102
103 //------------------------------------------------
104 // HISTOGRAMS
105 // --- Filled on an Event-by-event basis
106 //------------------------------------------------
107    fHistV0MultiplicityBeforeTrigSel(0),
108    fHistV0MultiplicityForTrigEvt(0), 
109    fHistV0MultiplicityForSelEvt(0),
110    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
111    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
112    fHistMultiplicityBeforeTrigSel(0),
113    fHistMultiplicityForTrigEvt(0),
114    fHistMultiplicity(0),
115    fHistMultiplicityNoTPCOnly(0),
116    fHistMultiplicityNoTPCOnlyNoPileup(0),
117
118 //------------------------------------------------
119 // PARTICLE HISTOGRAMS
120 // --- Filled on a Particle-by-Particle basis
121 //------------------------------------------------
122    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
123    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
124    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
125    f3dHistPrimRawPtVsYVsMultLambda(0),
126    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
127    f3dHistPrimRawPtVsYVsMultK0Short(0),
128    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
129    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
130    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
131    f3dHistGenPtVsYVsMultXiMinus(0),
132    f3dHistGenPtVsYVsMultXiPlus(0),
133    fHistPVx(0),
134    fHistPVy(0),
135    fHistPVz(0),
136    fHistPVxAnalysis(0),
137    fHistPVyAnalysis(0),
138    fHistPVzAnalysis(0),
139    fHistPVxAnalysisHasHighPtLambda(0),
140    fHistPVyAnalysisHasHighPtLambda(0),
141    fHistPVzAnalysisHasHighPtLambda(0),
142    fHistSwappedV0Counter(0)
143 {
144   // Dummy Constructor
145 }
146
147 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0(const char *name) 
148   : AliAnalysisTaskSE(name), fListHistV0(0), fTree(0), fPIDResponse(0),
149    fkIsNuclear   ( kFALSE ), 
150    fkLowEnergyPP ( kFALSE ),
151    fkUseOnTheFly ( kFALSE ),
152      
153 //------------------------------------------------
154 // HISTOGRAMS
155 // --- Filled on an Event-by-event basis
156 //------------------------------------------------
157    fHistV0MultiplicityBeforeTrigSel(0),
158    fHistV0MultiplicityForTrigEvt(0), 
159    fHistV0MultiplicityForSelEvt(0),
160    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
161    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
162    fHistMultiplicityBeforeTrigSel(0),
163    fHistMultiplicityForTrigEvt(0),
164    fHistMultiplicity(0),
165    fHistMultiplicityNoTPCOnly(0),
166    fHistMultiplicityNoTPCOnlyNoPileup(0),
167
168
169 //------------------------------------------------
170 // PARTICLE HISTOGRAMS
171 // --- Filled on a Particle-by-Particle basis
172 //------------------------------------------------
173    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
174    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
175    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
176    f3dHistPrimRawPtVsYVsMultLambda(0),
177    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
178    f3dHistPrimRawPtVsYVsMultK0Short(0),
179    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
180    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
181    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
182    f3dHistGenPtVsYVsMultXiMinus(0),
183    f3dHistGenPtVsYVsMultXiPlus(0),
184    fHistPVx(0),
185    fHistPVy(0),
186    fHistPVz(0),
187    fHistPVxAnalysis(0),
188    fHistPVyAnalysis(0),
189    fHistPVzAnalysis(0),
190    fHistPVxAnalysisHasHighPtLambda(0),
191    fHistPVyAnalysisHasHighPtLambda(0),
192    fHistPVzAnalysisHasHighPtLambda(0),
193    fHistSwappedV0Counter(0)
194 {
195    // Constructor
196    // Output slot #0 writes into a TList container (Cascade)
197    DefineOutput(1, TList::Class());
198    DefineOutput(2, TTree::Class());
199 }
200
201
202 AliAnalysisTaskExtractPerformanceV0::~AliAnalysisTaskExtractPerformanceV0()
203 {
204 //------------------------------------------------
205 // DESTRUCTOR
206 //------------------------------------------------
207
208    if (fListHistV0){
209       delete fListHistV0;
210       fListHistV0 = 0x0;
211    }
212    if (fTree){
213       delete fTree;
214       fTree = 0x0;
215    }
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskExtractPerformanceV0::UserCreateOutputObjects()
220 {
221
222    OpenFile(2); 
223    // Called once
224    fTree = new TTree("fTree","V0Candidates");
225
226 //------------------------------------------------
227 // fTree Branch definitions
228 //------------------------------------------------
229
230 //-----------BASIC-INFO---------------------------
231 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatus",&fTreeVariablePrimaryStatus,"fTreeVariablePrimaryStatus/I");        
232 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatusMother",&fTreeVariablePrimaryStatusMother,"fTreeVariablePrimaryStatusMother/I");      
233 /* 2*/   fTree->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"Chi2V0/F");
234 /* 3*/   fTree->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
235 /* 4*/   fTree->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
236 /* 5*/   fTree->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
237 /* 6*/   fTree->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
238 /* 7*/   fTree->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
239 /* 7*/   fTree->Branch("fTreeVariablePtMC",&fTreeVariablePtMC,"fTreeVariablePtMC/F");
240 /* 8*/   fTree->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
241 /* 9*/   fTree->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
242 /*10*/   fTree->Branch("fTreeVariableRapMC",&fTreeVariableRapMC,"fTreeVariableRapMC/F");
243 /*11*/   fTree->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
244 /*12*/   fTree->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
245 /*13*/   fTree->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
246 /*14*/   fTree->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
247 /*15*/   fTree->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
248 /*16*/   fTree->Branch("fTreeVariableNegTransvMomentum",&fTreeVariableNegTransvMomentum,"fTreeVariableNegTransvMomentum/F");
249 /*17*/   fTree->Branch("fTreeVariablePosTransvMomentum",&fTreeVariablePosTransvMomentum,"fTreeVariablePosTransvMomentum/F");
250 /*18*/   fTree->Branch("fTreeVariableNegTransvMomentumMC",&fTreeVariableNegTransvMomentumMC,"fTreeVariableNegTransvMomentumMC/F");
251 /*19*/   fTree->Branch("fTreeVariablePosTransvMomentumMC",&fTreeVariablePosTransvMomentumMC,"fTreeVariablePosTransvMomentumMC/F");
252 /*20*/   fTree->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
253 /*21*/   fTree->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
254 /*22*/   fTree->Branch("fTreeVariablePID",&fTreeVariablePID,"fTreeVariablePID/I");
255 /*23*/   fTree->Branch("fTreeVariablePIDPositive",&fTreeVariablePIDPositive,"fTreeVariablePIDPositive/I");
256 /*24*/   fTree->Branch("fTreeVariablePIDNegative",&fTreeVariablePIDNegative,"fTreeVariablePIDNegative/I");
257 /*25*/   fTree->Branch("fTreeVariablePIDMother",&fTreeVariablePIDMother,"fTreeVariablePIDMother/I");
258 /*26*/   fTree->Branch("fTreeVariablePtXiMother",&fTreeVariablePtMother,"fTreeVariablePtMother/F");
259 /*27*/   fTree->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
260 //-----------MULTIPLICITY-INFO--------------------
261 /*28*/   fTree->Branch("fTreeVariableMultiplicity",&fTreeVariableMultiplicity,"fTreeVariableMultiplicity/I");
262 //------------------------------------------------
263 /*29*/   fTree->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
264 /*30*/   fTree->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
265 /*31*/   fTree->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
266 /*32*/   fTree->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
267 /*33*/   fTree->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
268 //------------------------------------------------
269 /*34*/   fTree->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
270 /*35*/   fTree->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
271 /*36*/   fTree->Branch("fTreeVariableV0CreationRadius",&fTreeVariableV0CreationRadius,"fTreeVariableV0CreationRadius/F");
272 /*37*/   fTree->Branch("fTreeVariableIndexStatus",&fTreeVariableIndexStatus,"fTreeVariableIndexStatus/I");
273 /*38*/   fTree->Branch("fTreeVariableIndexStatusMother",&fTreeVariableIndexStatusMother,"fTreeVariableIndexStatusMother/I");
274
275 //------------------------------------------------
276 // Particle Identification Setup
277 //------------------------------------------------
278
279    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
280    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
281    fPIDResponse = inputHandler->GetPIDResponse();
282
283 //------------------------------------------------
284 // V0 Multiplicity Histograms
285 //------------------------------------------------
286
287    // Create histograms
288    OpenFile(1);
289    fListHistV0 = new TList();
290    fListHistV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
291
292
293    if(! fHistV0MultiplicityBeforeTrigSel) {
294       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", 
295          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", 
296          25, 0, 25);
297       fListHistV0->Add(fHistV0MultiplicityBeforeTrigSel);
298    }
299            
300    if(! fHistV0MultiplicityForTrigEvt) {
301       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", 
302          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", 
303          25, 0, 25);
304       fListHistV0->Add(fHistV0MultiplicityForTrigEvt);
305    }
306
307    if(! fHistV0MultiplicityForSelEvt) {
308       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", 
309          "V0s per event;Nbr of V0s/Evt;Events", 
310          25, 0, 25);
311       fListHistV0->Add(fHistV0MultiplicityForSelEvt);
312    }
313
314    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
315       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", 
316          "V0s per event;Nbr of V0s/Evt;Events", 
317          25, 0, 25);
318       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
319    }
320    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
321       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", 
322          "V0s per event;Nbr of V0s/Evt;Events", 
323          25, 0, 25);
324       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
325    }
326
327 //------------------------------------------------
328 // Track Multiplicity Histograms
329 //------------------------------------------------
330
331    if(! fHistMultiplicityBeforeTrigSel) {
332       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", 
333          "Tracks per event;Nbr of Tracks;Events", 
334          200, 0, 200);          
335       fListHistV0->Add(fHistMultiplicityBeforeTrigSel);
336    }
337    if(! fHistMultiplicityForTrigEvt) {
338       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", 
339          "Tracks per event;Nbr of Tracks;Events", 
340          200, 0, 200);          
341       fListHistV0->Add(fHistMultiplicityForTrigEvt);
342    }
343    if(! fHistMultiplicity) {
344       fHistMultiplicity = new TH1F("fHistMultiplicity", 
345          "Tracks per event;Nbr of Tracks;Events", 
346          200, 0, 200);          
347       fListHistV0->Add(fHistMultiplicity);
348    }
349    if(! fHistMultiplicityNoTPCOnly) {
350       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", 
351          "Tracks per event;Nbr of Tracks;Events", 
352          200, 0, 200);          
353       fListHistV0->Add(fHistMultiplicityNoTPCOnly);
354    }
355    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
356       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", 
357          "Tracks per event;Nbr of Tracks;Events", 
358          200, 0, 200);          
359       fListHistV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
360    }
361
362 //------------------------------------------------
363 // Generated Particle Histograms
364 //------------------------------------------------
365
366    Int_t lCustomNBins = 200; 
367    Double_t lCustomPtUpperLimit = 20; 
368    Int_t lCustomNBinsMultiplicity = 100;
369
370 //----------------------------------
371 // Raw Generated (Pre-physics-selection)
372 //----------------------------------
373
374 //--- 3D Histo (Pt, Y, Multiplicity)  
375
376    if(! f3dHistPrimRawPtVsYVsMultLambda) {
377       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);
378       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultLambda);
379    }
380    if(! f3dHistPrimRawPtVsYVsMultAntiLambda) {
381       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);
382       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultAntiLambda);
383    }
384    if(! f3dHistPrimRawPtVsYVsMultK0Short) {
385       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);
386       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultK0Short);
387    }
388
389 //--- 3D Histo (Pt, Y, Proper Decay Length)
390
391    if(! f3dHistPrimRawPtVsYVsDecayLengthLambda) {
392       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);
393       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthLambda);
394    }
395    if(! f3dHistPrimRawPtVsYVsDecayLengthAntiLambda) {
396       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);
397       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthAntiLambda);
398    }
399    if(! f3dHistPrimRawPtVsYVsDecayLengthK0Short) {
400       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);
401       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthK0Short);
402    }
403
404 //--- 3D Histo (Pt, Y, Multiplicity) for generated charged Xi (feeddown)
405
406    if(! f3dHistGenPtVsYVsMultXiMinus) {
407       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);
408       fListHistV0->Add(f3dHistGenPtVsYVsMultXiMinus);
409    }
410    if(! f3dHistGenPtVsYVsMultXiPlus) {
411       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);
412       fListHistV0->Add(f3dHistGenPtVsYVsMultXiPlus);
413    }
414
415 //----------------------------------
416 // Histos at analysis level 
417 //----------------------------------
418
419    if(! f3dHistPrimAnalysisPtVsYVsMultLambda) {
420       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);
421       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultLambda);
422    }
423    if(! f3dHistPrimAnalysisPtVsYVsMultAntiLambda) {
424       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);
425       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultAntiLambda);
426    }
427    if(! f3dHistPrimAnalysisPtVsYVsMultK0Short) {
428       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);
429       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultK0Short);
430    }
431
432 //----------------------------------
433 // Primary Vertex Position Histos
434 //----------------------------------
435
436    if(! fHistPVx) {
437          fHistPVx = new TH1F("fHistPVx", 
438             "PV x position;Nbr of Evts;x", 
439             2000, -0.5, 0.5);       
440       fListHistV0->Add(fHistPVx);
441    }
442    if(! fHistPVy) {
443          fHistPVy = new TH1F("fHistPVy", 
444             "PV y position;Nbr of Evts;y", 
445             2000, -0.5, 0.5);       
446       fListHistV0->Add(fHistPVy);
447    }
448    if(! fHistPVz) {
449          fHistPVz = new TH1F("fHistPVz", 
450             "PV z position;Nbr of Evts;z", 
451             400, -20, 20);       
452       fListHistV0->Add(fHistPVz);
453    }
454
455    if(! fHistPVxAnalysis) {
456          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis", 
457             "PV x position;Nbr of Evts;x", 
458             2000, -0.5, 0.5);       
459       fListHistV0->Add(fHistPVxAnalysis);
460    }
461    if(! fHistPVyAnalysis) {
462          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis", 
463             "PV y position;Nbr of Evts;y", 
464             2000, -0.5, 0.5);       
465       fListHistV0->Add(fHistPVyAnalysis);
466    }
467    if(! fHistPVzAnalysis) {
468          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis", 
469             "PV z position;Nbr of Evts;z", 
470             400, -20, 20);       
471       fListHistV0->Add(fHistPVzAnalysis);
472    }
473
474    if(! fHistPVxAnalysisHasHighPtLambda) {
475          fHistPVxAnalysisHasHighPtLambda = new TH1F("fHistPVxAnalysisHasHighPtLambda", 
476             "PV x position;Nbr of Evts;x", 
477             2000, -0.5, 0.5);       
478       fListHistV0->Add(fHistPVxAnalysisHasHighPtLambda);
479    }
480    if(! fHistPVyAnalysisHasHighPtLambda) {
481          fHistPVyAnalysisHasHighPtLambda = new TH1F("fHistPVyAnalysisHasHighPtLambda", 
482             "PV y position;Nbr of Evts;y", 
483             2000, -0.5, 0.5);       
484       fListHistV0->Add(fHistPVyAnalysisHasHighPtLambda);
485    }
486    if(! fHistPVzAnalysisHasHighPtLambda) {
487          fHistPVzAnalysisHasHighPtLambda = new TH1F("fHistPVzAnalysisHasHighPtLambda", 
488             "PV z position;Nbr of Evts;z", 
489             400, -20, 20);       
490       fListHistV0->Add(fHistPVzAnalysisHasHighPtLambda);
491    }
492    if(! fHistSwappedV0Counter) {
493       fHistSwappedV0Counter = new TH1F("fHistSwappedV0Counter", 
494          "Swap or not histo;Swapped (1) or not (0); count", 
495          2, 0, 2);              
496       fListHistV0->Add(fHistSwappedV0Counter);
497    }
498
499    //List of Histograms: Normal
500    PostData(1, fListHistV0);
501
502    //TTree Object: Saved to base directory. Should cache to disk while saving. 
503    //(Important to avoid excessive memory usage, particularly when merging)
504    PostData(2, fTree);
505
506 }// end UserCreateOutputObjects
507
508
509 //________________________________________________________________________
510 void AliAnalysisTaskExtractPerformanceV0::UserExec(Option_t *) 
511 {
512   // Main loop
513   // Called for each event
514
515    AliESDEvent *lESDevent = 0x0;
516    AliMCEvent  *lMCevent  = 0x0; 
517    AliStack    *lMCstack  = 0x0; 
518
519    Int_t    lNumberOfV0s                      = -1;
520    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
521    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
522    Double_t lMagneticField                 = -10.;
523    
524   // Connect to the InputEvent   
525   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
526
527    // Appropriate for ESD analysis! 
528       
529    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
530    if (!lESDevent) {
531       AliWarning("ERROR: lESDevent not available \n");
532       return;
533    }
534         
535    lMCevent = MCEvent();
536    if (!lMCevent) {
537       Printf("ERROR: Could not retrieve MC event \n");
538       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;   
539       return;
540    }
541
542    lMCstack = lMCevent->Stack();
543    if (!lMCstack) {
544       Printf("ERROR: Could not retrieve MC stack \n");
545       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
546       return;
547    }
548         
549 //------------------------------------------------
550 // Multiplicity Information Acquistion
551 //------------------------------------------------
552
553    //REVISED multiplicity estimator after 'multiplicity day' (2011)
554    Int_t lMultiplicity = -100; 
555
556    if(fkIsNuclear == kFALSE) lMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity( lESDevent );
557
558    //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
559    //---> Warning: Experimental
560    if(fkIsNuclear == kTRUE){ 
561       AliCentrality* centrality;
562       centrality = lESDevent->GetCentrality();
563       lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
564       if (centrality->GetQuality()>1) {
565         PostData(1, fListHistV0);
566         PostData(2, fTree);
567         return;
568       }
569    }
570   
571    //Set variable for filling tree afterwards!
572    //---> pp case......: GetReferenceMultiplicity
573    //---> Pb-Pb case...: Centrality by V0M
574    fTreeVariableMultiplicity = lMultiplicity;
575
576    fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
577    fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
578         
579 //------------------------------------------------
580 // MC Information Acquistion
581 //------------------------------------------------
582
583    Int_t iNumberOfPrimaries = -1;
584    iNumberOfPrimaries = lMCstack->GetNprimary();
585    if(iNumberOfPrimaries < 1) return; 
586    Bool_t lHasHighPtLambda = kFALSE;
587
588 //------------------------------------------------
589 // Variable Definition
590 //------------------------------------------------
591
592    Int_t lNbMCPrimary        = 0;
593
594    Int_t lPdgcodeCurrentPart = 0;
595    Double_t lRapCurrentPart  = 0;
596    Double_t lPtCurrentPart   = 0;
597   
598    //Int_t lComeFromSigma      = 0;
599
600    // current mc particle 's mother
601    //Int_t iCurrentMother  = 0;
602    lNbMCPrimary = lMCstack->GetNprimary();
603
604 //------------------------------------------------
605 // Pre-Physics Selection
606 //------------------------------------------------
607
608 //----- Loop on primary Xi, Omega --------------------------------------------------------------
609    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) 
610    {// This is the begining of the loop on primaries
611       
612       TParticle* lCurrentParticlePrimary = 0x0; 
613       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
614       if(!lCurrentParticlePrimary){
615          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
616          continue;
617       }
618       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 ){ 
619          Double_t lRapXiMCPrimary = 0.5*TMath::Log((lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13));
620
621          //=================================================================================
622          // Xi Histograms for Feeddown - Primary Charged Xis
623          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){ 
624             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
625             f3dHistGenPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
626          }
627          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){ 
628             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
629             f3dHistGenPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
630          }
631       } 
632    }
633 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
634
635 //----- Loop on Lambda, K0Short ----------------------------------------------------------------
636    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++) 
637    {// This is the begining of the loop on tracks
638       
639       TParticle* lCurrentParticleForLambdaCheck = 0x0; 
640       lCurrentParticleForLambdaCheck = lMCstack->Particle( iCurrentLabelStack );
641       if(!lCurrentParticleForLambdaCheck){
642          Printf("V0s loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
643          continue;
644       }
645
646       //=================================================================================
647       //Single-Strange checks
648       // Keep only K0s, Lambda and AntiLambda:
649       lPdgcodeCurrentPart = lCurrentParticleForLambdaCheck->GetPdgCode();             
650
651       if ( (lCurrentParticleForLambdaCheck->GetPdgCode() == 310   ) ||
652            (lCurrentParticleForLambdaCheck->GetPdgCode() == 3122  ) ||
653            (lCurrentParticleForLambdaCheck->GetPdgCode() == -3122 ) )
654            {
655          lRapCurrentPart   = MyRapidity(lCurrentParticleForLambdaCheck->Energy(),lCurrentParticleForLambdaCheck->Pz());
656          lPtCurrentPart    = lCurrentParticleForLambdaCheck->Pt();
657
658          //Use Physical Primaries only for filling PrimRaw Histograms!
659          if ( lMCstack->IsPhysicalPrimary(iCurrentLabelStack)!=kTRUE ) continue;
660
661          if( lPdgcodeCurrentPart == 3122 ){
662             f3dHistPrimRawPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
663             if( TMath::Abs( lCurrentParticleForLambdaCheck->Eta() )<1.2 && lPtCurrentPart>2 ){
664                lHasHighPtLambda = kTRUE; //Keep track of events with Lambda within |eta|<1.2 and pt>2
665             }
666          }
667          if( lPdgcodeCurrentPart == -3122 ){
668             f3dHistPrimRawPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
669          }
670          if( lPdgcodeCurrentPart == 310 ){
671             f3dHistPrimRawPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
672          }
673          //Decay Length Acquisition=====================================================
674          Double_t decaylength = -1; 
675          Double_t lV0Mass = -1; 
676           
677          if( !(lCurrentParticleForLambdaCheck->GetDaughter(0) < 0) ) {
678             TParticle* lDght0ofV0 = lMCstack->Particle(  lCurrentParticleForLambdaCheck->GetDaughter(0) ); //get first daughter
679             if(lDght0ofV0){ // skip if not defined. 
680                decaylength = TMath::Sqrt(
681                                         TMath::Power( lCurrentParticleForLambdaCheck->Vx() - lDght0ofV0->Vx() , 2) + 
682                                         TMath::Power( lCurrentParticleForLambdaCheck->Vy() - lDght0ofV0->Vy() , 2) + 
683                                         TMath::Power( lCurrentParticleForLambdaCheck->Vz() - lDght0ofV0->Vz() , 2)
684                );
685                //Need to correct for relativitity! Involves multiplying by mass and dividing by momentum. 
686                if(TMath::Abs( lPdgcodeCurrentPart ) == 3122 ) { lV0Mass = 1.115683; }
687                if(TMath::Abs( lPdgcodeCurrentPart ) == 310 ) { lV0Mass = 0.497614; }
688                decaylength = ( lV0Mass * decaylength ) / ( lCurrentParticleForLambdaCheck->P() + 1e-10 );
689             }
690          }
691          if( lPdgcodeCurrentPart == 3122) f3dHistPrimRawPtVsYVsDecayLengthLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
692          if( lPdgcodeCurrentPart == -3122) f3dHistPrimRawPtVsYVsDecayLengthAntiLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
693          if( lPdgcodeCurrentPart == 310) f3dHistPrimRawPtVsYVsDecayLengthK0Short ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength ); 
694       }
695    }//End of loop on tracks
696 //----- End Loop on Lambda, K0Short ------------------------------------------------------------
697
698 // ---> Set Variables to Zero again
699 // ---> Variable Definition
700
701    lPdgcodeCurrentPart = 0;
702    lRapCurrentPart  = 0;
703    lPtCurrentPart   = 0;
704
705 //------------------------------------------------
706 // Physics Selection
707 //------------------------------------------------
708
709    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
710    Bool_t isSelected = 0;
711    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
712
713    //pp at 2.76TeV: special case, ignore FastOnly
714    if ( (fkLowEnergyPP == kTRUE) && (maskIsSelected& AliVEvent::kFastOnly) ){
715       PostData(1, fListHistV0);
716       PostData(2, fTree);
717       return;
718    } 
719    //Standard Min-Bias Selection
720    if ( ! isSelected ) { 
721       PostData(1, fListHistV0);
722       PostData(2, fTree);
723       return;
724    }
725
726 //------------------------------------------------
727 // After Trigger Selection
728 //------------------------------------------------
729
730    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
731   
732    //Set variable for filling tree afterwards!
733    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
734    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
735
736 //------------------------------------------------
737 // Getting: Primary Vertex + MagField Info
738 //------------------------------------------------
739
740    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
741    // get the vtx stored in ESD found with tracks
742    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
743         
744    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();      
745    // get the best primary vertex available for the event
746    // As done in AliCascadeVertexer, we keep the one which is the best one available.
747    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
748    // This one will be used for next calculations (DCA essentially)
749    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
750
751    Double_t lPrimaryVtxPosition[3];
752    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
753    lPrimaryVtxPosition[0] = primaryVtx->GetX();
754    lPrimaryVtxPosition[1] = primaryVtx->GetY();
755    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
756    fHistPVx->Fill( lPrimaryVtxPosition[0] );
757    fHistPVy->Fill( lPrimaryVtxPosition[1] );
758    fHistPVz->Fill( lPrimaryVtxPosition[2] );
759
760 //------------------------------------------------
761 // Primary Vertex Z position: SKIP
762 //------------------------------------------------
763
764    if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) { 
765       AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !"); 
766       PostData(1, fListHistV0);
767       PostData(2, fTree);
768       return; 
769    }
770
771    lMagneticField = lESDevent->GetMagneticField( );
772    fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
773    fHistMultiplicity->Fill(lMultiplicity);
774
775 //------------------------------------------------
776 // SKIP: Events with well-established PVtx
777 //------------------------------------------------
778         
779    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
780    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
781    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
782       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
783       PostData(1, fListHistV0);
784       PostData(2, fTree);
785       return;
786    }
787    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
788    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
789
790 //------------------------------------------------
791 // Pileup Rejection Studies
792 //------------------------------------------------
793
794    // FIXME : quality selection regarding pile-up rejection 
795    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
796       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
797       PostData(1, fListHistV0);
798       PostData(2, fTree);
799       return;
800    }
801    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
802    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
803
804    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
805    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){ 
806       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
807       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
808       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
809       if ( lHasHighPtLambda == kTRUE ){ 
810          fHistPVxAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[0] );
811          fHistPVyAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[1] );
812          fHistPVzAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[2] );
813       }
814    }
815
816 //------------------------------------------------
817 // stack loop starts here
818 //------------------------------------------------
819
820 //---> Loop over ALL PARTICLES
821  
822    for (Int_t iMc = 0; iMc < (lMCstack->GetNtrack()); iMc++) {  
823       TParticle *p0 = lMCstack->Particle(iMc); 
824       if (!p0) {
825          //Printf("ERROR: particle with label %d not found in lMCstack (mc loop)", iMc);
826          continue;
827       }
828       lPdgcodeCurrentPart = p0->GetPdgCode();
829
830       // Keep only K0s, Lambda and AntiLambda:
831       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
832         
833       lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
834       lPtCurrentPart    = p0->Pt();
835
836         //Use Physical Primaries only for filling PrimRaw Histograms!
837       if ( lMCstack->IsPhysicalPrimary(iMc)!=kTRUE ) continue;
838
839       if( lPdgcodeCurrentPart == 3122 ){
840          f3dHistPrimAnalysisPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
841       }
842       if( lPdgcodeCurrentPart == -3122 ){
843          f3dHistPrimAnalysisPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
844       }
845       if( lPdgcodeCurrentPart == 310 ){
846          f3dHistPrimAnalysisPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
847       }
848    }
849
850 //------------------------------------------------
851 // MAIN LAMBDA LOOP STARTS HERE
852 //------------------------------------------------
853
854    //Variable definition
855    Int_t    lOnFlyStatus = 0;
856    Double_t lChi2V0 = 0;
857    Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
858    Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
859    Double_t lV0CosineOfPointingAngle = 0;
860    Double_t lV0Radius = 0, lPt = 0;
861    Double_t lRapK0Short = 0, lRapLambda = 0;
862    Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
863    Double_t lAlphaV0 = 0, lPtArmV0 = 0;
864    Double_t fMinV0Pt = 0; 
865    Double_t fMaxV0Pt = 100; 
866
867    Int_t nv0s = 0;
868    nv0s = lESDevent->GetNumberOfV0s();
869    
870    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) 
871         {// This is the begining of the V0 loop
872       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
873       if (!v0) continue;
874
875       //---> Fix On-the-Fly candidates
876       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
877         fHistSwappedV0Counter -> Fill( 1 );
878       }else{
879         fHistSwappedV0Counter -> Fill( 0 ); 
880       }
881       if ( fkUseOnTheFly ) CheckChargeV0(v0); 
882
883
884       Double_t tV0mom[3];
885       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] ); 
886       Double_t lV0TotalMomentum = TMath::Sqrt(
887          tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
888
889       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]); 
890       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
891       lPt = v0->Pt();
892       lRapK0Short = v0->RapK0Short();
893       lRapLambda  = v0->RapLambda();
894       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
895
896       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
897       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
898
899       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
900       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
901
902       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
903       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
904       if (!pTrack || !nTrack) {
905          Printf("ERROR: Could not retreive one of the daughter track");
906          continue;
907       }
908
909       fTreeVariableNegEta = nTrack->Eta();
910       fTreeVariablePosEta = pTrack->Eta();
911
912       // Filter like-sign V0 (next: add counter and distribution)
913       if ( pTrack->GetSign() == nTrack->GetSign()){
914          continue;
915       } 
916
917       //________________________________________________________________________
918       // Track quality cuts 
919       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
920       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
921       fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
922       if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
923          fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
924
925       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
926       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;      
927       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
928
929       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
930         
931       //GetKinkIndex condition
932       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
933
934       //Findable clusters > 0 condition
935       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
936
937       //Compute ratio Crossed Rows / Findable clusters
938       //Note: above test avoids division by zero! 
939       Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); 
940       Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); 
941
942       fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
943       if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
944          fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
945
946       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
947       if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
948
949       //End track Quality Cuts
950       //________________________________________________________________________
951
952       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(lPrimaryVtxPosition[0],
953                                                         lPrimaryVtxPosition[1],
954                                                         lMagneticField) );
955
956       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(lPrimaryVtxPosition[0],
957                                                         lPrimaryVtxPosition[1],
958                                                         lMagneticField) );
959
960       lOnFlyStatus = v0->GetOnFlyStatus();
961       lChi2V0 = v0->GetChi2V0();
962       lDcaV0Daughters = v0->GetDcaV0Daughters();
963       lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
964       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
965       fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
966
967       // Getting invariant mass infos directly from ESD
968       v0->ChangeMassHypothesis(310);
969       lInvMassK0s = v0->GetEffMass();
970       v0->ChangeMassHypothesis(3122);
971       lInvMassLambda = v0->GetEffMass();
972       v0->ChangeMassHypothesis(-3122);
973       lInvMassAntiLambda = v0->GetEffMass();
974       lAlphaV0 = v0->AlphaV0();
975       lPtArmV0 = v0->PtArmV0();
976
977       //fTreeVariableOnFlyStatus = lOnFlyStatus;
978       //fHistV0OnFlyStatus->Fill(lOnFlyStatus);
979
980 //===============================================
981 // Monte Carlo Association starts here
982 //===============================================
983
984       //---> Set Everything to "I don't know" before starting
985
986       fTreeVariablePIDPositive = 0;
987       fTreeVariablePIDNegative = 0;
988
989       fTreeVariableIndexStatus = 0;
990       fTreeVariableIndexStatusMother = 0;
991
992       fTreeVariablePtMother = -1;
993       fTreeVariablePtMC = -1;
994       fTreeVariableRapMC = -100;
995
996       fTreeVariablePID = -1; 
997       fTreeVariablePIDMother = -1;
998
999       fTreeVariablePrimaryStatus = 0; 
1000       fTreeVariablePrimaryStatusMother = 0; 
1001       fTreeVariableV0CreationRadius = -1;
1002
1003       Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrack->GetLabel() );  
1004       Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrack->GetLabel() );
1005                 
1006       TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1007       TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1008             
1009       fTreeVariablePosTransvMomentumMC = mcPosV0Dghter->Pt();
1010       fTreeVariableNegTransvMomentumMC = mcNegV0Dghter->Pt();
1011
1012       Int_t lPIDPositive = mcPosV0Dghter -> GetPdgCode();
1013       Int_t lPIDNegative = mcNegV0Dghter -> GetPdgCode();
1014
1015       fTreeVariablePIDPositive = lPIDPositive;
1016       fTreeVariablePIDNegative = lPIDNegative;
1017
1018       Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
1019       Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
1020
1021       if( lblMotherPosV0Dghter == lblMotherNegV0Dghter && lblMotherPosV0Dghter > -1 ){
1022          //either label is fine, they're equal at this stage
1023          TParticle* pThisV0 = lMCstack->Particle( lblMotherPosV0Dghter ); 
1024          //Set tree variables
1025          fTreeVariablePID   = pThisV0->GetPdgCode(); //PDG Code
1026          fTreeVariablePtMC  = pThisV0->Pt(); //Perfect Pt
1027          //Only Interested if it's a Lambda, AntiLambda or K0s 
1028          //Avoid the Junction Bug! PYTHIA has particles with Px=Py=Pz=E=0 occasionally, 
1029          //having particle code 88 (unrecognized by PDG), for documentation purposes.
1030          //Even ROOT's TParticle::Y() is not prepared to deal with that exception!
1031          //Note that TParticle::Pt() is immune (that would just return 0)...
1032          //Though granted that that should be extremely rare in this precise condition...
1033          if( TMath::Abs(fTreeVariablePID) == 3122 || fTreeVariablePID==310 ){
1034             fTreeVariableRapMC = pThisV0->Y(); //Perfect Y
1035          }
1036          fTreeVariableV0CreationRadius = pThisV0->R(); // Creation Radius
1037          if( lblMotherPosV0Dghter  < lNbMCPrimary ) fTreeVariableIndexStatus = 1; //looks primary
1038          if( lblMotherPosV0Dghter >= lNbMCPrimary ) fTreeVariableIndexStatus = 2; //looks secondary
1039          if( lMCstack->IsPhysicalPrimary       (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 1; //Is Primary!
1040          if( lMCstack->IsSecondaryFromWeakDecay(lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 2; //Weak Decay!
1041          if( lMCstack->IsSecondaryFromMaterial (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 3; //Material Int!
1042          
1043          //Now we try to acquire the V0 parent particle, if possible
1044          Int_t lblThisV0Parent = pThisV0->GetFirstMother();
1045          if ( lblThisV0Parent > -1 ){ //if it has a parent, get it and store specs
1046             TParticle* pThisV0Parent = lMCstack->Particle( lblThisV0Parent );
1047             fTreeVariablePIDMother   = pThisV0Parent->GetPdgCode(); //V0 Mother PDG
1048             fTreeVariablePtMother    = pThisV0Parent->Pt();         //V0 Mother Pt
1049             //Primary Status for the V0 Mother particle 
1050             if( lblThisV0Parent  < lNbMCPrimary ) fTreeVariableIndexStatusMother = 1; //looks primary
1051             if( lblThisV0Parent >= lNbMCPrimary ) fTreeVariableIndexStatusMother = 2; //looks secondary
1052             if( lMCstack->IsPhysicalPrimary       (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 1; //Is Primary!
1053             if( lMCstack->IsSecondaryFromWeakDecay(lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 2; //Weak Decay!
1054             if( lMCstack->IsSecondaryFromMaterial (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 3; //Material Int!
1055          }
1056       }
1057
1058       fTreeVariablePt = v0->Pt();
1059       fTreeVariableChi2V0 = lChi2V0; 
1060       fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
1061       fTreeVariableDcaV0Daughters = lDcaV0Daughters;
1062       fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle; 
1063       fTreeVariableV0Radius = lV0Radius;
1064       fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
1065       fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
1066       fTreeVariableInvMassK0s = lInvMassK0s;
1067       fTreeVariableInvMassLambda = lInvMassLambda;
1068       fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
1069       fTreeVariableRapK0Short = lRapK0Short;
1070
1071       fTreeVariableRapLambda = lRapLambda;
1072       fTreeVariableAlphaV0 = lAlphaV0;
1073       fTreeVariablePtArmV0 = lPtArmV0;
1074
1075       //Official means of acquiring N-sigmas 
1076       fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
1077       fTreeVariableNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
1078       fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
1079       fTreeVariableNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
1080
1081 //tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]
1082       fTreeVariableDistOverTotMom = TMath::Sqrt(
1083                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
1084                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
1085                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
1086                                         );
1087       fTreeVariableDistOverTotMom /= (lV0TotalMomentum + 1e-10); //avoid division by zero, to be sure
1088
1089       Double_t lMomentumPosTemp[3];
1090       pTrack->GetPxPyPz(lMomentumPosTemp);
1091       Double_t lPtPosTemporary = sqrt(pow(lMomentumPosTemp[0],2) + pow(lMomentumPosTemp[1],2));
1092
1093       Double_t lMomentumNegTemp[3];
1094       nTrack->GetPxPyPz(lMomentumNegTemp);
1095       Double_t lPtNegTemporary = sqrt(pow(lMomentumNegTemp[0],2) + pow(lMomentumNegTemp[1],2));
1096
1097       fTreeVariablePosTransvMomentum = lPtPosTemporary;
1098       fTreeVariableNegTransvMomentum = lPtNegTemporary;
1099
1100
1101 //------------------------------------------------
1102 // Fill Tree! 
1103 //------------------------------------------------
1104
1105       // The conditionals are meant to decrease excessive
1106       // memory usage! 
1107
1108       //Modified version: Keep only OnFlyStatus == 0
1109       //Keep only if included in a parametric InvMass Region 20 sigmas away from peak
1110
1111       //First Selection: Reject OnFly
1112       if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){
1113          //Second Selection: rough 20-sigma band, parametric. 
1114          //K0Short: Enough to parametrize peak broadening with linear function.    
1115          Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt; 
1116          Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
1117          //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
1118          //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
1119          Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt); 
1120          Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
1121          //Do Selection      
1122          if( (fTreeVariableInvMassLambda     < lUpperLimitLambda  && fTreeVariableInvMassLambda     > lLowerLimitLambda     ) || 
1123              (fTreeVariableInvMassAntiLambda < lUpperLimitLambda  && fTreeVariableInvMassAntiLambda > lLowerLimitLambda     ) || 
1124              (fTreeVariableInvMassK0s        < lUpperLimitK0Short && fTreeVariableInvMassK0s        > lLowerLimitK0Short    ) ){
1125              //Pre-selection in case this is AA...
1126              if( fkIsNuclear == kFALSE ) fTree->Fill();
1127              if( fkIsNuclear == kTRUE){ 
1128              //If this is a nuclear collision___________________
1129              // ... pre-filter with daughter eta selection only (not TPC)
1130                if ( TMath::Abs(fTreeVariableNegEta)<0.8 && TMath::Abs(fTreeVariablePosEta)<0.8 ) fTree->Fill();
1131              }//end nuclear_____________________________________
1132          }
1133       }
1134
1135 //------------------------------------------------
1136 // Fill tree over.
1137 //------------------------------------------------
1138
1139
1140    }// This is the end of the V0 loop
1141
1142    // Post output data.
1143    PostData(1, fListHistV0);
1144    PostData(2, fTree);
1145 }
1146
1147 //________________________________________________________________________
1148 void AliAnalysisTaskExtractPerformanceV0::Terminate(Option_t *)
1149 {
1150    // Draw result to the screen
1151    // Called once at the end of the query
1152
1153    TList *cRetrievedList = 0x0;
1154    cRetrievedList = (TList*)GetOutputData(1);
1155    if(!cRetrievedList){
1156       Printf("ERROR - AliAnalysisTaskExtractV0 : ouput data container list not available\n");
1157       return;
1158    }    
1159         
1160    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1161    if (!fHistV0MultiplicityForTrigEvt) {
1162       Printf("ERROR - AliAnalysisTaskExtractV0 : fHistV0MultiplicityForTrigEvt not available");
1163       return;
1164    }
1165   
1166    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractV0","V0 Multiplicity",10,10,510,510);
1167    canCheck->cd(1)->SetLogy();
1168
1169    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1170    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1171 }
1172
1173 //----------------------------------------------------------------------------
1174
1175 Double_t AliAnalysisTaskExtractPerformanceV0::MyRapidity(Double_t rE, Double_t rPz) const
1176 {
1177    // Local calculation for rapidity
1178    return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1179
1180
1181 //________________________________________________________________________
1182 void AliAnalysisTaskExtractPerformanceV0::CheckChargeV0(AliESDv0 *v0)
1183 {
1184    // This function checks charge of negative and positive daughter tracks. 
1185    // If incorrectly defined (onfly vertexer), swaps out. 
1186    if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
1187       //V0 daughter track swapping is required! Note: everything is swapped here... P->N, N->P
1188       Long_t lCorrectNidx = v0->GetPindex();
1189       Long_t lCorrectPidx = v0->GetNindex();
1190       Double32_t        lCorrectNmom[3];
1191       Double32_t        lCorrectPmom[3];
1192       v0->GetPPxPyPz( lCorrectNmom[0], lCorrectNmom[1], lCorrectNmom[2] );
1193       v0->GetNPxPyPz( lCorrectPmom[0], lCorrectPmom[1], lCorrectPmom[2] );
1194
1195       AliExternalTrackParam     lCorrectParamN(
1196         v0->GetParamP()->GetX() , 
1197         v0->GetParamP()->GetAlpha() , 
1198         v0->GetParamP()->GetParameter() , 
1199         v0->GetParamP()->GetCovariance() 
1200       );
1201       AliExternalTrackParam     lCorrectParamP(
1202         v0->GetParamN()->GetX() , 
1203         v0->GetParamN()->GetAlpha() , 
1204         v0->GetParamN()->GetParameter() , 
1205         v0->GetParamN()->GetCovariance() 
1206       );
1207       lCorrectParamN.SetMostProbablePt( v0->GetParamP()->GetMostProbablePt() );
1208       lCorrectParamP.SetMostProbablePt( v0->GetParamN()->GetMostProbablePt() );
1209
1210       //Get Variables___________________________________________________
1211       Double_t lDcaV0Daughters = v0 -> GetDcaV0Daughters();
1212       Double_t lCosPALocal     = v0 -> GetV0CosineOfPointingAngle(); 
1213       Bool_t lOnFlyStatusLocal = v0 -> GetOnFlyStatus();
1214
1215       //Create Replacement Object_______________________________________
1216       AliESDv0 *v0correct = new AliESDv0(lCorrectParamN,lCorrectNidx,lCorrectParamP,lCorrectPidx);
1217       v0correct->SetDcaV0Daughters          ( lDcaV0Daughters   );
1218       v0correct->SetV0CosineOfPointingAngle ( lCosPALocal       );
1219       v0correct->ChangeMassHypothesis       ( kK0Short          );
1220       v0correct->SetOnFlyStatus             ( lOnFlyStatusLocal );
1221
1222       //Reverse Cluster info..._________________________________________
1223       v0correct->SetClusters( v0->GetClusters( 1 ), v0->GetClusters ( 0 ) );
1224
1225       *v0 = *v0correct;
1226       //Proper cleanup..._______________________________________________
1227       v0correct->Delete();
1228       v0correct = 0x0;
1229
1230       //Just another cross-check and output_____________________________
1231       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ) {
1232         AliWarning("Found Swapped Charges, tried to correct but something FAILED!");
1233       }else{
1234         //AliWarning("Found Swapped Charges and fixed.");
1235       }
1236       //________________________________________________________________
1237    }else{
1238       //Don't touch it! ---
1239       //Printf("Ah, nice. Charges are already ordered...");
1240    }
1241    return;
1242
1243