]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/STRANGENESS/Cascades/AliAnalysisTaskStrangenessVsMultiplicityMC.cxx
another slight revision of calibration maps
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskStrangenessVsMultiplicityMC.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 // This task is meant to explore the possibility of using a VZERO amplitude
19 // based multiplicity estimator for proton-proton collisions. For this, two 
20 // main operation methods for this task are foreseen: 
21 // 
22 //  1) (under development) it should act as an auxiliary task and provide a 
23 //     calibrated estimator 
24 //
25 //  2) "Debug mode" which will also create a ROOT TTree object with event 
26 //     by event info potentially used for exploration / calibration. This 
27 //     includes the following info: 
28 //    
29 //      --- All VZERO Amplitudes (saved as Float_t) 
30 //      --- (optional) time for each channel
31 //      --- (optional) time width for each channel 
32 //      --- GetReferenceMultiplicity Estimator, |eta|<0.5 
33 //      --- GetReferenceMultiplicity Estimator, |eta|<0.8 
34 //      --- (if MC) True Multiplicity, |eta|<0.5
35 //      --- (if MC) True Multiplicity,  2.8 < eta < 5.8 (VZEROA region)
36 //      --- (if MC) True Multiplicity, -3.7 < eta <-1.7 (VZEROC region)
37 //      --- Run Number
38 //
39 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
40
41 class TTree;
42 class TParticle;
43 class TVector3;
44
45 //class AliMCEventHandler;
46 //class AliMCEvent;
47 //class AliStack;
48
49 class AliESDVertex;
50 class AliAODVertex;
51 class AliESDv0;
52 class AliAODv0;
53
54 #include <Riostream.h>
55 #include "TList.h"
56 #include "TH1.h"
57 #include "TH2.h"
58 #include "TH3.h"
59 #include "TFile.h"
60 #include "THnSparse.h"
61 #include "TVector3.h"
62 #include "TCanvas.h"
63 #include "TMath.h"
64 #include "TLegend.h"
65 //#include "AliLog.h"
66
67 #include "AliESDEvent.h"
68 #include "AliAODEvent.h"
69 #include "AliV0vertexer.h"
70 #include "AliCascadeVertexer.h"
71 #include "AliESDpid.h"
72 #include "AliESDtrack.h"
73 #include "AliESDtrackCuts.h"
74 #include "AliInputEventHandler.h"
75 #include "AliAnalysisManager.h"
76 #include "AliMCEventHandler.h"
77 #include "AliMCEvent.h"
78 #include "AliStack.h"
79 #include "AliCentrality.h"
80
81 #include "AliCFContainer.h"
82 #include "AliMultiplicity.h"
83 #include "AliAODMCParticle.h"
84 #include "AliESDcascade.h"
85 #include "AliAODcascade.h"
86 #include "AliESDUtils.h"
87 #include "AliGenEventHeader.h"
88 #include "AliAnalysisTaskSE.h"
89 #include "AliAnalysisUtils.h"
90 #include "AliAnalysisTaskStrangenessVsMultiplicityMC.h"
91
92 using std::cout;
93 using std::endl;
94
95 ClassImp(AliAnalysisTaskStrangenessVsMultiplicityMC)
96
97 AliAnalysisTaskStrangenessVsMultiplicityMC::AliAnalysisTaskStrangenessVsMultiplicityMC()
98   : AliAnalysisTaskSE(), fListHist(0), fTreeEvent(0), fTreeV0(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), 
99   fkSaveV0Tree      ( kFALSE ),
100   fkSaveCascadeTree ( kTRUE  ),
101   fkRunVertexers    ( kTRUE  ),
102   fkSkipEventSelection( kFALSE ),
103   //---> Variables for fTreeEvent
104   fAmplitude_V0A   (0),   
105   fAmplitude_V0C   (0),   
106   fAmplitude_V0AEq (0),   
107   fAmplitude_V0CEq (0),  
108   fCentrality_V0A(0), 
109   fCentrality_V0C(0), 
110   fCentrality_V0M(0), 
111   fCentrality_V0AEq(0), 
112   fCentrality_V0CEq(0), 
113   fCentrality_V0MEq(0), 
114   fRefMultEta5(0),
115   fRefMultEta8(0),
116   fTrueMultEta5(0),
117   fTrueMultEta8(0),
118   fTrueMultVZEROA(0),
119   fTrueMultVZEROC(0),
120   fRunNumber(0),
121     fEvSel_HasAtLeastSPDVertex(0),
122     fEvSel_VtxZCut(0),
123     fEvSel_IsNotPileup(0),
124     fEvSel_IsNotPileupMV(0),
125     fEvSel_IsNotPileupInMultBins(0),
126     fEvSel_HasVtxContributor(0),
127     fEvSel_Triggered(0),
128   //---> Variables for fTreeV0
129         fTreeVariableChi2V0(0),
130         fTreeVariableDcaV0Daughters(0),
131         fTreeVariableDcaV0ToPrimVertex(0),
132         fTreeVariableDcaPosToPrimVertex(0),
133         fTreeVariableDcaNegToPrimVertex(0),
134         fTreeVariableV0CosineOfPointingAngle(0),
135         fTreeVariableV0Radius(0),
136         fTreeVariablePt(0),
137         fTreeVariablePtMC(0),
138         fTreeVariableRapK0Short(0),
139         fTreeVariableRapLambda(0),
140         fTreeVariableRapMC(0),
141         fTreeVariableInvMassK0s(0),
142         fTreeVariableInvMassLambda(0),
143         fTreeVariableInvMassAntiLambda(0),
144         fTreeVariableAlphaV0(0),
145         fTreeVariablePtArmV0(0),
146         fTreeVariableNegEta(0),
147         fTreeVariablePosEta(0),
148
149         fTreeVariableNSigmasPosProton(0),
150         fTreeVariableNSigmasPosPion(0),
151         fTreeVariableNSigmasNegProton(0),
152         fTreeVariableNSigmasNegPion(0),
153
154         fTreeVariableDistOverTotMom(0),
155         fTreeVariableLeastNbrCrossedRows(0),
156         fTreeVariableLeastRatioCrossedRowsOverFindable(0),
157         
158   fTreeVariableCentV0A(0),
159   fTreeVariableCentV0C(0),
160   fTreeVariableCentV0M(0),
161   fTreeVariableCentV0AEq(0),
162   fTreeVariableCentV0CEq(0),
163   fTreeVariableCentV0MEq(0),
164   fTreeVariableAmpV0A(0),
165   fTreeVariableAmpV0C(0),
166   fTreeVariableAmpV0AEq(0),
167   fTreeVariableAmpV0CEq(0),
168   fTreeVariableRefMultEta8(0),
169   fTreeVariableRefMultEta5(0),
170   fTreeVariableRunNumber(0),
171
172   fTreeVariablePtMother(0),
173   fTreeVariablePID(0),
174   fTreeVariablePIDPositive(0),
175   fTreeVariablePIDNegative(0),
176   fTreeVariablePIDMother(0),
177   fTreeVariablePrimaryStatus(0),
178   fTreeVariablePrimaryStatusMother(0),
179   //---> Variables for fTreeCascade
180   fTreeCascVarCharge(0), 
181         fTreeCascVarMassAsXi(0),
182         fTreeCascVarMassAsOmega(0),
183         fTreeCascVarPt(0),
184         fTreeCascVarPtMC(0),
185         fTreeCascVarRapXi(0),
186         fTreeCascVarRapOmega(0),
187         fTreeCascVarRapMC(0),
188         fTreeCascVarNegEta(0),
189         fTreeCascVarPosEta(0),
190         fTreeCascVarBachEta(0),
191         fTreeCascVarDCACascDaughters(0),
192         fTreeCascVarDCABachToPrimVtx(0),
193         fTreeCascVarDCAV0Daughters(0),
194         fTreeCascVarDCAV0ToPrimVtx(0),
195         fTreeCascVarDCAPosToPrimVtx(0),
196         fTreeCascVarDCANegToPrimVtx(0),
197         fTreeCascVarCascCosPointingAngle(0),
198         fTreeCascVarCascRadius(0),
199         fTreeCascVarV0Mass(0),
200         fTreeCascVarV0CosPointingAngle(0),
201         fTreeCascVarV0CosPointingAngleSpecial(0),
202         fTreeCascVarV0Radius(0),
203   fTreeCascVarLeastNbrClusters(0),
204         fTreeCascVarDistOverTotMom(0),
205         fTreeCascVarNegNSigmaPion(0),
206         fTreeCascVarNegNSigmaProton(0),
207         fTreeCascVarPosNSigmaPion(0),
208         fTreeCascVarPosNSigmaProton(0),
209         fTreeCascVarBachNSigmaPion(0),
210         fTreeCascVarBachNSigmaKaon(0),
211         fTreeCascVarCentV0A(0),
212         fTreeCascVarCentV0C(0),
213         fTreeCascVarCentV0M(0),
214         fTreeCascVarCentV0AEq(0),
215         fTreeCascVarCentV0CEq(0),
216         fTreeCascVarCentV0MEq(0),
217         fTreeCascVarAmpV0A(0),
218         fTreeCascVarAmpV0C(0),
219         fTreeCascVarAmpV0AEq(0),
220         fTreeCascVarAmpV0CEq(0),
221   fTreeCascVarRefMultEta8(0),
222   fTreeCascVarRefMultEta5(0),
223   fTreeCascVarTrueMultEta5(0),
224   fTreeCascVarTrueMultEta8(0),
225   fTreeCascVarTrueMultVZEROA(0),
226   fTreeCascVarTrueMultVZEROC(0),
227   fTreeCascVarIsPhysicalPrimary(0), 
228   fTreeCascVarPID(0), 
229   fTreeCascVarRunNumber(0), 
230   //---> Histograms
231   fHistEventCounter(0), 
232   //---> MC Generated Histo (analysis level) 
233         fHistPt_GenK0Short(0),
234         fHistPt_GenLambda(0),
235         fHistPt_GenAntiLambda(0),
236         fHistPt_GenXiMinus(0),
237         fHistPt_GenXiPlus(0),
238         fHistPt_GenOmegaMinus(0),
239         fHistPt_GenOmegaPlus(0),
240
241   //VsRefMult
242         fHistPtVsRefMultEta5_GenXiMinus(0),
243         fHistPtVsRefMultEta5_GenXiPlus(0),
244         fHistPtVsRefMultEta5_GenOmegaMinus(0),
245         fHistPtVsRefMultEta5_GenOmegaPlus(0),
246         fHistPtVsRefMultEta8_GenXiMinus(0),
247         fHistPtVsRefMultEta8_GenXiPlus(0),
248         fHistPtVsRefMultEta8_GenOmegaMinus(0),
249         fHistPtVsRefMultEta8_GenOmegaPlus(0),
250
251   //VsCentralities
252         fHistPtVsCentV0A_GenXiMinus(0),
253         fHistPtVsCentV0A_GenXiPlus(0),
254         fHistPtVsCentV0A_GenOmegaMinus(0),
255         fHistPtVsCentV0A_GenOmegaPlus(0),
256         fHistPtVsCentV0C_GenXiMinus(0),
257         fHistPtVsCentV0C_GenXiPlus(0),
258         fHistPtVsCentV0C_GenOmegaMinus(0),
259         fHistPtVsCentV0C_GenOmegaPlus(0),
260         fHistPtVsCentV0M_GenXiMinus(0),
261         fHistPtVsCentV0M_GenXiPlus(0),
262         fHistPtVsCentV0M_GenOmegaMinus(0),
263         fHistPtVsCentV0M_GenOmegaPlus(0),
264
265   //Equalized
266         fHistPtVsCentV0AEq_GenXiMinus(0),
267         fHistPtVsCentV0AEq_GenXiPlus(0),
268         fHistPtVsCentV0AEq_GenOmegaMinus(0),
269         fHistPtVsCentV0AEq_GenOmegaPlus(0),
270         fHistPtVsCentV0CEq_GenXiMinus(0),
271         fHistPtVsCentV0CEq_GenXiPlus(0),
272         fHistPtVsCentV0CEq_GenOmegaMinus(0),
273         fHistPtVsCentV0CEq_GenOmegaPlus(0),
274         fHistPtVsCentV0MEq_GenXiMinus(0),
275         fHistPtVsCentV0MEq_GenXiPlus(0),
276         fHistPtVsCentV0MEq_GenOmegaMinus(0),
277         fHistPtVsCentV0MEq_GenOmegaPlus(0),
278
279   //VsAmp
280         fHistPtVsAmpV0A_GenXiMinus(0),
281         fHistPtVsAmpV0A_GenXiPlus(0),
282         fHistPtVsAmpV0A_GenOmegaMinus(0),
283         fHistPtVsAmpV0A_GenOmegaPlus(0),
284         fHistPtVsAmpV0C_GenXiMinus(0),
285         fHistPtVsAmpV0C_GenXiPlus(0),
286         fHistPtVsAmpV0C_GenOmegaMinus(0),
287         fHistPtVsAmpV0C_GenOmegaPlus(0),
288         fHistPtVsAmpV0M_GenXiMinus(0),
289         fHistPtVsAmpV0M_GenXiPlus(0),
290         fHistPtVsAmpV0M_GenOmegaMinus(0),
291         fHistPtVsAmpV0M_GenOmegaPlus(0),
292   //Equalized Amps
293         fHistPtVsAmpV0AEq_GenXiMinus(0),
294         fHistPtVsAmpV0AEq_GenXiPlus(0),
295         fHistPtVsAmpV0AEq_GenOmegaMinus(0),
296         fHistPtVsAmpV0AEq_GenOmegaPlus(0),
297         fHistPtVsAmpV0CEq_GenXiMinus(0),
298         fHistPtVsAmpV0CEq_GenXiPlus(0),
299         fHistPtVsAmpV0CEq_GenOmegaMinus(0),
300         fHistPtVsAmpV0CEq_GenOmegaPlus(0),
301         fHistPtVsAmpV0MEq_GenXiMinus(0),
302         fHistPtVsAmpV0MEq_GenXiPlus(0),
303         fHistPtVsAmpV0MEq_GenOmegaMinus(0),
304         fHistPtVsAmpV0MEq_GenOmegaPlus(0)  
305
306 //------------------------------------------------
307 // Tree Variables 
308 {
309
310 }
311
312 AliAnalysisTaskStrangenessVsMultiplicityMC::AliAnalysisTaskStrangenessVsMultiplicityMC(const char *name) 
313   : AliAnalysisTaskSE(name), fListHist(0), fTreeEvent(0), fTreeV0(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), 
314   fkSaveV0Tree      ( kFALSE ),
315   fkSaveCascadeTree ( kTRUE  ), 
316   fkRunVertexers    ( kTRUE  ),
317   fkSkipEventSelection( kFALSE ),
318   //---> Variables for fTreeEvent
319   fAmplitude_V0A (0),   
320   fAmplitude_V0C (0), 
321   fAmplitude_V0AEq (0),   
322   fAmplitude_V0CEq (0), 
323   fCentrality_V0A(0), 
324   fCentrality_V0C(0), 
325   fCentrality_V0M(0), 
326   fCentrality_V0AEq(0), 
327   fCentrality_V0CEq(0), 
328   fCentrality_V0MEq(0), 
329   fRefMultEta5(0),
330   fRefMultEta8(0),
331   fTrueMultEta5(0),
332   fTrueMultEta8(0),
333   fTrueMultVZEROA(0),
334   fTrueMultVZEROC(0),
335   fRunNumber(0),
336     fEvSel_HasAtLeastSPDVertex(0),
337     fEvSel_VtxZCut(0),
338     fEvSel_IsNotPileup(0),
339     fEvSel_IsNotPileupMV(0),
340     fEvSel_IsNotPileupInMultBins(0),
341     fEvSel_HasVtxContributor(0),
342     fEvSel_Triggered(0),
343   //---> Variables for fTreeV0
344         fTreeVariableChi2V0(0),
345         fTreeVariableDcaV0Daughters(0),
346         fTreeVariableDcaV0ToPrimVertex(0),
347         fTreeVariableDcaPosToPrimVertex(0),
348         fTreeVariableDcaNegToPrimVertex(0),
349         fTreeVariableV0CosineOfPointingAngle(0),
350         fTreeVariableV0Radius(0),
351         fTreeVariablePt(0),
352         fTreeVariablePtMC(0),
353         fTreeVariableRapK0Short(0),
354         fTreeVariableRapLambda(0),
355         fTreeVariableRapMC(0),
356         fTreeVariableInvMassK0s(0),
357         fTreeVariableInvMassLambda(0),
358         fTreeVariableInvMassAntiLambda(0),
359         fTreeVariableAlphaV0(0),
360         fTreeVariablePtArmV0(0),
361         fTreeVariableNegEta(0),
362         fTreeVariablePosEta(0),
363
364         fTreeVariableNSigmasPosProton(0),
365         fTreeVariableNSigmasPosPion(0),
366         fTreeVariableNSigmasNegProton(0),
367         fTreeVariableNSigmasNegPion(0),
368
369         fTreeVariableDistOverTotMom(0),
370         fTreeVariableLeastNbrCrossedRows(0),
371         fTreeVariableLeastRatioCrossedRowsOverFindable(0),
372         
373   fTreeVariableCentV0A(0),
374   fTreeVariableCentV0C(0),
375   fTreeVariableCentV0M(0),
376   fTreeVariableCentV0AEq(0),
377   fTreeVariableCentV0CEq(0),
378   fTreeVariableCentV0MEq(0),
379   fTreeVariableAmpV0A(0),
380   fTreeVariableAmpV0C(0),
381   fTreeVariableAmpV0AEq(0),
382   fTreeVariableAmpV0CEq(0),
383   fTreeVariableRefMultEta8(0),
384   fTreeVariableRefMultEta5(0),
385   fTreeVariableRunNumber(0),
386
387   fTreeVariablePtMother(0),
388   fTreeVariablePID(0),
389   fTreeVariablePIDPositive(0),
390   fTreeVariablePIDNegative(0),
391   fTreeVariablePIDMother(0),
392   fTreeVariablePrimaryStatus(0),
393   fTreeVariablePrimaryStatusMother(0),
394   //---> Variables for fTreeCascade
395   fTreeCascVarCharge(0), 
396         fTreeCascVarMassAsXi(0),
397         fTreeCascVarMassAsOmega(0),
398         fTreeCascVarPt(0),
399         fTreeCascVarPtMC(0),
400         fTreeCascVarRapXi(0),
401         fTreeCascVarRapOmega(0),
402         fTreeCascVarRapMC(0),
403         fTreeCascVarNegEta(0),
404         fTreeCascVarPosEta(0),
405         fTreeCascVarBachEta(0),
406         fTreeCascVarDCACascDaughters(0),
407         fTreeCascVarDCABachToPrimVtx(0),
408         fTreeCascVarDCAV0Daughters(0),
409         fTreeCascVarDCAV0ToPrimVtx(0),
410         fTreeCascVarDCAPosToPrimVtx(0),
411         fTreeCascVarDCANegToPrimVtx(0),
412         fTreeCascVarCascCosPointingAngle(0),
413         fTreeCascVarCascRadius(0),
414         fTreeCascVarV0Mass(0),
415         fTreeCascVarV0CosPointingAngle(0),
416         fTreeCascVarV0CosPointingAngleSpecial(0),
417         fTreeCascVarV0Radius(0),
418   fTreeCascVarLeastNbrClusters(0),
419         fTreeCascVarDistOverTotMom(0),
420         fTreeCascVarNegNSigmaPion(0),
421         fTreeCascVarNegNSigmaProton(0),
422         fTreeCascVarPosNSigmaPion(0),
423         fTreeCascVarPosNSigmaProton(0),
424         fTreeCascVarBachNSigmaPion(0),
425         fTreeCascVarBachNSigmaKaon(0),
426         fTreeCascVarCentV0A(0),
427         fTreeCascVarCentV0C(0),
428         fTreeCascVarCentV0M(0),
429         fTreeCascVarCentV0AEq(0),
430         fTreeCascVarCentV0CEq(0),
431         fTreeCascVarCentV0MEq(0),
432         fTreeCascVarAmpV0A(0),
433         fTreeCascVarAmpV0C(0),
434         fTreeCascVarAmpV0AEq(0),
435         fTreeCascVarAmpV0CEq(0),
436   fTreeCascVarRefMultEta8(0),
437   fTreeCascVarRefMultEta5(0),
438   fTreeCascVarTrueMultEta5(0),
439   fTreeCascVarTrueMultEta8(0),
440   fTreeCascVarTrueMultVZEROA(0),
441   fTreeCascVarTrueMultVZEROC(0),
442   fTreeCascVarIsPhysicalPrimary(0), 
443   fTreeCascVarPID(0), 
444   fTreeCascVarRunNumber(0), 
445   //---> Histograms
446   fHistEventCounter(0), 
447   //---> MC Generated Histo (analysis level) 
448         fHistPt_GenK0Short(0),
449         fHistPt_GenLambda(0),
450         fHistPt_GenAntiLambda(0),
451         fHistPt_GenXiMinus(0),
452         fHistPt_GenXiPlus(0),
453         fHistPt_GenOmegaMinus(0),
454         fHistPt_GenOmegaPlus(0),
455
456   //VsRefMult
457         fHistPtVsRefMultEta5_GenXiMinus(0),
458         fHistPtVsRefMultEta5_GenXiPlus(0),
459         fHistPtVsRefMultEta5_GenOmegaMinus(0),
460         fHistPtVsRefMultEta5_GenOmegaPlus(0),
461         fHistPtVsRefMultEta8_GenXiMinus(0),
462         fHistPtVsRefMultEta8_GenXiPlus(0),
463         fHistPtVsRefMultEta8_GenOmegaMinus(0),
464         fHistPtVsRefMultEta8_GenOmegaPlus(0),
465
466   //VsCentralities
467         fHistPtVsCentV0A_GenXiMinus(0),
468         fHistPtVsCentV0A_GenXiPlus(0),
469         fHistPtVsCentV0A_GenOmegaMinus(0),
470         fHistPtVsCentV0A_GenOmegaPlus(0),
471         fHistPtVsCentV0C_GenXiMinus(0),
472         fHistPtVsCentV0C_GenXiPlus(0),
473         fHistPtVsCentV0C_GenOmegaMinus(0),
474         fHistPtVsCentV0C_GenOmegaPlus(0),
475         fHistPtVsCentV0M_GenXiMinus(0),
476         fHistPtVsCentV0M_GenXiPlus(0),
477         fHistPtVsCentV0M_GenOmegaMinus(0),
478         fHistPtVsCentV0M_GenOmegaPlus(0),
479
480   //Equalized
481         fHistPtVsCentV0AEq_GenXiMinus(0),
482         fHistPtVsCentV0AEq_GenXiPlus(0),
483         fHistPtVsCentV0AEq_GenOmegaMinus(0),
484         fHistPtVsCentV0AEq_GenOmegaPlus(0),
485         fHistPtVsCentV0CEq_GenXiMinus(0),
486         fHistPtVsCentV0CEq_GenXiPlus(0),
487         fHistPtVsCentV0CEq_GenOmegaMinus(0),
488         fHistPtVsCentV0CEq_GenOmegaPlus(0),
489         fHistPtVsCentV0MEq_GenXiMinus(0),
490         fHistPtVsCentV0MEq_GenXiPlus(0),
491         fHistPtVsCentV0MEq_GenOmegaMinus(0),
492         fHistPtVsCentV0MEq_GenOmegaPlus(0),
493
494   //VsAmp
495         fHistPtVsAmpV0A_GenXiMinus(0),
496         fHistPtVsAmpV0A_GenXiPlus(0),
497         fHistPtVsAmpV0A_GenOmegaMinus(0),
498         fHistPtVsAmpV0A_GenOmegaPlus(0),
499         fHistPtVsAmpV0C_GenXiMinus(0),
500         fHistPtVsAmpV0C_GenXiPlus(0),
501         fHistPtVsAmpV0C_GenOmegaMinus(0),
502         fHistPtVsAmpV0C_GenOmegaPlus(0),
503         fHistPtVsAmpV0M_GenXiMinus(0),
504         fHistPtVsAmpV0M_GenXiPlus(0),
505         fHistPtVsAmpV0M_GenOmegaMinus(0),
506         fHistPtVsAmpV0M_GenOmegaPlus(0),
507   //Equalized Amps
508         fHistPtVsAmpV0AEq_GenXiMinus(0),
509         fHistPtVsAmpV0AEq_GenXiPlus(0),
510         fHistPtVsAmpV0AEq_GenOmegaMinus(0),
511         fHistPtVsAmpV0AEq_GenOmegaPlus(0),
512         fHistPtVsAmpV0CEq_GenXiMinus(0),
513         fHistPtVsAmpV0CEq_GenXiPlus(0),
514         fHistPtVsAmpV0CEq_GenOmegaMinus(0),
515         fHistPtVsAmpV0CEq_GenOmegaPlus(0),
516         fHistPtVsAmpV0MEq_GenXiMinus(0),
517         fHistPtVsAmpV0MEq_GenXiPlus(0),
518         fHistPtVsAmpV0MEq_GenOmegaMinus(0),
519         fHistPtVsAmpV0MEq_GenOmegaPlus(0)
520 {
521
522   //Re-vertex: Will only apply for cascade candidates
523
524   fV0VertexerSels[0] =  33.  ;  // max allowed chi2
525   fV0VertexerSels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
526   fV0VertexerSels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
527   fV0VertexerSels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
528   fV0VertexerSels[5] =   1.0 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
529   fV0VertexerSels[6] = 200.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
530         
531   fCascadeVertexerSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
532   fCascadeVertexerSels[1] =   0.05 ;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
533   fCascadeVertexerSels[2] =   0.010;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
534   fCascadeVertexerSels[3] =   0.03 ;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
535   fCascadeVertexerSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
536   fCascadeVertexerSels[5] =   0.95 ;  // min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
537   fCascadeVertexerSels[6] =   0.4  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
538   fCascadeVertexerSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
539         
540
541   DefineOutput(1, TList::Class()); // Event Counter Histo
542   DefineOutput(2, TTree::Class()); // Event Tree
543   DefineOutput(3, TTree::Class()); // V0 Tree
544   DefineOutput(4, TTree::Class()); // Cascade Tree
545 }
546
547
548 AliAnalysisTaskStrangenessVsMultiplicityMC::~AliAnalysisTaskStrangenessVsMultiplicityMC()
549 {
550 //------------------------------------------------
551 // DESTRUCTOR
552 //------------------------------------------------
553
554    if (fListHist){
555       delete fListHist;
556       fListHist = 0x0;
557    }
558    if (fTreeEvent){
559       delete fTreeEvent;
560       fTreeEvent = 0x0;
561    }
562    if (fTreeV0){
563       delete fTreeV0;
564       fTreeV0 = 0x0;
565    }
566    if (fTreeCascade){
567       delete fTreeCascade;
568       fTreeCascade = 0x0;
569    }
570 }
571
572 //________________________________________________________________________
573 void AliAnalysisTaskStrangenessVsMultiplicityMC::UserCreateOutputObjects()
574 {
575
576    OpenFile(2); 
577    // Called once
578
579 //------------------------------------------------
580
581    fTreeEvent = new TTree("fTreeEvent","Event");
582
583 //------------------------------------------------
584 // fTree Branch definitions - Event by Event info
585 //------------------------------------------------
586
587 //-----------BASIC-INFO---------------------------
588
589   //--- VZERO Data (Integrated)
590   fTreeEvent->Branch("fAmplitude_V0A",&fAmplitude_V0A,"fAmplitude_V0A/F");
591   fTreeEvent->Branch("fAmplitude_V0C",&fAmplitude_V0C,"fAmplitude_V0C/F");
592   fTreeEvent->Branch("fAmplitude_V0AEq",&fAmplitude_V0AEq,"fAmplitude_V0AEq/F");
593   fTreeEvent->Branch("fAmplitude_V0CEq",&fAmplitude_V0CEq,"fAmplitude_V0CEq/F");
594
595   //Info from AliCentrality (not necessarily 'centrality' per se) 
596   fTreeEvent->Branch("fCentrality_V0A",&fCentrality_V0A,"fCentrality_V0A/F");
597   fTreeEvent->Branch("fCentrality_V0C",&fCentrality_V0C,"fCentrality_V0C/F");
598   fTreeEvent->Branch("fCentrality_V0M",&fCentrality_V0M,"fCentrality_V0M/F");
599   fTreeEvent->Branch("fCentrality_V0AEq",&fCentrality_V0AEq,"fCentrality_V0AEq/F");
600   fTreeEvent->Branch("fCentrality_V0CEq",&fCentrality_V0CEq,"fCentrality_V0CEq/F");
601   fTreeEvent->Branch("fCentrality_V0MEq",&fCentrality_V0MEq,"fCentrality_V0MEq/F");
602   
603   //Official GetReferenceMultiplicity
604   fTreeEvent->Branch("fRefMultEta5",&fRefMultEta5,"fRefMultEta5/I");
605   fTreeEvent->Branch("fRefMultEta8",&fRefMultEta8,"fRefMultEta8/I");
606
607   fTreeEvent->Branch("fTrueMultEta5",&fTrueMultEta5,"fTrueMultEta5/I");
608   fTreeEvent->Branch("fTrueMultEta8",&fTrueMultEta8,"fTrueMultEta8/I");
609   fTreeEvent->Branch("fTrueMultVZEROA",&fTrueMultVZEROA,"fTrueMultVZEROA/I");
610   fTreeEvent->Branch("fTrueMultVZEROC",&fTrueMultVZEROC,"fTrueMultVZEROC/I");
611
612   //Run Number
613   fTreeEvent->Branch("fRunNumber", &fRunNumber, "fRunNumber/I");
614     
615     //Booleans for Event Selection
616     fTreeEvent->Branch("fEvSel_HasAtLeastSPDVertex", &fEvSel_HasAtLeastSPDVertex, "fEvSel_HasAtLeastSPDVertex/O");
617     fTreeEvent->Branch("fEvSel_VtxZCut", &fEvSel_VtxZCut, "fEvSel_VtxZCut/O");
618     fTreeEvent->Branch("fEvSel_IsNotPileup", &fEvSel_IsNotPileup, "fEvSel_IsNotPileup/O");
619     fTreeEvent->Branch("fEvSel_IsNotPileupMV", &fEvSel_IsNotPileupMV, "fEvSel_IsNotPileupMV/O");
620     fTreeEvent->Branch("fEvSel_IsNotPileupInMultBins", &fEvSel_IsNotPileupInMultBins, "fEvSel_IsNotPileupInMultBins/O");
621     fTreeEvent->Branch("fEvSel_HasVtxContributor", &fEvSel_HasVtxContributor, "fEvSel_HasVtxContributor/O");
622     fTreeEvent->Branch("fEvSel_Triggered", &fEvSel_Triggered, "fEvSel_Triggered/O");
623     
624
625   //Create Basic V0 Output Tree
626   fTreeV0 = new TTree( "fTreeV0", "V0 Candidates");
627
628 //------------------------------------------------
629 // fTreeV0 Branch definitions
630 //------------------------------------------------
631
632 //-----------BASIC-INFO---------------------------
633   fTreeV0->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"fTreeVariableChi2V0/F");
634   fTreeV0->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
635   fTreeV0->Branch("fTreeVariableDcaV0ToPrimVertex",&fTreeVariableDcaV0ToPrimVertex,"fTreeVariableDcaV0ToPrimVertex/F");
636   fTreeV0->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
637   fTreeV0->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
638   fTreeV0->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
639   fTreeV0->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
640   fTreeV0->Branch("fTreeVariablePtMC",&fTreeVariablePtMC,"fTreeVariablePtMC/F");
641   fTreeV0->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
642   fTreeV0->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
643   fTreeV0->Branch("fTreeVariableRapMC",&fTreeVariableRapMC,"fTreeVariableRapMC/F");
644   fTreeV0->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
645   fTreeV0->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
646   fTreeV0->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
647   fTreeV0->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
648   fTreeV0->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
649   fTreeV0->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
650   fTreeV0->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
651   fTreeV0->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
652   fTreeV0->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
653   fTreeV0->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
654   fTreeV0->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
655   fTreeV0->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
656   fTreeV0->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
657   fTreeV0->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
658   fTreeV0->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
659 //-----------MULTIPLICITY-INFO--------------------
660   fTreeV0->Branch("fTreeVariableCentV0A",&fTreeVariableCentV0A,"fTreeVariableCentV0A/F");
661   fTreeV0->Branch("fTreeVariableCentV0C",&fTreeVariableCentV0C,"fTreeVariableCentV0C/F");
662   fTreeV0->Branch("fTreeVariableCentV0M",&fTreeVariableCentV0M,"fTreeVariableCentV0M/F");
663   fTreeV0->Branch("fTreeVariableCentV0AEq",&fTreeVariableCentV0AEq,"fTreeVariableCentV0AEq/F");
664   fTreeV0->Branch("fTreeVariableCentV0CEq",&fTreeVariableCentV0CEq,"fTreeVariableCentV0CEq/F");
665   fTreeV0->Branch("fTreeVariableCentV0MEq",&fTreeVariableCentV0MEq,"fTreeVariableCentV0MEq/F");
666   fTreeV0->Branch("fTreeVariableAmpV0A",&fTreeVariableAmpV0A,"fTreeVariableAmpV0A/F");
667   fTreeV0->Branch("fTreeVariableAmpV0C",&fTreeVariableAmpV0C,"fTreeVariableAmpV0C/F");
668   fTreeV0->Branch("fTreeVariableAmpV0AEq",&fTreeVariableAmpV0AEq,"fTreeVariableAmpV0AEq/F");
669   fTreeV0->Branch("fTreeVariableAmpV0CEq",&fTreeVariableAmpV0CEq,"fTreeVariableAmpV0CEq/F");
670   fTreeV0->Branch("fTreeVariableRefMultEta8",&fTreeVariableRefMultEta8,"fTreeVariableRefMultEta8/I");
671   fTreeV0->Branch("fTreeVariableRefMultEta5",&fTreeVariableRefMultEta5,"fTreeVariableRefMultEta5/I");
672   fTreeV0->Branch("fTreeVariableRunNumber",&fTreeVariableRunNumber,"fTreeVariableRunNumber/I");
673 //-----------MC Exclusive info--------------------
674   fTreeV0->Branch("fTreeVariablePtMother",&fTreeVariablePtMother,"fTreeVariablePtMother/F");
675   fTreeV0->Branch("fTreeVariablePID",&fTreeVariablePID,"fTreeVariablePID/I");
676   fTreeV0->Branch("fTreeVariablePIDPositive",&fTreeVariablePIDPositive,"fTreeVariablePIDPositive/I");
677   fTreeV0->Branch("fTreeVariablePIDNegative",&fTreeVariablePIDNegative,"fTreeVariablePIDNegative/I");
678   fTreeV0->Branch("fTreeVariablePIDMother",&fTreeVariablePIDMother,"fTreeVariablePIDMother/I");
679   fTreeV0->Branch("fTreeVariablePrimaryStatus",&fTreeVariablePrimaryStatus,"fTreeVariablePrimaryStatus/I");
680   fTreeV0->Branch("fTreeVariablePrimaryStatusMother",&fTreeVariablePrimaryStatusMother,"fTreeVariablePrimaryStatusMother/I");
681 //------------------------------------------------
682
683   //Create Cascade output tree
684   fTreeCascade = new TTree("fTreeCascade","CascadeCandidates");
685
686 //------------------------------------------------
687 // fTreeCascade Branch definitions - Cascade Tree
688 //------------------------------------------------
689
690 //-----------BASIC-INFO---------------------------
691         fTreeCascade->Branch("fTreeCascVarCharge",&fTreeCascVarCharge,"fTreeCascVarCharge/I");  
692   fTreeCascade->Branch("fTreeCascVarMassAsXi",&fTreeCascVarMassAsXi,"fTreeCascVarMassAsXi/F");
693   fTreeCascade->Branch("fTreeCascVarMassAsOmega",&fTreeCascVarMassAsOmega,"fTreeCascVarMassAsOmega/F");
694   fTreeCascade->Branch("fTreeCascVarPt",&fTreeCascVarPt,"fTreeCascVarPt/F");
695   fTreeCascade->Branch("fTreeCascVarPtMC",&fTreeCascVarPtMC,"fTreeCascVarPtMC/F");
696   fTreeCascade->Branch("fTreeCascVarRapXi",&fTreeCascVarRapXi,"fTreeCascVarRapXi/F");
697   fTreeCascade->Branch("fTreeCascVarRapOmega",&fTreeCascVarRapOmega,"fTreeCascVarRapOmega/F");
698   fTreeCascade->Branch("fTreeCascVarRapMC",&fTreeCascVarRapMC,"fTreeCascVarRapMC/F");
699   fTreeCascade->Branch("fTreeCascVarNegEta",&fTreeCascVarNegEta,"fTreeCascVarNegEta/F");
700   fTreeCascade->Branch("fTreeCascVarPosEta",&fTreeCascVarPosEta,"fTreeCascVarPosEta/F");
701   fTreeCascade->Branch("fTreeCascVarBachEta",&fTreeCascVarBachEta,"fTreeCascVarBachEta/F");
702 //-----------INFO-FOR-CUTS------------------------
703   fTreeCascade->Branch("fTreeCascVarDCACascDaughters",&fTreeCascVarDCACascDaughters,"fTreeCascVarDCACascDaughters/F");
704   fTreeCascade->Branch("fTreeCascVarDCABachToPrimVtx",&fTreeCascVarDCABachToPrimVtx,"fTreeCascVarDCABachToPrimVtx/F");
705   fTreeCascade->Branch("fTreeCascVarDCAV0Daughters",&fTreeCascVarDCAV0Daughters,"fTreeCascVarDCAV0Daughters/F");
706   fTreeCascade->Branch("fTreeCascVarDCAV0ToPrimVtx",&fTreeCascVarDCAV0ToPrimVtx,"fTreeCascVarDCAV0ToPrimVtx/F");
707   fTreeCascade->Branch("fTreeCascVarDCAPosToPrimVtx",&fTreeCascVarDCAPosToPrimVtx,"fTreeCascVarDCAPosToPrimVtx/F");
708   fTreeCascade->Branch("fTreeCascVarDCANegToPrimVtx",&fTreeCascVarDCANegToPrimVtx,"fTreeCascVarDCANegToPrimVtx/F");
709   fTreeCascade->Branch("fTreeCascVarCascCosPointingAngle",&fTreeCascVarCascCosPointingAngle,"fTreeCascVarCascCosPointingAngle/F");
710   fTreeCascade->Branch("fTreeCascVarCascRadius",&fTreeCascVarCascRadius,"fTreeCascVarCascRadius/F");
711   fTreeCascade->Branch("fTreeCascVarV0Mass",&fTreeCascVarV0Mass,"fTreeCascVarV0Mass/F");
712   fTreeCascade->Branch("fTreeCascVarV0CosPointingAngle",&fTreeCascVarV0CosPointingAngle,"fTreeCascVarV0CosPointingAngle/F");
713   fTreeCascade->Branch("fTreeCascVarV0CosPointingAngleSpecial",&fTreeCascVarV0CosPointingAngleSpecial,"fTreeCascVarV0CosPointingAngleSpecial/F");
714   fTreeCascade->Branch("fTreeCascVarV0Radius",&fTreeCascVarV0Radius,"fTreeCascVarV0Radius/F");
715   fTreeCascade->Branch("fTreeCascVarLeastNbrClusters",&fTreeCascVarLeastNbrClusters,"fTreeCascVarLeastNbrClusters/I");
716 //-----------MULTIPLICITY-INFO--------------------
717   fTreeCascade->Branch("fTreeCascVarCentV0A",&fTreeCascVarCentV0A,"fTreeCascVarCentV0A/F");
718   fTreeCascade->Branch("fTreeCascVarCentV0C",&fTreeCascVarCentV0C,"fTreeCascVarCentV0C/F");
719   fTreeCascade->Branch("fTreeCascVarCentV0M",&fTreeCascVarCentV0M,"fTreeCascVarCentV0M/F");
720   fTreeCascade->Branch("fTreeCascVarCentV0AEq",&fTreeCascVarCentV0AEq,"fTreeCascVarCentV0AEq/F");
721   fTreeCascade->Branch("fTreeCascVarCentV0CEq",&fTreeCascVarCentV0CEq,"fTreeCascVarCentV0CEq/F");
722   fTreeCascade->Branch("fTreeCascVarCentV0MEq",&fTreeCascVarCentV0MEq,"fTreeCascVarCentV0MEq/F");
723   fTreeCascade->Branch("fTreeCascVarAmpV0A",&fTreeCascVarAmpV0A,"fTreeCascVarAmpV0A/F");
724   fTreeCascade->Branch("fTreeCascVarAmpV0C",&fTreeCascVarAmpV0C,"fTreeCascVarAmpV0C/F");
725   fTreeCascade->Branch("fTreeCascVarAmpV0AEq",&fTreeCascVarAmpV0AEq,"fTreeCascVarAmpV0AEq/F");
726   fTreeCascade->Branch("fTreeCascVarAmpV0CEq",&fTreeCascVarAmpV0CEq,"fTreeCascVarAmpV0CEq/F");
727   fTreeCascade->Branch("fTreeCascVarRefMultEta8",&fTreeCascVarRefMultEta8,"fTreeCascVarRefMultEta8/I");
728   fTreeCascade->Branch("fTreeCascVarRefMultEta5",&fTreeCascVarRefMultEta5,"fTreeCascVarRefMultEta5/I");
729   fTreeCascade->Branch("fTreeCascVarTrueMultEta5",&fTreeCascVarTrueMultEta5,"fTreeCascVarTrueMultEta5/I");
730   fTreeCascade->Branch("fTreeCascVarTrueMultEta8",&fTreeCascVarTrueMultEta8,"fTreeCascVarTrueMultEta8/I");
731   fTreeCascade->Branch("fTreeCascVarTrueMultVZEROA",&fTreeCascVarTrueMultVZEROA,"fTreeCascVarTrueMultVZEROA/I");
732   fTreeCascade->Branch("fTreeCascVarTrueMultVZEROC",&fTreeCascVarTrueMultVZEROC,"fTreeCascVarTrueMultVZEROC/I");
733   fTreeCascade->Branch("fTreeCascVarIsPhysicalPrimary",&fTreeCascVarIsPhysicalPrimary,"fTreeCascVarIsPhysicalPrimary/I");
734   fTreeCascade->Branch("fTreeCascVarPID",&fTreeCascVarPID,"fTreeCascVarPID/I");
735   fTreeCascade->Branch("fTreeCascVarRunNumber",&fTreeCascVarRunNumber,"fTreeCascVarRunNumber/I");
736 //-----------DECAY-LENGTH-INFO--------------------
737   fTreeCascade->Branch("fTreeCascVarDistOverTotMom",&fTreeCascVarDistOverTotMom,"fTreeCascVarDistOverTotMom/F");
738 //------------------------------------------------
739   fTreeCascade->Branch("fTreeCascVarNegNSigmaPion",&fTreeCascVarNegNSigmaPion,"fTreeCascVarNegNSigmaPion/F");
740   fTreeCascade->Branch("fTreeCascVarNegNSigmaProton",&fTreeCascVarNegNSigmaProton,"fTreeCascVarNegNSigmaProton/F");
741   fTreeCascade->Branch("fTreeCascVarPosNSigmaPion",&fTreeCascVarPosNSigmaPion,"fTreeCascVarPosNSigmaPion/F");
742   fTreeCascade->Branch("fTreeCascVarPosNSigmaProton",&fTreeCascVarPosNSigmaProton,"fTreeCascVarPosNSigmaProton/F");
743   fTreeCascade->Branch("fTreeCascVarBachNSigmaPion",&fTreeCascVarBachNSigmaPion,"fTreeCascVarBachNSigmaPion/F");
744   fTreeCascade->Branch("fTreeCascVarBachNSigmaKaon",&fTreeCascVarBachNSigmaKaon,"fTreeCascVarBachNSigmaKaon/F");
745
746 //------------------------------------------------
747 // Particle Identification Setup
748 //------------------------------------------------
749   
750    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
751    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
752    fPIDResponse = inputHandler->GetPIDResponse();
753
754   // Multiplicity
755   if(! fESDtrackCuts ){
756     fESDtrackCuts = new AliESDtrackCuts();
757   }
758
759 //------------------------------------------------
760 // V0 Multiplicity Histograms
761 //------------------------------------------------
762
763    // Create histograms
764    OpenFile(1);
765    fListHist = new TList();
766    fListHist->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
767
768    if(! fHistEventCounter ) {
769     //Histogram Output: Event-by-Event
770     fHistEventCounter = new TH1D( "fHistEventCounter", ";Evt. Sel. Step;Count",5,0,5); 
771     fHistEventCounter->GetXaxis()->SetBinLabel(1, "Processed");
772     fHistEventCounter->GetXaxis()->SetBinLabel(2, "Phys-Sel");  
773     fHistEventCounter->GetXaxis()->SetBinLabel(3, "Has Vtx");  
774     fHistEventCounter->GetXaxis()->SetBinLabel(4, "Vtx |z|<10cm");  
775     fHistEventCounter->GetXaxis()->SetBinLabel(5, "Isn't Pileup");
776     fListHist->Add(fHistEventCounter); 
777    }
778
779   //Histograms for Efficiency corrections... a bunch of them 
780   //1D Histograms - Fine if efficiency doesn't change vs mult (expected) 
781   //---> Always filled for |y|<0.5 
782   //V0s: basic Histos
783   if(! fHistPt_GenK0Short ) {
784     fHistPt_GenK0Short    = new TH1D( "fHistPt_GenK0Short",    "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenK0Short);    }
785   if(! fHistPt_GenLambda ) {
786     fHistPt_GenLambda     = new TH1D( "fHistPt_GenLambda",     "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenLambda);     }
787   if(! fHistPt_GenAntiLambda ) {
788     fHistPt_GenAntiLambda = new TH1D( "fHistPt_GenAntiLambda", "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenAntiLambda); }
789   //Cascades: basic Histos
790   if(! fHistPt_GenXiMinus ) {
791     fHistPt_GenXiMinus    = new TH1D( "fHistPt_GenXiMinus",    "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenXiMinus);    }
792   if(! fHistPt_GenXiPlus ) {
793      fHistPt_GenXiPlus     = new TH1D( "fHistPt_GenXiPlus",    "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenXiPlus);     }
794   if(! fHistPt_GenOmegaMinus ) {
795     fHistPt_GenOmegaMinus = new TH1D( "fHistPt_GenOmegaMinus", "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenOmegaMinus); }
796   if(! fHistPt_GenOmegaPlus ) {
797     fHistPt_GenOmegaPlus  = new TH1D( "fHistPt_GenOmegaPlus",  "Generated;p_{T} (GeV/c)",200,0,20);   fListHist->Add(fHistPt_GenOmegaPlus);  }
798   //2D Histos for vs Mult calculation 
799   if(! fHistPtVsRefMultEta5_GenXiMinus ) {
800     fHistPtVsRefMultEta5_GenXiMinus    = new TH2D( "fHistPtVsRefMultEta5_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta5_GenXiMinus);    }
801   if(! fHistPtVsRefMultEta5_GenXiPlus ) {
802     fHistPtVsRefMultEta5_GenXiPlus     = new TH2D( "fHistPtVsRefMultEta5_GenXiPlus",        "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta5_GenXiPlus);    }
803   if(! fHistPtVsRefMultEta5_GenOmegaMinus ) {
804     fHistPtVsRefMultEta5_GenOmegaMinus    = new TH2D( "fHistPtVsRefMultEta5_GenOmegaMinus", "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta5_GenOmegaMinus);    }
805   if(! fHistPtVsRefMultEta5_GenOmegaPlus ) {
806     fHistPtVsRefMultEta5_GenOmegaPlus     = new TH2D( "fHistPtVsRefMultEta5_GenOmegaPlus",  "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta5_GenOmegaPlus);    }
807   if(! fHistPtVsRefMultEta8_GenXiMinus ) {
808     fHistPtVsRefMultEta8_GenXiMinus    = new TH2D( "fHistPtVsRefMultEta8_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenXiMinus);    }
809   if(! fHistPtVsRefMultEta8_GenXiPlus ) {
810     fHistPtVsRefMultEta8_GenXiPlus     = new TH2D( "fHistPtVsRefMultEta8_GenXiPlus",        "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenXiPlus);    }
811   if(! fHistPtVsRefMultEta8_GenOmegaMinus ) {
812     fHistPtVsRefMultEta8_GenOmegaMinus    = new TH2D( "fHistPtVsRefMultEta8_GenOmegaMinus", "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenOmegaMinus);    }
813   if(! fHistPtVsRefMultEta8_GenOmegaPlus ) {
814     fHistPtVsRefMultEta8_GenOmegaPlus     = new TH2D( "fHistPtVsRefMultEta8_GenOmegaPlus",  "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   fListHist->Add(fHistPtVsRefMultEta8_GenOmegaPlus);    }
815
816   //Centralities: V0A, V0C, V0M, +Eq
817   if(! fHistPtVsCentV0A_GenXiMinus ) {
818     fHistPtVsCentV0A_GenXiMinus    = new TH2D( 
819     "fHistPtVsCentV0A_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
820     fListHist->Add(fHistPtVsCentV0A_GenXiMinus); }
821   if(! fHistPtVsCentV0A_GenXiPlus ) {
822     fHistPtVsCentV0A_GenXiPlus    = new TH2D( 
823     "fHistPtVsCentV0A_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
824     fListHist->Add(fHistPtVsCentV0A_GenXiPlus); }
825   if(! fHistPtVsCentV0A_GenOmegaMinus ) {
826     fHistPtVsCentV0A_GenOmegaMinus    = new TH2D( 
827     "fHistPtVsCentV0A_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
828     fListHist->Add(fHistPtVsCentV0A_GenOmegaMinus); }
829   if(! fHistPtVsCentV0A_GenOmegaPlus ) {
830     fHistPtVsCentV0A_GenOmegaPlus    = new TH2D( 
831     "fHistPtVsCentV0A_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
832     fListHist->Add(fHistPtVsCentV0A_GenOmegaPlus); }  
833
834   if(! fHistPtVsCentV0C_GenXiMinus ) {
835     fHistPtVsCentV0C_GenXiMinus    = new TH2D( 
836     "fHistPtVsCentV0C_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
837     fListHist->Add(fHistPtVsCentV0C_GenXiMinus); }
838   if(! fHistPtVsCentV0C_GenXiPlus ) {
839     fHistPtVsCentV0C_GenXiPlus    = new TH2D( 
840     "fHistPtVsCentV0C_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
841     fListHist->Add(fHistPtVsCentV0C_GenXiPlus); }
842   if(! fHistPtVsCentV0C_GenOmegaMinus ) {
843     fHistPtVsCentV0C_GenOmegaMinus    = new TH2D( 
844     "fHistPtVsCentV0C_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
845     fListHist->Add(fHistPtVsCentV0C_GenOmegaMinus); }
846   if(! fHistPtVsCentV0C_GenOmegaPlus ) {
847     fHistPtVsCentV0C_GenOmegaPlus    = new TH2D( 
848     "fHistPtVsCentV0C_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
849     fListHist->Add(fHistPtVsCentV0C_GenOmegaPlus); }  
850
851   if(! fHistPtVsCentV0M_GenXiMinus ) {
852     fHistPtVsCentV0M_GenXiMinus    = new TH2D( 
853     "fHistPtVsCentV0M_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
854     fListHist->Add(fHistPtVsCentV0M_GenXiMinus); }
855   if(! fHistPtVsCentV0M_GenXiPlus ) {
856     fHistPtVsCentV0M_GenXiPlus    = new TH2D( 
857     "fHistPtVsCentV0M_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
858     fListHist->Add(fHistPtVsCentV0M_GenXiPlus); }
859   if(! fHistPtVsCentV0M_GenOmegaMinus ) {
860     fHistPtVsCentV0M_GenOmegaMinus    = new TH2D( 
861     "fHistPtVsCentV0M_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
862     fListHist->Add(fHistPtVsCentV0M_GenOmegaMinus); }
863   if(! fHistPtVsCentV0M_GenOmegaPlus ) {
864     fHistPtVsCentV0M_GenOmegaPlus    = new TH2D( 
865     "fHistPtVsCentV0M_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
866     fListHist->Add(fHistPtVsCentV0M_GenOmegaPlus); }  
867
868   //Equalized
869   if(! fHistPtVsCentV0AEq_GenXiMinus ) {
870     fHistPtVsCentV0AEq_GenXiMinus    = new TH2D( 
871     "fHistPtVsCentV0AEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
872     fListHist->Add(fHistPtVsCentV0AEq_GenXiMinus); }
873   if(! fHistPtVsCentV0AEq_GenXiPlus ) {
874     fHistPtVsCentV0AEq_GenXiPlus    = new TH2D( 
875     "fHistPtVsCentV0AEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
876     fListHist->Add(fHistPtVsCentV0AEq_GenXiPlus); }
877   if(! fHistPtVsCentV0AEq_GenOmegaMinus ) {
878     fHistPtVsCentV0AEq_GenOmegaMinus    = new TH2D( 
879     "fHistPtVsCentV0AEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
880     fListHist->Add(fHistPtVsCentV0AEq_GenOmegaMinus); }
881   if(! fHistPtVsCentV0AEq_GenOmegaPlus ) {
882     fHistPtVsCentV0AEq_GenOmegaPlus    = new TH2D( 
883     "fHistPtVsCentV0AEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
884     fListHist->Add(fHistPtVsCentV0AEq_GenOmegaPlus); }  
885
886   if(! fHistPtVsCentV0CEq_GenXiMinus ) {
887     fHistPtVsCentV0CEq_GenXiMinus    = new TH2D( 
888     "fHistPtVsCentV0CEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
889     fListHist->Add(fHistPtVsCentV0CEq_GenXiMinus); }
890   if(! fHistPtVsCentV0CEq_GenXiPlus ) {
891     fHistPtVsCentV0CEq_GenXiPlus    = new TH2D( 
892     "fHistPtVsCentV0CEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
893     fListHist->Add(fHistPtVsCentV0CEq_GenXiPlus); }
894   if(! fHistPtVsCentV0CEq_GenOmegaMinus ) {
895     fHistPtVsCentV0CEq_GenOmegaMinus    = new TH2D( 
896     "fHistPtVsCentV0CEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
897     fListHist->Add(fHistPtVsCentV0CEq_GenOmegaMinus); }
898   if(! fHistPtVsCentV0CEq_GenOmegaPlus ) {
899     fHistPtVsCentV0CEq_GenOmegaPlus    = new TH2D( 
900     "fHistPtVsCentV0CEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
901     fListHist->Add(fHistPtVsCentV0CEq_GenOmegaPlus); }  
902
903   if(! fHistPtVsCentV0MEq_GenXiMinus ) {
904     fHistPtVsCentV0MEq_GenXiMinus    = new TH2D( 
905     "fHistPtVsCentV0MEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
906     fListHist->Add(fHistPtVsCentV0MEq_GenXiMinus); }
907   if(! fHistPtVsCentV0MEq_GenXiPlus ) {
908     fHistPtVsCentV0MEq_GenXiPlus    = new TH2D( 
909     "fHistPtVsCentV0MEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
910     fListHist->Add(fHistPtVsCentV0MEq_GenXiPlus); }
911   if(! fHistPtVsCentV0MEq_GenOmegaMinus ) {
912     fHistPtVsCentV0MEq_GenOmegaMinus    = new TH2D( 
913     "fHistPtVsCentV0MEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
914     fListHist->Add(fHistPtVsCentV0MEq_GenOmegaMinus); }
915   if(! fHistPtVsCentV0MEq_GenOmegaPlus ) {
916     fHistPtVsCentV0MEq_GenOmegaPlus    = new TH2D( 
917     "fHistPtVsCentV0MEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,100,0,100);   
918     fListHist->Add(fHistPtVsCentV0MEq_GenOmegaPlus); }  
919
920   //AMPLITUDES: V0A, V0C, V0M, +Eq
921   Double_t lMaxAmplitude = 2500; 
922   Long_t lAmplitudeBins = 10000;
923   if(! fHistPtVsAmpV0A_GenXiMinus ) {
924     fHistPtVsAmpV0A_GenXiMinus    = new TH2D( 
925     "fHistPtVsAmpV0A_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
926     fListHist->Add(fHistPtVsAmpV0A_GenXiMinus); }
927   if(! fHistPtVsAmpV0A_GenXiPlus ) {
928     fHistPtVsAmpV0A_GenXiPlus    = new TH2D( 
929     "fHistPtVsAmpV0A_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
930     fListHist->Add(fHistPtVsAmpV0A_GenXiPlus); }
931   if(! fHistPtVsAmpV0A_GenOmegaMinus ) {
932     fHistPtVsAmpV0A_GenOmegaMinus    = new TH2D( 
933     "fHistPtVsAmpV0A_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
934     fListHist->Add(fHistPtVsAmpV0A_GenOmegaMinus); }
935   if(! fHistPtVsAmpV0A_GenOmegaPlus ) {
936     fHistPtVsAmpV0A_GenOmegaPlus    = new TH2D( 
937     "fHistPtVsAmpV0A_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
938     fListHist->Add(fHistPtVsAmpV0A_GenOmegaPlus); }  
939
940   if(! fHistPtVsAmpV0C_GenXiMinus ) {
941     fHistPtVsAmpV0C_GenXiMinus    = new TH2D( 
942     "fHistPtVsAmpV0C_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
943     fListHist->Add(fHistPtVsAmpV0C_GenXiMinus); }
944   if(! fHistPtVsAmpV0C_GenXiPlus ) {
945     fHistPtVsAmpV0C_GenXiPlus    = new TH2D( 
946     "fHistPtVsAmpV0C_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
947     fListHist->Add(fHistPtVsAmpV0C_GenXiPlus); }
948   if(! fHistPtVsAmpV0C_GenOmegaMinus ) {
949     fHistPtVsAmpV0C_GenOmegaMinus    = new TH2D( 
950     "fHistPtVsAmpV0C_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
951     fListHist->Add(fHistPtVsAmpV0C_GenOmegaMinus); }
952   if(! fHistPtVsAmpV0C_GenOmegaPlus ) {
953     fHistPtVsAmpV0C_GenOmegaPlus    = new TH2D( 
954     "fHistPtVsAmpV0C_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
955     fListHist->Add(fHistPtVsAmpV0C_GenOmegaPlus); }  
956
957   if(! fHistPtVsAmpV0M_GenXiMinus ) {
958     fHistPtVsAmpV0M_GenXiMinus    = new TH2D( 
959     "fHistPtVsAmpV0M_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
960     fListHist->Add(fHistPtVsAmpV0M_GenXiMinus); }
961   if(! fHistPtVsAmpV0M_GenXiPlus ) {
962     fHistPtVsAmpV0M_GenXiPlus    = new TH2D( 
963     "fHistPtVsAmpV0M_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
964     fListHist->Add(fHistPtVsAmpV0M_GenXiPlus); }
965   if(! fHistPtVsAmpV0M_GenOmegaMinus ) {
966     fHistPtVsAmpV0M_GenOmegaMinus    = new TH2D( 
967     "fHistPtVsAmpV0M_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
968     fListHist->Add(fHistPtVsAmpV0M_GenOmegaMinus); }
969   if(! fHistPtVsAmpV0M_GenOmegaPlus ) {
970     fHistPtVsAmpV0M_GenOmegaPlus    = new TH2D( 
971     "fHistPtVsAmpV0M_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
972     fListHist->Add(fHistPtVsAmpV0M_GenOmegaPlus); }  
973
974   //Equalized
975   if(! fHistPtVsAmpV0AEq_GenXiMinus ) {
976     fHistPtVsAmpV0AEq_GenXiMinus    = new TH2D( 
977     "fHistPtVsAmpV0AEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
978     fListHist->Add(fHistPtVsAmpV0AEq_GenXiMinus); }
979   if(! fHistPtVsAmpV0AEq_GenXiPlus ) {
980     fHistPtVsAmpV0AEq_GenXiPlus    = new TH2D( 
981     "fHistPtVsAmpV0AEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
982     fListHist->Add(fHistPtVsAmpV0AEq_GenXiPlus); }
983   if(! fHistPtVsAmpV0AEq_GenOmegaMinus ) {
984     fHistPtVsAmpV0AEq_GenOmegaMinus    = new TH2D( 
985     "fHistPtVsAmpV0AEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
986     fListHist->Add(fHistPtVsAmpV0AEq_GenOmegaMinus); }
987   if(! fHistPtVsAmpV0AEq_GenOmegaPlus ) {
988     fHistPtVsAmpV0AEq_GenOmegaPlus    = new TH2D( 
989     "fHistPtVsAmpV0AEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
990     fListHist->Add(fHistPtVsAmpV0AEq_GenOmegaPlus); }  
991
992   if(! fHistPtVsAmpV0CEq_GenXiMinus ) {
993     fHistPtVsAmpV0CEq_GenXiMinus    = new TH2D( 
994     "fHistPtVsAmpV0CEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
995     fListHist->Add(fHistPtVsAmpV0CEq_GenXiMinus); }
996   if(! fHistPtVsAmpV0CEq_GenXiPlus ) {
997     fHistPtVsAmpV0CEq_GenXiPlus    = new TH2D( 
998     "fHistPtVsAmpV0CEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
999     fListHist->Add(fHistPtVsAmpV0CEq_GenXiPlus); }
1000   if(! fHistPtVsAmpV0CEq_GenOmegaMinus ) {
1001     fHistPtVsAmpV0CEq_GenOmegaMinus    = new TH2D( 
1002     "fHistPtVsAmpV0CEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1003     fListHist->Add(fHistPtVsAmpV0CEq_GenOmegaMinus); }
1004   if(! fHistPtVsAmpV0CEq_GenOmegaPlus ) {
1005     fHistPtVsAmpV0CEq_GenOmegaPlus    = new TH2D( 
1006     "fHistPtVsAmpV0CEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1007     fListHist->Add(fHistPtVsAmpV0CEq_GenOmegaPlus); }  
1008
1009   if(! fHistPtVsAmpV0MEq_GenXiMinus ) {
1010     fHistPtVsAmpV0MEq_GenXiMinus    = new TH2D( 
1011     "fHistPtVsAmpV0MEq_GenXiMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1012     fListHist->Add(fHistPtVsAmpV0MEq_GenXiMinus); }
1013   if(! fHistPtVsAmpV0MEq_GenXiPlus ) {
1014     fHistPtVsAmpV0MEq_GenXiPlus    = new TH2D( 
1015     "fHistPtVsAmpV0MEq_GenXiPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1016     fListHist->Add(fHistPtVsAmpV0MEq_GenXiPlus); }
1017   if(! fHistPtVsAmpV0MEq_GenOmegaMinus ) {
1018     fHistPtVsAmpV0MEq_GenOmegaMinus    = new TH2D( 
1019     "fHistPtVsAmpV0MEq_GenOmegaMinus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1020     fListHist->Add(fHistPtVsAmpV0MEq_GenOmegaMinus); }
1021   if(! fHistPtVsAmpV0MEq_GenOmegaPlus ) {
1022     fHistPtVsAmpV0MEq_GenOmegaPlus    = new TH2D( 
1023     "fHistPtVsAmpV0MEq_GenOmegaPlus",       "Generated;p_{T} (GeV/c); Mult",200,0,20,lAmplitudeBins,0,lMaxAmplitude);   
1024     fListHist->Add(fHistPtVsAmpV0MEq_GenOmegaPlus); }  
1025
1026    //List of Histograms: Normal
1027    PostData(1, fListHist);
1028
1029    //TTree Object: Saved to base directory. Should cache to disk while saving. 
1030    //(Important to avoid excessive memory usage, particularly when merging)
1031    PostData(2, fTreeEvent);
1032    PostData(3, fTreeV0);
1033    PostData(4, fTreeCascade);
1034
1035 }// end UserCreateOutputObjects
1036
1037
1038 //________________________________________________________________________
1039 void AliAnalysisTaskStrangenessVsMultiplicityMC::UserExec(Option_t *) 
1040 {
1041   // Main loop
1042   // Called for each event
1043
1044    AliESDEvent *lESDevent = 0x0;
1045    AliMCEvent  *lMCevent  = 0x0; 
1046    AliStack    *lMCstack  = 0x0; 
1047   
1048     //Zero all booleans, etc
1049     fEvSel_HasAtLeastSPDVertex    = kFALSE;
1050     fEvSel_VtxZCut                = kFALSE;
1051     fEvSel_IsNotPileup            = kFALSE;
1052     fEvSel_IsNotPileupInMultBins  = kFALSE;
1053     fEvSel_HasVtxContributor      = kFALSE;
1054     fEvSel_Triggered              = kFALSE;
1055   // Connect to the InputEvent   
1056   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
1057
1058    // Appropriate for ESD analysis! 
1059       
1060    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1061    if (!lESDevent) {
1062       AliWarning("ERROR: lESDevent not available \n");
1063       return;
1064    }
1065
1066   //Get VZERO Information for multiplicity later
1067   AliVVZERO* esdV0 = lESDevent->GetVZEROData();
1068   if (!esdV0) {
1069     AliError("AliVVZERO not available");
1070     return;
1071   }
1072         
1073   lMCevent = MCEvent();
1074   if (!lMCevent) {
1075     Printf("ERROR: Could not retrieve MC event \n");
1076     cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;   
1077     return;
1078   }
1079
1080   lMCstack = lMCevent->Stack();
1081   if (!lMCstack) {
1082     Printf("ERROR: Could not retrieve MC stack \n");
1083     cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1084     return;
1085   }
1086
1087    fRunNumber = lESDevent->GetRunNumber();
1088
1089   Double_t lMagneticField = -10; 
1090   lMagneticField = lESDevent->GetMagneticField( );
1091
1092 //------------------------------------------------
1093 // Variable Definition
1094 //------------------------------------------------
1095
1096 //------------------------------------------------
1097 // Physics Selection
1098 //------------------------------------------------
1099   
1100   fHistEventCounter->Fill(0.5); 
1101
1102   UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1103   Bool_t isSelected = 0;
1104   isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1105   fEvSel_Triggered = isSelected;
1106   
1107   //Standard Min-Bias Selection
1108   if ( (! isSelected) && (! fkSkipEventSelection ) ) {
1109     PostData(1, fListHist);
1110     PostData(2, fTreeEvent);
1111     PostData(3, fTreeV0);
1112     PostData(4, fTreeCascade);
1113     return;
1114   }
1115
1116   fHistEventCounter->Fill(1.5);
1117  
1118   //------------------------------------------------
1119   // Primary Vertex Requirements Section:
1120   //  ---> pp: has vertex, |z|<10cm
1121   //------------------------------------------------
1122   
1123   //classical Proton-proton like selection 
1124   const AliESDVertex *lPrimaryBestESDVtx     = lESDevent->GetPrimaryVertex();   
1125   const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
1126   const AliESDVertex *lPrimarySPDVtx         = lESDevent->GetPrimaryVertexSPD();
1127
1128   Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
1129   lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
1130
1131   //Only accept if Tracking or SPD vertex is fine 
1132   if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtx->GetStatus() && !fkSkipEventSelection ){
1133     AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1134     PostData(1, fListHist); 
1135     PostData(2, fTreeEvent);
1136     PostData(3, fTreeV0);
1137     PostData(4, fTreeCascade);
1138     return;
1139   }
1140     
1141     if(! (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtx->GetStatus()) ){
1142         //Passed selection!
1143         fEvSel_HasAtLeastSPDVertex = kTRUE;
1144     }
1145     
1146   //Has SPD or Tracking Vertex
1147   fHistEventCounter -> Fill(2.5); 
1148
1149   //Always do Primary Vertex Selection 
1150   if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 && !fkSkipEventSelection ) {
1151     AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1152     PostData(1, fListHist); 
1153     PostData(2, fTreeEvent);
1154     PostData(3, fTreeV0);
1155     PostData(4, fTreeCascade);
1156     return;
1157   }
1158     
1159     if(TMath::Abs(lBestPrimaryVtxPos[2]) <= 10.0 ){
1160         //Passed selection!
1161         fEvSel_VtxZCut = kTRUE;
1162     }
1163
1164   //Fill Event selected counter
1165   fHistEventCounter -> Fill(3.5);
1166
1167   //------------------------------------------------
1168   // Check if this isn't pileup
1169   //------------------------------------------------
1170
1171   if(lESDevent->IsPileupFromSPDInMultBins() && !fkSkipEventSelection ){
1172     // minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  
1173     //-> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1174     AliWarning("Pb / Event tagged as pile-up by SPD... return !"); 
1175     PostData(1, fListHist); 
1176     PostData(2, fTreeEvent);
1177     PostData(3, fTreeV0);
1178     PostData(4, fTreeCascade);
1179     return; 
1180   }
1181     
1182     if( !lESDevent->IsPileupFromSPD()           ) fEvSel_IsNotPileup           = kTRUE;
1183     if( !lESDevent->IsPileupFromSPDInMultBins() ) fEvSel_IsNotPileupInMultBins = kTRUE;
1184     
1185     //First implementation of pileup from multi-vertexer (simple use of analysis utils)
1186     //if ( !fUtils->IsPileUpMV( lESDevent ) ) fEvSel_IsNotPileupMV = kTRUE;
1187     fEvSel_IsNotPileupMV = kFALSE ; //dummy
1188     
1189   //Fill Event isn't pileup counter
1190   fHistEventCounter -> Fill(4.5);
1191
1192 //------------------------------------------------
1193 // Multiplicity Information Acquistion
1194 //------------------------------------------------
1195
1196   //Monte Carlo Level information ! 
1197   //--------- GENERATED NUMBER OF CHARGED PARTICLES
1198   // ---> Variable Definition
1199
1200   Long_t lNchEta5   = 0; 
1201   Long_t lNchEta8   = 0; 
1202   Long_t lNchVZEROA = 0; 
1203   Long_t lNchVZEROC = 0; 
1204
1205   //----- Loop on Stack ----------------------------------------------------------------
1206   for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++) 
1207   {// This is the begining of the loop on tracks
1208       TParticle* particleOne = lMCstack->Particle(iCurrentLabelStack);
1209       if(!particleOne) continue;
1210       if(!particleOne->GetPDG()) continue;
1211       Double_t lThisCharge = particleOne->GetPDG()->Charge()/3.;
1212       if(TMath::Abs(lThisCharge)<0.001) continue;
1213       if(! (lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) ) continue;
1214      
1215       //Double_t gpt = particleOne -> Pt();
1216       Double_t geta = particleOne -> Eta(); 
1217
1218       if( TMath::Abs(geta) < 0.5 ) lNchEta5++; 
1219       if( TMath::Abs(geta) < 0.8 ) lNchEta8++; 
1220       if( 2.8 < geta && geta < 5.1 ) lNchVZEROA++; 
1221       if(-3.7 < geta && geta <-1.7 ) lNchVZEROC++; 
1222   }//End of loop on tracks
1223
1224   //Attribution 
1225   fTrueMultEta5 = lNchEta5; 
1226   fTrueMultEta8 = lNchEta8; 
1227   fTrueMultVZEROA = lNchVZEROA; 
1228   fTrueMultVZEROC = lNchVZEROC; 
1229   //----- End Loop on Stack ------------------------------------------------------------
1230
1231   //Standard GetReferenceMultiplicity Estimator (0.5 and 0.8)
1232   fRefMultEta5 = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1233   fRefMultEta8 = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.8);
1234
1235   // VZERO PART
1236   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
1237   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
1238   Float_t  multV0AEq  = 0;          //  multiplicity from V0 reco side A
1239   Float_t  multV0CEq  = 0;          //  multiplicity from V0 reco side C
1240   Float_t  multV0ACorr  = 0;            //  multiplicity from V0 reco side A
1241   Float_t  multV0CCorr  = 0;            //  multiplicity from V0 reco side C
1242
1243   //Non-Equalized Signal: copy of multV0ACorr and multV0CCorr from AliCentralitySelectionTask
1244   //Getters for uncorrected multiplicity  
1245   multV0A=esdV0->GetMTotV0A();
1246   multV0C=esdV0->GetMTotV0C();
1247
1248   //Get Z vertex position of SPD vertex (why not Tracking if available?) 
1249   Float_t zvtx = lPrimarySPDVtx->GetZ(); 
1250
1251   //Acquire Corrected multV0A 
1252   multV0ACorr = AliESDUtils::GetCorrV0A(multV0A,zvtx);    
1253   multV0CCorr = AliESDUtils::GetCorrV0C(multV0C,zvtx);   
1254     
1255   //Copy to Event Tree for extra information 
1256   fAmplitude_V0A = multV0ACorr; 
1257   fAmplitude_V0C = multV0CCorr; 
1258
1259   // Equalized signals // From AliCentralitySelectionTask
1260   for(Int_t iCh = 4; iCh < 7; ++iCh) {
1261     Double_t mult = lESDevent->GetVZEROEqMultiplicity(iCh);
1262     multV0AEq += mult;
1263   }
1264   for(Int_t iCh = 0; iCh < 3; ++iCh) {
1265     Double_t mult = lESDevent->GetVZEROEqMultiplicity(iCh);
1266     multV0CEq += mult;
1267   }
1268   fAmplitude_V0AEq = multV0AEq; 
1269   fAmplitude_V0CEq = multV0CEq; 
1270
1271   fCentrality_V0A   = -100; 
1272   fCentrality_V0C   = -100; 
1273   fCentrality_V0M   = -100; 
1274   fCentrality_V0AEq = -100; 
1275   fCentrality_V0CEq = -100; 
1276   fCentrality_V0MEq = -100; 
1277
1278   //AliCentrality... Check if working? 
1279   AliCentrality* centrality;
1280   centrality = lESDevent->GetCentrality();
1281   if ( !(centrality->GetQuality()>1) ){ 
1282     fCentrality_V0A   = centrality->GetCentralityPercentile( "V0A"   ); 
1283     fCentrality_V0C   = centrality->GetCentralityPercentile( "V0C"   ); 
1284     fCentrality_V0M   = centrality->GetCentralityPercentile( "V0M"   ); 
1285     fCentrality_V0AEq = centrality->GetCentralityPercentile( "V0AEq" ); 
1286     fCentrality_V0CEq = centrality->GetCentralityPercentile( "V0CEq" ); 
1287     fCentrality_V0MEq = centrality->GetCentralityPercentile( "V0MEq" ); 
1288   }
1289   
1290   //Event-level fill 
1291   fTreeEvent->Fill() ;
1292   
1293 //------------------------------------------------
1294
1295 //------------------------------------------------
1296 // Fill Efficiency Denominators, please 
1297 //------------------------------------------------
1298
1299    Int_t    lThisPDG  = 0;
1300    Double_t lThisRap  = 0;
1301    Double_t lThisPt   = 0;
1302
1303 //----- Loop on Generated CASCADES ---------------
1304    for (Int_t ilab = 0;  ilab < (lMCstack->GetNtrack()); ilab++) 
1305    {// This is the begining of the loop on tracks
1306       
1307       TParticle* lPart = 0x0; 
1308       lPart = lMCstack->Particle( ilab );
1309       if(!lPart){
1310          Printf("Generated loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", ilab );
1311          continue;
1312       }
1313
1314       lThisPDG = lPart->GetPdgCode();         
1315
1316       if ( (TMath::Abs(lThisPDG) == 3312) || (TMath::Abs(lThisPDG) == 3334) || (TMath::Abs(lThisPDG) == 3122) || lThisPDG == 310 ) 
1317       {
1318         lThisRap   = MyRapidity(lPart->Energy(),lPart->Pz());
1319         lThisPt    = lPart->Pt();
1320
1321         //Use Physical Primaries only for filling These Histos
1322         if ( lMCstack->IsPhysicalPrimary(ilab)!=kTRUE ) continue;
1323
1324         if( lThisPDG ==   310 && TMath::Abs(lThisRap) < 0.5 ) fHistPt_GenK0Short       -> Fill ( lThisPt ); 
1325         if( lThisPDG ==  3122 && TMath::Abs(lThisRap) < 0.5 ) fHistPt_GenLambda     -> Fill ( lThisPt ); 
1326         if( lThisPDG == -3122 && TMath::Abs(lThisRap) < 0.5 ) fHistPt_GenAntiLambda -> Fill ( lThisPt ); 
1327
1328         if( lThisPDG ==  3312 && TMath::Abs(lThisRap) < 0.5 ){
1329           fHistPt_GenXiMinus                -> Fill (lThisPt);          
1330           fHistPtVsRefMultEta5_GenXiMinus   -> Fill (lThisPt, fRefMultEta5); 
1331           fHistPtVsRefMultEta8_GenXiMinus   -> Fill (lThisPt, fRefMultEta8);
1332           //Centralities 
1333           fHistPtVsCentV0A_GenXiMinus       -> Fill (lThisPt, fCentrality_V0A);  
1334           fHistPtVsCentV0C_GenXiMinus       -> Fill (lThisPt, fCentrality_V0C);  
1335           fHistPtVsCentV0M_GenXiMinus       -> Fill (lThisPt, fCentrality_V0M);  
1336           fHistPtVsCentV0AEq_GenXiMinus       -> Fill (lThisPt, fCentrality_V0AEq);  
1337           fHistPtVsCentV0CEq_GenXiMinus       -> Fill (lThisPt, fCentrality_V0CEq);  
1338           fHistPtVsCentV0MEq_GenXiMinus       -> Fill (lThisPt, fCentrality_V0MEq);  
1339           //Amplitudes 
1340           fHistPtVsAmpV0A_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0A);  
1341           fHistPtVsAmpV0C_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0C);  
1342           fHistPtVsAmpV0M_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0A + fAmplitude_V0C);  
1343           fHistPtVsAmpV0AEq_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0AEq);  
1344           fHistPtVsAmpV0CEq_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0CEq);  
1345           fHistPtVsAmpV0MEq_GenXiMinus       -> Fill (lThisPt, fAmplitude_V0AEq + fAmplitude_V0CEq);  
1346         }
1347         if( lThisPDG == -3312 && TMath::Abs(lThisRap) < 0.5 ){
1348           fHistPt_GenXiPlus                -> Fill (lThisPt);          
1349           fHistPtVsRefMultEta5_GenXiPlus   -> Fill (lThisPt, fRefMultEta5); 
1350           fHistPtVsRefMultEta8_GenXiPlus   -> Fill (lThisPt, fRefMultEta8);
1351           //Centralities 
1352           fHistPtVsCentV0A_GenXiPlus       -> Fill (lThisPt, fCentrality_V0A);  
1353           fHistPtVsCentV0C_GenXiPlus       -> Fill (lThisPt, fCentrality_V0C);  
1354           fHistPtVsCentV0M_GenXiPlus       -> Fill (lThisPt, fCentrality_V0M);  
1355           fHistPtVsCentV0AEq_GenXiPlus       -> Fill (lThisPt, fCentrality_V0AEq);  
1356           fHistPtVsCentV0CEq_GenXiPlus       -> Fill (lThisPt, fCentrality_V0CEq);  
1357           fHistPtVsCentV0MEq_GenXiPlus       -> Fill (lThisPt, fCentrality_V0MEq);  
1358           //Amplitudes 
1359           fHistPtVsAmpV0A_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0A);  
1360           fHistPtVsAmpV0C_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0C);  
1361           fHistPtVsAmpV0M_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0A + fAmplitude_V0C);  
1362           fHistPtVsAmpV0AEq_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0AEq);  
1363           fHistPtVsAmpV0CEq_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0CEq);  
1364           fHistPtVsAmpV0MEq_GenXiPlus       -> Fill (lThisPt, fAmplitude_V0AEq + fAmplitude_V0CEq);  
1365         }
1366         if( lThisPDG ==  3334 && TMath::Abs(lThisRap) < 0.5 ){
1367           fHistPt_GenOmegaMinus                -> Fill (lThisPt);          
1368           fHistPtVsRefMultEta5_GenOmegaMinus   -> Fill (lThisPt, fRefMultEta5); 
1369           fHistPtVsRefMultEta8_GenOmegaMinus   -> Fill (lThisPt, fRefMultEta8);
1370           //Centralities 
1371           fHistPtVsCentV0A_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0A);  
1372           fHistPtVsCentV0C_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0C);  
1373           fHistPtVsCentV0M_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0M);  
1374           fHistPtVsCentV0AEq_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0AEq);  
1375           fHistPtVsCentV0CEq_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0CEq);  
1376           fHistPtVsCentV0MEq_GenOmegaMinus       -> Fill (lThisPt, fCentrality_V0MEq);  
1377           //Amplitudes 
1378           fHistPtVsAmpV0A_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0A);  
1379           fHistPtVsAmpV0C_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0C);  
1380           fHistPtVsAmpV0M_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0A + fAmplitude_V0C);  
1381           fHistPtVsAmpV0AEq_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0AEq);  
1382           fHistPtVsAmpV0CEq_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0CEq);  
1383           fHistPtVsAmpV0MEq_GenOmegaMinus       -> Fill (lThisPt, fAmplitude_V0AEq + fAmplitude_V0CEq);  
1384         }
1385         if( lThisPDG == -3334 && TMath::Abs(lThisRap) < 0.5 ){
1386           fHistPt_GenOmegaPlus                -> Fill (lThisPt);          
1387           fHistPtVsRefMultEta5_GenOmegaPlus   -> Fill (lThisPt, fRefMultEta5); 
1388           fHistPtVsRefMultEta8_GenOmegaPlus   -> Fill (lThisPt, fRefMultEta8);
1389           //Centralities 
1390           fHistPtVsCentV0A_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0A);  
1391           fHistPtVsCentV0C_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0C);  
1392           fHistPtVsCentV0M_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0M);  
1393           fHistPtVsCentV0AEq_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0AEq);  
1394           fHistPtVsCentV0CEq_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0CEq);  
1395           fHistPtVsCentV0MEq_GenOmegaPlus       -> Fill (lThisPt, fCentrality_V0MEq);  
1396           //Amplitudes 
1397           fHistPtVsAmpV0A_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0A);  
1398           fHistPtVsAmpV0C_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0C);  
1399           fHistPtVsAmpV0M_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0A + fAmplitude_V0C);  
1400           fHistPtVsAmpV0AEq_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0AEq);  
1401           fHistPtVsAmpV0CEq_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0CEq);  
1402           fHistPtVsAmpV0MEq_GenOmegaPlus       -> Fill (lThisPt, fAmplitude_V0AEq + fAmplitude_V0CEq);  
1403         }
1404       }
1405    }//End of loop on tracks
1406 //----- End Loop on Cascades ------------------------------------------------------------
1407
1408 //------------------------------------------------
1409 // Fill V0 Tree as needed
1410 //------------------------------------------------
1411
1412 //Variable definition
1413    Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;
1414    Double_t lChi2V0 = 0;
1415    Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
1416    Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
1417    Double_t lV0CosineOfPointingAngle = 0;
1418    Double_t lV0Radius = 0, lPt = 0;
1419    Double_t lRapK0Short = 0, lRapLambda = 0;
1420    Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
1421    Double_t lAlphaV0 = 0, lPtArmV0 = 0;
1422
1423    Double_t fMinV0Pt = 0; 
1424    Double_t fMaxV0Pt = 100; 
1425
1426    Int_t nv0s = 0;
1427    nv0s = lESDevent->GetNumberOfV0s();
1428
1429    for (Int_t iV0 = 0; iV0 < nv0s; iV0++) //extra-crazy test
1430    {// This is the begining of the V0 loop
1431       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
1432       if (!v0) continue;
1433
1434       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]); 
1435
1436       Double_t tV0mom[3];
1437       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] ); 
1438       Double_t lV0TotalMomentum = TMath::Sqrt(
1439       tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
1440
1441       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
1442
1443       lPt = v0->Pt();
1444       lRapK0Short = v0->RapK0Short();
1445       lRapLambda  = v0->RapLambda();
1446       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
1447
1448       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
1449       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
1450
1451       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1452       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1453
1454       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
1455       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
1456       if (!pTrack || !nTrack) {
1457          Printf("ERROR: Could not retreive one of the daughter track");
1458          continue;
1459       }
1460
1461       //Daughter Eta for Eta selection, afterwards
1462       fTreeVariableNegEta = nTrack->Eta();
1463       fTreeVariablePosEta = pTrack->Eta();
1464
1465       // Filter like-sign V0 (next: add counter and distribution)
1466       if ( pTrack->GetSign() == nTrack->GetSign()){
1467          continue;
1468       } 
1469
1470       //________________________________________________________________________
1471       // Track quality cuts 
1472       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
1473       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
1474       fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
1475       if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
1476          fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
1477
1478       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
1479       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
1480       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
1481      
1482   
1483       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
1484         
1485       //GetKinkIndex condition
1486       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
1487
1488       //Findable clusters > 0 condition
1489       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
1490
1491       //Compute ratio Crossed Rows / Findable clusters
1492       //Note: above test avoids division by zero! 
1493       Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); 
1494       Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); 
1495
1496       fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
1497       if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
1498          fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
1499
1500       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
1501       if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
1502
1503       //End track Quality Cuts
1504       //________________________________________________________________________
1505
1506       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(lBestPrimaryVtxPos[0],
1507                                                         lBestPrimaryVtxPos[1],
1508                                                         lMagneticField) );
1509
1510       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(lBestPrimaryVtxPos[0],
1511                                                         lBestPrimaryVtxPos[1],
1512                                                         lMagneticField) );
1513
1514       lOnFlyStatus = v0->GetOnFlyStatus();
1515       lChi2V0 = v0->GetChi2V0();
1516       lDcaV0Daughters = v0->GetDcaV0Daughters();
1517       lDcaV0ToPrimVertex = v0->GetD(lBestPrimaryVtxPos[0],lBestPrimaryVtxPos[1],lBestPrimaryVtxPos[2]);
1518       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(lBestPrimaryVtxPos[0],lBestPrimaryVtxPos[1],lBestPrimaryVtxPos[2]);
1519       fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
1520
1521       // Getting invariant mass infos directly from ESD
1522       v0->ChangeMassHypothesis(310);
1523       lInvMassK0s = v0->GetEffMass();
1524       v0->ChangeMassHypothesis(3122);
1525       lInvMassLambda = v0->GetEffMass();
1526       v0->ChangeMassHypothesis(-3122);
1527       lInvMassAntiLambda = v0->GetEffMass();
1528       lAlphaV0 = v0->AlphaV0();
1529       lPtArmV0 = v0->PtArmV0();
1530
1531
1532 //===============================================
1533 // Monte Carlo Association starts here
1534 //===============================================
1535
1536       //---> Set Everything to "I don't know" before starting
1537
1538       fTreeVariablePIDPositive = 0;
1539       fTreeVariablePIDNegative = 0;
1540
1541       fTreeVariablePtMother = -1;
1542       fTreeVariablePtMC = -1;
1543       fTreeVariableRapMC = -100;
1544
1545       fTreeVariablePID = -1; 
1546       fTreeVariablePIDMother = -1;
1547
1548       fTreeVariablePrimaryStatus = 0; 
1549       fTreeVariablePrimaryStatusMother = 0; 
1550     
1551       Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrack->GetLabel() );
1552       Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrack->GetLabel() );
1553                 
1554       TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1555       TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1556             
1557       Int_t lPIDPositive = mcPosV0Dghter -> GetPdgCode();
1558       Int_t lPIDNegative = mcNegV0Dghter -> GetPdgCode();
1559
1560       fTreeVariablePIDPositive = lPIDPositive;
1561       fTreeVariablePIDNegative = lPIDNegative;
1562
1563       Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
1564       Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
1565     
1566       if( lblMotherPosV0Dghter == lblMotherNegV0Dghter && lblMotherPosV0Dghter > -1 ){
1567          //either label is fine, they're equal at this stage
1568          TParticle* pThisV0 = lMCstack->Particle( lblMotherPosV0Dghter ); 
1569          //Set tree variables
1570          fTreeVariablePID   = pThisV0->GetPdgCode(); //PDG Code
1571          fTreeVariablePtMC  = pThisV0->Pt(); //Perfect Pt
1572
1573          //Only Interested if it's a Lambda, AntiLambda or K0s 
1574          //Avoid the Junction Bug! PYTHIA has particles with Px=Py=Pz=E=0 occasionally, 
1575          //having particle code 88 (unrecognized by PDG), for documentation purposes.
1576          //Even ROOT's TParticle::Y() is not prepared to deal with that exception!
1577          //Note that TParticle::Pt() is immune (that would just return 0)...
1578          //Though granted that that should be extremely rare in this precise condition...
1579          if( TMath::Abs(fTreeVariablePID) == 3122 || fTreeVariablePID==310 ){
1580             fTreeVariableRapMC = pThisV0->Y(); //Perfect Y
1581          }
1582          if( lMCstack->IsPhysicalPrimary       (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 1; //Is Primary!
1583          if( lMCstack->IsSecondaryFromWeakDecay(lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 2; //Weak Decay!
1584          if( lMCstack->IsSecondaryFromMaterial (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 3; //Material Int!
1585          
1586          //Now we try to acquire the V0 parent particle, if possible
1587          Int_t lblThisV0Parent = pThisV0->GetFirstMother();
1588          if ( lblThisV0Parent > -1 ){ //if it has a parent, get it and store specs
1589             TParticle* pThisV0Parent = lMCstack->Particle( lblThisV0Parent );
1590             fTreeVariablePIDMother   = pThisV0Parent->GetPdgCode(); //V0 Mother PDG
1591             fTreeVariablePtMother    = pThisV0Parent->Pt();         //V0 Mother Pt
1592             //Primary Status for the V0 Mother particle 
1593             if( lMCstack->IsPhysicalPrimary       (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 1; //Is Primary!
1594             if( lMCstack->IsSecondaryFromWeakDecay(lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 2; //Weak Decay!
1595             if( lMCstack->IsSecondaryFromMaterial (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 3; //Material Int!
1596          }
1597       }
1598
1599       fTreeVariablePt = v0->Pt();
1600       fTreeVariableChi2V0 = lChi2V0; 
1601       fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
1602       fTreeVariableDcaV0Daughters = lDcaV0Daughters;
1603       fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle; 
1604       fTreeVariableV0Radius = lV0Radius;
1605       fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
1606       fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
1607       fTreeVariableInvMassK0s = lInvMassK0s;
1608       fTreeVariableInvMassLambda = lInvMassLambda;
1609       fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
1610       fTreeVariableRapK0Short = lRapK0Short;
1611       fTreeVariableRapLambda = lRapLambda;
1612       fTreeVariableAlphaV0 = lAlphaV0;
1613       fTreeVariablePtArmV0 = lPtArmV0;
1614
1615       //Official means of acquiring N-sigmas 
1616       fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
1617       fTreeVariableNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
1618       fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
1619       fTreeVariableNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
1620
1621 //This requires an Invariant Mass Hypothesis afterwards
1622       fTreeVariableDistOverTotMom = TMath::Sqrt(
1623                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
1624                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
1625                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
1626                                         );
1627       fTreeVariableDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure
1628
1629       //Copy Multiplicity information 
1630       fTreeVariableCentV0A = fCentrality_V0A; 
1631       fTreeVariableCentV0C = fCentrality_V0C; 
1632       fTreeVariableCentV0M = fCentrality_V0M; 
1633       fTreeVariableCentV0AEq = fCentrality_V0AEq; 
1634       fTreeVariableCentV0CEq = fCentrality_V0CEq; 
1635       fTreeVariableCentV0MEq = fCentrality_V0MEq; 
1636       fTreeVariableAmpV0A = fAmplitude_V0A; 
1637       fTreeVariableAmpV0C = fAmplitude_V0C; 
1638       fTreeVariableAmpV0AEq = fAmplitude_V0AEq; 
1639       fTreeVariableAmpV0CEq = fAmplitude_V0CEq; 
1640       fTreeVariableRefMultEta8 = fRefMultEta8;
1641       fTreeVariableRefMultEta5 = fRefMultEta5;
1642       fTreeVariableRunNumber = fRunNumber; 
1643
1644 //------------------------------------------------
1645 // Fill Tree! 
1646 //------------------------------------------------
1647      
1648      // The conditionals are meant to decrease excessive
1649      // memory usage!
1650      
1651      //First Selection: Reject OnFly
1652      if( lOnFlyStatus == 0 ){
1653        //Second Selection: rough 20-sigma band, parametric.
1654        //K0Short: Enough to parametrize peak broadening with linear function.
1655        Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt;
1656        Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
1657        //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
1658        //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
1659        Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt);
1660        Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
1661        //Do Selection
1662        if( (fTreeVariableInvMassLambda    < lUpperLimitLambda  && fTreeVariableInvMassLambda     > lLowerLimitLambda     ) ||
1663           (fTreeVariableInvMassAntiLambda < lUpperLimitLambda  && fTreeVariableInvMassAntiLambda > lLowerLimitLambda     ) ||
1664           (fTreeVariableInvMassK0s        < lUpperLimitK0Short && fTreeVariableInvMassK0s        > lLowerLimitK0Short    ) ){
1665          //Pre-selection in case this is AA...
1666          if ( TMath::Abs(fTreeVariableNegEta)<0.8 && TMath::Abs(fTreeVariablePosEta)<0.8 && fkSaveV0Tree ) fTreeV0->Fill();
1667        }
1668      }
1669    }// This is the end of the V0 loop
1670
1671 //------------------------------------------------
1672 // Fill V0 tree over.
1673 //------------------------------------------------
1674
1675 //------------------------------------------------
1676 // Rerun cascade vertexer! 
1677 //------------------------------------------------
1678     
1679   if( fkRunVertexers ){ 
1680     lESDevent->ResetCascades();
1681     lESDevent->ResetV0s();
1682
1683     AliV0vertexer lV0vtxer;
1684     AliCascadeVertexer lCascVtxer;
1685                   
1686     lV0vtxer.SetDefaultCuts(fV0VertexerSels);
1687     lCascVtxer.SetDefaultCuts(fCascadeVertexerSels);
1688
1689     lV0vtxer.Tracks2V0vertices(lESDevent);
1690     lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
1691   }
1692
1693 //------------------------------------------------
1694 // MAIN CASCADE LOOP STARTS HERE
1695 //------------------------------------------------
1696 // Code Credit: Antonin Maire (thanks^100)
1697 // ---> This is an adaptation
1698
1699   Long_t ncascades = 0;
1700         ncascades = lESDevent->GetNumberOfCascades();
1701   
1702   for (Int_t iXi = 0; iXi < ncascades; iXi++){
1703     //------------------------------------------------
1704     // Initializations
1705     //------------------------------------------------  
1706           //Double_t lTrkgPrimaryVtxRadius3D = -500.0;
1707           //Double_t lBestPrimaryVtxRadius3D = -500.0;
1708
1709           // - 1st part of initialisation : variables needed to store AliESDCascade data members
1710           Double_t lEffMassXi      = 0. ;
1711           //Double_t lChi2Xi         = -1. ;
1712           Double_t lDcaXiDaughters = -1. ;
1713           Double_t lXiCosineOfPointingAngle = -1. ;
1714           Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
1715           Double_t lXiRadius = -1000. ;
1716           
1717           // - 2nd part of initialisation : Nbr of clusters within TPC for the 3 daughter cascade tracks
1718           Int_t    lPosTPCClusters    = -1; // For ESD only ...//FIXME : wait for availability in AOD
1719           Int_t    lNegTPCClusters    = -1; // For ESD only ...
1720           Int_t    lBachTPCClusters   = -1; // For ESD only ...
1721                         
1722           // - 3rd part of initialisation : about V0 part in cascades
1723           Double_t lInvMassLambdaAsCascDghter = 0.;
1724           //Double_t lV0Chi2Xi         = -1. ;
1725           Double_t lDcaV0DaughtersXi = -1.;
1726                 
1727           Double_t lDcaBachToPrimVertexXi = -1., lDcaV0ToPrimVertexXi = -1.;
1728           Double_t lDcaPosToPrimVertexXi  = -1.;
1729           Double_t lDcaNegToPrimVertexXi  = -1.;
1730           Double_t lV0CosineOfPointingAngleXi = -1. ;
1731           Double_t lV0CosineOfPointingAngleXiSpecial = -1. ;
1732           Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
1733           Double_t lV0RadiusXi = -1000.0;
1734           Double_t lV0quality  = 0.;
1735         
1736           // - 4th part of initialisation : Effective masses
1737           Double_t lInvMassXiMinus    = 0.;
1738           Double_t lInvMassXiPlus     = 0.;
1739           Double_t lInvMassOmegaMinus = 0.;
1740           Double_t lInvMassOmegaPlus  = 0.;
1741     
1742           // - 6th part of initialisation : extra info for QA
1743           Double_t lXiMomX       = 0. , lXiMomY = 0., lXiMomZ = 0.;
1744           Double_t lXiTransvMom  = 0. ;
1745           //Double_t lXiTransvMomMC= 0. ;
1746           Double_t lXiTotMom     = 0. ;
1747                 
1748           Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
1749           //Double_t lBachTransvMom  = 0.;
1750           //Double_t lBachTotMom     = 0.;
1751
1752     fTreeCascVarNegNSigmaPion   = -100;
1753     fTreeCascVarNegNSigmaProton = -100;
1754     fTreeCascVarPosNSigmaPion   = -100;
1755     fTreeCascVarPosNSigmaProton = -100;
1756     fTreeCascVarBachNSigmaPion  = -100;
1757     fTreeCascVarBachNSigmaKaon  = -100;
1758         
1759           Short_t  lChargeXi = -2;
1760           //Double_t lV0toXiCosineOfPointingAngle = 0. ;
1761         
1762           Double_t lRapXi   = -20.0, lRapOmega = -20.0, lRapMC = -20;//  lEta = -20.0, lTheta = 360., lPhi = 720. ;
1763           //Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
1764             
1765     // -------------------------------------
1766     // II.ESD - Calculation Part dedicated to Xi vertices (ESD)
1767     
1768           AliESDcascade *xi = lESDevent->GetCascade(iXi);
1769           if (!xi) continue;
1770         
1771                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)        
1772                 //-------------
1773           lV0quality = 0.;
1774           xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
1775
1776           lEffMassXi                    = xi->GetEffMassXi();
1777           //lChi2Xi                         = xi->GetChi2Xi();
1778           lDcaXiDaughters       = xi->GetDcaXiDaughters();
1779           lXiCosineOfPointingAngle                  = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
1780                                                                                  lBestPrimaryVtxPos[1],
1781                                                                                  lBestPrimaryVtxPos[2] );
1782                   // Take care : the best available vertex should be used (like in AliCascadeVertexer)
1783         
1784           xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] ); 
1785           lXiRadius                     = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );           
1786
1787                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
1788                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
1789                 //-------------
1790                 
1791         UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
1792         UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
1793         UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
1794                 // Care track label can be negative in MC production (linked with the track quality)
1795                 // However = normally, not the case for track index ...
1796           
1797           // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
1798           if(lBachIdx == lIdxNegXi) {
1799                   AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
1800           }
1801     if(lBachIdx == lIdxPosXi) {
1802         AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
1803           }
1804           
1805           AliESDtrack *pTrackXi         = lESDevent->GetTrack( lIdxPosXi );
1806           AliESDtrack *nTrackXi         = lESDevent->GetTrack( lIdxNegXi );
1807           AliESDtrack *bachTrackXi      = lESDevent->GetTrack( lBachIdx );
1808
1809           if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1810                   AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
1811                   continue;
1812           }
1813
1814       fTreeCascVarPosEta = pTrackXi->Eta();
1815       fTreeCascVarNegEta = nTrackXi->Eta();
1816       fTreeCascVarBachEta = bachTrackXi->Eta();
1817       
1818       Double_t lBMom[3], lNMom[3], lPMom[3];
1819       xi->GetBPxPyPz( lBMom[0], lBMom[1], lBMom[2] );
1820       xi->GetPPxPyPz( lPMom[0], lPMom[1], lPMom[2] );
1821       xi->GetNPxPyPz( lNMom[0], lNMom[1], lNMom[2] );
1822       
1823       //fTreeCascVarBachTransMom = TMath::Sqrt( lBMom[0]*lBMom[0] + lBMom[1]*lBMom[1] );
1824       //fTreeCascVarPosTransMom  = TMath::Sqrt( lPMom[0]*lPMom[0] + lPMom[1]*lPMom[1] );
1825       //fTreeCascVarNegTransMom  = TMath::Sqrt( lNMom[0]*lNMom[0] + lNMom[1]*lNMom[1] );
1826   
1827     //------------------------------------------------
1828     // TPC dEdx information 
1829     //------------------------------------------------
1830     fTreeCascVarNegNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kPion   );
1831     fTreeCascVarNegNSigmaProton = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kProton );
1832     fTreeCascVarPosNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kPion );
1833     fTreeCascVarPosNSigmaProton = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kProton );
1834     fTreeCascVarBachNSigmaPion  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kPion );
1835     fTreeCascVarBachNSigmaKaon  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kKaon );
1836
1837     //------------------------------------------------
1838     // TPC Number of clusters info
1839     // --- modified to save the smallest number 
1840     // --- of TPC clusters for the 3 tracks
1841     //------------------------------------------------
1842               
1843           lPosTPCClusters   = pTrackXi->GetTPCNcls();
1844           lNegTPCClusters   = nTrackXi->GetTPCNcls();
1845           lBachTPCClusters  = bachTrackXi->GetTPCNcls(); 
1846
1847     // 1 - Poor quality related to TPCrefit
1848           ULong_t pStatus    = pTrackXi->GetStatus();
1849           ULong_t nStatus    = nTrackXi->GetStatus();
1850           ULong_t bachStatus = bachTrackXi->GetStatus();
1851
1852     //fTreeCascVarkITSRefitBachelor = kTRUE; 
1853     //fTreeCascVarkITSRefitNegative = kTRUE; 
1854     //fTreeCascVarkITSRefitPositive = kTRUE; 
1855
1856     if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
1857     if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
1858     if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
1859
1860     //Extra Debug Information: booleans for ITS refit
1861     //if ((pStatus&AliESDtrack::kITSrefit)    == 0) { fTreeCascVarkITSRefitPositive = kFALSE; }
1862     //if ((nStatus&AliESDtrack::kITSrefit)    == 0) { fTreeCascVarkITSRefitNegative = kFALSE; }
1863     //if ((bachStatus&AliESDtrack::kITSrefit) == 0) { fTreeCascVarkITSRefitBachelor = kFALSE; }
1864
1865           // 2 - Poor quality related to TPC clusters: lowest cut of 70 clusters
1866     if(lPosTPCClusters  < 70) { AliWarning("Pb / V0 Pos. track has less than 70 TPC clusters ... continue!"); continue; }
1867           if(lNegTPCClusters  < 70) { AliWarning("Pb / V0 Neg. track has less than 70 TPC clusters ... continue!"); continue; }
1868           if(lBachTPCClusters < 70) { AliWarning("Pb / Bach.   track has less than 70 TPC clusters ... continue!"); continue; }
1869           Int_t leastnumberofclusters = 1000; 
1870           if( lPosTPCClusters < leastnumberofclusters ) leastnumberofclusters = lPosTPCClusters;
1871           if( lNegTPCClusters < leastnumberofclusters ) leastnumberofclusters = lNegTPCClusters;
1872           if( lBachTPCClusters < leastnumberofclusters ) leastnumberofclusters = lBachTPCClusters;
1873
1874           lInvMassLambdaAsCascDghter    = xi->GetEffMass();
1875           // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
1876           lDcaV0DaughtersXi             = xi->GetDcaV0Daughters(); 
1877           //lV0Chi2Xi                   = xi->GetChi2V0();
1878         
1879           lV0CosineOfPointingAngleXi    = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
1880                                                                             lBestPrimaryVtxPos[1],
1881                                                                             lBestPrimaryVtxPos[2] );
1882     //Modification: V0 CosPA wrt to Cascade decay vertex
1883           lV0CosineOfPointingAngleXiSpecial     = xi->GetV0CosineOfPointingAngle( lPosXi[0],
1884                                                                             lPosXi[1],
1885                                                                             lPosXi[2] );
1886
1887           lDcaV0ToPrimVertexXi          = xi->GetD( lBestPrimaryVtxPos[0], 
1888                                                       lBestPrimaryVtxPos[1], 
1889                                                       lBestPrimaryVtxPos[2] );
1890                 
1891           lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD(       lBestPrimaryVtxPos[0], 
1892                                                                 lBestPrimaryVtxPos[1], 
1893                                                                 lMagneticField  ) ); 
1894                                           // Note : AliExternalTrackParam::GetD returns an algebraic value ...
1895                 
1896           xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] ); 
1897           lV0RadiusXi           = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1898         
1899           lDcaPosToPrimVertexXi         = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0], 
1900                                                                 lBestPrimaryVtxPos[1], 
1901                                                                 lMagneticField  )     ); 
1902         
1903           lDcaNegToPrimVertexXi         = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0], 
1904                                                                 lBestPrimaryVtxPos[1], 
1905                                                                 lMagneticField  )     ); 
1906                 
1907           // - II.Step 4 : around effective masses (ESD)
1908           // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1909                 
1910           if( bachTrackXi->Charge() < 0 )       {
1911                   lV0quality = 0.;
1912                   xi->ChangeMassHypothesis(lV0quality , 3312);  
1913                           // Calculate the effective mass of the Xi- candidate. 
1914                           // pdg code 3312 = Xi-
1915                   lInvMassXiMinus = xi->GetEffMassXi();
1916                 
1917                   lV0quality = 0.;
1918                   xi->ChangeMassHypothesis(lV0quality , 3334);  
1919                           // Calculate the effective mass of the Xi- candidate. 
1920                           // pdg code 3334 = Omega-
1921                   lInvMassOmegaMinus = xi->GetEffMassXi();
1922                                         
1923                   lV0quality = 0.;
1924                   xi->ChangeMassHypothesis(lV0quality , 3312);  // Back to default hyp.
1925           }// end if negative bachelor
1926         
1927         
1928           if( bachTrackXi->Charge() >  0 ){
1929                   lV0quality = 0.;
1930                   xi->ChangeMassHypothesis(lV0quality , -3312);         
1931                           // Calculate the effective mass of the Xi+ candidate. 
1932                           // pdg code -3312 = Xi+
1933                   lInvMassXiPlus = xi->GetEffMassXi();
1934                 
1935                   lV0quality = 0.;
1936                   xi->ChangeMassHypothesis(lV0quality , -3334);         
1937                           // Calculate the effective mass of the Xi+ candidate. 
1938                           // pdg code -3334  = Omega+
1939                   lInvMassOmegaPlus = xi->GetEffMassXi();
1940                 
1941                   lV0quality = 0.;
1942                   xi->ChangeMassHypothesis(lV0quality , -3312);         // Back to "default" hyp.
1943           }// end if positive bachelor
1944                   // - II.Step 6 : extra info for QA (ESD)
1945                   // miscellaneous pieces of info that may help regarding data quality assessment.
1946                   //-------------
1947
1948           xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1949                   lXiTransvMom          = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1950                   lXiTotMom     = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1951                 
1952           xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
1953                   //lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1954                   //lBachTotMom         = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1955
1956           lChargeXi = xi->Charge();
1957
1958           //lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1959         
1960           lRapXi    = xi->RapXi();
1961           lRapOmega = xi->RapOmega();
1962           //lEta      = xi->Eta();
1963           //lTheta    = xi->Theta() *180.0/TMath::Pi();
1964           //lPhi      = xi->Phi()   *180.0/TMath::Pi();
1965           //lAlphaXi  = xi->AlphaXi();
1966           //lPtArmXi  = xi->PtArmXi();
1967
1968 //------------------------------------------------
1969 // Associate Cascade Candidates to Monte Carlo!
1970 //------------------------------------------------      
1971
1972 //Warning: Not using Continues... Need to fill tree later!
1973
1974   Double_t lXiTransvMomMC= 0. ; 
1975         Int_t lPDGCodeCascade = 0;      
1976         Int_t lPID_BachMother = 0;
1977         Int_t lPID_NegMother = 0;
1978         Int_t lPID_PosMother = 0;
1979   fTreeCascVarIsPhysicalPrimary = 0; // 0: not defined, any candidate may have this
1980
1981         if(fDebug > 5)
1982                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent() 
1983                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
1984         
1985 //----------------------------------------
1986 // Regular MC ASSOCIATION STARTS HERE
1987 //----------------------------------------
1988
1989           Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );  
1990                   // Abs value = needed ! question of quality track association ...
1991           Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
1992           Int_t lblBach        = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
1993
1994           TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1995           TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1996           TParticle* mcBach        = lMCstack->Particle( lblBach );
1997       
1998     //fTreeCascVarPosTransMomMC = mcPosV0Dghter->Pt();
1999     //fTreeCascVarNegTransMomMC = mcNegV0Dghter->Pt();
2000
2001           //fTreeCascVarPIDPositive = mcPosV0Dghter -> GetPdgCode();
2002           //fTreeCascVarPIDNegative = mcNegV0Dghter -> GetPdgCode();
2003           //fTreeCascVarPIDBachelor = mcBach->GetPdgCode();
2004
2005           // - Step 4.2 : level of the Xi daughters
2006                 
2007           Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
2008           Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
2009         
2010           //Rather uncivilized: Open brackets for each 'continue'
2011           if(! (lblMotherPosV0Dghter != lblMotherNegV0Dghter) ) { // same mother
2012           if(! (lblMotherPosV0Dghter < 0) ) { // mother != primary (!= -1)
2013           if(! (lblMotherNegV0Dghter < 0) ) {
2014                                         
2015                 // mothers = Lambda candidate ... a priori
2016         
2017           TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
2018           TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );
2019                         
2020           // - Step 4.3 : level of Xi candidate
2021         
2022           Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetFirstMother() ;
2023           Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetFirstMother() ;
2024                                 
2025                 if(! (lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter) ) {
2026                 if(! (lblGdMotherPosV0Dghter < 0) ) { // primary lambda ...
2027                 if(! (lblGdMotherNegV0Dghter < 0) ) { // primary lambda ...
2028
2029                   // Gd mothers = Xi candidate ... a priori
2030         
2031           TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
2032           TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
2033                                         
2034           Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );
2035         
2036   //            if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
2037                   if(!(lblMotherBach != lblGdMotherPosV0Dghter)) { //same mother for bach and V0 daughters
2038         
2039           TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
2040         
2041     lPID_BachMother = mcMotherBach->GetPdgCode();
2042           lPID_NegMother = mcGdMotherPosV0Dghter->GetPdgCode();
2043           lPID_PosMother = mcGdMotherNegV0Dghter->GetPdgCode();
2044    
2045           if(lPID_BachMother==lPID_NegMother && lPID_BachMother==lPID_PosMother){ 
2046                   lPDGCodeCascade = lPID_BachMother; 
2047       lXiTransvMomMC = mcMotherBach->Pt();
2048       if( lMCstack->IsPhysicalPrimary       (lblMotherBach) ) fTreeCascVarIsPhysicalPrimary = 1; //Is Primary!
2049       if( lMCstack->IsSecondaryFromWeakDecay(lblMotherBach) ) fTreeCascVarIsPhysicalPrimary = 2; //Weak Decay!
2050       if( lMCstack->IsSecondaryFromMaterial (lblMotherBach) ) fTreeCascVarIsPhysicalPrimary = 3; //From Material!
2051       if ( (mcMotherBach->Energy() + mcMotherBach->Pz()) / (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) !=0 ){
2052         lRapMC = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) / (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
2053       }
2054           }
2055
2056   }}}}}}} //Ends all conditionals above...
2057
2058 //----------------------------------------
2059 // Regular MC ASSOCIATION ENDS HERE
2060 //----------------------------------------
2061
2062
2063   //------------------------------------------------
2064   // Set Variables for adding to tree
2065   //------------------------------------------------            
2066         
2067           fTreeCascVarCharge    = lChargeXi;
2068           fTreeCascVarPID = lPDGCodeCascade; 
2069           if(lInvMassXiMinus!=0)    fTreeCascVarMassAsXi = lInvMassXiMinus;
2070           if(lInvMassXiPlus!=0)     fTreeCascVarMassAsXi = lInvMassXiPlus;
2071           if(lInvMassOmegaMinus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaMinus;
2072           if(lInvMassOmegaPlus!=0)  fTreeCascVarMassAsOmega = lInvMassOmegaPlus;
2073           fTreeCascVarPt = lXiTransvMom;
2074           fTreeCascVarPtMC = lXiTransvMomMC;
2075           fTreeCascVarRapXi = lRapXi ;
2076           fTreeCascVarRapMC = lRapMC ;
2077           fTreeCascVarRapOmega = lRapOmega ;
2078           fTreeCascVarDCACascDaughters = lDcaXiDaughters;
2079           fTreeCascVarDCABachToPrimVtx = lDcaBachToPrimVertexXi;
2080           fTreeCascVarDCAV0Daughters = lDcaV0DaughtersXi;
2081           fTreeCascVarDCAV0ToPrimVtx = lDcaV0ToPrimVertexXi;
2082           fTreeCascVarDCAPosToPrimVtx = lDcaPosToPrimVertexXi;
2083           fTreeCascVarDCANegToPrimVtx = lDcaNegToPrimVertexXi;
2084           fTreeCascVarCascCosPointingAngle = lXiCosineOfPointingAngle;
2085           fTreeCascVarCascRadius = lXiRadius;
2086           fTreeCascVarV0Mass = lInvMassLambdaAsCascDghter;
2087           fTreeCascVarV0CosPointingAngle = lV0CosineOfPointingAngleXi;
2088           fTreeCascVarV0CosPointingAngleSpecial = lV0CosineOfPointingAngleXiSpecial;
2089           fTreeCascVarV0Radius = lV0RadiusXi;
2090           fTreeCascVarLeastNbrClusters = leastnumberofclusters;
2091
2092           //Copy Multiplicity information 
2093           fTreeCascVarCentV0A = fCentrality_V0A; 
2094           fTreeCascVarCentV0C = fCentrality_V0C; 
2095           fTreeCascVarCentV0M = fCentrality_V0M; 
2096           fTreeCascVarCentV0AEq = fCentrality_V0AEq; 
2097           fTreeCascVarCentV0CEq = fCentrality_V0CEq; 
2098           fTreeCascVarCentV0MEq = fCentrality_V0MEq; 
2099           fTreeCascVarAmpV0A = fAmplitude_V0A; 
2100           fTreeCascVarAmpV0C = fAmplitude_V0C; 
2101           fTreeCascVarAmpV0AEq = fAmplitude_V0AEq; 
2102           fTreeCascVarAmpV0CEq = fAmplitude_V0CEq; 
2103           fTreeCascVarRefMultEta8 = fRefMultEta8;
2104           fTreeCascVarRefMultEta5 = fRefMultEta5;
2105           fTreeCascVarRunNumber = fRunNumber; 
2106
2107           fTreeCascVarDistOverTotMom = TMath::Sqrt(
2108                                                 TMath::Power( lPosXi[0] - lBestPrimaryVtxPos[0] , 2) +
2109                                                 TMath::Power( lPosXi[1] - lBestPrimaryVtxPos[1] , 2) +
2110                                                 TMath::Power( lPosXi[2] - lBestPrimaryVtxPos[2] , 2)
2111                                         );
2112           fTreeCascVarDistOverTotMom /= (lXiTotMom+1e-13);
2113
2114 //All vars not specified here: specified elsewhere!
2115
2116 //------------------------------------------------
2117 // Fill Tree! 
2118 //------------------------------------------------
2119
2120 // The conditional is meant to decrease excessive
2121 // memory usage! Be careful when loosening the 
2122 // cut!
2123
2124   //Xi    Mass window: 150MeV wide
2125   //Omega mass window: 150MeV wide
2126
2127   if( fkSaveCascadeTree && ( (fTreeCascVarMassAsXi<1.32+0.075&&fTreeCascVarMassAsXi>1.32-0.075) ||
2128       (fTreeCascVarMassAsOmega<1.68+0.075&&fTreeCascVarMassAsOmega>1.68-0.075) ) ){
2129       fTreeCascade->Fill();
2130   }
2131
2132 //------------------------------------------------
2133 // Fill tree over.
2134 //------------------------------------------------
2135
2136         }// end of the Cascade loop (ESD or AOD)
2137
2138   // Post output data.
2139   PostData(1, fListHist); 
2140   PostData(2, fTreeEvent);
2141   PostData(3, fTreeV0);
2142   PostData(4, fTreeCascade);
2143 }
2144
2145 //________________________________________________________________________
2146 void AliAnalysisTaskStrangenessVsMultiplicityMC::Terminate(Option_t *)
2147 {
2148    // Draw result to the screen
2149    // Called once at the end of the query
2150
2151    TList *cRetrievedList = 0x0;
2152    cRetrievedList = (TList*)GetOutputData(1);
2153    if(!cRetrievedList){
2154       Printf("ERROR - AliAnalysisTaskStrangenessVsMultiplicityMC : ouput data container list not available\n");
2155       return;
2156    }    
2157         
2158    fHistEventCounter = dynamic_cast<TH1D*> (  cRetrievedList->FindObject("fHistEventCounter")  );
2159    if (!fHistEventCounter) {
2160       Printf("ERROR - AliAnalysisTaskStrangenessVsMultiplicityMC : fHistEventCounter not available");
2161       return;
2162    }
2163   
2164    TCanvas *canCheck = new TCanvas("AliAnalysisTaskStrangenessVsMultiplicityMC","V0 Multiplicity",10,10,510,510);
2165    canCheck->cd(1)->SetLogy();
2166
2167    fHistEventCounter->SetMarkerStyle(22);
2168    fHistEventCounter->DrawCopy("E");
2169 }
2170
2171 //----------------------------------------------------------------------------
2172
2173 Double_t AliAnalysisTaskStrangenessVsMultiplicityMC::MyRapidity(Double_t rE, Double_t rPz) const
2174 {
2175    // Local calculation for rapidity
2176    Double_t ReturnValue = -100;
2177    if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){ 
2178       ReturnValue =  0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
2179    }
2180    return ReturnValue;
2181