]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/STRANGENESS/Cascades/AliAnalysisTaskCheckPerformanceCascadepp276.cxx
Correction to calibration maps (obj names fixed)
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskCheckPerformanceCascadepp276.cxx
1 /***************************************************************          *
2  *  Authors : Antonin Maire, Boris Hippolyte 
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 //-----------------------------------------------------------------
15 //            AliAnalysisTaskCheckPerformanceCascadePbPb class
16 //            This task is for a performance study of cascade identification in PbPb.
17 //            It works with MC info and ESD/AOD.
18 //            Origin   : AliAnalysisTaskCheckPerformanceCascade class by A. Maire Nov2010, antonin.maire@ires.in2p3.fr
19 //            Modified for PbPb analysis: M. Nicassio Feb2011, maria.nicassio@ba.infn.it:
20 //                        - physics selection moved to the run.C macro
21 //                        - added centrality selection and possibility to select events in nTracks ranges 
22 //                        - added new histograms 
23 //                        - modified binning of some histograms and containers 
24 //                        - flag to enable CF container usage 
25 //                        - check in the destructor for CAF usage
26 //                        - flag for acceptance cut in the MC part
27 //                        - in the MC particle selection IsPhysicalPrimary added and number of particles taken as appropriate for HIJING 
28 //                          (however for cascades one gets the same if runs on Nprimaries in the stack and does not check IsPhysicalPrimary)
29 //                        - automatic settings for PID 
30 //                        - selection of injected cascades and HIJING cascades (kind of "bug" in method IsFromBGEvent())
31 //                        - added proper time histograms for cascades and lambdas
32 //                        - cos of PA V0 wrt Xi vertex and not primary vertex  
33 //                        - distance xi-V0 added in the container
34 //                        - AOD analysis developed (January 2012)
35 //
36 //
37 //
38 //              Adapted to pp 2.76 analysis: D. Colella, domenico.colella@ba.infn.it (Nov. 2012):
39 //                        - added new and removed other histograms 
40 //                        - Physics selection moved here (mainly for normalization in the efficiency calcuation)
41 //                        - Centrality selection deleted
42 //                        - 3DHisto denominator moved before any event selection for Normalization
43 //                        - injected and natural part of MC selection removed
44 // 
45 //
46 //
47 //-----------------------------------------------------------------
48
49
50 #include <Riostream.h>
51
52 #include "TList.h"
53 #include "TFile.h"
54 #include "TH1F.h"
55 #include "TH2F.h"
56 #include "TH3F.h"
57 #include "TVector3.h"
58 #include "TCanvas.h"
59 #include "TParticle.h"
60 #include "TMath.h"
61
62 #include "AliLog.h"
63 #include "AliHeader.h"
64 #include "AliMCEvent.h"
65 #include "AliStack.h"
66 #include "AliMultiplicity.h"
67 #include "AliInputEventHandler.h"
68 #include "AliAnalysisManager.h"
69
70 #include "AliCFContainer.h"
71
72 #include "AliESDVZERO.h"
73
74 #include "AliGenEventHeader.h"
75 #include "AliGenCocktailEventHeader.h"
76 #include "AliGenHijingEventHeader.h"
77 #include "AliESDtrackCuts.h"
78 #include "AliPIDResponse.h"
79 #include "AliV0vertexer.h"
80 #include "AliCascadeVertexer.h"
81 #include "AliESDEvent.h"
82 #include "AliESDcascade.h"
83 #include "AliAODEvent.h"
84 #include "AliAODMCParticle.h" 
85 #include "AliAnalysisTaskCheckPerformanceCascadepp276.h"
86
87 using std::cout;
88 using std::endl;
89
90 ClassImp(AliAnalysisTaskCheckPerformanceCascadepp276)
91
92
93
94 //________________________________________________________________________________________
95 AliAnalysisTaskCheckPerformanceCascadepp276::AliAnalysisTaskCheckPerformanceCascadepp276() 
96 : AliAnalysisTaskSE(), // <- take care to AliAnalysisTask( empty )
97     fAnalysisType                  ("ESD"), 
98     fESDtrackCuts                  (0), 
99     fPIDResponse                   (0),
100     fkRerunV0CascVertexers         (0),
101     fkSDDselectionOn               (kTRUE),
102     fkQualityCutZprimVtxPos        (kTRUE),
103     fkRejectEventPileUp            (kTRUE),
104     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
105     fkQualityCutTPCrefit           (kTRUE),
106     fkQualityCutnTPCcls            (kTRUE),
107     fwithSDD                       (kTRUE),
108     fMinnTPCcls                    (0),
109     fkExtraSelections              (0),
110     fVtxRange                      (0),
111     fVtxRangeMin                   (0),
112     fApplyAccCut                   (0),
113     fMinPtCutOnDaughterTracks      (0),
114     fEtaCutOnDaughterTracks        (0),
115     
116     // - Plots initialisation
117     fListHistCascade(0),
118
119     // - General Plots
120     // Cascade multiplicity plots
121     fHistCascadeMultiplicityBeforeAnySel(0),
122     fHistCascadeMultiplicityAfterSDDSel(0),
123     fHistCascadeMultiplicityAfterPhysicsSel(0),
124     fHistCascadeMultiplicityForSelEvtNoTPCOnly(0),
125     fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup(0),
126     fHistCascadeMultiplicityAfterVertexCutSel(0),
127     fHistnXiPlusPerEvTot(0),                  // After any event selections, in all the eta and pt range
128     fHistnXiMinusPerEvTot(0),                 // After any event selections, in all the eta and pt range
129     fHistnOmegaPlusPerEvTot(0),               // After any event selections, in all the eta and pt range
130     fHistnOmegaMinusPerEvTot(0),              // After any event selections, in all the eta and pt range
131     fHistnXiPlusPerEv(0),                     // After any event selections, in the detector acceptance and over a pt minimum
132     fHistnXiMinusPerEv(0),                    // After any event selections, in the detector acceptance and over a pt minimum
133     fHistnOmegaPlusPerEv(0),                  // After any event selections, in the detector acceptance and over a pt minimum
134     fHistnOmegaMinusPerEv(0),                 // After any event selections, in the detector acceptance and over a pt minimum
135     fHistnAssoXiMinus(0),                     // For the Reconstructed-Associated cascades 
136     fHistnAssoXiPlus(0),                      // For the Reconstructed-Associated cascades 
137     fHistnAssoOmegaMinus(0),                  // For the Reconstructed-Associated cascades 
138     fHistnAssoOmegaPlus(0),                   // For the Reconstructed-Associated cascades 
139     // Tracks multiplicity plots
140     fHistTrackMultiplicityBeforeAnySel(0),
141     fHistTrackMultiplicityAfterSDDSel(0),
142     fHistTrackMultiplicityAfterPhysicsSel(0),
143     fHistTrackMultiplicityForSelEvtNoTPCOnly(0),
144     fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup(0),
145     fHistTrackMultiplicityAfterVertexCutSel(0),
146     // Vertex position plots (BestVertex)
147     fHistPVx(0),                              // After any selections but before |Z| < 10 cm
148     fHistPVy(0),                              // After any selections but before |Z| < 10 cm
149     fHistPVz(0),                              // After any selections but before |Z| < 10 cm
150     fHistPVxAnalysis(0),                      // After any event selections
151     fHistPVyAnalysis(0),                      // After any event selections
152     fHistPVzAnalysis(0),                      // After any event selections
153     // - Plots before Physics Selection
154     f3dHistGenPtVsGenYvsNtracksXiMinus(0),    // After the SDD event selection (For efficinecy calculation)
155     f3dHistGenPtVsGenctauvsYXiMinus(0),       // After the SDD event selection (For efficinecy calculation)
156     f3dHistGenPtVsGenYvsNtracksXiPlus(0),     // After the SDD event selection (For efficinecy calculation)
157     f3dHistGenPtVsGenctauvsYXiPlus(0),        // After the SDD event selection (For efficinecy calculation)
158     f3dHistGenPtVsGenYvsNtracksOmegaMinus(0), // After the SDD event selection (For efficinecy calculation)
159     f3dHistGenPtVsGenctauvsYOmegaMinus(0),    // After the SDD event selection (For efficinecy calculation)
160     f3dHistGenPtVsGenYvsNtracksOmegaPlus(0),  // After the SDD event selection (For efficinecy calculation)
161     f3dHistGenPtVsGenctauvsYOmegaPlus(0),     // After the SDD event selection (For efficinecy calculation)
162     // - Generated cascade plots
163     // After all the event selections 
164     //Xi-
165     fHistEtaGenCascXiMinus(0),                // In all the eta and pt range (as they are generated)
166     fHistThetaGenCascXiMinus(0),              // In all the eta and pt range (as they are generated)
167     f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff(0),    // 
168     f3dHistGenPtVsGenctauvsYXiMinusPhysEff(0),       // 
169     f2dHistGenPtVsGenYFdblXiMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
170     fHistThetaLambdaXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
171     fHistThetaBachXiMinus(0),                 // In the detector acceptance and over a pt minimum (Findable particle)
172     fHistThetaMesDghterXiMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
173     fHistThetaBarDghterXiMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
174     fHistPtBachXiMinus(0),                    // In the detector acceptance and over a pt minimum (Findable particle)
175     fHistPtMesDghterXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
176     fHistPtBarDghterXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
177     //Xi+
178     fHistEtaGenCascXiPlus(0),                 // In all the eta and pt range (as they are generated)
179     fHistThetaGenCascXiPlus(0),               // In all the eta and pt range (as they are generated)
180     f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff(0),    // 
181     f3dHistGenPtVsGenctauvsYXiPlusPhysEff(0),       //
182     f2dHistGenPtVsGenYFdblXiPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
183     fHistThetaLambdaXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
184     fHistThetaBachXiPlus(0),                  // In the detector acceptance and over a pt minimum (Findable particle)
185     fHistThetaMesDghterXiPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
186     fHistThetaBarDghterXiPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
187     fHistPtBachXiPlus(0),                     // In the detector acceptance and over a pt minimum (Findable particle)
188     fHistPtMesDghterXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
189     fHistPtBarDghterXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
190     //Omega-
191     fHistEtaGenCascOmegaMinus(0),             // In all the eta and pt range (as they are generated)
192     fHistThetaGenCascOmegaMinus(0),           // In all the eta and pt range (as they are generated)
193     f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff(0),    // 
194     f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff(0),       //
195     f2dHistGenPtVsGenYFdblOmegaMinus(0),      // In the detector acceptance and over a pt minimum (Findable particle)
196     fHistThetaLambdaOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
197     fHistThetaBachOmegaMinus(0),              // In the detector acceptance and over a pt minimum (Findable particle)
198     fHistThetaMesDghterOmegaMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
199     fHistThetaBarDghterOmegaMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
200     fHistPtBachOmegaMinus(0),                 // In the detector acceptance and over a pt minimum (Findable particle)
201     fHistPtMesDghterOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
202     fHistPtBarDghterOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
203     //Omega+      
204     fHistEtaGenCascOmegaPlus(0),              // In all the eta and pt range (as they are generated)
205     fHistThetaGenCascOmegaPlus(0),            // In all the eta and pt range (as they are generated)
206     f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff(0),    // 
207     f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff(0),       //
208     f2dHistGenPtVsGenYFdblOmegaPlus(0),       // In the detector acceptance and over a pt minimum (Findable particle)
209     fHistThetaLambdaOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
210     fHistThetaBachOmegaPlus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
211     fHistThetaMesDghterOmegaPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
212     fHistThetaBarDghterOmegaPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
213     fHistPtBachOmegaPlus(0),                  // In the detector acceptance and over a pt minimum (Findable particle)
214     fHistPtMesDghterOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
215     fHistPtBarDghterOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
216
217     // - Associated to MC cascade plots
218     fHistMassXiMinus(0),                      // For the Reconstructed-Associated cascades
219     fHistMassXiPlus(0),                       // For the Reconstructed-Associated cascades
220     fHistMassOmegaMinus(0),                   // For the Reconstructed-Associated cascades
221     fHistMassOmegaPlus(0),                    // For the Reconstructed-Associated cascades
222     // Effective mass histos with combined PID
223     fHistMassWithCombPIDXiMinus(0),           
224     fHistMassWithCombPIDXiPlus(0),
225     fHistMassWithCombPIDOmegaMinus(0), 
226     fHistMassWithCombPIDOmegaPlus(0),   
227     // PID Probability versus MC Pt(bachelor track)
228     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
229     // Effective mass histos with perfect MC PID on the bachelor
230     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
231     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
232     // Effective mass histos for the cascade candidates associated with MC
233     fHistAsMCMassXiMinus(0),            
234     fHistAsMCMassXiPlus(0),             
235     fHistAsMCMassOmegaMinus(0),
236     fHistAsMCMassOmegaPlus(0),
237     // Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
238     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
239     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
240     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
241     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
242     // Generated Pt Vs generated y, for the cascade candidates associated with MC
243     f2dHistAsMCGenPtVsGenYXiMinus(0),
244     f2dHistAsMCGenPtVsGenYXiPlus(0),
245     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
246     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
247     // Generated Eta of the the cascade candidates associated with MC
248     fHistAsMCGenEtaXiMinus(0),
249     fHistAsMCGenEtaXiPlus(0),
250     fHistAsMCGenEtaOmegaMinus(0),
251     fHistAsMCGenEtaOmegaPlus(0),
252     // Resolution in Pt as function of generated Pt
253     f2dHistAsMCResPtXiMinus(0),         
254     f2dHistAsMCResPtXiPlus(0),          
255     f2dHistAsMCResPtOmegaMinus(0),
256     f2dHistAsMCResPtOmegaPlus(0),       
257     // Resolution in R(2D) as function of generated R
258     f2dHistAsMCResRXiMinus(0),          
259     f2dHistAsMCResRXiPlus(0),           
260     f2dHistAsMCResROmegaMinus(0),
261     f2dHistAsMCResROmegaPlus(0),
262     // Resolution in phi as function of generated Pt
263     f2dHistAsMCResPhiXiMinus(0),
264     f2dHistAsMCResPhiXiPlus(0),
265     f2dHistAsMCResPhiOmegaMinus(0),
266     f2dHistAsMCResPhiOmegaPlus(0),
267     // Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geat/Fluka correction)
268     f2dHistAsMCptProtonMCptXiMinus(0),
269     f2dHistAsMCptAntiprotonMCptXiPlus(0),
270     f2dHistAsMCptProtonMCptOmegaMinus(0),
271     f2dHistAsMCptAntiprotonMCptOmegaPlus(0),
272     // QA plots
273     fHistV0toXiCosineOfPointingAngle(0),
274     fHistV0CosineOfPointingAnglevsPtXi(0),
275     fHistV0CosineOfPointingAnglevsPtOmega(0), 
276     
277     // Containers                       
278     fCFContCascadePIDAsXiMinus(0),
279     fCFContCascadePIDAsXiPlus(0),
280     fCFContCascadePIDAsOmegaMinus(0),
281     fCFContCascadePIDAsOmegaPlus(0),
282     fCFContAsCascadeCuts(0)
283
284     //____Dummy costructor____
285     {
286         for(Int_t iV0selIdx   = 0; iV0selIdx   < 7; iV0selIdx++   ) { fV0Sels          [iV0selIdx   ] = -1.; }
287         for(Int_t iCascSelIdx = 0; iCascSelIdx < 8; iCascSelIdx++ ) { fCascSels        [iCascSelIdx ] = -1.; }
288     }
289      
290         
291      
292 //_____Non-default Constructor________________________________________________________________
293 AliAnalysisTaskCheckPerformanceCascadepp276::AliAnalysisTaskCheckPerformanceCascadepp276(const char *name) 
294   : AliAnalysisTaskSE(name),
295     fAnalysisType                  ("ESD"), 
296     fESDtrackCuts                  (0),
297     fPIDResponse                   (0),
298     fkRerunV0CascVertexers         (0),
299     fkSDDselectionOn               (kTRUE),
300     fkQualityCutZprimVtxPos        (kTRUE),
301     fkRejectEventPileUp            (kTRUE),
302     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
303     fkQualityCutTPCrefit           (kTRUE),
304     fkQualityCutnTPCcls            (kTRUE),
305     fwithSDD                       (kTRUE),
306     fMinnTPCcls                    (0),
307     fkExtraSelections              (0),
308     fVtxRange                      (0),
309     fVtxRangeMin                   (0),
310     fApplyAccCut                   (0),
311     fMinPtCutOnDaughterTracks      (0),
312     fEtaCutOnDaughterTracks        (0),
313
314     // - Plots initialisation
315     fListHistCascade(0),
316
317     // - General Plots
318     // Cascade multiplicity plots
319     fHistCascadeMultiplicityBeforeAnySel(0),
320     fHistCascadeMultiplicityAfterSDDSel(0),
321     fHistCascadeMultiplicityAfterPhysicsSel(0),
322     fHistCascadeMultiplicityForSelEvtNoTPCOnly(0),
323     fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup(0),
324     fHistCascadeMultiplicityAfterVertexCutSel(0),
325     fHistnXiPlusPerEvTot(0),                  // After any event selections, in all the eta and pt range
326     fHistnXiMinusPerEvTot(0),                 // After any event selections, in all the eta and pt range
327     fHistnOmegaPlusPerEvTot(0),               // After any event selections, in all the eta and pt range
328     fHistnOmegaMinusPerEvTot(0),              // After any event selections, in all the eta and pt range
329     fHistnXiPlusPerEv(0),                     // After any event selections, in the detector acceptance and over a pt minimum
330     fHistnXiMinusPerEv(0),                    // After any event selections, in the detector acceptance and over a pt minimum
331     fHistnOmegaPlusPerEv(0),                  // After any event selections, in the detector acceptance and over a pt minimum
332     fHistnOmegaMinusPerEv(0),                 // After any event selections, in the detector acceptance and over a pt minimum
333     fHistnAssoXiMinus(0),                     // For the Reconstructed-Associated cascades 
334     fHistnAssoXiPlus(0),                      // For the Reconstructed-Associated cascades 
335     fHistnAssoOmegaMinus(0),                  // For the Reconstructed-Associated cascades 
336     fHistnAssoOmegaPlus(0),                   // For the Reconstructed-Associated cascades 
337     // Tracks multiplicity plots
338     fHistTrackMultiplicityBeforeAnySel(0),
339     fHistTrackMultiplicityAfterSDDSel(0),
340     fHistTrackMultiplicityAfterPhysicsSel(0),
341     fHistTrackMultiplicityForSelEvtNoTPCOnly(0),
342     fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup(0),
343     fHistTrackMultiplicityAfterVertexCutSel(0),
344     // Vertex position plots (BestVertex)
345     fHistPVx(0),                              // After any selections but before |Z| < 10 cm
346     fHistPVy(0),                              // After any selections but before |Z| < 10 cm
347     fHistPVz(0),                              // After any selections but before |Z| < 10 cm
348     fHistPVxAnalysis(0),                      // After any event selections
349     fHistPVyAnalysis(0),                      // After any event selections
350     fHistPVzAnalysis(0),                      // After any event selections
351     // - Plots before Physics Selection
352     f3dHistGenPtVsGenYvsNtracksXiMinus(0),    // After the SDD event selection (For efficiency calculation)
353     f3dHistGenPtVsGenctauvsYXiMinus(0),       // After the SDD event selection (For efficiency calculation)
354     f3dHistGenPtVsGenYvsNtracksXiPlus(0),     // After the SDD event selection (For efficiency calculation)
355     f3dHistGenPtVsGenctauvsYXiPlus(0),        // After the SDD event selection (For efficiency calculation)
356     f3dHistGenPtVsGenYvsNtracksOmegaMinus(0), // After the SDD event selection (For efficiency calculation)
357     f3dHistGenPtVsGenctauvsYOmegaMinus(0),    // After the SDD event selection (For efficiency calculation)
358     f3dHistGenPtVsGenYvsNtracksOmegaPlus(0),  // After the SDD event selection (For efficiency calculation)
359     f3dHistGenPtVsGenctauvsYOmegaPlus(0),     // After the SDD event selection (For efficiency calculation)
360     // - Generated cascade plots
361     // After all the event selections 
362     //Xi-
363     fHistEtaGenCascXiMinus(0),                // In all the eta and pt range (as they are generated)
364     fHistThetaGenCascXiMinus(0),              // In all the eta and pt range (as they are generated)
365     f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff(0),    // 
366     f3dHistGenPtVsGenctauvsYXiMinusPhysEff(0),       //
367     f2dHistGenPtVsGenYFdblXiMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
368     fHistThetaLambdaXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
369     fHistThetaBachXiMinus(0),                 // In the detector acceptance and over a pt minimum (Findable particle)
370     fHistThetaMesDghterXiMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
371     fHistThetaBarDghterXiMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
372     fHistPtBachXiMinus(0),                    // In the detector acceptance and over a pt minimum (Findable particle)
373     fHistPtMesDghterXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
374     fHistPtBarDghterXiMinus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
375     //Xi+
376     fHistEtaGenCascXiPlus(0),                 // In all the eta and pt range (as they are generated)
377     fHistThetaGenCascXiPlus(0),               // In all the eta and pt range (as they are generated)
378     f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff(0),    // 
379     f3dHistGenPtVsGenctauvsYXiPlusPhysEff(0),       //
380     f2dHistGenPtVsGenYFdblXiPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
381     fHistThetaLambdaXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
382     fHistThetaBachXiPlus(0),                  // In the detector acceptance and over a pt minimum (Findable particle)
383     fHistThetaMesDghterXiPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
384     fHistThetaBarDghterXiPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
385     fHistPtBachXiPlus(0),                     // In the detector acceptance and over a pt minimum (Findable particle)
386     fHistPtMesDghterXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
387     fHistPtBarDghterXiPlus(0),                // In the detector acceptance and over a pt minimum (Findable particle)
388     //Omega-
389     fHistEtaGenCascOmegaMinus(0),             // In all the eta and pt range (as they are generated)
390     fHistThetaGenCascOmegaMinus(0),           // In all the eta and pt range (as they are generated)
391     f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff(0),    // 
392     f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff(0),       //
393     f2dHistGenPtVsGenYFdblOmegaMinus(0),      // In the detector acceptance and over a pt minimum (Findable particle)
394     fHistThetaLambdaOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
395     fHistThetaBachOmegaMinus(0),              // In the detector acceptance and over a pt minimum (Findable particle)
396     fHistThetaMesDghterOmegaMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
397     fHistThetaBarDghterOmegaMinus(0),         // In the detector acceptance and over a pt minimum (Findable particle)
398     fHistPtBachOmegaMinus(0),                 // In the detector acceptance and over a pt minimum (Findable particle)
399     fHistPtMesDghterOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
400     fHistPtBarDghterOmegaMinus(0),            // In the detector acceptance and over a pt minimum (Findable particle)
401     //Omega+      
402     fHistEtaGenCascOmegaPlus(0),              // In all the eta and pt range (as they are generated)
403     fHistThetaGenCascOmegaPlus(0),            // In all the eta and pt range (as they are generated)
404     f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff(0),    // 
405     f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff(0),       //
406     f2dHistGenPtVsGenYFdblOmegaPlus(0),       // In the detector acceptance and over a pt minimum (Findable particle)
407     fHistThetaLambdaOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
408     fHistThetaBachOmegaPlus(0),               // In the detector acceptance and over a pt minimum (Findable particle)
409     fHistThetaMesDghterOmegaPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
410     fHistThetaBarDghterOmegaPlus(0),          // In the detector acceptance and over a pt minimum (Findable particle)
411     fHistPtBachOmegaPlus(0),                  // In the detector acceptance and over a pt minimum (Findable particle)
412     fHistPtMesDghterOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
413     fHistPtBarDghterOmegaPlus(0),             // In the detector acceptance and over a pt minimum (Findable particle)
414
415     // - Associated to MC cascade plots
416     fHistMassXiMinus(0),                      // For the Reconstructed-Associated cascades
417     fHistMassXiPlus(0),                       // For the Reconstructed-Associated cascades
418     fHistMassOmegaMinus(0),                   // For the Reconstructed-Associated cascades
419     fHistMassOmegaPlus(0),                    // For the Reconstructed-Associated cascades
420     // Effective mass histos with combined PID
421     fHistMassWithCombPIDXiMinus(0),
422     fHistMassWithCombPIDXiPlus(0),
423     fHistMassWithCombPIDOmegaMinus(0),
424     fHistMassWithCombPIDOmegaPlus(0),
425     // PID Probability versus MC Pt(bachelor track)
426     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
427     // Effective mass histos with perfect MC PID on the bachelor
428     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
429     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
430     // Effective mass histos for the cascade candidates associated with MC
431     fHistAsMCMassXiMinus(0),
432     fHistAsMCMassXiPlus(0),
433     fHistAsMCMassOmegaMinus(0),
434     fHistAsMCMassOmegaPlus(0),
435     // Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
436     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
437     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
438     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
439     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
440     // Generated Pt Vs generated y, for the cascade candidates associated with MC
441     f2dHistAsMCGenPtVsGenYXiMinus(0),
442     f2dHistAsMCGenPtVsGenYXiPlus(0),
443     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
444     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
445     // Generated Eta of the the cascade candidates associated with MC
446     fHistAsMCGenEtaXiMinus(0),
447     fHistAsMCGenEtaXiPlus(0),
448     fHistAsMCGenEtaOmegaMinus(0),
449     fHistAsMCGenEtaOmegaPlus(0),
450     // Resolution in Pt as function of generated Pt
451     f2dHistAsMCResPtXiMinus(0),
452     f2dHistAsMCResPtXiPlus(0),
453     f2dHistAsMCResPtOmegaMinus(0),
454     f2dHistAsMCResPtOmegaPlus(0),
455     // Resolution in R(2D) as function of generated R
456     f2dHistAsMCResRXiMinus(0),
457     f2dHistAsMCResRXiPlus(0),
458     f2dHistAsMCResROmegaMinus(0),
459     f2dHistAsMCResROmegaPlus(0),
460     // Resolution in phi as function of generated Pt
461     f2dHistAsMCResPhiXiMinus(0),
462     f2dHistAsMCResPhiXiPlus(0),
463     f2dHistAsMCResPhiOmegaMinus(0),
464     f2dHistAsMCResPhiOmegaPlus(0),
465     // Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geat/Fluka correction)
466     f2dHistAsMCptProtonMCptXiMinus(0),
467     f2dHistAsMCptAntiprotonMCptXiPlus(0),
468     f2dHistAsMCptProtonMCptOmegaMinus(0),
469     f2dHistAsMCptAntiprotonMCptOmegaPlus(0),
470     // QA plots
471     fHistV0toXiCosineOfPointingAngle(0),
472     fHistV0CosineOfPointingAnglevsPtXi(0),
473     fHistV0CosineOfPointingAnglevsPtOmega(0),
474
475     // Containers                       
476     fCFContCascadePIDAsXiMinus(0),
477     fCFContCascadePIDAsXiPlus(0),
478     fCFContCascadePIDAsOmegaMinus(0),
479     fCFContCascadePIDAsOmegaPlus(0),
480     fCFContAsCascadeCuts(0)
481
482     //____Costructor____
483     {
484       // Define input and output slots here
485       // Input slot #0 works with a TChain
486       // Output slot #1 writes into a TList container (cascade)
487         
488         // PbPb default cuts  
489         fV0Sels[0] =  33.;       // max allowed chi2
490         fV0Sels[1] =   0.1;      // min allowed impact parameter for the 1st daughter 
491         fV0Sels[2] =   0.1;      // min allowed impact parameter for the 2nd daughter 
492         fV0Sels[3] =   1.5 ;     // max allowed DCA between the daughter tracks       
493         fV0Sels[4] =   0.9 ;     // min allowed cosine of V0's pointing angle         
494         fV0Sels[5] =   0.2;      // min radius of the fiducial volume                 
495         fV0Sels[6] = 200.  ;     // max radius of the fiducial volume                 
496         fCascSels[0] =  33.;     // max allowed chi2 
497         fCascSels[1] =   0.01;   // min allowed V0 impact parameter                    
498         fCascSels[2] =   0.008;  // "window" around the Lambda mass                    
499         fCascSels[3] =   0.01;   // min allowed bachelor's impact parameter            
500         fCascSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    
501         fCascSels[5] =   0.95;   // min allowed cosine of the cascade pointing angle   
502         fCascSels[6] =   0.2  ;  // min radius of the fiducial volume                  
503         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  
504               
505         DefineOutput(1, TList::Class());
506         DefineOutput(2, AliCFContainer::Class());
507         DefineOutput(3, AliCFContainer::Class());
508         DefineOutput(4, AliCFContainer::Class());
509         DefineOutput(5, AliCFContainer::Class());
510         DefineOutput(6, AliCFContainer::Class());
511     }
512
513     //____Destructor____
514     AliAnalysisTaskCheckPerformanceCascadepp276::~AliAnalysisTaskCheckPerformanceCascadepp276()
515     {
516       // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
517       // They will be deleted when fListCascade is deleted by the TSelector dtor
518       // Because of TList::SetOwner()
519       if (fListHistCascade && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())              {delete fListHistCascade;              fListHistCascade = 0x0;}  
520       if (fCFContCascadePIDAsXiMinus && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())    {delete fCFContCascadePIDAsXiMinus;    fCFContCascadePIDAsXiMinus = 0x0;}
521       if (fCFContCascadePIDAsXiPlus && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())     {delete fCFContCascadePIDAsXiPlus;     fCFContCascadePIDAsXiPlus = 0x0;}
522       if (fCFContCascadePIDAsOmegaMinus && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {delete fCFContCascadePIDAsOmegaMinus; fCFContCascadePIDAsOmegaMinus = 0x0;}
523       if (fCFContCascadePIDAsOmegaPlus && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())  {delete fCFContCascadePIDAsOmegaPlus;  fCFContCascadePIDAsOmegaPlus = 0x0;}
524       if (fCFContAsCascadeCuts && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())          {delete fCFContAsCascadeCuts;          fCFContAsCascadeCuts = 0x0;}
525       if (fESDtrackCuts)                                                                             {delete fESDtrackCuts;                 fESDtrackCuts = 0x0;}
526     }
527
528
529 //________________________________________________________________________
530 void AliAnalysisTaskCheckPerformanceCascadepp276::UserCreateOutputObjects() {
531   // Create histograms
532   // Called once
533
534  // - Option for AliLog: to suppress the extensive info prompted by a run with MC
535  AliLog::SetGlobalLogLevel(AliLog::kError); 
536
537  // - Definition of the output datamembers      
538  fListHistCascade = new TList();
539  fListHistCascade->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
540
541  //-----------------------------------------------
542  // Particle Identification Setup (new PID object)
543  //-----------------------------------------------
544  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
545  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
546  fPIDResponse = inputHandler->GetPIDResponse();
547         
548  // - Only used to get the number of primary reconstructed tracks
549  if (! fESDtrackCuts ){
550       fESDtrackCuts = new AliESDtrackCuts();
551  }
552  
553  //----------------------
554  // Initialize the histos
555  //----------------------
556
557  //----------------------------------
558  // - Same general binning definition
559  Double_t ptBinLimits[101];
560  for (Int_t iptbin = 0; iptbin<101; ++iptbin) ptBinLimits[iptbin]=iptbin*0.1;
561  Double_t yBinLimits[111];
562  for (Int_t iybin = 0; iybin<111; ++iybin) yBinLimits[iybin]=-1.1+iybin*0.02;
563  Double_t ctauBinLimits[112];
564  for (Int_t ict = 0; ict<112; ++ict) ctauBinLimits[ict] = (Double_t) (ict-1.); 
565  
566  //------------------
567  // - General plots
568    // - Cascades multiplicity plots 
569    if(! fHistCascadeMultiplicityBeforeAnySel) {
570         fHistCascadeMultiplicityBeforeAnySel = new TH1F("fHistCascadeMultiplicityBeforeAnySel",
571                         "Cascades per event (before any selections);Nbr of Cascades/Evt;Events", 50, 0, 50);
572         fListHistCascade->Add(fHistCascadeMultiplicityBeforeAnySel);
573    }
574    if(! fHistCascadeMultiplicityAfterSDDSel) {
575         fHistCascadeMultiplicityAfterSDDSel = new TH1F("fHistCascadeMultiplicityAfterSDDSel",
576                         "Cascades per event (after only the SDD selection);Nbr of Cascades/Evt;Events", 50, 0, 50);
577         fListHistCascade->Add(fHistCascadeMultiplicityAfterSDDSel);
578    }
579    if(! fHistCascadeMultiplicityAfterPhysicsSel) {
580         fHistCascadeMultiplicityAfterPhysicsSel = new TH1F("fHistCascadeMultiplicityAfterPhysicsSel",
581                         "Cascades per event (after physics selection);Nbr of Cascades/Evt;Events", 50, 0, 50);
582         fListHistCascade->Add(fHistCascadeMultiplicityAfterPhysicsSel);
583    }
584    if(! fHistCascadeMultiplicityForSelEvtNoTPCOnly) {
585         fHistCascadeMultiplicityForSelEvtNoTPCOnly = new TH1F("fHistCascadeMultiplicityForSelEvtNoTPCOnly",
586                         "Cascades per event (for selected events with well-established PV);Nbr of Cascades/Evt;Events", 50, 0, 50);
587         fListHistCascade->Add(fHistCascadeMultiplicityForSelEvtNoTPCOnly);
588    }
589    if(! fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup) {
590         fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup",
591                         "Cascades per event (for selected events with well-establisched PV and no pile-up);Nbr of Cascades/Evt;Events", 50, 0, 50);
592         fListHistCascade->Add(fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup);
593    }
594    if(! fHistCascadeMultiplicityAfterVertexCutSel) {
595         fHistCascadeMultiplicityAfterVertexCutSel = new TH1F("fHistCascadeMultiplicityAfterVertexCutSel",
596                                                              "Cascades per event (after vertex cut selection);Nbr of Cascades/Evt;Events", 50, 0, 50);
597         fListHistCascade->Add(fHistCascadeMultiplicityAfterVertexCutSel);
598    }
599    // - Tracks multiplicity plots 
600    if(! fHistTrackMultiplicityBeforeAnySel) {
601         fHistTrackMultiplicityBeforeAnySel = new TH1F("fHistTrackMultiplicityBeforeAnySel",
602                         "Tracks per event (before any selections);Nbr of Tracks/Evt;Events", 200, 0, 200);
603         fListHistCascade->Add(fHistTrackMultiplicityBeforeAnySel);
604    }
605    if(! fHistTrackMultiplicityAfterSDDSel) {
606         fHistTrackMultiplicityAfterSDDSel = new TH1F("fHistTrackMultiplicityAfterSDDSel",
607                         "Tracks per event (after only the SDD selection);Nbr of Tracks/Evt;Events", 200, 0, 200);
608         fListHistCascade->Add(fHistTrackMultiplicityAfterSDDSel);
609    }
610    if(! fHistTrackMultiplicityAfterPhysicsSel) {
611         fHistTrackMultiplicityAfterPhysicsSel = new TH1F("fHistTrackMultiplicityAfterPhysicsSel",
612                         "Tracks per event (after physics selection);Nbr of Tracks/Evt;Events", 200, 0, 200);
613         fListHistCascade->Add(fHistTrackMultiplicityAfterPhysicsSel);
614    }
615    if(! fHistTrackMultiplicityForSelEvtNoTPCOnly) {
616         fHistTrackMultiplicityForSelEvtNoTPCOnly = new TH1F("fHistTrackMultiplicityForSelEvtNoTPCOnly",
617                         "Tracks per event (for selected events with well-established PV);Nbr of Tracks/Evt;Events", 200, 0, 200);
618         fListHistCascade->Add(fHistTrackMultiplicityForSelEvtNoTPCOnly);
619    }
620    if(! fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup) {
621         fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup",
622                         "Tracks per event (for selected events with well-establisched PV and no pile-up);Nbr of Tracks/Evt;Events", 200, 0, 200);
623         fListHistCascade->Add(fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup);
624    }
625    if(! fHistTrackMultiplicityAfterVertexCutSel) {
626         fHistTrackMultiplicityAfterVertexCutSel = new TH1F("fHistTrackMultiplicityAfterVertexCutSel",
627                                                            "Tracks per event (after vertex cut selection);Nbr of Tracks/Evt;Events", 200, 0, 200);
628         fListHistCascade->Add(fHistTrackMultiplicityAfterVertexCutSel);
629    }
630    // - Vertex position plots
631    if(! fHistPVx ){
632         fHistPVx = new TH1F("fHistPVx", "Best PV position in x; x (cm); Events", 2000, -0.5, 0.5);
633         fListHistCascade->Add(fHistPVx);
634    }
635    if(! fHistPVy ){
636         fHistPVy = new TH1F("fHistPVy", "Best PV position in y; y (cm); Events", 2000, -0.5, 0.5);
637         fListHistCascade->Add(fHistPVy);
638    }
639    if(! fHistPVz ){
640         fHistPVz = new TH1F("fHistPVz", "Best PV position in z; z (cm); Events", 400, -20, 20);
641         fListHistCascade->Add(fHistPVz);
642    }
643    if(! fHistPVxAnalysis ){
644         fHistPVxAnalysis = new TH1F("fHistPVxAnalysis", "Best PV position in x (after events selections); x (cm); Events", 2000, -0.5, 0.5);
645         fListHistCascade->Add(fHistPVxAnalysis);
646    }
647    if(! fHistPVyAnalysis ){
648         fHistPVyAnalysis = new TH1F("fHistPVyAnalysis", "Best PV position in y (after events selections); y (cm); Events" , 2000, -0.5, 0.5);
649         fListHistCascade->Add(fHistPVyAnalysis);
650    }
651    if(! fHistPVzAnalysis ){
652         fHistPVzAnalysis = new TH1F("fHistPVzAnalysis", "Best PV position in z (after events selections); z (cm); Events", 400, -20, 20);
653         fListHistCascade->Add(fHistPVzAnalysis);
654    }
655
656  //--------------------------
657  // - Generated cascade plots
658    // - Generated Cascade multiplicity distributions (for singol cascade)
659    fHistnXiPlusPerEvTot = new TH1F("fHistnXiPlusPerEvTot", "", 25, 0, 25);
660    fListHistCascade->Add(fHistnXiPlusPerEvTot);
661    fHistnXiMinusPerEvTot = new TH1F("fHistnXiMinusPerEvTot", "", 25, 0, 25);
662    fListHistCascade->Add(fHistnXiMinusPerEvTot);
663    fHistnOmegaPlusPerEvTot = new TH1F("fHistnOmegaPlusPerEvTot", "", 25, 0, 25);
664    fListHistCascade->Add(fHistnOmegaPlusPerEvTot);
665    fHistnOmegaMinusPerEvTot = new TH1F("fHistnOmegaMinusPerEvTot", "", 25, 0, 25);
666    fListHistCascade->Add(fHistnOmegaMinusPerEvTot);   
667    fHistnXiPlusPerEv = new TH1F("fHistnXiPlusPerEv", "", 25, 0, 25);
668    fListHistCascade->Add(fHistnXiPlusPerEv);
669    fHistnXiMinusPerEv = new TH1F("fHistnXiMinusPerEv", "", 25, 0, 25);
670    fListHistCascade->Add(fHistnXiMinusPerEv);
671    fHistnOmegaPlusPerEv = new TH1F("fHistnOmegaPlusPerEv", "", 25, 0, 25);
672    fListHistCascade->Add(fHistnOmegaPlusPerEv);
673    fHistnOmegaMinusPerEv = new TH1F("fHistnOmegaMinusPerEv", "", 25, 0, 25);
674    fListHistCascade->Add(fHistnOmegaMinusPerEv);
675    // - Xi- 
676    // - Pseudo-Rapidity distribution
677    if (!fHistEtaGenCascXiMinus) {
678      fHistEtaGenCascXiMinus = new TH1F("fHistEtaGenCascXiMinus", "#eta of any gen. #Xi^{-}; #eta; Number of Casc", 200, -10, 10);
679      fListHistCascade->Add(fHistEtaGenCascXiMinus);
680    }
681    if (!f3dHistGenPtVsGenYvsNtracksXiMinus) {
682      f3dHistGenPtVsGenYvsNtracksXiMinus = new TH3D("f3dHistGenPtVsGenYvsNtracksXiMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
683      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiMinus);
684    }
685    if (!f3dHistGenPtVsGenctauvsYXiMinus) {
686       f3dHistGenPtVsGenctauvsYXiMinus = new TH3D("f3dHistGenPtVsGenctauvsYXiMinus", "MC P_{t} Vs MC ctau Vs Y of Gen #Xi^{-}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
687       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiMinus);
688    }
689    if (!f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff) {
690      f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff = new TH3D("f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff", "MC P_{t} Vs MC Y of Gen #Xi^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
691      fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff);
692    }
693    if (!f3dHistGenPtVsGenctauvsYXiMinusPhysEff) {
694       f3dHistGenPtVsGenctauvsYXiMinusPhysEff = new TH3D("f3dHistGenPtVsGenctauvsYXiMinusPhysEff", "MC P_{t} Vs MC ctau Vs Y of Gen #Xi^{-}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
695       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiMinusPhysEff);
696    }
697    // - Info at the generation level of multi-strange particle
698    if (!fHistThetaGenCascXiMinus) {
699       fHistThetaGenCascXiMinus = new TH1F("fHistThetaGenCascXiMinus", "#theta of gen. #Xi^{-}; #theta; Number of Casc.", 200, -10, 190);
700       fListHistCascade->Add(fHistThetaGenCascXiMinus);
701    }
702    if (!f2dHistGenPtVsGenYFdblXiMinus) {
703       f2dHistGenPtVsGenYFdblXiMinus = new TH2D("f2dHistGenPtVsGenYFdblXiMinus", "MC P_{t} Vs MC Y of findable Gen #Xi^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
704       fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiMinus);
705    }
706    // - Theta distribution the daughters (control plots)
707    if (!fHistThetaLambdaXiMinus) {
708       fHistThetaLambdaXiMinus = new TH1F("fHistThetaLambdaXiMinus", "#theta of gen. #Lambda (Xi dghter); #theta_{#Lambda}; Number of #Lambda^0", 200, -10, 190);
709       fListHistCascade->Add(fHistThetaLambdaXiMinus);
710    }
711    if (!fHistThetaBachXiMinus) {
712       fHistThetaBachXiMinus = new TH1F("fHistThetaBachXiMinus", "#theta of gen. Bach.; #theta_{Bach}; Number of Bach.", 200, -10, 190);
713       fListHistCascade->Add(fHistThetaBachXiMinus);
714    }
715    if (!fHistThetaMesDghterXiMinus) {
716       fHistThetaMesDghterXiMinus = new TH1F("fHistThetaMesDghterXiMinus", "#theta of gen. Meson #Lambda dghter; #theta_{MesDght}; Number of Mes.", 200, -10, 190);
717       fListHistCascade->Add(fHistThetaMesDghterXiMinus);
718    }
719    if (!fHistThetaBarDghterXiMinus) {
720       fHistThetaBarDghterXiMinus = new TH1F("fHistThetaBarDghterXiMinus", "#theta of gen. Baryon #Lambda dghter; #theta_{BarDght}; Number of Bar.", 200, -10, 190);
721       fListHistCascade->Add(fHistThetaBarDghterXiMinus);
722    }
723    // - Pt distribution (control plots)
724    if (!fHistPtBachXiMinus) {
725       fHistPtBachXiMinus = new TH1F("fHistPtBachXiMinus", "p_{t} of gen. Bach.; pt_{Bach}; Number of Bach.", 200, 0, 10);
726       fListHistCascade->Add(fHistPtBachXiMinus);
727    }
728    if (!fHistPtMesDghterXiMinus) {
729       fHistPtMesDghterXiMinus = new TH1F("fHistPtMesDghterXiMinus", "p_{t} of gen. Meson #Lambda dghter; pt_{MesDght}; Number of Mes.", 200, 0, 10);
730       fListHistCascade->Add(fHistPtMesDghterXiMinus);
731    }
732    if (!fHistPtBarDghterXiMinus) {
733       fHistPtBarDghterXiMinus = new TH1F("fHistPtBarDghterXiMinus", "p_{t} of gen. Baryon #Lambda dghter; pt_{BarDght}; Number of Bar.", 200, 0, 10);
734       fListHistCascade->Add(fHistPtBarDghterXiMinus);
735    }
736    // - Xi+ 
737    // - Pseudo-Rapidity distribution
738    if (!fHistEtaGenCascXiPlus) {
739       fHistEtaGenCascXiPlus = new TH1F("fHistEtaGenCascXiPlus", "#eta of any gen. #Xi^{+}; #eta; Number of Casc", 200, -10, 10);
740       fListHistCascade->Add(fHistEtaGenCascXiPlus);
741    }
742    if (!f3dHistGenPtVsGenYvsNtracksXiPlus) {
743       f3dHistGenPtVsGenYvsNtracksXiPlus = new TH3D("f3dHistGenPtVsGenYvsNtracksXiPlus", "MC P_{t} Vs MC Y of Gen #Xi^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
744       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiPlus);
745    }
746    if (!f3dHistGenPtVsGenctauvsYXiPlus) {
747       f3dHistGenPtVsGenctauvsYXiPlus = new TH3D("f3dHistGenPtVsGenctauvsYXiPlus", "MC P_{t} Vs MC ctau Vs Yof Gen #Xi^{+}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
748       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiPlus);
749    }
750    if (!f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff) {
751       f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff = new TH3D("f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff", "MC P_{t} Vs MC Y of Gen #Xi^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
752       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff);
753    }
754    if (!f3dHistGenPtVsGenctauvsYXiPlusPhysEff) {
755       f3dHistGenPtVsGenctauvsYXiPlusPhysEff = new TH3D("f3dHistGenPtVsGenctauvsYXiPlusPhysEff", "MC P_{t} Vs MC ctau Vs Yof Gen #Xi^{+}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
756       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYXiPlusPhysEff);
757    }
758    // - Info at the generation level of multi-strange particle
759    if (!fHistThetaGenCascXiPlus) {
760       fHistThetaGenCascXiPlus = new TH1F("fHistThetaGenCascXiPlus", "#theta of gen. #Xi^{+}; #theta; Number of Casc.", 200, -10, 190);
761       fListHistCascade->Add(fHistThetaGenCascXiPlus);
762    }
763    if (!f2dHistGenPtVsGenYFdblXiPlus) {
764       f2dHistGenPtVsGenYFdblXiPlus = new TH2D("f2dHistGenPtVsGenYFdblXiPlus", "MC P_{t} Vs MC Y of findable Gen #Xi^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
765       fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiPlus);
766    }
767    // - Theta distribution the daughters (control plots)
768    if (!fHistThetaLambdaXiPlus) {
769       fHistThetaLambdaXiPlus = new TH1F("fHistThetaLambdaXiPlus", "#theta of gen. #Lambda (Xi dghter); #theta_{#Lambda}; Number of #Lambda", 200, -10, 190);
770       fListHistCascade->Add(fHistThetaLambdaXiPlus);
771    }
772    if (!fHistThetaBachXiPlus) {
773       fHistThetaBachXiPlus = new TH1F("fHistThetaBachXiPlus", "#theta of gen. Bach.; #theta_{Bach}; Number of Bach.", 200, -10, 190);
774       fListHistCascade->Add(fHistThetaBachXiPlus);
775    }
776    if (!fHistThetaMesDghterXiPlus) {
777       fHistThetaMesDghterXiPlus = new TH1F("fHistThetaMesDghterXiPlus", "#theta of gen. Meson #Lambda dghter; #theta_{MesDght}; Number of Mes.", 200, -10, 190);
778       fListHistCascade->Add(fHistThetaMesDghterXiPlus);
779    }
780    if (!fHistThetaBarDghterXiPlus) {
781       fHistThetaBarDghterXiPlus = new TH1F("fHistThetaBarDghterXiPlus", "#theta of gen. Baryon #Lambda dghter; #theta_{BarDght}; Number of Bar.", 200, -10, 190);
782       fListHistCascade->Add(fHistThetaBarDghterXiPlus);
783    }
784    // - Pt distribution (control plots)
785    if (!fHistPtBachXiPlus) {
786       fHistPtBachXiPlus = new TH1F("fHistPtBachXiPlus", "p_{t} of gen. Bach.; pt_{Bach}; Number of Bach.", 200, 0, 10);
787       fListHistCascade->Add(fHistPtBachXiPlus);
788    }
789    if (!fHistPtMesDghterXiPlus) {
790       fHistPtMesDghterXiPlus = new TH1F("fHistPtMesDghterXiPlus", "p_{t} of gen. Meson #Lambda dghter; pt_{MesDght}; Number of Mes.", 200, 0, 10);
791       fListHistCascade->Add(fHistPtMesDghterXiPlus);
792    }
793    if (!fHistPtBarDghterXiPlus) {
794       fHistPtBarDghterXiPlus = new TH1F("fHistPtBarDghterXiPlus", "p_{t} of gen. Baryon #Lambda dghter); pt_{BarDght}; Number of Bar.", 200, 0, 10);
795       fListHistCascade->Add(fHistPtBarDghterXiPlus);
796    }
797    // - Omega- 
798    // - Pseudo-Rapidity distribution
799    if (!fHistEtaGenCascOmegaMinus) {
800       fHistEtaGenCascOmegaMinus = new TH1F("fHistEtaGenCascOmegaMinus", "#eta of any gen. #Omega^{-}; #eta; Number of Casc", 200, -10, 10);
801       fListHistCascade->Add(fHistEtaGenCascOmegaMinus);
802    }
803    if (!f3dHistGenPtVsGenYvsNtracksOmegaMinus) {
804       f3dHistGenPtVsGenYvsNtracksOmegaMinus = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaMinus", "MC P_{t} Vs MC Y of Gen #Omega^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
805       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaMinus);
806    }
807    if (!f3dHistGenPtVsGenctauvsYOmegaMinus) {
808       f3dHistGenPtVsGenctauvsYOmegaMinus = new TH3D("f3dHistGenPtVsGenctauvsYOmegaMinus", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{-} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
809       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaMinus);
810    }
811    if (!f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff) {
812       f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff", "MC P_{t} Vs MC Y of Gen #Omega^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
813       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff);
814    }
815    if (!f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff) {
816       f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff = new TH3D("f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{-}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
817       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff);
818    }
819    // - Info at the generation level of multi-strange particle
820    if (!fHistThetaGenCascOmegaMinus) {
821       fHistThetaGenCascOmegaMinus = new TH1F("fHistThetaGenCascOmegaMinus", "#theta of gen. #Omega^{-}; #theta; Number of Casc.", 200, -10, 190);
822       fListHistCascade->Add(fHistThetaGenCascOmegaMinus);
823    }
824    if (!f2dHistGenPtVsGenYFdblOmegaMinus) {
825       f2dHistGenPtVsGenYFdblOmegaMinus = new TH2D("f2dHistGenPtVsGenYFdblOmegaMinus", "MC P_{t} Vs MC Y of findable Gen #Omega^{-}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
826       fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaMinus);
827    }
828    // - Theta distribution the daughters (control plots)
829    if (!fHistThetaLambdaOmegaMinus) {
830       fHistThetaLambdaOmegaMinus = new TH1F("fHistThetaLambdaOmegaMinus", "#theta of gen. #Lambda (Omega dghter); #theta_{#Lambda}; Number of #Lambda", 200, -10, 190);
831       fListHistCascade->Add(fHistThetaLambdaOmegaMinus);
832    }
833    if (!fHistThetaBachOmegaMinus) {
834       fHistThetaBachOmegaMinus = new TH1F("fHistThetaBachOmegaMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
835       fListHistCascade->Add(fHistThetaBachOmegaMinus);
836    }
837    if (!fHistThetaMesDghterOmegaMinus) {
838       fHistThetaMesDghterOmegaMinus = new TH1F("fHistThetaMesDghterOmegaMinus", "#theta of gen. Meson #Lambda dghter; #theta_{MesDght}; Number of Mes.", 200, -10, 190);
839       fListHistCascade->Add(fHistThetaMesDghterOmegaMinus);
840    }
841    if (!fHistThetaBarDghterOmegaMinus) {
842       fHistThetaBarDghterOmegaMinus = new TH1F("fHistThetaBarDghterOmegaMinus", "#theta of gen. Baryon #Lambda dghter; #theta_{BarDght}; Number of Bar.", 200, -10, 190);
843       fListHistCascade->Add(fHistThetaBarDghterOmegaMinus);
844    }
845    // - Pt distribution (control plots)
846    if (!fHistPtBachOmegaMinus) {
847       fHistPtBachOmegaMinus = new TH1F("fHistPtBachOmegaMinus", "p_{t} of gen. Bach.; pt_{Bach}; Number of Bach.", 200, 0, 10);
848       fListHistCascade->Add(fHistPtBachOmegaMinus);
849    }
850    if (!fHistPtMesDghterOmegaMinus) {
851       fHistPtMesDghterOmegaMinus = new TH1F("fHistPtMesDghterOmegaMinus", "p_{t} of gen. Meson #Lambda dghter); pt_{MesDght}; Number of Mes.", 200, 0, 10);
852       fListHistCascade->Add(fHistPtMesDghterOmegaMinus);
853    }
854    if (!fHistPtBarDghterOmegaMinus) {
855       fHistPtBarDghterOmegaMinus = new TH1F("fHistPtBarDghterOmegaMinus", "p_{t} of gen. Baryon #Lambda dghter); pt_{BarDght}; Number of Bar.", 200, 0, 10);
856       fListHistCascade->Add(fHistPtBarDghterOmegaMinus);
857    }
858    // - Omega+ 
859    // - Pseudo-Rapidity distribution
860    if (!fHistEtaGenCascOmegaPlus) {
861       fHistEtaGenCascOmegaPlus = new TH1F("fHistEtaGenCascOmegaPlus", "#eta of any gen. #Omega^{+}; #eta; Number of Casc", 200, -10, 10);
862       fListHistCascade->Add(fHistEtaGenCascOmegaPlus);
863    }
864    if (!f3dHistGenPtVsGenYvsNtracksOmegaPlus) {
865       f3dHistGenPtVsGenYvsNtracksOmegaPlus = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaPlus", "MC P_{t} Vs MC Y of Gen #Omega^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
866       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaPlus);
867    }
868    if (!f3dHistGenPtVsGenctauvsYOmegaPlus) {
869       f3dHistGenPtVsGenctauvsYOmegaPlus = new TH3D("f3dHistGenPtVsGenctauvsYOmegaPlus", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{+} ", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
870       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaPlus);
871    }
872    if (!f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff) {
873       f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff = new TH3D("f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff", "MC P_{t} Vs MC Y of Gen #Omega^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 110, -1.1, 1.1, 200, 0., 200.);
874       fListHistCascade->Add(f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff);
875    }
876    if (!f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff) {
877       f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff = new TH3D("f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff", "MC P_{t} Vs MC ctau Vs Y of Gen #Omega^{+}", 100, ptBinLimits, 111, ctauBinLimits, 110, yBinLimits);
878       fListHistCascade->Add(f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff);
879    }
880    // - Info at the generation level of multi-strange particle
881    if (!fHistThetaGenCascOmegaPlus) {
882       fHistThetaGenCascOmegaPlus = new TH1F("fHistThetaGenCascOmegaPlus", "#theta of gen. #Omega^{+}; #theta; Number of Casc.", 200, -10, 190);
883       fListHistCascade->Add(fHistThetaGenCascOmegaPlus);
884    }
885    if (!f2dHistGenPtVsGenYFdblOmegaPlus) {
886       f2dHistGenPtVsGenYFdblOmegaPlus = new TH2D("f2dHistGenPtVsGenYFdblOmegaPlus", "MC P_{t} Vs MC Y of findable Gen #Omega^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
887       fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaPlus);
888    }
889    // - Info at the generation level of multi-strange particle
890    if (!fHistThetaGenCascOmegaPlus) {
891       fHistThetaGenCascOmegaPlus = new TH1F("fHistThetaGenCascOmegaPlus", "#theta of gen. #Omega^{+}; #theta; Number of Casc.", 200, -10, 190);
892       fListHistCascade->Add(fHistThetaGenCascOmegaPlus);
893    }
894    if (!f2dHistGenPtVsGenYFdblOmegaPlus) {
895       f2dHistGenPtVsGenYFdblOmegaPlus = new TH2D("f2dHistGenPtVsGenYFdblOmegaPlus", "MC P_{t} Vs MC Y of findable Gen #Omega^{+}; Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
896       fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaPlus);
897    }
898    // - Theta distribution the daughters (control plots)
899    if (!fHistThetaLambdaOmegaPlus) {
900       fHistThetaLambdaOmegaPlus = new TH1F("fHistThetaLambdaOmegaPlus", "#theta of gen. #Lambda (Omega dghter); #theta_{#Lambda}; Number of #Lambda", 200, -10, 190);
901       fListHistCascade->Add(fHistThetaLambdaOmegaPlus);
902    }
903    if (!fHistThetaBachOmegaPlus) {
904       fHistThetaBachOmegaPlus = new TH1F("fHistThetaBachOmegaPlus", "#theta of gen. Bach.; #theta_{Bach}; Number of Bach.", 200, -10, 190);
905       fListHistCascade->Add(fHistThetaBachOmegaPlus);
906    }
907    if (!fHistThetaMesDghterOmegaPlus) {
908       fHistThetaMesDghterOmegaPlus = new TH1F("fHistThetaMesDghterOmegaPlus", "#theta of gen. Meson #Lambda dghter; #theta_{MesDght}; Number of Mes.", 200, -10, 190);
909       fListHistCascade->Add(fHistThetaMesDghterOmegaPlus);
910    }
911    if (!fHistThetaBarDghterOmegaPlus) {
912       fHistThetaBarDghterOmegaPlus = new TH1F("fHistThetaBarDghterOmegaPlus", "#theta of gen. Baryon #Lambda dghter; #theta_{BarDght}; Number of Bar.", 200, -10, 190);
913       fListHistCascade->Add(fHistThetaBarDghterOmegaPlus);
914    }
915    // - Pt distribution (control plots)
916    if (!fHistPtBachOmegaPlus) {
917       fHistPtBachOmegaPlus = new TH1F("fHistPtBachOmegaPlus", "p_{t} of gen. Bach.; pt_{Bach}; Number of Bach.", 200, 0, 10);
918       fListHistCascade->Add(fHistPtBachOmegaPlus);
919    }
920    if (!fHistPtMesDghterOmegaPlus) {
921       fHistPtMesDghterOmegaPlus = new TH1F("fHistPtMesDghterOmegaPlus", "p_{t} of gen. Meson #Lambda dghter; pt_{MesDght}; Number of Mes.", 200, 0, 10);
922       fListHistCascade->Add(fHistPtMesDghterOmegaPlus);
923    }
924    if (!fHistPtBarDghterOmegaPlus) {
925       fHistPtBarDghterOmegaPlus = new TH1F("fHistPtBarDghterOmegaPlus", "p_{t} of gen. Baryon #Lambda dghter); pt_{BarDght}; Number of Bar.", 200, 0, 10);
926       fListHistCascade->Add(fHistPtBarDghterOmegaPlus);
927    }
928  
929  //-------------------------------------------------------------------------
930  // - Any reconstructed cascades + reconstructed cascades associated with MC
931   
932    // - Multiplicity cascde plots
933    fHistnAssoXiMinus= new TH1F("fHistnAssoXiMinus", "", 25, 0, 25);
934    fListHistCascade->Add(fHistnAssoXiMinus);
935    fHistnAssoXiPlus= new TH1F("fHistnAssoXiPlus", "", 25, 0, 25);
936    fListHistCascade->Add(fHistnAssoXiPlus);
937    fHistnAssoOmegaMinus= new TH1F("fHistnAssoOmegaMinus", "", 25, 0, 25);
938    fListHistCascade->Add(fHistnAssoOmegaMinus);
939    fHistnAssoOmegaPlus= new TH1F("fHistnAssoOmegaPlus", "", 25, 0, 25);
940    fListHistCascade->Add(fHistnAssoOmegaPlus);
941    // - Effective mass histos for cascades candidates. 
942    if (! fHistMassXiMinus) {
943           fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates; M( #Lambda , #pi^{-} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
944           fListHistCascade->Add(fHistMassXiMinus);
945    }
946    if (! fHistMassXiPlus) {
947           fHistMassXiPlus = new TH1F("fHistMassXiPlus","#Xi^{+} candidates; M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
948           fListHistCascade->Add(fHistMassXiPlus);
949    }
950    if (! fHistMassOmegaMinus) {
951           fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates; M( #Lambda , K^{-} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
952           fListHistCascade->Add(fHistMassOmegaMinus);
953    } 
954    if (! fHistMassOmegaPlus) {
955           fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#Omega^{+} candidates; M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
956           fListHistCascade->Add(fHistMassOmegaPlus);
957    }
958    // - Effective mass histos with combined PID
959    if (! fHistMassWithCombPIDXiMinus) {
960       fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID; M( #Lambda , #pi^{-} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
961       fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
962    }
963    if (! fHistMassWithCombPIDXiPlus) {
964       fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#Xi^{+} candidates, with Bach. comb. PID; M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
965       fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
966    }
967    if (! fHistMassWithCombPIDOmegaMinus) {
968       fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID; M( #Lambda , K^{-} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
969       fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
970    }
971    if (! fHistMassWithCombPIDOmegaPlus) {
972       fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#Omega^{+} candidates, with Bach. comb. PID; M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
973       fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
974    }
975    // - PID Probability versus MC Pt(bachelor track)
976    if (! f2dHistPIDprobaKaonVsMCPtBach ){
977       f2dHistPIDprobaKaonVsMCPtBach  = new TH2F("f2dHistPIDprobaKaonVsMCPtBach", "Comb. PID proba to be K^{#pm} Vs MC Bach. Pt; Pt_{MC}(Bach.) (GeV/c); Comb. PID Proba (Bach. = K^{#pm})", 100, 0.0, 5.0, 110, 0.0, 1.10);
978       fListHistCascade->Add(f2dHistPIDprobaKaonVsMCPtBach);
979    }
980    if(! f2dHistPIDprobaPionVsMCPtBach ){
981       f2dHistPIDprobaPionVsMCPtBach  = new TH2F("f2dHistPIDprobaPionVsMCPtBach", "Comb. PID proba to be #pi^{#pm} Vs MC Bach. Pt; Pt_{MC}(Bach.) (GeV/c); Comb. PID Proba (Bach. = #pi^{#pm})", 100, 0.0, 5.0, 110, 0.0, 1.10);
982       fListHistCascade->Add(f2dHistPIDprobaPionVsMCPtBach);
983    }
984    // - Effective mass histos with perfect MC PID on the bachelor
985    if (! fHistMassWithMcPIDXiMinus) {
986       fHistMassWithMcPIDXiMinus = new TH1F("fHistMassWithMcPIDXiMinus", "#Xi^{-} candidates, with Bach. MC PID; M( #Lambda , #pi^{-} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
987       fListHistCascade->Add(fHistMassWithMcPIDXiMinus);
988    }
989    if (! fHistMassWithMcPIDXiPlus) {
990       fHistMassWithMcPIDXiPlus = new TH1F("fHistMassWithMcPIDXiPlus", "#Xi^{+} candidates, with Bach. MC PID; M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
991       fListHistCascade->Add(fHistMassWithMcPIDXiPlus);
992    }
993    if (! fHistMassWithMcPIDOmegaMinus) {
994       fHistMassWithMcPIDOmegaMinus = new TH1F("fHistMassWithMcPIDOmegaMinus", "#Omega^{-} candidates, with Bach. MC PID; M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500, 1.5, 2.5);
995       fListHistCascade->Add(fHistMassWithMcPIDOmegaMinus);
996    }
997    if (! fHistMassWithMcPIDOmegaPlus) {
998       fHistMassWithMcPIDOmegaPlus = new TH1F("fHistMassWithMcPIDOmegaPlus", "#Omega^{+} candidates, with Bach. MC PID; M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
999       fListHistCascade->Add(fHistMassWithMcPIDOmegaPlus);
1000    }
1001    // - Effective mass histos for cascades candidates ASSOCIATED with MC.
1002    if (! fHistAsMCMassXiMinus) {
1003       fHistAsMCMassXiMinus = new TH1F("fHistAsMCMassXiMinus", "#Xi^{-} candidates associated to MC; M( #Lambda , #pi^{-} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
1004       fListHistCascade->Add(fHistAsMCMassXiMinus);
1005    }
1006    if (! fHistAsMCMassXiPlus) {
1007       fHistAsMCMassXiPlus = new TH1F("fHistAsMCMassXiPlus", "#Xi^{+} candidates associated to MC; M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2}); Counts", 400, 1.2, 2.0);
1008       fListHistCascade->Add(fHistAsMCMassXiPlus);
1009    }
1010    if (! fHistAsMCMassOmegaMinus) {
1011       fHistAsMCMassOmegaMinus = new TH1F("fHistAsMCMassOmegaMinus", "#Omega^{-} candidates associated to MC; M( #Lambda , K^{-} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
1012       fListHistCascade->Add(fHistAsMCMassOmegaMinus);
1013    }
1014    if (! fHistAsMCMassOmegaPlus) {
1015       fHistAsMCMassOmegaPlus = new TH1F("fHistAsMCMassOmegaPlus", "#Omega^{+} candidates associated to MC; M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2}); Counts", 500, 1.5, 2.5);
1016       fListHistCascade->Add(fHistAsMCMassOmegaPlus);
1017    }
1018    // -  Generated Pt Vs generated Y of the cascade candidates associated with MC + having the proper maximum proba of combined PID for the bachelor
1019    if (!f2dHistAsMCandCombPIDGenPtVsGenYXiMinus) {
1020       f2dHistAsMCandCombPIDGenPtVsGenYXiMinus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYXiMinus", "MC P_{t} Vs MC Y of #Xi^{-} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1021       fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiMinus);
1022    }
1023    if (!f2dHistAsMCandCombPIDGenPtVsGenYXiPlus) {
1024       f2dHistAsMCandCombPIDGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of #Xi^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 100, 0., 10., 220, -1.1, 1.1);
1025       fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiPlus);
1026    } 
1027    if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus) {
1028       f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus", "MC P_{t} Vs MC Y of #Omega^{-} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1029       fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus);
1030    }
1031    if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus) {
1032       f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of #Omega^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1033       fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus);
1034    }
1035    // - Generated Pt Vs Generated Y, for the cascade candidates associated with MC
1036    if (!f2dHistAsMCGenPtVsGenYXiMinus) {
1037       f2dHistAsMCGenPtVsGenYXiMinus = new TH2F("f2dHistAsMCGenPtVsGenYXiMinus", "MC P_{t} Vs MC Y of gen. #Xi^{-} (associated); Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1038       fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiMinus );
1039    }
1040    if (!f2dHistAsMCGenPtVsGenYXiPlus) {
1041       f2dHistAsMCGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of gen. #Xi^{+} (associated); Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1042       fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiPlus );
1043    }
1044    if (!f2dHistAsMCGenPtVsGenYOmegaMinus) {
1045       f2dHistAsMCGenPtVsGenYOmegaMinus = new TH2F("f2dHistAsMCGenPtVsGenYOmegaMinus", "MC P_{t} Vs MC Y of gen. #Omega^{-} (associated); Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1046       fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaMinus );
1047    }
1048    if (!f2dHistAsMCGenPtVsGenYOmegaPlus) {
1049       f2dHistAsMCGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of gen. #Omega^{+} (associated); Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1050       fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaPlus );
1051    } 
1052    // - Generated Eta of the the cascade candidates associated with MC
1053    if (!fHistAsMCGenEtaXiMinus) {
1054       fHistAsMCGenEtaXiMinus = new TH1F("fHistAsMCGenEtaXiMinus", "#eta of gen. #Xi^{-} (associated); #eta; Count", 100, -5, 5);
1055       fListHistCascade->Add( fHistAsMCGenEtaXiMinus );
1056    }
1057    if (!fHistAsMCGenEtaXiPlus) {
1058       fHistAsMCGenEtaXiPlus = new TH1F("fHistAsMCGenEtaXiPlus", "#eta of gen. #Xi^{+} (associated); #eta; Count", 100, -5, 5);
1059       fListHistCascade->Add( fHistAsMCGenEtaXiPlus );
1060    }
1061    if (!fHistAsMCGenEtaOmegaMinus) {
1062       fHistAsMCGenEtaOmegaMinus = new TH1F("fHistAsMCGenEtaOmegaMinus", "#eta of gen. #Omega^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1063       fListHistCascade->Add( fHistAsMCGenEtaOmegaMinus );
1064    }
1065    if (!fHistAsMCGenEtaOmegaPlus) {
1066       fHistAsMCGenEtaOmegaPlus = new TH1F("fHistAsMCGenEtaOmegaPlus", "#eta of gen. #Omega^{+} (associated); #eta; Count", 100, -5, 5);
1067       fListHistCascade->Add( fHistAsMCGenEtaOmegaPlus );
1068    }
1069    // - Resolution in Pt as function of generated Pt
1070    if (! f2dHistAsMCResPtXiMinus) {
1071       f2dHistAsMCResPtXiMinus = new TH2F("f2dHistAsMCResPtXiMinus", "Resolution in Pt reconstruction for #Xi^{-}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1072       fListHistCascade->Add(f2dHistAsMCResPtXiMinus);
1073    }
1074    if (! f2dHistAsMCResPtXiPlus) {
1075       f2dHistAsMCResPtXiPlus = new TH2F("f2dHistAsMCResPtXiPlus", "Resolution in Pt reconstruction for #Xi^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1076       fListHistCascade->Add(f2dHistAsMCResPtXiPlus);
1077    }
1078    if (! f2dHistAsMCResPtOmegaMinus) {
1079       f2dHistAsMCResPtOmegaMinus = new TH2F("f2dHistAsMCResPtOmegaMinus", "Resolution in Pt reconstruction for #Omega^{-}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1080       fListHistCascade->Add(f2dHistAsMCResPtOmegaMinus);
1081    }
1082    if (! f2dHistAsMCResPtOmegaPlus) {
1083       f2dHistAsMCResPtOmegaPlus = new TH2F("f2dHistAsMCResPtOmegaPlus", "Resolution in Pt reconstruction for #Omega^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1084       fListHistCascade->Add(f2dHistAsMCResPtOmegaPlus);
1085    }
1086    // - Resolution in R(2D) as function of generated R
1087    if (! f2dHistAsMCResRXiMinus) {
1088       f2dHistAsMCResRXiMinus = new TH2F("f2dHistAsMCResRXiMinus", "Resolution in transv. position for #Xi^{-}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1089       fListHistCascade->Add(f2dHistAsMCResRXiMinus);
1090    }
1091    if (! f2dHistAsMCResRXiPlus) {
1092       f2dHistAsMCResRXiPlus = new TH2F("f2dHistAsMCResRXiPlus", "Resolution in transv. position for #Xi^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1093       fListHistCascade->Add(f2dHistAsMCResRXiPlus);
1094    } 
1095    if (! f2dHistAsMCResROmegaMinus) {
1096       f2dHistAsMCResROmegaMinus = new TH2F("f2dHistAsMCResROmegaMinus", "Resolution in transv. position for #Omega^{-}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1097       fListHistCascade->Add(f2dHistAsMCResROmegaMinus);
1098    }
1099    if (! f2dHistAsMCResROmegaPlus) {
1100       f2dHistAsMCResROmegaPlus = new TH2F("f2dHistAsMCResROmegaPlus", "Resolution in transv. position for #Omega^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1101       fListHistCascade->Add(f2dHistAsMCResROmegaPlus);
1102    }
1103    // - Resolution in phi as function of generated Pt 
1104    if (! f2dHistAsMCResPhiXiMinus) {
1105       f2dHistAsMCResPhiXiMinus = new TH2F("f2dHistAsMCResPhiXiMinus", "Resolution in #phi for #Xi^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1106       fListHistCascade->Add(f2dHistAsMCResPhiXiMinus);
1107    }
1108    if (! f2dHistAsMCResPhiXiPlus) {
1109       f2dHistAsMCResPhiXiPlus = new TH2F("f2dHistAsMCResPhiXiPlus", "Resolution in #phi for #Xi^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1110       fListHistCascade->Add(f2dHistAsMCResPhiXiPlus);
1111    }
1112    if (! f2dHistAsMCResPhiOmegaMinus) {
1113       f2dHistAsMCResPhiOmegaMinus = new TH2F("f2dHistAsMCResPhiOmegaMinus", "Resolution in #phi for #Omega^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);  
1114       fListHistCascade->Add(f2dHistAsMCResPhiOmegaMinus);
1115    }
1116    if (! f2dHistAsMCResPhiOmegaPlus) {
1117       f2dHistAsMCResPhiOmegaPlus = new TH2F("f2dHistAsMCResPhiOmegaPlus", "Resolution in #phi for #Omega^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1118       fListHistCascade->Add(f2dHistAsMCResPhiOmegaPlus);
1119    }
1120    //  - Correlation between proton (antiproton) daughter MC pt and Xi/Omega MC pt (to apply Geant/Fluka correction)
1121    if (!f2dHistAsMCptProtonMCptXiMinus) {
1122       f2dHistAsMCptProtonMCptXiMinus = new TH2F("f2dHistAsMCptProtonMCptXiMinus", "Proton MC pt vs Xi- MC pt", 100, 0., 10., 100, 0., 10.); 
1123       fListHistCascade->Add(f2dHistAsMCptProtonMCptXiMinus);
1124    }
1125    if (!f2dHistAsMCptAntiprotonMCptXiPlus) {
1126       f2dHistAsMCptAntiprotonMCptXiPlus = new TH2F("f2dHistAsMCptAntiprotonMCptXiPlus", "Antiproton MC pt vs Xi+ MC pt", 100, 0., 10., 100, 0., 10.);
1127       fListHistCascade->Add(f2dHistAsMCptAntiprotonMCptXiPlus);
1128    }
1129    if (!f2dHistAsMCptProtonMCptOmegaMinus) {
1130       f2dHistAsMCptProtonMCptOmegaMinus = new TH2F("f2dHistAsMCptProtonMCptOmegaMinus", "Proton MC pt vs Omega- MC pt", 100, 0., 10., 100, 0., 10.);
1131       fListHistCascade->Add(f2dHistAsMCptProtonMCptOmegaMinus);
1132    }
1133    if (!f2dHistAsMCptAntiprotonMCptOmegaPlus) {
1134       f2dHistAsMCptAntiprotonMCptOmegaPlus = new TH2F("f2dHistAsMCptAntiprotonMCptOmegaPlus", "Antiproton MC pt vs Omega+ MC pt", 100, 0., 10., 100, 0., 10.);
1135       fListHistCascade->Add(f2dHistAsMCptAntiprotonMCptOmegaPlus);
1136    }
1137    // - Cosine of Pointing angle
1138    if (! fHistV0toXiCosineOfPointingAngle) {
1139       fHistV0toXiCosineOfPointingAngle = new TH1F("fHistV0toXiCosineOfPointingAngle", "Cos. of V0 Ptng Angl / Xi vtx ; Cos(V0 Point. Angl / Xi vtx); Counts", 200, 0.95, 1.0001);
1140       fListHistCascade->Add(fHistV0toXiCosineOfPointingAngle);
1141    }
1142    if (! fHistV0CosineOfPointingAnglevsPtXi) {
1143       fHistV0CosineOfPointingAnglevsPtXi = new TH2F("fHistV0CosineOfPointingAnglevsPtXi", "Cos. of V0 Ptng Angl vs cascade Pt; Cos(V0 Point. Angl); Counts", 100, 0., 10., 200, 0.95, 1.0001);
1144       fListHistCascade->Add(fHistV0CosineOfPointingAnglevsPtXi);
1145    }
1146    if (! fHistV0CosineOfPointingAnglevsPtOmega) {
1147       fHistV0CosineOfPointingAnglevsPtOmega = new TH2F("fHistV0CosineOfPointingAnglevsPtOmega", "Cos. of V0 Ptng Angl vs cascade Pt; Cos(V0 Point. Angl); Counts", 100, 0., 10., 200, 0.95, 1.0001);
1148       fListHistCascade->Add(fHistV0CosineOfPointingAnglevsPtOmega);
1149    }
1150
1151   //--------------
1152   // - CFContainer
1153   // PID container Xi-
1154   if(! fCFContCascadePIDAsXiMinus)  {
1155      const Int_t  lNbSteps      =  7;
1156      const Int_t  lNbVariables  =  3;
1157        //Array for the number of bins in each dimension:
1158      Int_t lNbBinsPerVar[3] = {0};
1159      lNbBinsPerVar[0] = 100;
1160      lNbBinsPerVar[1] = 800;
1161      lNbBinsPerVar[2] = 22;
1162      fCFContCascadePIDAsXiMinus = new AliCFContainer(Form("fCFContCascadePIDAsXiMinus_minnTPCcls%i_vtxlim%.1f-%.1f_minptdghtrk%.1f_etacutdghtrk%.1f",fMinnTPCcls,fVtxRange,fVtxRangeMin,fMinPtCutOnDaughterTracks,fEtaCutOnDaughterTracks),"Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1163        //Setting the bin limits 
1164      fCFContCascadePIDAsXiMinus->SetBinLimits(0,   0.0  ,  10.0 );      // Pt(Cascade)
1165      fCFContCascadePIDAsXiMinus->SetBinLimits(1,   1.2  ,   2.0 );      // Xi Effective mass
1166      fCFContCascadePIDAsXiMinus->SetBinLimits(2,  -1.1  ,   1.1 );      // Rapidity
1167        //Setting the step title : one per PID case
1168      fCFContCascadePIDAsXiMinus->SetStepTitle(0, "No PID");
1169      fCFContCascadePIDAsXiMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1170      fCFContCascadePIDAsXiMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1171      fCFContCascadePIDAsXiMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1172      fCFContCascadePIDAsXiMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1173      fCFContCascadePIDAsXiMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1174      fCFContCascadePIDAsXiMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1175        //Setting the variable title, per axis
1176      fCFContCascadePIDAsXiMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1177      fCFContCascadePIDAsXiMinus->SetVarTitle(1, "M( #Lambda , #pi^{-} ) (GeV/c^{2})");
1178      fCFContCascadePIDAsXiMinus->SetVarTitle(2, "Y_{#Xi}");
1179      fListHistCascade->Add(fCFContCascadePIDAsXiMinus);  
1180   }
1181   // PID container Xi+
1182   if(! fCFContCascadePIDAsXiPlus)  {
1183      const Int_t  lNbSteps      =  7;
1184      const Int_t  lNbVariables  =  3;
1185        //Array for the number of bins in each dimension :
1186      Int_t lNbBinsPerVar[3] = {0};
1187      lNbBinsPerVar[0] = 100;
1188      lNbBinsPerVar[1] = 800;
1189      lNbBinsPerVar[2] = 22;
1190      fCFContCascadePIDAsXiPlus = new AliCFContainer(Form("fCFContCascadePIDAsXiPlus_minnTPCcls%i_vtxlim%.1f-%.1f_minptdghtrk%.1f_etacutdghtrk%.1f",fMinnTPCcls,fVtxRange,fVtxRangeMin,fMinPtCutOnDaughterTracks,fEtaCutOnDaughterTracks),"Pt_{cascade} Vs M_{#Xi^{+} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1191        //Setting the bin limits (valid  for v4-18-10-AN)
1192      fCFContCascadePIDAsXiPlus->SetBinLimits(0,   0.0  ,  10.0 );       // Pt(Cascade)
1193      fCFContCascadePIDAsXiPlus->SetBinLimits(1,   1.2  ,   2.0 );       // Xi Effective mass
1194      fCFContCascadePIDAsXiPlus->SetBinLimits(2,  -1.1  ,   1.1 );       // Rapidity 
1195        //Setting the step title : one per PID case
1196      fCFContCascadePIDAsXiPlus->SetStepTitle(0, "No PID");
1197      fCFContCascadePIDAsXiPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1198      fCFContCascadePIDAsXiPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1199      fCFContCascadePIDAsXiPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1200      fCFContCascadePIDAsXiPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1201      fCFContCascadePIDAsXiPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1202      fCFContCascadePIDAsXiPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");  
1203        //Setting the variable title, per axis
1204      fCFContCascadePIDAsXiPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1205      fCFContCascadePIDAsXiPlus->SetVarTitle(1, "M( #Lambda , #pi^{+} ) (GeV/c^{2})");
1206      fCFContCascadePIDAsXiPlus->SetVarTitle(2, "Y_{#Xi}");
1207      fListHistCascade->Add(fCFContCascadePIDAsXiPlus);
1208   }
1209   // PID container Omega-
1210   if(! fCFContCascadePIDAsOmegaMinus)  {
1211      const Int_t  lNbSteps      =  7;
1212      const Int_t  lNbVariables  =  3;
1213        //Array for the number of bins in each dimension :
1214      Int_t lNbBinsPerVar[3] = {0};
1215      lNbBinsPerVar[0] = 100;
1216      lNbBinsPerVar[1] = 1000;
1217      lNbBinsPerVar[2] = 22;
1218      fCFContCascadePIDAsOmegaMinus = new AliCFContainer(Form("fCFContCascadePIDAsOmegaMinus_minnTPCcls%i_vtxlim%.1f-%.1f_minptdghtrk%.1f_etacutdghtrk%.1f",fMinnTPCcls,fVtxRange,fVtxRangeMin,fMinPtCutOnDaughterTracks,fEtaCutOnDaughterTracks),"Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1219        //Setting the bin limits 
1220      fCFContCascadePIDAsOmegaMinus->SetBinLimits(0,   0.0  ,  10.0 );   // Pt(Cascade)
1221      fCFContCascadePIDAsOmegaMinus->SetBinLimits(1,   1.5  ,   2.5 );   // Omega Effective mass
1222      fCFContCascadePIDAsOmegaMinus->SetBinLimits(2,  -1.1  ,   1.1 );   // Rapidity
1223        //Setting the step title : one per PID case
1224      fCFContCascadePIDAsOmegaMinus->SetStepTitle(0, "No PID");
1225      fCFContCascadePIDAsOmegaMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1226      fCFContCascadePIDAsOmegaMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1227      fCFContCascadePIDAsOmegaMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1228      fCFContCascadePIDAsOmegaMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1229      fCFContCascadePIDAsOmegaMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1230      fCFContCascadePIDAsOmegaMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1231        //Setting the variable title, per axis
1232      fCFContCascadePIDAsOmegaMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1233      fCFContCascadePIDAsOmegaMinus->SetVarTitle(1, "M( #Lambda , K^{-} ) (GeV/c^{2})");
1234      fCFContCascadePIDAsOmegaMinus->SetVarTitle(2, "Y_{#Omega}");
1235      fListHistCascade->Add(fCFContCascadePIDAsOmegaMinus);
1236   }
1237   // PID container Omega+
1238   if(! fCFContCascadePIDAsOmegaPlus)  {
1239      const Int_t  lNbSteps      =  7;
1240      const Int_t  lNbVariables  =  3;
1241        //Array for the number of bins in each dimension :
1242      Int_t lNbBinsPerVar[3]= {0};
1243      lNbBinsPerVar[0] = 100;
1244      lNbBinsPerVar[1] = 1000;
1245      lNbBinsPerVar[2] = 22;  
1246      fCFContCascadePIDAsOmegaPlus = new AliCFContainer(Form("fCFContCascadePIDAsOmegaPlus_minnTPCcls%i_vtxlim%.1f-%.1f_minptdghtrk%.1f_etacutdghtrk%.1f",fMinnTPCcls,fVtxRange,fVtxRangeMin,fMinPtCutOnDaughterTracks,fEtaCutOnDaughterTracks),"Pt_{cascade} Vs M_{#Omega^{+} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1247        //Setting the bin limits 
1248      fCFContCascadePIDAsOmegaPlus->SetBinLimits(0,   0.0  ,  10.0 );    // Pt(Cascade)
1249      fCFContCascadePIDAsOmegaPlus->SetBinLimits(1,   1.5  ,   2.5 );    // Omega Effective mass
1250      fCFContCascadePIDAsOmegaPlus->SetBinLimits(2,  -1.1  ,   1.1 );    // Rapidity
1251        //Setting the step title : one per PID case
1252      fCFContCascadePIDAsOmegaPlus->SetStepTitle(0, "No PID");
1253      fCFContCascadePIDAsOmegaPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1254      fCFContCascadePIDAsOmegaPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1255      fCFContCascadePIDAsOmegaPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1256      fCFContCascadePIDAsOmegaPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1257      fCFContCascadePIDAsOmegaPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1258      fCFContCascadePIDAsOmegaPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1259        //Setting the variable title, per axis
1260      fCFContCascadePIDAsOmegaPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1261      fCFContCascadePIDAsOmegaPlus->SetVarTitle(1, "M( #Lambda , K^{+} ) (GeV/c^{2})");
1262      fCFContCascadePIDAsOmegaPlus->SetVarTitle(2, "Y_{#Omega}");
1263      fListHistCascade->Add(fCFContCascadePIDAsOmegaPlus);
1264   }
1265   // Container for optimisation of topological selections 
1266   if(! fCFContAsCascadeCuts){
1267         // Container meant to store all the relevant distributions corresponding to the cut variables.
1268         //          - NB overflow/underflow of variables on which we want to cut later should be 0!!!
1269      const Int_t  lNbSteps      =  4;
1270      const Int_t  lNbVariables  =  19;
1271        //Array for the number of bins in each dimension :
1272      Int_t lNbBinsPerVar[lNbVariables] = {0};
1273      lNbBinsPerVar[0]  = 25;   //DcaCascDaughters                : [0.0,2.,3.0]        -> Rec.Cut = 2.0; 
1274      lNbBinsPerVar[1]  = 25;   //DcaBachToPrimVertex             : [0.0,0.24,100.0]    -> Rec.Cur = 0.01;
1275      lNbBinsPerVar[2]  = 60;     //CascCosineOfPointingAngle    :  [0.94,1.0]          -> Rec.Cut = 0.95;
1276      //lNbBinsPerVar[2]  = 30;   //CascCosineOfPointingAngle       : [0.97,1.]           -> Rec.Cut = 0.98;
1277      lNbBinsPerVar[3]  = 40;   //CascRadius                      : [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
1278      lNbBinsPerVar[4]  = 30;   //InvMassLambdaAsCascDghter       : [1.1,1.3]           -> Rec.Cut = 0.008;
1279      lNbBinsPerVar[5]  = 20;   //DcaV0Daughters                  : [0.0,2.0]           -> Rec.Cut = 1.5;
1280      lNbBinsPerVar[6]  = 201;  //V0CosineOfPointingAngle         : [0.89,1.0]          -> Rec.Cut = 0.9;
1281      lNbBinsPerVar[7]  = 40;   //V0Radius                        : [0.0,3.9,1000.0]    -> Rec.Cut = 0.2;
1282      lNbBinsPerVar[8]  = 40;   //DcaV0ToPrimVertex               : [0.0,0.39,110.0]    -> Rec.Cut = 0.01;
1283      lNbBinsPerVar[9]  = 25;   //DcaPosToPrimVertex              : [0.0,0.24,100.0]    -> Rec.Cut = 0.05;
1284      lNbBinsPerVar[10] = 25;   //DcaNegToPrimVertex              : [0.0,0.24,100.0]    -> Rec.Cut = 0.05;
1285      lNbBinsPerVar[11] = 150;  //InvMassXi                       :  2-MeV/c2 bins
1286      lNbBinsPerVar[12] = 120;  //InvMassOmega                    :  2-MeV/c2 bins
1287      lNbBinsPerVar[13] = 100;  //CascTransvMom                   : [0.0,10.0]
1288      lNbBinsPerVar[14] = 110;  //Y(Xi)                           :  0.02 unit of y per bin 
1289      lNbBinsPerVar[15] = 110;  //Y(Omega)                        :  0.02 unit of y per bin
1290      lNbBinsPerVar[16] = 112;  //Proper lenght of cascade
1291      lNbBinsPerVar[17] = 112;  //Proper lenght of V0 
1292      lNbBinsPerVar[18] = 112;  //Distance V0-Xi in the transverse plane  
1293      fCFContAsCascadeCuts = new AliCFContainer(Form("fCFContAsCascadeCuts_minnTPCcls%i_vtxlim%.1f-%.1f_minptdghtrk%.1f_etacutdghtrk%.1f",fMinnTPCcls,fVtxRange,fVtxRangeMin,fMinPtCutOnDaughterTracks,fEtaCutOnDaughterTracks),"Cut Container for Asso. Cascades", lNbSteps, lNbVariables, lNbBinsPerVar );
1294        //Setting the bin limits 
1295        //0 - DcaCascDaughters
1296      Double_t *lBinLim0 = new Double_t[ lNbBinsPerVar[0]+1 ];
1297      for(Int_t i=0; i<lNbBinsPerVar[0]; i++) lBinLim0[i] = (Double_t)0.0 + (2.4 -0.0)/(lNbBinsPerVar[0] - 1) * (Double_t)i;
1298      lBinLim0[ lNbBinsPerVar[0] ] = 3.0;
1299      fCFContAsCascadeCuts -> SetBinLimits(0, lBinLim0);
1300      delete[] lBinLim0;
1301        //1 - DcaBachToPrimVertex
1302      Double_t *lBinLim1 = new Double_t[ lNbBinsPerVar[1]+1 ];
1303      for(Int_t i=0; i<lNbBinsPerVar[1]; i++) lBinLim1[i] = (Double_t)0.0 + (0.24 - 0.0)/(lNbBinsPerVar[1] - 1) * (Double_t)i;
1304      lBinLim1[ lNbBinsPerVar[1] ] = 100.0;
1305      fCFContAsCascadeCuts -> SetBinLimits(1, lBinLim1);
1306      delete [] lBinLim1;
1307        //2 - CascCosineOfPointingAngle
1308      fCFContAsCascadeCuts -> SetBinLimits(2, .94, 1.);        
1309        //3 - CascRadius
1310      Double_t *lBinLim3 = new Double_t[ lNbBinsPerVar[3]+1 ];
1311      for(Int_t i=0; i<lNbBinsPerVar[3]; i++) lBinLim3[i] = (Double_t)0.0 + (3.9 -0.0)/(lNbBinsPerVar[3] - 1) * (Double_t)i;
1312      lBinLim3[ lNbBinsPerVar[3] ] = 1000.0;
1313      fCFContAsCascadeCuts -> SetBinLimits(3, lBinLim3);
1314      delete[] lBinLim3;
1315        //4 - InvMassLambdaAsCascDghter
1316      fCFContAsCascadeCuts->SetBinLimits(4, 1.1, 1.13); 
1317        //5 - DcaV0Daughters
1318      fCFContAsCascadeCuts->SetBinLimits(5, 0., 2.);        
1319        //6 - V0CosineOfPointingAngle
1320      fCFContAsCascadeCuts->SetBinLimits(6, 0.8, 1.001);
1321        //7 - V0Radius
1322      Double_t *lBinLim7 = new Double_t[ lNbBinsPerVar[7]+1 ];
1323      for(Int_t i=0; i<lNbBinsPerVar[7]; i++) lBinLim7[i] = (Double_t)0.0 + (3.9 - 0.0)/(lNbBinsPerVar[7] - 1) * (Double_t)i ;
1324      lBinLim7[ lNbBinsPerVar[7] ] = 1000.0;
1325      fCFContAsCascadeCuts -> SetBinLimits(7, lBinLim7);
1326      delete [] lBinLim7;      
1327        //8 - DcaV0ToPrimVertexXi : 0. to 0.4 
1328      Double_t *lBinLim8 = new Double_t[ lNbBinsPerVar[8]+1 ];
1329      for(Int_t i=0; i<lNbBinsPerVar[8]; i++) lBinLim8[i] = (Double_t)0.0 + (0.39 - 0.0)/(lNbBinsPerVar[8] - 1) * (Double_t)i ;
1330      lBinLim8[ lNbBinsPerVar[8] ] = 100.0;
1331      fCFContAsCascadeCuts -> SetBinLimits(8, lBinLim8);
1332      delete [] lBinLim8;      
1333        //9 - DcaPosToPrimVertexXi
1334      Double_t *lBinLim9 = new Double_t[ lNbBinsPerVar[9]+1 ];
1335      for(Int_t i=0; i<lNbBinsPerVar[9]; i++) lBinLim9[i] = (Double_t)0.0 + (0.24 - 0.0)/(lNbBinsPerVar[9] - 1) * (Double_t)i;
1336      lBinLim9[ lNbBinsPerVar[9] ] = 100.0;
1337      fCFContAsCascadeCuts -> SetBinLimits(9, lBinLim9);
1338      delete [] lBinLim9;   
1339        //10 - DcaNegToPrimVertexXi
1340      Double_t *lBinLim10 = new Double_t[ lNbBinsPerVar[10]+1 ];
1341      for(Int_t i=0; i<lNbBinsPerVar[10]; i++) lBinLim10[i] = (Double_t)0.0 + (0.24 - 0.0 )/(lNbBinsPerVar[10] - 1) * (Double_t)i;
1342      lBinLim10[ lNbBinsPerVar[10] ] = 100.0;
1343      fCFContAsCascadeCuts -> SetBinLimits(10, lBinLim10);
1344      delete [] lBinLim10;  
1345        //11 - InvMassXi
1346      fCFContAsCascadeCuts -> SetBinLimits(11, 1.25, 1.40);
1347        //12 - InvMassOmega
1348      fCFContAsCascadeCuts -> SetBinLimits(12, 1.62, 1.74);
1349        //13 - XiTransvMom 
1350      fCFContAsCascadeCuts -> SetBinLimits(13, 0.0, 10.0);
1351        //14 - Y(Xi) 
1352      fCFContAsCascadeCuts -> SetBinLimits(14, -1.1, 1.1);
1353        //15 - Y(Omega)
1354      fCFContAsCascadeCuts -> SetBinLimits(15, -1.1, 1.1); 
1355        //16 - Proper time cascade 
1356      Double_t *lBinLim16 = new Double_t[ lNbBinsPerVar[16]+1 ];
1357      for(Int_t i=0; i<lNbBinsPerVar[16]; i++) lBinLim16[i] = (Double_t)-1. + (110. + 1.0 )/(lNbBinsPerVar[16] - 1) * (Double_t)i;
1358      lBinLim16[ lNbBinsPerVar[16] ] = 2000.0;
1359      fCFContAsCascadeCuts -> SetBinLimits(16, lBinLim16);
1360        //17 - Proper time V0 
1361      fCFContAsCascadeCuts -> SetBinLimits(17, lBinLim16);
1362        //18 - Distance V0-Xi in the transverse plane
1363      fCFContAsCascadeCuts -> SetBinLimits(18, lBinLim16);
1364      delete [] lBinLim16;
1365        // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
1366      fCFContAsCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates associated to MC");
1367      fCFContAsCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates associated to MC");
1368      fCFContAsCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates associated to MC");
1369      fCFContAsCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates associated to MC");
1370        // Setting the variable title, per axis
1371      fCFContAsCascadeCuts->SetVarTitle(0,  "DCA(cascade daughters) (cm)");
1372      fCFContAsCascadeCuts->SetVarTitle(1,  "ImpactParamToPV(bachelor) (cm)");
1373      fCFContAsCascadeCuts->SetVarTitle(2,  "cos(cascade PA)");
1374      fCFContAsCascadeCuts->SetVarTitle(3,  "R_{2d}(cascade decay) (cm)");
1375      fCFContAsCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(as casc dghter) (GeV/c^{2})");
1376      fCFContAsCascadeCuts->SetVarTitle(5,  "DCA(V0 daughters) (cm)");
1377      fCFContAsCascadeCuts->SetVarTitle(6,  "cos(V0 PA) in cascade");
1378      fCFContAsCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
1379      fCFContAsCascadeCuts->SetVarTitle(8,  "ImpactParamToPV(V0) (cm)");
1380      fCFContAsCascadeCuts->SetVarTitle(9,  "ImpactParamToPV(Pos) (cm)");
1381      fCFContAsCascadeCuts->SetVarTitle(10, "ImpactParamToPV(Neg) (cm)");
1382      fCFContAsCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
1383      fCFContAsCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
1384      fCFContAsCascadeCuts->SetVarTitle(13, "Pt_{MC}(cascade) (GeV/c)");
1385      fCFContAsCascadeCuts->SetVarTitle(14, "Y_{MC}(Xi)");
1386      fCFContAsCascadeCuts->SetVarTitle(15, "Y_{MC}(Omega)");
1387      fCFContAsCascadeCuts->SetVarTitle(16, "mL/p cascade (cm)");
1388      fCFContAsCascadeCuts->SetVarTitle(17, "mL/p V0 (cm)"); 
1389      fCFContAsCascadeCuts->SetVarTitle(18, "Distance V0-Cascade in the transverse plane (cm)");
1390      fListHistCascade->Add(fCFContAsCascadeCuts);
1391   }
1392
1393  PostData(1, fListHistCascade); 
1394  PostData(2, fCFContCascadePIDAsXiMinus);
1395  PostData(3, fCFContCascadePIDAsXiPlus);
1396  PostData(4, fCFContCascadePIDAsOmegaMinus);
1397  PostData(5, fCFContCascadePIDAsOmegaPlus);
1398  PostData(6, fCFContAsCascadeCuts);
1399
1400 }// end CreateOutputObjects
1401
1402
1403 //________________________________________________________________________
1404 void AliAnalysisTaskCheckPerformanceCascadepp276::UserExec(Option_t *) {
1405         
1406   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1407   // Main loop (called for each event)
1408   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1409    
1410    //------------------
1411    // - Define variables
1412    AliESDEvent  *lESDevent = 0x0;
1413    AliAODEvent  *lAODevent = 0x0;
1414    AliMCEvent   *lMCevent  = 0x0; 
1415    AliStack     *lMCstack  = 0x0; 
1416    TClonesArray *arrayMC = 0;
1417
1418    //-------------------------
1419    // - Check the PID response
1420    if (!fPIDResponse) {
1421         AliError("Cannot get pid response");
1422         return;
1423    }
1424
1425
1426   //////////////////    
1427   // Event selection    
1428   //////////////////
1429   // In order:
1430   // 1) SDD selection
1431   // 2) Physics selection
1432   // 3) Select only looking at events with well-established PV
1433   // 4) Pileup selection
1434   // 5) |Z| < 10 cm
1435     
1436    //---------------------------------------------------------
1437    // Load the InputEvent and check it (for the ESD and AOD)
1438    if (fAnalysisType == "ESD") {
1439        lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1440        if (!lESDevent) {
1441            Printf("ERROR: lESDevent not available \n");
1442            cout << "Name of the file with pb :" <<  CurrentFileName() << endl;  
1443            return;
1444        }
1445        lMCevent = MCEvent();
1446        if (!lMCevent) {
1447            Printf("ERROR: Could not retrieve MC event \n");
1448            cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1449            return;
1450        }
1451        lMCstack = lMCevent->Stack();
1452        if (!lMCstack) {
1453            Printf("ERROR: Could not retrieve MC stack \n");
1454            cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1455            return;
1456        }
1457        // - Cascade vertexer (ESD)
1458        // Relaunch V0 and Cascade vertexer
1459        if (fkRerunV0CascVertexers) { 
1460            lESDevent->ResetCascades();
1461            lESDevent->ResetV0s();
1462            AliV0vertexer *lV0vtxer = new AliV0vertexer();
1463            AliCascadeVertexer *lCascVtxer = new AliCascadeVertexer();
1464            //lV0vtxer->GetCuts(fV0Sels);
1465            //lCascVtxer->GetCuts(fCascSels);
1466            lV0vtxer->SetCuts(fV0Sels);      // NB don't use SetDefaultCuts!! because it acts on static variables 
1467            lCascVtxer->SetCuts(fCascSels);
1468            lV0vtxer->Tracks2V0vertices(lESDevent);
1469            lCascVtxer->V0sTracks2CascadeVertices(lESDevent);
1470            //delete lV0vtxer;
1471            //delete lCascVtxer; 
1472            //---
1473            //lESDevent->ResetCascades();
1474            //lESDevent->ResetV0s();
1475            //AliV0vertexer lV0vtxer;
1476            //AliCascadeVertexer lCascVtxer;
1477            //lV0vtxer.SetCuts(fV0Sels);
1478            //lCascVtxer.SetCuts(fCascSels);
1479            //lV0vtxer.Tracks2V0vertices(lESDevent);
1480            //lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
1481        }
1482    } else if (fAnalysisType == "AOD") {  
1483        lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() ); 
1484        if (!lAODevent) {
1485            Printf("ERROR: lAODevent not available \n");
1486            cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1487            return;
1488        }
1489        arrayMC = (TClonesArray*) lAODevent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1490        if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
1491    } else {
1492        Printf("Analysis type (ESD or AOD) not specified \n");
1493        return;
1494    }
1495
1496    //------------------------------
1497    // - Plots Before any selections
1498    //------------------------------
1499    // - Define variables
1500    Int_t  ncascadesBeforeAnySel          = -1; //number of cascades before any selections
1501    Int_t  nTrackMultiplicityBeforeAnySel = -1; //number of tracks before any selections
1502    if (fAnalysisType == "ESD") {
1503        //Multiplicity
1504        Int_t lMultiplicity = -100;
1505        lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC, 0.5);
1506        nTrackMultiplicityBeforeAnySel = lMultiplicity;
1507        ncascadesBeforeAnySel = lESDevent->GetNumberOfCascades();
1508    } else if (fAnalysisType == "AOD") {
1509        //Multiplicity
1510        Int_t lMultiplicity = -100;
1511        nTrackMultiplicityBeforeAnySel = lMultiplicity;
1512        ncascadesBeforeAnySel  = lAODevent->GetNumberOfCascades();
1513    }
1514    fHistTrackMultiplicityBeforeAnySel->Fill(nTrackMultiplicityBeforeAnySel);
1515    fHistCascadeMultiplicityBeforeAnySel->Fill(ncascadesBeforeAnySel);
1516
1517    //----------------
1518    // - SDD selection
1519    //----------------
1520    // - Define variables
1521    Int_t  ncascadesAfterSDDSel          = -1; //number of cascades after SDD selection
1522    Int_t  nTrackMultiplicityAfterSDDSel = -1; //number of tracks after SDD selection
1523    if (fkSDDselectionOn) {
1524         TString trcl = " ";
1525         trcl = lESDevent->GetFiredTriggerClasses();
1526         if (fAnalysisType == "ESD") trcl = lESDevent->GetFiredTriggerClasses();
1527         else if (fAnalysisType == "AOD") trcl = lAODevent->GetFiredTriggerClasses();
1528         if (fwithSDD){   // ---> Select event with SDD ON
1529             if(!(trcl.Contains("ALLNOTRD"))) {
1530                  PostData(1, fListHistCascade);
1531                  PostData(2, fCFContCascadePIDAsXiMinus);
1532                  PostData(3, fCFContCascadePIDAsXiPlus);
1533                  PostData(4, fCFContCascadePIDAsOmegaMinus);
1534                  PostData(5, fCFContCascadePIDAsOmegaPlus);
1535                  PostData(6, fCFContAsCascadeCuts);
1536                  cout<<"Bad event: SDD turn OFF =>  RETURN!! (Exclude it)..."<<endl;
1537                  return;
1538             } else {
1539                  cout<<"Good event: SDD turn ON."<<endl;
1540             }
1541         } else if (!fwithSDD){  // ---> Select event with SDD OFF
1542             if((trcl.Contains("ALLNOTRD"))) {
1543                  PostData(1, fListHistCascade);
1544                  PostData(2, fCFContCascadePIDAsXiMinus);
1545                  PostData(3, fCFContCascadePIDAsXiPlus);
1546                  PostData(4, fCFContCascadePIDAsOmegaMinus);
1547                  PostData(5, fCFContCascadePIDAsOmegaPlus);
1548                  PostData(6, fCFContAsCascadeCuts);
1549                  cout<<"Bad event:  SDD turn ON =>  RETURN!! (Exclude it)..."<<endl;
1550                  return;
1551             } else {
1552                  cout<<"Good event: SDD turn OFF."<<endl;
1553             }
1554         }
1555    }
1556    // - Take the number of cascades and tracks after the SDD selection
1557    if (fAnalysisType == "ESD") {
1558        Int_t lMultiplicity = -100;
1559        lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC, 0.5);
1560        ncascadesAfterSDDSel = lESDevent->GetNumberOfCascades();
1561        nTrackMultiplicityAfterSDDSel = lMultiplicity;
1562    } else if (fAnalysisType == "AOD") {
1563        Int_t lMultiplicity = -100;
1564        ncascadesAfterSDDSel = lAODevent->GetNumberOfCascades();
1565        nTrackMultiplicityAfterSDDSel = lMultiplicity;
1566    }
1567    // - Fill the plots
1568    fHistTrackMultiplicityAfterSDDSel->Fill(nTrackMultiplicityAfterSDDSel);
1569    fHistCascadeMultiplicityAfterSDDSel->Fill(ncascadesAfterSDDSel);
1570
1571    //------------------------------
1572    // - Plots pre-physics selection
1573    //------------------------------
1574    // - Produce the 3Dhisto for the efficiency denominator
1575    Int_t lNbMCPrimary = 0;
1576    lNbMCPrimary = lMCstack->GetNprimary();
1577
1578    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) {
1579
1580      Double_t partEnergy = 0.;
1581      Double_t partPz     = 0.;
1582      Double_t partP      = 0.;
1583      Double_t partPt     = 0.;
1584      Double_t partVx     = 0.;
1585      Double_t partVy     = 0.;
1586      Double_t partVz     = 0.;
1587      Double_t bacVx      = 0.;
1588      Double_t bacVy      = 0.;
1589      Double_t bacVz      = 0.;
1590      Double_t partMass   = 0.;
1591      Int_t    PDGcode    = 0;
1592      Int_t lPrimaryTrackMultiplicity = nTrackMultiplicityAfterSDDSel;
1593
1594        if ( fAnalysisType == "ESD" ) {
1595             TParticle* lCurrentParticlePrimary = 0x0;
1596             lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );        
1597             if (!lCurrentParticlePrimary) {
1598                   Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1599                   continue;
1600             }
1601             if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue;
1602             TParticle* xiMC = 0x0;
1603             xiMC = lCurrentParticlePrimary;
1604             if (!xiMC) {
1605                   Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
1606                   continue;
1607             }
1608             partEnergy = xiMC->Energy();
1609             partPz     = xiMC->Pz();
1610             partPt     = xiMC->Pt();
1611             partP      = xiMC->P();
1612             partMass   = xiMC->GetMass();
1613             partVx     = xiMC->Vx();
1614             partVy     = xiMC->Vy();
1615             partVz     = xiMC->Vz();
1616             if (xiMC->GetDaughter(0)>=0) {    
1617                  TParticle *mcBach = lMCstack->Particle(xiMC->GetDaughter(0));
1618                  if (mcBach) {
1619                      bacVx  = mcBach->Vx();
1620                      bacVy  = mcBach->Vy();
1621                      bacVz  = mcBach->Vz();
1622                  }
1623             }
1624             PDGcode = lCurrentParticlePrimary->GetPdgCode();
1625        } else if ( fAnalysisType == "AOD" ) {
1626             AliAODMCParticle *lCurrentParticleaod = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
1627             if (!lCurrentParticleaod) {
1628                   Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1629                   continue;
1630             }
1631             if (!lCurrentParticleaod->IsPhysicalPrimary()) continue;
1632             partEnergy = lCurrentParticleaod->E();
1633             partPz     = lCurrentParticleaod->Pz();
1634             partP      = lCurrentParticleaod->P();
1635             partPt     = lCurrentParticleaod->Pt();
1636             partMass   = lCurrentParticleaod->M(); 
1637             partVx     = lCurrentParticleaod->Xv();
1638             partVy     = lCurrentParticleaod->Yv();
1639             partVz     = lCurrentParticleaod->Zv();
1640             if (lCurrentParticleaod->GetDaughter(0)>=0) {
1641                  AliAODMCParticle *mcBach = (AliAODMCParticle*) arrayMC->At(lCurrentParticleaod->GetDaughter(0));
1642                  if (mcBach) {
1643                      bacVx  = mcBach->Xv();
1644                      bacVy  = mcBach->Yv();
1645                      bacVz  = mcBach->Zv();
1646                  }
1647             }     
1648             PDGcode = lCurrentParticleaod->GetPdgCode();
1649        }
1650
1651        // - Calculate rapidity
1652        Double_t lRapXiMC = 0.5*TMath::Log((partEnergy + partPz) / (partEnergy - partPz + 1.e-13));
1653        // - Calculate proper lenght
1654        Double_t lctau = TMath::Sqrt((partVx-bacVx)*(partVx-bacVx)+(partVy-bacVy)*(partVy-bacVy)+(partVz-bacVz)*(partVz-bacVz));
1655        if (partP != 0.) lctau = lctau*partMass/partP;
1656        else lctau = -1.;
1657        // - Fill Histograms
1658        if (PDGcode ==  3312) {
1659            f3dHistGenPtVsGenYvsNtracksXiMinus->Fill(partPt, lRapXiMC, lPrimaryTrackMultiplicity);
1660            f3dHistGenPtVsGenctauvsYXiMinus->Fill(partPt, lctau, lRapXiMC);
1661        }
1662        if (PDGcode == -3312) {
1663            f3dHistGenPtVsGenYvsNtracksXiPlus->Fill(partPt, lRapXiMC, lPrimaryTrackMultiplicity);
1664            f3dHistGenPtVsGenctauvsYXiPlus->Fill(partPt, lctau, lRapXiMC);
1665        }
1666        if (PDGcode ==  3334) {
1667            f3dHistGenPtVsGenYvsNtracksOmegaMinus->Fill(partPt, lRapXiMC, lPrimaryTrackMultiplicity);
1668            f3dHistGenPtVsGenctauvsYOmegaMinus->Fill(partPt, lctau, lRapXiMC);
1669        }
1670        if (PDGcode == -3334) {
1671            f3dHistGenPtVsGenYvsNtracksOmegaPlus->Fill(partPt, lRapXiMC, lPrimaryTrackMultiplicity);
1672            f3dHistGenPtVsGenctauvsYOmegaPlus->Fill(partPt, lctau, lRapXiMC);
1673        }
1674    }
1675
1676  
1677    //--------------------
1678    // - Physics selection
1679    //--------------------
1680    // - Define new variables
1681    Int_t    ncascadesAfterPhysicsSel          = -1; //number of cascades after physics selection
1682    Int_t    nTrackMultiplicityAfterPhysicsSel = -1; //number of tracks after physics selection
1683    // - Selection for ESD and AOD
1684    if (fAnalysisType == "ESD") {
1685        UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1686        Bool_t isSelected = 0;
1687        isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1688        if(!isSelected){
1689            PostData(1, fListHistCascade);
1690            PostData(2, fCFContCascadePIDAsXiMinus);
1691            PostData(3, fCFContCascadePIDAsXiPlus);
1692            PostData(4, fCFContCascadePIDAsOmegaMinus);
1693            PostData(5, fCFContCascadePIDAsOmegaPlus);
1694            PostData(6, fCFContAsCascadeCuts);
1695            return;
1696        }
1697        // - Take the number of cascades and tracks after physics selection
1698        ncascadesAfterPhysicsSel = lESDevent->GetNumberOfCascades();
1699        nTrackMultiplicityAfterPhysicsSel = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC, 0.5);
1700    } else if (fAnalysisType == "AOD") {
1701        UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1702        Bool_t isSelected = 0;
1703        isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1704        if(!isSelected){
1705            PostData(1, fListHistCascade);
1706            PostData(2, fCFContCascadePIDAsXiMinus);
1707            PostData(3, fCFContCascadePIDAsXiPlus);
1708            PostData(4, fCFContCascadePIDAsOmegaMinus);
1709            PostData(5, fCFContCascadePIDAsOmegaPlus);
1710            PostData(6, fCFContAsCascadeCuts);
1711            return;
1712        }
1713        // - Take the number of cascades and tracks after physics selection
1714        ncascadesAfterPhysicsSel = lAODevent->GetNumberOfCascades();
1715        nTrackMultiplicityAfterPhysicsSel = -100;
1716    }
1717    fHistCascadeMultiplicityAfterPhysicsSel->Fill(ncascadesAfterPhysicsSel);
1718    fHistTrackMultiplicityAfterPhysicsSel->Fill(nTrackMultiplicityAfterPhysicsSel);
1719
1720    //-------------------------------------------------------
1721    // Select only looking at events with well-established PV
1722    //-------------------------------------------------------
1723    Int_t    ncascadesForSelEvtNoTPCOnly                  = -1; //number of cascades after the TPConly selection
1724    Int_t    nTrackMultiplicityForSelEvtNoTPCOnly         = -1; //number of tracks after the TPConly selection
1725    if (fAnalysisType == "ESD" ) {
1726        // - Select only looking at events with well-established PV
1727        if (fkQualityCutNoTPConlyPrimVtx) {
1728            const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
1729            const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();
1730            if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingVtx->GetStatus() ){
1731                AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1732                PostData(1, fListHistCascade);
1733                PostData(2, fCFContCascadePIDAsXiMinus);
1734                PostData(3, fCFContCascadePIDAsXiPlus);
1735                PostData(4, fCFContCascadePIDAsOmegaMinus);
1736                PostData(5, fCFContCascadePIDAsOmegaPlus);
1737                PostData(6, fCFContAsCascadeCuts);
1738                return;
1739            }
1740        }
1741        // - Take the number of cascades and tracks after TPConly selection
1742        ncascadesForSelEvtNoTPCOnly = lESDevent->GetNumberOfCascades();
1743        nTrackMultiplicityForSelEvtNoTPCOnly = fESDtrackCuts->GetReferenceMultiplicity(lESDevent,AliESDtrackCuts::kTrackletsITSTPC,0.5);
1744    } else if (fAnalysisType == "AOD") {
1745        // - Select only looking at events with well-established PV
1746        if (fkQualityCutNoTPConlyPrimVtx) {
1747            const AliAODVertex *lPrimarySPDVtx = lAODevent->GetPrimaryVertexSPD();
1748            const AliAODVertex *lPrimaryTrackingAODVtx = lAODevent->GetPrimaryVertex();
1749            if (!lPrimarySPDVtx && !lPrimaryTrackingAODVtx) {
1750                AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1751                PostData(1, fListHistCascade);
1752                PostData(2, fCFContCascadePIDAsXiMinus);
1753                PostData(3, fCFContCascadePIDAsXiPlus);
1754                PostData(4, fCFContCascadePIDAsOmegaMinus);
1755                PostData(5, fCFContCascadePIDAsOmegaPlus);
1756                PostData(6, fCFContAsCascadeCuts);
1757                return;
1758            }
1759        }
1760        // - Take the number of cascades and tracks after TPConly selection
1761        ncascadesForSelEvtNoTPCOnly = lAODevent->GetNumberOfCascades();
1762        nTrackMultiplicityForSelEvtNoTPCOnly = -100;  //FIXME
1763    }
1764    fHistCascadeMultiplicityForSelEvtNoTPCOnly->Fill(ncascadesForSelEvtNoTPCOnly);
1765    fHistTrackMultiplicityForSelEvtNoTPCOnly->Fill(nTrackMultiplicityForSelEvtNoTPCOnly);
1766     
1767    //-----------------
1768    // Pileup selection
1769    //-----------------
1770    Int_t    ncascadesForSelEvtNoTPCOnlyNoPileup          = -1; //number of cascades after the NoPileup selection
1771    Int_t    nTrackMultiplicityForSelEvtNoTPCOnlyNoPileup = -1; //number of tracks after the Pileup selection
1772    if (fAnalysisType == "ESD" ) {
1773        // - Selection for pile up
1774        if (fkRejectEventPileUp) {
1775            if(lESDevent->IsPileupFromSPD()){
1776                AliWarning("Pb / Pile-up event ... return!");
1777                PostData(1, fListHistCascade);
1778                PostData(2, fCFContCascadePIDAsXiMinus);
1779                PostData(3, fCFContCascadePIDAsXiPlus);
1780                PostData(4, fCFContCascadePIDAsOmegaMinus);
1781                PostData(5, fCFContCascadePIDAsOmegaPlus);
1782                PostData(6, fCFContAsCascadeCuts);
1783                return;
1784            }
1785        }
1786        // - Take the number of cascades and tracks after Pileup selection
1787        ncascadesForSelEvtNoTPCOnlyNoPileup = lESDevent->GetNumberOfCascades();
1788        nTrackMultiplicityForSelEvtNoTPCOnlyNoPileup = fESDtrackCuts->GetReferenceMultiplicity(lESDevent,AliESDtrackCuts::kTrackletsITSTPC,0.5);
1789    } else if (fAnalysisType == "AOD") {
1790        // - Selection for pile up
1791        if (fkRejectEventPileUp) {
1792            if(lAODevent->IsPileupFromSPD()){
1793                AliWarning("Pb / Pile-up event ... return!");
1794                PostData(1, fListHistCascade);
1795                PostData(2, fCFContCascadePIDAsXiMinus);
1796                PostData(3, fCFContCascadePIDAsXiPlus);
1797                PostData(4, fCFContCascadePIDAsOmegaMinus);
1798                PostData(5, fCFContCascadePIDAsOmegaPlus);
1799                PostData(6, fCFContAsCascadeCuts);
1800                return;
1801            }
1802        }
1803        // - Take the number of cascades and tracks after Pileup selection
1804        ncascadesForSelEvtNoTPCOnlyNoPileup = lAODevent->GetNumberOfCascades();
1805        nTrackMultiplicityForSelEvtNoTPCOnlyNoPileup = -100;
1806    }
1807    fHistCascadeMultiplicityForSelEvtNoTPCOnlyNoPileup->Fill(ncascadesForSelEvtNoTPCOnlyNoPileup);
1808    fHistTrackMultiplicityForSelEvtNoTPCOnlyNoPileup->Fill(nTrackMultiplicityForSelEvtNoTPCOnlyNoPileup);
1809     
1810    //-------------------
1811    // - Vertex selection 
1812    //-------------------
1813    Int_t    ncascadesAfterVertexSel                      = -1; //number of cascades after vertex selection
1814    Int_t    nTrackMultiplicityAfterVertexSel             = -1; //number of tracks after vertex selection
1815    Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
1816    Double_t tPrimaryVtxPosition[3] = {-100.0, -100.0, -100.0}; 
1817    Double_t lMagneticField = -10.;
1818    if (fAnalysisType == "ESD" ) {
1819        // - Primary vertex definition
1820        const AliESDVertex *lPrimaryBestVtx = lESDevent->GetPrimaryVertex();
1821        if (!lPrimaryBestVtx) {
1822             AliWarning("No prim. vertex in AOD... return!");
1823             PostData(1, fListHistCascade);
1824             PostData(2, fCFContCascadePIDAsXiMinus);
1825             PostData(3, fCFContCascadePIDAsXiPlus);
1826             PostData(4, fCFContCascadePIDAsOmegaMinus);
1827             PostData(5, fCFContCascadePIDAsOmegaPlus);
1828             PostData(6, fCFContAsCascadeCuts);
1829             return;
1830        }
1831        lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1832        // - Vertex position before any event selection on vertex position
1833        const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
1834        tPrimaryVtxPosition[0] = primaryVtx->GetX();
1835        tPrimaryVtxPosition[1] = primaryVtx->GetY();
1836        tPrimaryVtxPosition[2] = primaryVtx->GetZ();
1837        fHistPVx->Fill( tPrimaryVtxPosition[0] );
1838        fHistPVy->Fill( tPrimaryVtxPosition[1] );
1839        fHistPVz->Fill( tPrimaryVtxPosition[2] );
1840        // - Get magnetic filed info
1841        lMagneticField = lESDevent->GetMagneticField();
1842        // - Selection on the primary vertex Z position 
1843        if (fkQualityCutZprimVtxPos) {
1844            if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange || TMath::Abs(lBestPrimaryVtxPos[2]) < fVtxRangeMin) {
1845                 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1846                 PostData(1, fListHistCascade);
1847                 PostData(2, fCFContCascadePIDAsXiMinus);
1848                 PostData(3, fCFContCascadePIDAsXiPlus);
1849                 PostData(4, fCFContCascadePIDAsOmegaMinus);
1850                 PostData(5, fCFContCascadePIDAsOmegaPlus);
1851                 PostData(6, fCFContAsCascadeCuts);
1852                 return;
1853           }
1854        }
1855        // - Take the number of cascades and tracks after vertex Z position selection
1856        ncascadesAfterVertexSel = lESDevent->GetNumberOfCascades();
1857        nTrackMultiplicityAfterVertexSel = fESDtrackCuts->GetReferenceMultiplicity(lESDevent,AliESDtrackCuts::kTrackletsITSTPC,0.5);
1858    } else if (fAnalysisType == "AOD") {
1859        // - Primary vertex definition
1860        const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex(); // get the best primary vertex available for the event GetVertex(0)
1861        if (!lPrimaryBestAODVtx) {
1862             AliWarning("No prim. vertex in AOD... return!");
1863             PostData(1, fListHistCascade);
1864             PostData(2, fCFContCascadePIDAsXiMinus);
1865             PostData(3, fCFContCascadePIDAsXiPlus);
1866             PostData(4, fCFContCascadePIDAsOmegaMinus);
1867             PostData(5, fCFContCascadePIDAsOmegaPlus);
1868             PostData(6, fCFContAsCascadeCuts);
1869             return;
1870        }
1871        lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
1872        // - Vertex position before any event selection on vertex position
1873        const AliVVertex *primaryVtx = lAODevent->GetPrimaryVertex();
1874        tPrimaryVtxPosition[0] = primaryVtx->GetX();
1875        tPrimaryVtxPosition[1] = primaryVtx->GetY();
1876        tPrimaryVtxPosition[2] = primaryVtx->GetZ();
1877        fHistPVx->Fill( tPrimaryVtxPosition[0] );
1878        fHistPVy->Fill( tPrimaryVtxPosition[1] );
1879        fHistPVz->Fill( tPrimaryVtxPosition[2] );
1880        // - Get magnetic filed info
1881        lMagneticField = lAODevent->GetMagneticField();
1882        // - Selection on the primary vertex Z position 
1883        if (fkQualityCutZprimVtxPos) {
1884            if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange && TMath::Abs(lBestPrimaryVtxPos[2]) < fVtxRangeMin) {
1885                 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1886                 PostData(1, fListHistCascade);
1887                 PostData(2, fCFContCascadePIDAsXiMinus);
1888                 PostData(3, fCFContCascadePIDAsXiPlus);
1889                 PostData(4, fCFContCascadePIDAsOmegaMinus);
1890                 PostData(5, fCFContCascadePIDAsOmegaPlus);
1891                 PostData(6, fCFContAsCascadeCuts);
1892                 return;
1893            }
1894        }
1895        // - Take the number of cascades and tracks after vertex Z position selection
1896        ncascadesAfterVertexSel = lAODevent->GetNumberOfCascades();
1897        nTrackMultiplicityAfterVertexSel = -100;
1898    }
1899    // - Fill the plots
1900    fHistCascadeMultiplicityAfterVertexCutSel->Fill(ncascadesAfterVertexSel);
1901    fHistTrackMultiplicityAfterVertexCutSel->Fill(nTrackMultiplicityAfterVertexSel);
1902
1903    // - Vertex position plots: after any event selections
1904    tPrimaryVtxPosition[0] = 0;
1905    tPrimaryVtxPosition[1] = 0;
1906    tPrimaryVtxPosition[2] = 0;
1907    if (fAnalysisType == "ESD" ) {
1908        const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
1909        tPrimaryVtxPosition[0] = primaryVtx->GetX();
1910        tPrimaryVtxPosition[1] = primaryVtx->GetY();
1911        tPrimaryVtxPosition[2] = primaryVtx->GetZ();
1912    } else if (fAnalysisType == "AOD") {
1913        const AliVVertex *primaryVtx = lAODevent->GetPrimaryVertex();
1914        tPrimaryVtxPosition[0] = primaryVtx->GetX();
1915        tPrimaryVtxPosition[1] = primaryVtx->GetY();
1916        tPrimaryVtxPosition[2] = primaryVtx->GetZ();
1917    }
1918    fHistPVxAnalysis->Fill( tPrimaryVtxPosition[0] );
1919    fHistPVyAnalysis->Fill( tPrimaryVtxPosition[1] );
1920    fHistPVzAnalysis->Fill( tPrimaryVtxPosition[2] );
1921     
1922
1923    //----------------------------------------------------------------------     
1924    // - Loop over the different types of GENERATED cascades (Xi-+, Omega-+)     
1925    //----------------------------------------------------------------------
1926    // - Initialisation of useful local variables                
1927    Int_t lPdgCodeCasc            = 0;
1928    Int_t lPdgCodeBach            = 0;
1929    Int_t lPdgCodeLambda          = 0;
1930    Int_t lPdgCodeDghtMesV0       = 0;
1931    Int_t lPdgCodeDghtBarV0       = 0;   
1932    TH1F *lHistEtaGenCasc         = 0;   
1933    TH3D *l3dHistGenPtVsGenYvsNtracksPhysEff = 0;
1934    TH3D *l3dHistGenPtVsGenctauvsYPhysEff    = 0;
1935    TH1F *lHistThetaGenCasc       = 0;
1936    TH2D *l2dHistGenPtVsGenYFdbl  = 0;
1937    TH1F *lHistThetaLambda        = 0;
1938    TH1F *lHistThetaBach          = 0;
1939    TH1F *lHistThetaBarDghter     = 0;
1940    TH1F *lHistThetaMesDghter     = 0;
1941    TH1F *lHistPtBach             = 0;
1942    TH1F *lHistPtBarDghter        = 0;
1943    TH1F *lHistPtMesDghter        = 0;
1944    Int_t ncascperev = 0; 
1945    Int_t ncascperevtot = 0;
1946
1947    for (Int_t iCascType = 1; iCascType < 5; iCascType++) { 
1948          ncascperev = 0;
1949          ncascperevtot = 0;
1950          Int_t lPrimaryTrackMultiplicity = nTrackMultiplicityAfterSDDSel;
1951
1952          switch (iCascType) {
1953            case 1: // Xi-
1954                lPdgCodeCasc       =   3312;  //Xi-
1955                lPdgCodeBach       =   -211;  //Pi-
1956                lPdgCodeLambda     =   3122;  //Lambda0
1957                lPdgCodeDghtMesV0  =   -211;  //Pi-
1958                lPdgCodeDghtBarV0  =   2212;  //Proton   
1959                lHistEtaGenCasc        = fHistEtaGenCascXiMinus;         // this plot for any Xi- 
1960                lHistThetaGenCasc      = fHistThetaGenCascXiMinus;       // cascades generated within acceptance (cut in pt + theta)
1961                l3dHistGenPtVsGenYvsNtracksPhysEff = f3dHistGenPtVsGenYvsNtracksXiMinusPhysEff;
1962                l3dHistGenPtVsGenctauvsYPhysEff    = f3dHistGenPtVsGenctauvsYXiMinusPhysEff;
1963                l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiMinus;
1964                lHistThetaLambda       = fHistThetaLambdaXiMinus;
1965                lHistThetaBach         = fHistThetaBachXiMinus;
1966                lHistThetaBarDghter    = fHistThetaBarDghterXiMinus;
1967                lHistThetaMesDghter    = fHistThetaMesDghterXiMinus;
1968                lHistPtBach            = fHistPtBachXiMinus;
1969                lHistPtBarDghter       = fHistPtBarDghterXiMinus;
1970                lHistPtMesDghter       = fHistPtMesDghterXiMinus;
1971                break; 
1972            case 2: // Xi+
1973                lPdgCodeCasc        =  -3312;  //Xi+
1974                lPdgCodeBach        =    211;  //Pi+
1975                lPdgCodeLambda      =  -3122;  //AntiLambda0
1976                lPdgCodeDghtMesV0   =    211;  //Pi+
1977                lPdgCodeDghtBarV0   =  -2212;  //AntiProton  
1978                lHistEtaGenCasc        = fHistEtaGenCascXiPlus;       // this plot for any Xi+
1979                lHistThetaGenCasc      = fHistThetaGenCascXiPlus;     // cascades generated within acceptance (cut in pt + theta)
1980                l3dHistGenPtVsGenYvsNtracksPhysEff = f3dHistGenPtVsGenYvsNtracksXiPlusPhysEff;
1981                l3dHistGenPtVsGenctauvsYPhysEff    = f3dHistGenPtVsGenctauvsYXiPlusPhysEff;
1982                l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiPlus;
1983                lHistThetaLambda       = fHistThetaLambdaXiPlus;
1984                lHistThetaBach         = fHistThetaBachXiPlus;
1985                lHistThetaBarDghter    = fHistThetaBarDghterXiPlus;
1986                lHistThetaMesDghter    = fHistThetaMesDghterXiPlus;
1987                lHistPtBach            = fHistPtBachXiPlus;
1988                lHistPtBarDghter       = fHistPtBarDghterXiPlus;
1989                lHistPtMesDghter       = fHistPtMesDghterXiPlus;  
1990                break;
1991            case 3: // Omega-
1992                lPdgCodeCasc       =   3334;  //Omega-
1993                lPdgCodeBach       =   -321;  //K-
1994                lPdgCodeLambda     =   3122;  //Lambda0
1995                lPdgCodeDghtMesV0  =   -211;  //Pi-
1996                lPdgCodeDghtBarV0  =   2212;  //Proton
1997                lHistEtaGenCasc        = fHistEtaGenCascOmegaMinus;        // this plot for any Omega+           
1998                lHistThetaGenCasc      = fHistThetaGenCascOmegaMinus;      // cascades generated within acceptance (cut in pt + theta)
1999                l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaMinus;
2000                l3dHistGenPtVsGenYvsNtracksPhysEff = f3dHistGenPtVsGenYvsNtracksOmegaMinusPhysEff;
2001                l3dHistGenPtVsGenctauvsYPhysEff    = f3dHistGenPtVsGenctauvsYOmegaMinusPhysEff;
2002                lHistThetaLambda       = fHistThetaLambdaOmegaMinus;
2003                lHistThetaBach         = fHistThetaBachOmegaMinus;
2004                lHistThetaBarDghter    = fHistThetaBarDghterOmegaMinus;
2005                lHistThetaMesDghter    = fHistThetaMesDghterOmegaMinus;
2006                lHistPtBach            = fHistPtBachOmegaMinus;
2007                lHistPtBarDghter       = fHistPtBarDghterOmegaMinus;
2008                lHistPtMesDghter       = fHistPtMesDghterOmegaMinus;   
2009                break;
2010            case 4:  // Omega+
2011                lPdgCodeCasc       =  -3334;  //Omega+
2012                lPdgCodeBach       =    321;  //K+
2013                lPdgCodeLambda     =  -3122;  //AntiLambda0
2014                lPdgCodeDghtMesV0  =    211;  //Pi+
2015                lPdgCodeDghtBarV0  =  -2212;  //AntiProton 
2016                lHistEtaGenCasc        = fHistEtaGenCascOmegaPlus;        // this plot for any Omega-
2017                lHistThetaGenCasc      = fHistThetaGenCascOmegaPlus;      // cascades generated within acceptance (cut in pt + theta)
2018                l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaPlus;
2019                l3dHistGenPtVsGenYvsNtracksPhysEff = f3dHistGenPtVsGenYvsNtracksOmegaPlusPhysEff;
2020                l3dHistGenPtVsGenctauvsYPhysEff    = f3dHistGenPtVsGenctauvsYOmegaPlusPhysEff;
2021                lHistThetaLambda       = fHistThetaLambdaOmegaPlus;
2022                lHistThetaBach         = fHistThetaBachOmegaPlus;
2023                lHistThetaBarDghter    = fHistThetaBarDghterOmegaPlus;
2024                lHistThetaMesDghter    = fHistThetaMesDghterOmegaPlus;
2025                lHistPtBach            = fHistPtBachOmegaPlus;
2026                lHistPtBarDghter       = fHistPtBarDghterOmegaPlus;
2027                lHistPtMesDghter       = fHistPtMesDghterOmegaPlus;  
2028                break;
2029          }
2030
2031          for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++) {
2032
2033                Double_t partEnergy = 0.;
2034                Double_t partPz     = 0.;
2035                Double_t partEta    = 0.;
2036                Double_t partTheta  = 0.;
2037                Double_t partP      = 0.;
2038                Double_t partPt     = 0.;
2039                Double_t partVx     = 0.;
2040                Double_t partVy     = 0.; 
2041                Double_t partVz     = 0.;
2042                Double_t bacVx      = 0.;
2043                Double_t bacVy      = 0.;
2044                Double_t bacVz      = 0.;    
2045                Double_t partMass   = 0.;
2046
2047                if ( fAnalysisType == "ESD" ) {      
2048                     TParticle* lCurrentParticle = 0x0; 
2049                     lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
2050                     if (!lCurrentParticle) {
2051                         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2052                         continue;
2053                     }
2054                     if (!lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) continue; 
2055                     if (lCurrentParticle->GetPdgCode() == lPdgCodeCasc) {  // Here !    
2056                         TParticle* xiMC = 0x0;
2057                         xiMC = lCurrentParticle;
2058                         if (!xiMC) {
2059                             Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
2060                             continue;
2061                         }
2062                         partEnergy = xiMC->Energy();
2063                         partPz     = xiMC->Pz();
2064                         partEta    = xiMC->Eta();
2065                         partPt     = xiMC->Pt();
2066                         partP      = xiMC->P();
2067                         partTheta  = xiMC->Theta();
2068                         partMass   = xiMC->GetMass();
2069                         partVx     = xiMC->Vx();
2070                         partVy     = xiMC->Vy();
2071                         partVz     = xiMC->Vz();
2072                         if (xiMC->GetDaughter(0)>=0) {
2073                              TParticle *mcBach = lMCstack->Particle(xiMC->GetDaughter(0));
2074                              if (mcBach) {
2075                                  bacVx  = mcBach->Vx();
2076                                  bacVy  = mcBach->Vy();
2077                                  bacVz  = mcBach->Vz();
2078                              }
2079                         }
2080                     } else continue;
2081                } else if ( fAnalysisType == "AOD" ) {
2082                     AliAODMCParticle *lCurrentParticleaod = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2083                     if (!lCurrentParticleaod) {
2084                         Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
2085                         continue;
2086                     }
2087                     if (!lCurrentParticleaod->IsPhysicalPrimary()) continue;  
2088                     if (!(lCurrentParticleaod->PdgCode() == lPdgCodeCasc)) continue;
2089                     partEnergy = lCurrentParticleaod->E();
2090                     partPz     = lCurrentParticleaod->Pz();
2091                     partEta    = lCurrentParticleaod->Eta();
2092                     partP      = lCurrentParticleaod->P();
2093                     partPt     = lCurrentParticleaod->Pt();
2094                     partTheta  = lCurrentParticleaod->Theta();
2095                     partMass   = lCurrentParticleaod->M();   //FIXME: not sure this works, seems not implemented
2096                     partVx     = lCurrentParticleaod->Xv();
2097                     partVy     = lCurrentParticleaod->Yv();
2098                     partVz     = lCurrentParticleaod->Zv();
2099                     if (lCurrentParticleaod->GetDaughter(0)>=0) {
2100                          AliAODMCParticle *mcBach = (AliAODMCParticle*) arrayMC->At(lCurrentParticleaod->GetDaughter(0));
2101                          if (mcBach) {
2102                               bacVx  = mcBach->Xv();
2103                               bacVy  = mcBach->Yv();
2104                               bacVz  = mcBach->Zv();
2105                          } 
2106                     }
2107                }
2108                ncascperevtot++; 
2109                // - Fill the first histos : = any generated Xi, not necessarily within the acceptance
2110                Double_t lRapXiMC = 0.5*TMath::Log((partEnergy + partPz) / (partEnergy - partPz +1.e-13));
2111                // - Calculate proper time
2112                Double_t lctau = TMath::Sqrt((partVx-bacVx)*(partVx-bacVx)+(partVy-bacVy)*(partVy-bacVy)+(partVz-bacVz)*(partVz-bacVz));
2113                if (partP!=0.)    lctau = lctau*partMass/partP;
2114                else lctau = -1.;
2115                Double_t lRadToDeg = 180.0/TMath::Pi();
2116                // - Fill the first histos : = any generated Xi, not necessarily within the acceptance           
2117                lHistEtaGenCasc->Fill( partEta );         
2118                l3dHistGenPtVsGenYvsNtracksPhysEff->Fill( partPt, lRapXiMC, lPrimaryTrackMultiplicity );
2119                l3dHistGenPtVsGenctauvsYPhysEff->Fill( partPt, lctau, lRapXiMC );
2120                lHistThetaGenCasc->Fill( lRadToDeg * partTheta );
2121
2122                //--------------------------------------------------------------------------------------------
2123                // - Check the emission of particle stays within the acceptance of the detector (cut in theta)
2124                if (fApplyAccCut) { if( partTheta < TMath::Pi()/4.0 || partTheta > 3.0*TMath::Pi()/4.0 ) continue;}      
2125
2126                Float_t lambdaTheta = 0.;
2127                Float_t bacTheta    = 0.;
2128                Float_t dghtBarV0Theta = 0.;
2129                Float_t dghtMesV0Theta = 0.;
2130                Float_t bacPt       = 0.;
2131                Float_t dghtBarV0Pt = 0.;
2132                Float_t dghtMesV0Pt = 0.;
2133
2134                if ( fAnalysisType == "ESD" ) { 
2135                     TParticle* xiMC = lMCstack->Particle( iCurrentLabelStack ); 
2136                     if ( xiMC->GetNDaughters() != 2) continue;
2137                     if ( xiMC->GetDaughter(0) < 0 )  continue;
2138                     if ( xiMC->GetDaughter(1) < 0 )  continue;  
2139                     TParticle* lDght0ofXi = lMCstack->Particle(  xiMC->GetDaughter(0) );
2140                     TParticle* lDght1ofXi = lMCstack->Particle(  xiMC->GetDaughter(1) );
2141                     TParticle* lLambda = 0;
2142                     TParticle* lBach   = 0;
2143
2144                     // Xi - Case 1
2145                     if ( lDght0ofXi->GetPdgCode() == lPdgCodeLambda && lDght1ofXi->GetPdgCode() == lPdgCodeBach ){                      
2146                            lLambda = lDght0ofXi;   // dghter0 = Lambda
2147                            lBach   = lDght1ofXi;   // dghter1 = Pi-
2148                     }           
2149                     // Xi - Case 2
2150                     else if ( lDght0ofXi->GetPdgCode() == lPdgCodeBach && lDght1ofXi->GetPdgCode() == lPdgCodeLambda ){ 
2151                            lBach   = lDght0ofXi;   // dghter0 = Pi-
2152                            lLambda = lDght1ofXi;   // dghter1 = Lambda
2153                     }
2154                     // Otherwise - Case 3       
2155                     else continue;
2156
2157                     // - Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2158                     if (fApplyAccCut) { 
2159                          if( lLambda->Theta() < TMath::Pi()/4.0 || lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2160                          if( lBach->Theta() < TMath::Pi()/4.0   || lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2161                          if( lBach->Pt() < 0.150 ) continue; //FIXME: maybe tuned for Xi but not for K- from Omega ...
2162                     }                           
2163  
2164                     //---------
2165                     // - V0 level
2166                     TParticle* lDghtBarV0 = 0;
2167                     TParticle* lDghtMesV0 = 0;
2168                     if( lLambda->GetNDaughters() != 2 )  continue;
2169                     if( lLambda->GetDaughter(0) < 0 )    continue;
2170                     if( lLambda->GetDaughter(1) < 0 )    continue;
2171                     TParticle* lDght0ofLambda = lMCstack->Particle(  lLambda->GetDaughter(0) );
2172                     TParticle* lDght1ofLambda = lMCstack->Particle(  lLambda->GetDaughter(1) );
2173
2174                     // V0 - Case 1
2175                     if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 && lDght1ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 ) {    // Here !                        
2176                            lDghtBarV0 = lDght0ofLambda;   // dghter0 = Proton
2177                            lDghtMesV0 = lDght1ofLambda;   // dghter1 = Pi-
2178                     }           
2179                     // V0 - Case 2
2180                     else if ( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 && lDght1ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 ) {      // Here !
2181                            lDghtMesV0 = lDght0ofLambda;   // dghter0 = Pi-
2182                            lDghtBarV0 = lDght1ofLambda;   // dghter1 = Proton
2183                     }           
2184                     // Otherwise - Case 3
2185                     else continue;
2186                                 
2187                     // - Check the emission of particle stays within the acceptance of the detector
2188                     if (fApplyAccCut) { 
2189                          if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2190                          if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2191                          if( lDghtBarV0->Pt() < 0.250 ) continue;
2192                          if( lDghtMesV0->Pt() < 0.150 ) continue;
2193                     }
2194  
2195                     lambdaTheta    = lLambda->Theta();
2196                     bacTheta       = lBach->Theta();
2197                     dghtBarV0Theta = lDghtBarV0->Theta();                       
2198                     dghtMesV0Theta = lDghtMesV0->Theta();
2199                     bacPt          = lBach->Pt();
2200                     dghtBarV0Pt    = lDghtBarV0->Pt();
2201                     dghtMesV0Pt    = lDghtMesV0->Pt();
2202                         
2203                } else if ( fAnalysisType == "AOD") {
2204
2205                     AliAODMCParticle *xiMC = (AliAODMCParticle*) arrayMC->At(iCurrentLabelStack);
2206                     if (xiMC->GetNDaughters() != 2) continue;
2207                     if (xiMC->GetDaughter(0) < 0 )  continue;
2208                     if (xiMC->GetDaughter(1) < 0 )  continue;
2209
2210                     AliAODMCParticle* lDght0ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(0) );
2211                     AliAODMCParticle* lDght1ofXi = (AliAODMCParticle*) arrayMC->At(  xiMC->GetDaughter(1) );
2212
2213                     AliAODMCParticle* lLambda = 0;
2214                     AliAODMCParticle* lBach   = 0;
2215
2216                     // Xi - Case 1
2217                     if ( lDght0ofXi->PdgCode() == lPdgCodeLambda  &&  lDght1ofXi->PdgCode() == lPdgCodeBach ){ 
2218                             lLambda = lDght0ofXi;   // dghter0 = Lambda
2219                             lBach   = lDght1ofXi;   // dghter1 = Pi-
2220                     }
2221                     // Xi - Case 2
2222                     else if ( lDght0ofXi->PdgCode() == lPdgCodeBach && lDght1ofXi->PdgCode() == lPdgCodeLambda ){
2223                             lBach   = lDght0ofXi;   // dghter0 = Pi
2224                             lLambda = lDght1ofXi;   //dghter1 = Lambda
2225                     }
2226                     // Otherwise - Case 3
2227                     else continue;
2228
2229                     // - Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2230                     if (fApplyAccCut) {
2231                           if ( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2232                           if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2233                           if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2234                     }
2235
2236                     //-----------
2237                     // - V0 level 
2238                     AliAODMCParticle* lDghtBarV0 = 0;
2239                     AliAODMCParticle* lDghtMesV0 = 0;
2240
2241                     if( lLambda->GetNDaughters() != 2 )  continue;
2242                     if( lLambda->GetDaughter(0) < 0 )    continue;
2243                     if( lLambda->GetDaughter(1) < 0 )    continue;
2244
2245                     AliAODMCParticle* lDght0ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(0) );
2246                     AliAODMCParticle* lDght1ofLambda = (AliAODMCParticle*) arrayMC->At(  lLambda->GetDaughter(1) );
2247
2248                     // V0 - Case 1
2249                     if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtBarV0 && lDght1ofLambda->PdgCode() == lPdgCodeDghtMesV0 ) { 
2250                             lDghtBarV0 = lDght0ofLambda;   // dghter0 = Proton
2251                             lDghtMesV0 = lDght1ofLambda;   // dghter1 = Pi-
2252                     } 
2253                     // V0 - Case 2
2254                     else if ( lDght0ofLambda->PdgCode() == lPdgCodeDghtMesV0 && lDght1ofLambda->PdgCode() == lPdgCodeDghtBarV0 ) { 
2255                             lDghtMesV0 = lDght0ofLambda;   // dghter0 = Pi-
2256                             lDghtBarV0 = lDght1ofLambda;   // dghter1 = proton
2257                     } 
2258                     // V0 otherwise - Case 3
2259                     else continue;
2260
2261                     // - Check the emission of particle stays within the acceptance of the detector
2262                     if (fApplyAccCut) {
2263                          if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2264                          if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2265                          if( lDghtBarV0->Pt() < 0.250 ) continue;
2266                          if( lDghtMesV0->Pt() < 0.150 ) continue;
2267                     }
2268
2269                     lambdaTheta    = lLambda->Theta();
2270                     bacTheta       = lBach->Theta();
2271                     dghtBarV0Theta = lDghtBarV0->Theta();
2272                     dghtMesV0Theta = lDghtMesV0->Theta();
2273                     bacPt          = lBach->Pt();
2274                     dghtBarV0Pt    = lDghtBarV0->Pt();
2275                     dghtMesV0Pt    = lDghtMesV0->Pt();
2276                }
2277
2278                //---------------------------------------
2279                // - Filling histos for findable cascades
2280                // - Fill theta histos 
2281                lHistThetaLambda->Fill( lRadToDeg * lambdaTheta );
2282                lHistThetaBach->Fill( lRadToDeg * bacTheta );
2283                lHistThetaBarDghter->Fill( lRadToDeg * dghtBarV0Theta );
2284                lHistThetaMesDghter->Fill( lRadToDeg * dghtMesV0Theta );
2285                // - Fill pt histos
2286                lHistPtBach             ->Fill( bacPt );
2287                lHistPtBarDghter        ->Fill( dghtBarV0Pt );
2288                lHistPtMesDghter        ->Fill( dghtMesV0Pt );
2289                l2dHistGenPtVsGenYFdbl  ->Fill( partPt, lRapXiMC );
2290
2291                ncascperev++;                    
2292                              
2293          }// This is the end of the loop on primaries
2294   
2295          if (iCascType == 1) {
2296               fHistnXiMinusPerEv->Fill(ncascperev);
2297               fHistnXiMinusPerEvTot->Fill(ncascperevtot);
2298          }
2299          if (iCascType == 2) {
2300               fHistnXiPlusPerEv->Fill(ncascperev);
2301               fHistnXiPlusPerEvTot->Fill(ncascperevtot);
2302          }
2303          if (iCascType == 3) {
2304               fHistnOmegaMinusPerEv->Fill(ncascperev);
2305               fHistnOmegaMinusPerEvTot->Fill(ncascperevtot);
2306          }
2307          if (iCascType == 4) {
2308               fHistnOmegaPlusPerEv->Fill(ncascperev);
2309               fHistnOmegaPlusPerEvTot->Fill(ncascperevtot);
2310          }
2311
2312          // - Re-initialisation of the local THF pointers
2313          lHistEtaGenCasc         = 0x0;
2314          lHistThetaGenCasc       = 0x0;
2315          l2dHistGenPtVsGenYFdbl  = 0x0;
2316          lHistThetaLambda        = 0x0;
2317          lHistThetaBach          = 0x0;
2318          lHistThetaBarDghter     = 0x0;
2319          lHistThetaMesDghter     = 0x0;
2320          lHistPtBach             = 0x0;
2321          lHistPtBarDghter        = 0x0;
2322          lHistPtMesDghter        = 0x0; 
2323
2324    } // end of loop over the different types of cascades (Xi-+, Omega-+)
2325         
2326  
2327  
2328    //-----------------------------------------  
2329    // - Loop over the reconstructed candidates
2330    //-----------------------------------------
2331    Int_t nAssoXiMinus     = 0;
2332    Int_t nAssoXiPlus      = 0;
2333    Int_t nAssoOmegaMinus  = 0;
2334    Int_t nAssoOmegaPlus   = 0;
2335    Int_t lPosTPCClusters  = 0;
2336    Int_t lNegTPCClusters  = 0;
2337    Int_t lBachTPCClusters = 0;
2338    Double_t lDcaXiDaughters            = -1. ;
2339    Double_t lDcaBachToPrimVertexXi     = -1. ;
2340    Double_t lXiCosineOfPointingAngle   = -1. ;
2341    Double_t lPosXi[3]                  = { -1000.0, -1000.0, -1000.0 };
2342    Double_t lXiRadius                  = -1000. ;
2343    Double_t lInvMassLambdaAsCascDghter = 0.;
2344    Double_t lDcaV0DaughtersXi          = -1.;
2345    Double_t lV0CosineOfPointingAngleXi = -1.;
2346    Double_t lV0CosineOfPointingAngle   = -1.;
2347    Double_t lPosV0Xi[3]                = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
2348    Double_t lV0RadiusXi                = -1000.;
2349    Double_t lDcaV0ToPrimVertexXi       = -1.;
2350    Double_t lDcaPosToPrimVertexXi      = -1.;
2351    Double_t lDcaNegToPrimVertexXi      = -1.;
2352    Double_t lChargeXi                  = -1.;
2353    Double_t lV0mom                     = -1000.;
2354    Double_t lmcPt                      = -1.;         
2355    Double_t lmcRapCasc                 = -1.; 
2356    Double_t lmcEta                     = -1000.;     
2357    Double_t lmcTransvRadius            = -1000.; 
2358    Double_t lrecoPt                    = -100.;     
2359    Double_t lrecoTransvRadius          = -1000.; 
2360    Double_t lDeltaPhiMcReco            = -1.;
2361    Double_t lBachTransvMom             = 0.;
2362    Double_t lpTrackTransvMom           = 0.;
2363    Double_t lnTrackTransvMom           = 0.;
2364    Double_t lmcPtPosV0Dghter           = -100.;
2365    Double_t lmcPtNegV0Dghter           = -100.;
2366    Double_t lrecoP                     = -100.;
2367    Double_t lmcPtBach                  = -100.;
2368    Double_t cascadeMass                = 0.;
2369
2370    // - Get the number of cascades
2371    Int_t ncascades = 0;
2372    if      ( fAnalysisType == "ESD" ) { ncascades = lESDevent->GetNumberOfCascades(); }
2373    else if ( fAnalysisType == "AOD" ) { ncascades = lAODevent->GetNumberOfCascades(); }
2374
2375    //-------------------------------
2376    // - Begining of the Cascade Loop
2377    for (Int_t iXi = 0; iXi < ncascades; iXi++) {
2378
2379         Bool_t   lIsPosInXiProton      = kFALSE;
2380         Bool_t   lIsPosInXiPion        = kFALSE;
2381         Bool_t   lIsPosInOmegaProton   = kFALSE;
2382         Bool_t   lIsPosInOmegaPion     = kFALSE;
2383         Bool_t   lIsNegInXiProton      = kFALSE;
2384         Bool_t   lIsNegInXiPion        = kFALSE;
2385         Bool_t   lIsNegInOmegaProton   = kFALSE;
2386         Bool_t   lIsNegInOmegaPion     = kFALSE;
2387         Bool_t   lIsBachelorKaon       = kFALSE;
2388         Bool_t   lIsBachelorPion       = kFALSE;
2389         Bool_t   lIsBachelorKaonForTPC = kFALSE;
2390         Bool_t   lIsBachelorPionForTPC = kFALSE;
2391         Bool_t   lIsNegPionForTPC      = kFALSE;
2392         Bool_t   lIsPosPionForTPC      = kFALSE;
2393         Bool_t   lIsNegProtonForTPC    = kFALSE;
2394         Bool_t   lIsPosProtonForTPC    = kFALSE;
2395
2396         // - Combined PID
2397         // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
2398         Double_t lPriorsGuessXi[14]    = {0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
2399         Double_t lPriorsGuessOmega[14] = {0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0};
2400         Double_t ppionBach = 0.0, pkaonBach = 0.0;
2401         Bool_t   lIsBachelorMCPiMinus  = kFALSE;
2402         Bool_t   lIsBachelorMCPiPlus   = kFALSE;
2403         Bool_t   lIsBachelorMCKMinus   = kFALSE;
2404         Bool_t   lIsBachelorMCKPlus    = kFALSE;
2405         Double_t lInvMassXiMinus    = 0.;
2406         Double_t lInvMassXiPlus     = 0.;
2407         Double_t lInvMassOmegaMinus = 0.;
2408         Double_t lInvMassOmegaPlus  = 0.;
2409         Bool_t lAssoXiMinus    = kFALSE;
2410         Bool_t lAssoXiPlus     = kFALSE;
2411         Bool_t lAssoOmegaMinus = kFALSE;
2412         Bool_t lAssoOmegaPlus  = kFALSE;
2413        
2414         Float_t  etaBach = 0.;
2415         Float_t  etaPos  = 0.;
2416         Float_t  etaNeg  = 0.;
2417
2418         if ( fAnalysisType == "ESD" ) {         
2419
2420              // - Load the cascade
2421              AliESDcascade *xiESD = lESDevent->GetCascade(iXi);
2422              if (!xiESD) continue;
2423         
2424              // - Connection to daughter tracks of the current cascade          
2425              UInt_t lIdxPosXi = (UInt_t) TMath::Abs( xiESD->GetPindex() );
2426              UInt_t lIdxNegXi = (UInt_t) TMath::Abs( xiESD->GetNindex() );
2427              UInt_t lBachIdx  = (UInt_t) TMath::Abs( xiESD->GetBindex() );
2428                    
2429              // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2430              if(lBachIdx == lIdxNegXi) {
2431                  AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2432              }
2433              if(lBachIdx == lIdxPosXi) {
2434                  AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2435              }
2436       
2437              // - Get the daughter tracks
2438              AliESDtrack *pTrackXi = lESDevent->GetTrack( lIdxPosXi );
2439              AliESDtrack *nTrackXi = lESDevent->GetTrack( lIdxNegXi );
2440              AliESDtrack *bachTrackXi = lESDevent->GetTrack( lBachIdx  );
2441              if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2442                  Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
2443                  continue;
2444              }
2445         
2446              // Get the number of TPC clusters
2447              lPosTPCClusters   = pTrackXi->GetTPCNcls();
2448              lNegTPCClusters   = nTrackXi->GetTPCNcls();
2449              lBachTPCClusters  = bachTrackXi->GetTPCNcls(); 
2450              // - Rejection of a poor quality tracks
2451              if(fkQualityCutTPCrefit){
2452                  // - Poor quality related to TPCrefit
2453                  ULong_t pStatus    = pTrackXi->GetStatus();
2454                  ULong_t nStatus    = nTrackXi->GetStatus();
2455                  ULong_t bachStatus = bachTrackXi->GetStatus();
2456                  if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2457                  if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2458                  if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2459              } 
2460              if(fkQualityCutnTPCcls){
2461                 // - Poor quality related to TPC clusters
2462                 if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than 80 TPC clusters ... continue!"); continue; }
2463                 if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than 80 TPC clusters ... continue!"); continue; }
2464                 if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than 80 TPC clusters ... continue!"); continue; }
2465              }
2466
2467              etaPos  = pTrackXi->Eta();
2468              etaNeg  = nTrackXi->Eta();
2469              etaBach = bachTrackXi->Eta();
2470         
2471              // - Info over reconstructed cascades
2472              Double_t lV0quality = 0.;
2473              if( bachTrackXi->Charge() < 0 ) {
2474                   //Calculate the effective mass of the Xi- candidate: Xi- hyp. (pdg code 3312
2475                   lV0quality = 0.;
2476                   xiESD->ChangeMassHypothesis(lV0quality , 3312);       
2477                   lInvMassXiMinus = xiESD->GetEffMassXi();
2478                   //Calculate the effective mass of the Xi- candidate: Omega- hyp. (pdg code 3334)
2479                   lV0quality = 0.;
2480                   xiESD->ChangeMassHypothesis(lV0quality , 3334);       
2481                   lInvMassOmegaMinus = xiESD->GetEffMassXi();
2482                   //Back to "default" hyp. (Xi-)        
2483                   lV0quality = 0.;
2484                   xiESD->ChangeMassHypothesis(lV0quality , 3312);
2485              }
2486              if( bachTrackXi->Charge() >  0 ){
2487                   //Calculate the effective mass of the Xi- candidate: Xi+ hyp. (pdg code -3312)
2488                   lV0quality = 0.;
2489                   xiESD->ChangeMassHypothesis(lV0quality , -3312);      
2490                   lInvMassXiPlus = xiESD->GetEffMassXi();
2491                   //Calculate the effective mass of the Xi- candidate: Omega+ hyp. (pdg code -3334)
2492                   lV0quality = 0.;
2493                   xiESD->ChangeMassHypothesis(lV0quality , -3334);      
2494                   lInvMassOmegaPlus = xiESD->GetEffMassXi();
2495                   //Back to "default" hyp. (Xi-)
2496                   lV0quality = 0.;
2497                   xiESD->ChangeMassHypothesis(lV0quality , -3312);
2498              }
2499              lDcaXiDaughters          = xiESD->GetDcaXiDaughters();
2500              lDcaBachToPrimVertexXi   = TMath::Abs( bachTrackXi->GetD(lBestPrimaryVtxPos[0],lBestPrimaryVtxPos[1],lMagneticField) );
2501              lXiCosineOfPointingAngle = xiESD->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
2502              xiESD->GetXYZcascade( lPosXi[0], lPosXi[1], lPosXi[2] );  
2503              lInvMassLambdaAsCascDghter = xiESD->GetEffMass();
2504              lDcaV0DaughtersXi          = xiESD->GetDcaV0Daughters();
2505              lV0CosineOfPointingAngleXi = xiESD->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
2506              lV0CosineOfPointingAngle   = xiESD->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2]); 
2507              xiESD->GetXYZ( lPosV0Xi[0], lPosV0Xi[1], lPosV0Xi[2] );
2508              lDcaV0ToPrimVertexXi       = xiESD->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
2509              lDcaPosToPrimVertexXi      = TMath::Abs( pTrackXi->GetD(lBestPrimaryVtxPos[0],lBestPrimaryVtxPos[1],lMagneticField) );
2510              lDcaNegToPrimVertexXi      = TMath::Abs( nTrackXi->GetD(lBestPrimaryVtxPos[0],lBestPrimaryVtxPos[1],lMagneticField) );
2511              lChargeXi = xiESD->Charge();
2512         
2513              //------------------
2514              // - PID Information
2515
2516              // - Combined VO-positive-daughter PID
2517              AliPID pPidXi;    pPidXi.SetPriors( lPriorsGuessXi );
2518              AliPID pPidOmega; pPidOmega.SetPriors( lPriorsGuessOmega );                
2519              if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  
2520                   Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
2521                   pPidXi.SetProbabilities(r);
2522                   pPidOmega.SetProbabilities(r);                
2523                   // Check if the V0 positive track is a proton (case for Xi-)
2524                   Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
2525                   if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
2526                       pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
2527                       pproton > pPidXi.GetProbability(AliPID::kPion)     &&
2528                       pproton > pPidXi.GetProbability(AliPID::kKaon)        ) lIsPosInXiProton = kTRUE;
2529                   // Check if the V0 positive track is a pi+ (case for Xi+)
2530                   Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
2531                   if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
2532                       ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
2533                       ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
2534                       ppion > pPidXi.GetProbability(AliPID::kProton)      ) lIsPosInXiPion = kTRUE;
2535                   // Check if the V0 positive track is a proton (case for Omega-)
2536                   pproton = 0.;
2537                   pproton = pPidOmega.GetProbability(AliPID::kProton);
2538                   if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
2539                       pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
2540                       pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
2541                       pproton > pPidOmega.GetProbability(AliPID::kKaon)       ) lIsPosInOmegaProton = kTRUE;
2542                   // Check if the V0 positive track is a pi+ (case for Omega+)
2543                   ppion = 0.;
2544                   ppion = pPidOmega.GetProbability(AliPID::kPion);
2545                   if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
2546                       ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
2547                       ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
2548                       ppion > pPidOmega.GetProbability(AliPID::kProton)      ) lIsPosInOmegaPion = kTRUE;
2549              }          
2550              // - Combined VO-negative-daughter PID
2551              AliPID nPidXi;    nPidXi.SetPriors( lPriorsGuessXi );
2552              AliPID nPidOmega; nPidOmega.SetPriors( lPriorsGuessOmega );                
2553              if( nTrackXi->IsOn(AliESDtrack::kESDpid) ) {  
2554                   Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
2555                   nPidXi.SetProbabilities(r);
2556                   nPidOmega.SetProbabilities(r);
2557                   // Check if the V0 negative track is a pi- (case for Xi-)
2558                   Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
2559                   if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
2560                       ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
2561                       ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
2562                       ppion > nPidXi.GetProbability(AliPID::kProton)      ) lIsNegInXiPion = kTRUE;
2563                   // Check if the V0 negative track is an anti-proton (case for Xi+)
2564                   Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
2565                   if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
2566                       pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
2567                       pproton > nPidXi.GetProbability(AliPID::kPion)     &&
2568                       pproton > nPidXi.GetProbability(AliPID::kKaon)       ) lIsNegInXiProton = kTRUE;
2569                   // Check if the V0 negative track is a pi- (case for Omega-)
2570                   ppion = 0.;
2571                   ppion = nPidOmega.GetProbability(AliPID::kPion);
2572                   if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
2573                       ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
2574                       ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
2575                       ppion > nPidOmega.GetProbability(AliPID::kProton)     ) lIsNegInOmegaPion = kTRUE;
2576                   // Check if the V0 negative track is an anti-proton (case for Omega+)
2577                   pproton = 0.;
2578                   pproton = nPidOmega.GetProbability(AliPID::kProton);
2579                   if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
2580                       pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
2581                       pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
2582                       pproton > nPidOmega.GetProbability(AliPID::kKaon)       ) lIsNegInOmegaProton = kTRUE;
2583              }
2584              // - Combined bachelor PID
2585              AliPID bachPidXi;    bachPidXi.SetPriors( lPriorsGuessXi );
2586              AliPID bachPidOmega; bachPidOmega.SetPriors( lPriorsGuessOmega );
2587              if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  
2588                    Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
2589                    bachPidXi.SetProbabilities(r);
2590                    bachPidOmega.SetProbabilities(r);
2591                    // Check if the bachelor track is a pion
2592                    ppionBach = bachPidXi.GetProbability(AliPID::kPion);
2593                    if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
2594                        ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
2595                        ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
2596                        ppionBach > bachPidXi.GetProbability(AliPID::kProton)     ) lIsBachelorPion = kTRUE;
2597                    // Check if the bachelor track is a kaon
2598                    pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
2599                    if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
2600                        pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
2601                        pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
2602                        pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)     ) lIsBachelorKaon = kTRUE;  
2603              }
2604              // - 4-sigma bands on Bethe-Bloch curve
2605              // Bachelor
2606              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
2607              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
2608              // Negative V0 daughter
2609              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
2610              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
2611              // Positive V0 daughter
2612              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
2613              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
2614              /*        
2615              const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
2616              const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
2617              const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
2618              if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){   
2619                  Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
2620                  Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
2621                  Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
2622                  // Bachelor
2623                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
2624                  if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
2625                  if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
2626                  // Negative V0 daughter
2627                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
2628                  if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
2629                  if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
2630                  // Positive V0 daughter
2631                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
2632                  if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
2633                  if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
2634              }*/
2635              // - PID proba Vs Pt(Bach)
2636              Int_t      lblBachForPID = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2637              TParticle* mcBachForPID  = lMCstack->Particle( lblBachForPID );
2638              lmcPtBach = mcBachForPID->Pt();
2639              // - MC perfect PID
2640              if( mcBachForPID->GetPdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
2641              if( mcBachForPID->GetPdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
2642              if( mcBachForPID->GetPdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
2643              if( mcBachForPID->GetPdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
2644         
2645         
2646              //---------------------------------------------------------------
2647              // -  MC association (care : lots of "continue;" below this line)
2648              if(fDebug > 5) cout<< "MC EventNumber: "<<lMCevent->Header()->GetEvent()<<" / MC event Number in Run : "<<lMCevent->Header()->GetEventNrInRun()<<endl;
2649              // - Level of the V0 daughters
2650              Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );  
2651              Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );                 
2652              TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
2653              TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
2654              // - Level of the Xi daughters     
2655              Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
2656              Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
2657              if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
2658              if( lblMotherPosV0Dghter < 0 ) continue;                    // this particle is primary, no mother   
2659              if( lblMotherNegV0Dghter < 0 ) continue;                    // this particle is primary, no mother
2660                                                                                  // mothers = Lambda candidate ... a priori
2661              TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
2662              TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );  // MN: redundant?? already checked that labels are the same...-->same part from stack
2663              Int_t lblBach = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2664              TParticle* mcBach = lMCstack->Particle( lblBach ); 
2665              // - Level of Xi candidate
2666              Int_t lblGdMotherPosV0Dghter = mcMotherPosV0Dghter->GetFirstMother() ;
2667              Int_t lblGdMotherNegV0Dghter = mcMotherNegV0Dghter->GetFirstMother() ;
2668              if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
2669              if( lblGdMotherPosV0Dghter < 0 ) continue;                  // primary lambda ...   
2670              if( lblGdMotherNegV0Dghter < 0 ) continue;                  // primary lambda ...                                  
2671                                                                               // Gd mothers = Xi candidate ... a priori
2672              TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
2673              TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );                           
2674              Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );  
2675              if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
2676              TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
2677             
2678              // - Check if cascade is primary
2679              if (!(lMCstack->IsPhysicalPrimary(lblMotherBach))) continue;  
2680
2681              // - Manage boolean for association
2682              if     ( mcMotherBach          ->GetPdgCode() == 3312 &&
2683                       mcGdMotherPosV0Dghter ->GetPdgCode() == 3312 &&
2684                       mcGdMotherNegV0Dghter ->GetPdgCode() == 3312   ) {lAssoXiMinus = kTRUE;
2685                                                                         cascadeMass = 1.321;
2686                                                                          nAssoXiMinus++; }
2687              else if( mcMotherBach           ->GetPdgCode() == -3312 &&
2688                       mcGdMotherPosV0Dghter  ->GetPdgCode() == -3312 &&
2689                       mcGdMotherNegV0Dghter  ->GetPdgCode() == -3312   ) {lAssoXiPlus = kTRUE;
2690                                                                           cascadeMass = 1.321;
2691                                                                           nAssoXiPlus++; }
2692              else if( mcMotherBach           ->GetPdgCode() == 3334 &&
2693                       mcGdMotherPosV0Dghter  ->GetPdgCode() == 3334 &&
2694                       mcGdMotherNegV0Dghter  ->GetPdgCode() == 3334    ) {lAssoOmegaMinus = kTRUE;
2695                                                                           cascadeMass = 1.672;
2696                                                                           nAssoOmegaMinus++; }
2697              else if( mcMotherBach           ->GetPdgCode() == -3334 &&
2698                       mcGdMotherPosV0Dghter  ->GetPdgCode() == -3334 &&
2699                       mcGdMotherNegV0Dghter  ->GetPdgCode() == -3334   ) {lAssoOmegaPlus = kTRUE;
2700                                                                              cascadeMass = 1.672;
2701                                                                              nAssoOmegaPlus++; }
2702              // If a proper association  exists ...     
2703              if(fDebug > 4){
2704                       cout<<"XiMinus    = "<<lAssoXiMinus   <<endl;
2705                       cout<<"XiPlus     = "<<lAssoXiPlus    <<endl;
2706                       cout<<"OmegaMinus = "<<lAssoOmegaMinus<<endl;
2707                       cout<<"OmegaPlus  = "<<lAssoOmegaPlus <<endl 
2708                           <<"----"            <<endl;   
2709              }
2710              if(fDebug > 5){
2711                          cout<<endl;
2712                       cout<<"- V0 daughters - "<<endl;
2713                       cout<<"     + V0 Pos. / Label : "<<lblPosV0Dghter<<" - Pdg Code : "<<mcPosV0Dghter->GetTitle()<<endl;
2714                       cout<<"     - V0 Neg. / Label : "<<lblNegV0Dghter<<" - Pdg Code : "<<mcNegV0Dghter->GetTitle()<<endl;
2715
2716                       cout<<"- Xi daughters - "<<endl;
2717                       cout<<"     + V0 Pos. mother / Label : "<<lblMotherPosV0Dghter<<" - Pdg Code : "<<mcMotherPosV0Dghter->GetTitle()<<endl;
2718                       cout<<"     - V0 Neg. mother / Label : "<<lblMotherNegV0Dghter<<" - Pdg Code : "<<mcMotherNegV0Dghter->GetTitle()<<endl;
2719                 
2720                       cout<<"     --  Bach. / Label :"<<lblBach<<" -  Pdg Code : "<<mcBach->GetTitle()<<endl;
2721                 
2722                       cout<<"- Xi candidate -"<<endl;
2723                       cout<<"    +  V0 Pos. Gd Mother / Label : "<<lblGdMotherPosV0Dghter<<" - Pdg Code : "<< mcGdMotherPosV0Dghter->GetTitle()<<endl;
2724                       cout<<"    -  V0 Neg. Gd Mother / Label : "<<lblGdMotherNegV0Dghter<<" - Pdg Code : "<< mcGdMotherNegV0Dghter->GetTitle()<<endl;
2725                 
2726                       cout<<"    --  Mother Bach. / Label : "<<lblMotherBach<<" - Pdg Code    : "<<mcMotherBach->GetTitle()<<endl;
2727                       cout<<endl;
2728              }
2729         
2730              lmcPt             = mcMotherBach->Pt();
2731              lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) / (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
2732              lmcEta            = mcMotherBach->Eta();
2733              lmcTransvRadius   = mcBach->R(); // to get the decay point of Xi, = the production vertex of Bachelor ...
2734              TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
2735              lrecoPt           = xiESD->Pt();
2736              lrecoTransvRadius = TMath::Sqrt( xiESD->Xv() * xiESD->Xv() + xiESD->Yv() * xiESD->Yv() );
2737              TVector3 lrecoTVect3Mom( xiESD->Px(), xiESD->Py(), xiESD->Pz() );
2738              lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
2739              lmcPtPosV0Dghter  = mcPosV0Dghter->Pt() ;
2740              lmcPtNegV0Dghter  = mcNegV0Dghter->Pt();
2741              lrecoP            = xiESD->P();
2742              Double_t nV0mom[3] = {0. ,0. ,0. };
2743              Double_t pV0mom[3] = {0. ,0. ,0. };
2744              xiESD->GetNPxPyPz(nV0mom[0],nV0mom[1],nV0mom[2]);    
2745              xiESD->GetPPxPyPz(pV0mom[0],pV0mom[1],pV0mom[2]);
2746              lV0mom = TMath::Sqrt(TMath::Power(nV0mom[0]+pV0mom[0],2)+TMath::Power(nV0mom[1]+pV0mom[1],2)+TMath::Power(nV0mom[2]+pV0mom[2],2));
2747              Double_t lBachMomX = 0.; Double_t lBachMomY = 0.; Double_t lBachMomZ = 0.;
2748              xiESD->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
2749              lBachTransvMom   = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
2750              lnTrackTransvMom = TMath::Sqrt( nV0mom[0]*nV0mom[0] + nV0mom[1]*nV0mom[1] );
2751              lpTrackTransvMom = TMath::Sqrt( pV0mom[0]*pV0mom[0] + pV0mom[1]*pV0mom[1] );
2752             
2753         } else if ( fAnalysisType == "AOD" ) {
2754
2755              // - Load the cascade
2756              const AliAODcascade *xiAOD = lAODevent->GetCascade(iXi);
2757              if (!xiAOD) continue;
2758
2759              // - Connection to daughter tracks of the current cascade
2760              AliAODTrack *pTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(0) );
2761              AliAODTrack *nTrackXi    = dynamic_cast<AliAODTrack*>( xiAOD->GetDaughter(1) );
2762              AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xiAOD->GetDecayVertexXi()->GetDaughter(0) );
2763              if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2764                   AliWarning("ERROR: Could not retrieve one of the 3 AOD daughter tracks of the cascade ...");
2765                   continue;
2766              }
2767              UInt_t lIdxPosXi  = (UInt_t) TMath::Abs( pTrackXi->GetID() );
2768              UInt_t lIdxNegXi  = (UInt_t) TMath::Abs( nTrackXi->GetID() );
2769              UInt_t lBachIdx   = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
2770
2771              // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2772              if(lBachIdx == lIdxNegXi) {
2773                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2774              }
2775              if(lBachIdx == lIdxPosXi) {
2776                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2777              }
2778              lPosTPCClusters   = pTrackXi->GetTPCNcls();
2779              lNegTPCClusters   = nTrackXi->GetTPCNcls();
2780              lBachTPCClusters  = bachTrackXi->GetTPCNcls();
2781
2782              // - Rejection of a poor quality tracks
2783              if (fkQualityCutTPCrefit) {
2784                  // - Poor quality related to TPCrefit
2785                  if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2786                  if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2787                  if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2788              }
2789              if (fkQualityCutnTPCcls) {
2790                  // - Poor quality related to TPC clusters
2791                  if(lPosTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Pos. track has less than 80 TPC clusters ... continue!"); continue; }
2792                  if(lNegTPCClusters  < fMinnTPCcls) { AliWarning("Pb / V0 Neg. track has less than 80 TPC clusters ... continue!"); continue; }
2793                  if(lBachTPCClusters < fMinnTPCcls) { AliWarning("Pb / Bach.   track has less than 80 TPC clusters ... continue!"); continue; }
2794              }
2795
2796              etaPos  = pTrackXi->Eta();
2797              etaNeg  = nTrackXi->Eta();
2798              etaBach = bachTrackXi->Eta();
2799
2800              // - Info over reconstructed cascades
2801              if( bachTrackXi->Charge() < 0 ) {
2802                   lInvMassXiMinus = xiAOD->MassXi();
2803                   lInvMassOmegaMinus = xiAOD->MassOmega();
2804              }
2805              if( bachTrackXi->Charge() >  0 ){
2806                   lInvMassXiPlus = xiAOD->MassXi();
2807                   lInvMassOmegaPlus = xiAOD->MassOmega();
2808              }
2809              lDcaXiDaughters            = xiAOD->DcaXiDaughters();
2810              lDcaBachToPrimVertexXi     = xiAOD->DcaBachToPrimVertex();
2811              lXiCosineOfPointingAngle   = xiAOD->CosPointingAngleXi( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
2812              lPosXi[0]                  = xiAOD->DecayVertexXiX();
2813              lPosXi[1]                  = xiAOD->DecayVertexXiY();
2814              lPosXi[2]                  = xiAOD->DecayVertexXiZ();
2815              lInvMassLambdaAsCascDghter = xiAOD->MassLambda();
2816              lDcaV0DaughtersXi          = xiAOD->DcaV0Daughters();
2817              lV0CosineOfPointingAngleXi = xiAOD->CosPointingAngle( lPosXi );
2818              lV0CosineOfPointingAngle   = xiAOD->CosPointingAngle( lBestPrimaryVtxPos );
2819              lPosV0Xi[0]                = xiAOD->DecayVertexV0X();
2820              lPosV0Xi[1]                = xiAOD->DecayVertexV0Y();
2821              lPosV0Xi[2]                = xiAOD->DecayVertexV0Z();
2822              lDcaV0ToPrimVertexXi       = xiAOD->DcaV0ToPrimVertex();
2823              lDcaPosToPrimVertexXi      = xiAOD->DcaPosToPrimVertex();
2824              lDcaNegToPrimVertexXi      = xiAOD->DcaNegToPrimVertex();
2825              lChargeXi                  = xiAOD->ChargeXi();
2826
2827              //------------------
2828              // - PID Information
2829              // Combined VO-positive-daughter PID
2830              // Combined bachelor PID
2831              /* 
2832              AliPID bachPidXi;       bachPidXi.SetPriors(    lPriorsGuessXi    );
2833              AliPID bachPidOmega;    bachPidOmega.SetPriors( lPriorsGuessOmega );
2834
2835              if ( bachTrackXi->IsOn(AliESDtrack::kESDpid) ) {  // Combined PID exists
2836                   Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
2837                   bachPidXi.SetProbabilities(r);
2838                   bachPidOmega.SetProbabilities(r);
2839                   // Check if the bachelor track is a pion
2840                   ppionBach = bachPidXi.GetProbability(AliPID::kPion);
2841                   if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
2842                       ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
2843                       ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
2844                       ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
2845                   // Check if the bachelor track is a kaon
2846                   pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
2847                   if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
2848                       pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
2849                       pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
2850                       pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;
2851              }// end if bachelor track with existing combined PID
2852              */
2853
2854              // - TPC PID: 4-sigma bands on Bethe-Bloch curve
2855              // Bachelor
2856              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
2857              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
2858              // Negative V0 daughter
2859              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
2860              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
2861              // Positive V0 daughter
2862              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
2863              if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
2864              /*
2865              const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
2866              const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
2867              const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
2868              if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
2869                  Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
2870                  Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
2871                  Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
2872                  // Bachelor
2873                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
2874                  if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
2875                  if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
2876                  // Negative V0 daughter
2877                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
2878                  if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
2879                  if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
2880                  // Positive V0 daughter
2881                  if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
2882                  if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
2883                  if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
2884              }*/
2885
2886              // - PID proba Vs Pt(Bach)
2887              Int_t lblBachForPID = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2888              AliAODMCParticle* mcBachForPID   = (AliAODMCParticle*) arrayMC->At( lblBachForPID );
2889              lmcPtBach = mcBachForPID->Pt();
2890
2891              // - MC perfect PID
2892              if( mcBachForPID->PdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
2893              if( mcBachForPID->PdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
2894              if( mcBachForPID->PdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
2895              if( mcBachForPID->PdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
2896
2897              //--------------------------------------------------------------
2898              // - MC association (care : lots of "continue;" below this line)
2899              if(fDebug > 5) cout<<"MC EventNumber : "<<lMCevent->Header()->GetEvent()<<" / MC event Number in Run : "<<lMCevent->Header()->GetEventNrInRun()<<endl;
2900              // - Level of the V0 daughters
2901              Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );
2902              Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
2903              AliAODMCParticle* mcPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblPosV0Dghter );
2904              AliAODMCParticle* mcNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblNegV0Dghter );
2905              // - Level of the Xi daughters
2906              Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetMother();   
2907              Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetMother();
2908              if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
2909              if( lblMotherPosV0Dghter < 0 ) continue;                    // this particle is primary, no mother
2910              if( lblMotherNegV0Dghter < 0 ) continue;                    // this particle is primary, no mother
2911                                                                          // mothers = Lambda candidate ... a priori
2912              AliAODMCParticle* mcMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherPosV0Dghter );
2913              AliAODMCParticle* mcMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblMotherNegV0Dghter );  
2914              Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2915              AliAODMCParticle* mcBach   = (AliAODMCParticle*) arrayMC->At( lblBach );
2916              // - Level of Xi candidate
2917              Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetMother() ;
2918              Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetMother() ;
2919              if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
2920              if( lblGdMotherPosV0Dghter < 0 ) continue;                    // primary lambda ...
2921              if( lblGdMotherNegV0Dghter < 0 ) continue;                    // primary lambda ...
2922                                                                            // Gd mothers = Xi candidate ... a priori
2923              AliAODMCParticle* mcGdMotherPosV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherPosV0Dghter );
2924              AliAODMCParticle* mcGdMotherNegV0Dghter = (AliAODMCParticle*) arrayMC->At( lblGdMotherNegV0Dghter );
2925              Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetMother() );
2926              if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
2927              AliAODMCParticle* mcMotherBach = (AliAODMCParticle*) arrayMC->At( lblMotherBach );
2928
2929              // - Check if cascade is primary
2930              if (!(mcMotherBach->IsPhysicalPrimary())) continue;
2931
2932              // - Manage boolean for association
2933              if     ( mcMotherBach          ->GetPdgCode() == 3312 &&
2934                       mcGdMotherPosV0Dghter ->GetPdgCode() == 3312 &&
2935                       mcGdMotherNegV0Dghter ->GetPdgCode() == 3312   ) {lAssoXiMinus = kTRUE;
2936                                                                         cascadeMass = 1.321;
2937                                                                         nAssoXiMinus++; }
2938              else if( mcMotherBach           ->GetPdgCode() == -3312 &&
2939                       mcGdMotherPosV0Dghter  ->GetPdgCode() == -3312 &&
2940                       mcGdMotherNegV0Dghter  ->GetPdgCode() == -3312   ) {lAssoXiPlus = kTRUE;
2941                                                                           cascadeMass = 1.321;
2942                                                                           nAssoXiPlus++; }
2943              else if( mcMotherBach           ->GetPdgCode() == 3334 &&
2944                       mcGdMotherPosV0Dghter  ->GetPdgCode() == 3334 &&
2945                       mcGdMotherNegV0Dghter  ->GetPdgCode() == 3334   ) {lAssoOmegaMinus = kTRUE;
2946                                                                          cascadeMass = 1.672;
2947                                                                          nAssoOmegaMinus++; }
2948              else if( mcMotherBach           ->GetPdgCode() == -3334 &&
2949                       mcGdMotherPosV0Dghter  ->GetPdgCode() == -3334 &&
2950                       mcGdMotherNegV0Dghter  ->GetPdgCode() == -3334   ) {lAssoOmegaPlus = kTRUE;
2951                                                                           cascadeMass = 1.672;
2952                                                                           nAssoOmegaPlus++; }
2953
2954              lmcPt              = mcMotherBach->Pt();
2955              lmcRapCasc         = 0.5*TMath::Log( (mcMotherBach->E() + mcMotherBach->Pz()) / (mcMotherBach->E() - mcMotherBach->Pz() +1.e-13) );
2956              lmcEta             = mcMotherBach->Eta();
2957              Float_t decayCascX = mcBach->Xv();
2958              Float_t decayCascY = mcBach->Yv();
2959              lmcTransvRadius    = TMath::Sqrt(decayCascX*decayCascX+decayCascY*decayCascY); // decay point of Xi, = the production vertex of Bachelor ...
2960              TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
2961              Double_t xiMomX    = xiAOD->MomXiX();
2962              Double_t xiMomY    = xiAOD->MomXiY();
2963              Double_t xiMomZ    = xiAOD->MomXiZ();
2964              lrecoPt            = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY ); 
2965              lrecoTransvRadius  = TMath::Sqrt( xiAOD->DecayVertexXiX() * xiAOD->DecayVertexXiX() + xiAOD->DecayVertexXiY() * xiAOD->DecayVertexXiY() );
2966              TVector3 lrecoTVect3Mom( xiMomX, xiMomY, xiMomZ );
2967              lDeltaPhiMcReco    = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
2968              lmcPtPosV0Dghter   = mcPosV0Dghter->Pt() ;
2969              lmcPtNegV0Dghter   = mcNegV0Dghter->Pt();
2970              lrecoP             = TMath::Sqrt( xiMomX*xiMomX   + xiMomY*xiMomY   + xiMomZ*xiMomZ );;
2971              Double_t lV0momX   = xiAOD->MomV0X();
2972              Double_t lV0momY   = xiAOD->MomV0Y();
2973              Double_t lV0momZ   = xiAOD->MomV0Z();
2974              lV0mom             = TMath::Sqrt(TMath::Power(lV0momX,2)+TMath::Power(lV0momY,2)+TMath::Power(lV0momZ,2));
2975              Double_t lBachMomX = xiAOD->MomBachX();
2976              Double_t lBachMomY = xiAOD->MomBachY();
2977              lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
2978              Double_t lV0NMomX = xiAOD->MomNegX();
2979              Double_t lV0NMomY = xiAOD->MomNegY();
2980              Double_t lV0PMomX = xiAOD->MomPosX();
2981              Double_t lV0PMomY = xiAOD->MomPosY();
2982              lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX   + lV0NMomY*lV0NMomY );
2983              lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX   + lV0PMomY*lV0PMomY );
2984             
2985         }
2986
2987         lXiRadius   = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
2988         lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] ); 
2989
2990         // - Cut on pt of the three daughter tracks
2991         if (lBachTransvMom<fMinPtCutOnDaughterTracks) continue;
2992         if (lpTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
2993         if (lnTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
2994        
2995         // - Cut on pseudorapidity of the three daughter tracks
2996         if (TMath::Abs(etaBach)>fEtaCutOnDaughterTracks) continue;
2997         if (TMath::Abs(etaPos)>fEtaCutOnDaughterTracks) continue;
2998         if (TMath::Abs(etaNeg)>fEtaCutOnDaughterTracks) continue;
2999        
3000         // - Extra-selection for cascade candidates
3001         if (fkExtraSelections) {
3002                 if (lDcaXiDaughters > 0.3) continue;              // in AliCascadeVertexer
3003                 if (lXiCosineOfPointingAngle < 0.999 ) continue;  // in AliCascadeVertexer
3004                 if (lDcaV0ToPrimVertexXi < 0.05) continue;        // in AliCascadeVertexer
3005                 if (lDcaBachToPrimVertexXi < 0.03) continue;      // in AliCascadeVertexer
3006                 if (lDcaV0DaughtersXi > 1.) continue;             // in AliV0vertexer
3007                 if (lV0CosineOfPointingAngleXi < 0.998) continue; // in AliV0vertexer
3008                 if (lDcaPosToPrimVertexXi < 0.1) continue;        // in AliV0vertexer
3009                 if (lDcaNegToPrimVertexXi < 0.1) continue;        // in AliV0vertexer
3010                 if(lXiRadius < .9) continue;                      // in AliCascadeVertexer
3011                 if(lV0RadiusXi < 0.9) continue;                   // in AliV0vertexer
3012         }
3013
3014         //-------------------------
3015         // - Fill combined PID TH1s
3016         if( lChargeXi < 0 && lIsBachelorPion )    fHistMassWithCombPIDXiMinus    ->Fill( lInvMassXiMinus    );
3017         if( lChargeXi > 0 && lIsBachelorPion )    fHistMassWithCombPIDXiPlus     ->Fill( lInvMassXiPlus     );
3018         if( lChargeXi < 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaMinus ->Fill( lInvMassOmegaMinus );
3019         if( lChargeXi > 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaPlus  ->Fill( lInvMassOmegaPlus  );
3020         if( lChargeXi < 0 )   fHistMassXiMinus    ->Fill( lInvMassXiMinus );
3021         if( lChargeXi > 0 )   fHistMassXiPlus     ->Fill( lInvMassXiPlus );
3022         if( lChargeXi < 0 )   fHistMassOmegaMinus ->Fill( lInvMassOmegaMinus );
3023         if( lChargeXi > 0 )   fHistMassOmegaPlus  ->Fill( lInvMassOmegaPlus );
3024         if(lIsBachelorPion)   f2dHistPIDprobaPionVsMCPtBach->Fill( lmcPtBach, ppionBach );
3025         if(lIsBachelorKaon)   f2dHistPIDprobaKaonVsMCPtBach->Fill( lmcPtBach, pkaonBach );
3026         if( lChargeXi < 0 && lIsBachelorMCPiMinus )    fHistMassWithMcPIDXiMinus     ->Fill( lInvMassXiMinus );
3027         if( lChargeXi > 0 && lIsBachelorMCPiPlus  )    fHistMassWithMcPIDXiPlus      ->Fill( lInvMassXiPlus );
3028         if( lChargeXi < 0 && lIsBachelorMCKMinus  )    fHistMassWithMcPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
3029         if( lChargeXi > 0 && lIsBachelorMCKPlus   )    fHistMassWithMcPIDOmegaPlus   ->Fill( lInvMassOmegaPlus );
3030
3031
3032         // - No association, skip the rest of the code
3033         if(!lAssoXiMinus && !lAssoXiPlus && !lAssoOmegaMinus && !lAssoOmegaPlus) continue; 
3034
3035         //--------------
3036         // - Proper time         
3037         // For cascade (reconstructed)   
3038         Double_t lctau = TMath::Sqrt(TMath::Power((lPosXi[0]-lBestPrimaryVtxPos[0]),2)+TMath::Power((lPosXi[1]-lBestPrimaryVtxPos[1]),2)+TMath::Power((lPosXi[2]-lBestPrimaryVtxPos[2]),2));
3039         if (lrecoP!=0) lctau = lctau*cascadeMass/lrecoP;   
3040         else lctau = -1.;
3041         // For Lambda (reconstructed)
3042         Float_t lambdaMass = 1.115683; // PDG mass 
3043         Float_t distV0Xi = TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2)+TMath::Power((lPosV0Xi[2]-lPosXi[2]),2)); 
3044         Float_t lctauV0 = -1.;
3045         if (lV0mom!=0) lctauV0 = distV0Xi*lambdaMass/lV0mom; 
3046         // Distance
3047         Float_t distTV0Xi = TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2));
3048        
3049         //------------------------------------------------------------
3050         // - Fill histos for the cascade candidates associated with MC
3051         if( lChargeXi < 0 && lAssoXiMinus){     
3052                 fHistAsMCMassXiMinus          ->Fill( lInvMassXiMinus  );
3053                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus->Fill( lmcPt, lmcRapCasc );
3054                 f2dHistAsMCGenPtVsGenYXiMinus ->Fill( lmcPt, lmcRapCasc);
3055                 fHistAsMCGenEtaXiMinus        ->Fill( lmcEta           );
3056                 f2dHistAsMCResPtXiMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3057                 f2dHistAsMCResRXiMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3058                 f2dHistAsMCResPhiXiMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3059                 f2dHistAsMCptProtonMCptXiMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3060                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3061         }       
3062         else if( lChargeXi > 0 && lAssoXiPlus){ 
3063                 fHistAsMCMassXiPlus           ->Fill( lInvMassXiPlus   );
3064                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus->Fill( lmcPt, lmcRapCasc );
3065                 f2dHistAsMCGenPtVsGenYXiPlus  ->Fill( lmcPt, lmcRapCasc);
3066                 fHistAsMCGenEtaXiPlus         ->Fill( lmcEta           );
3067                 f2dHistAsMCResPtXiPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3068                 f2dHistAsMCResRXiPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3069                 f2dHistAsMCResPhiXiPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3070                 f2dHistAsMCptAntiprotonMCptXiPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3071                 fHistV0CosineOfPointingAnglevsPtXi->Fill(lmcPt,lV0CosineOfPointingAngle);
3072         }
3073         else if( lChargeXi < 0 && lAssoOmegaMinus){     
3074                 fHistAsMCMassOmegaMinus          ->Fill( lInvMassOmegaMinus );
3075                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus->Fill( lmcPt, lmcRapCasc );
3076                 f2dHistAsMCGenPtVsGenYOmegaMinus ->Fill( lmcPt, lmcRapCasc  );
3077                 fHistAsMCGenEtaOmegaMinus        ->Fill( lmcEta             );
3078                 f2dHistAsMCResPtOmegaMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3079                 f2dHistAsMCResROmegaMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3080                 f2dHistAsMCResPhiOmegaMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
3081                 f2dHistAsMCptProtonMCptOmegaMinus->Fill(lmcPt,lmcPtPosV0Dghter);
3082                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3083         }       
3084         else if( lChargeXi > 0 && lAssoOmegaPlus){      
3085                 fHistAsMCMassOmegaPlus           ->Fill( lInvMassOmegaPlus );
3086                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus->Fill( lmcPt, lmcRapCasc );
3087                 f2dHistAsMCGenPtVsGenYOmegaPlus  ->Fill( lmcPt, lmcRapCasc   );
3088                 fHistAsMCGenEtaOmegaPlus         ->Fill( lmcEta            );
3089                 f2dHistAsMCResPtOmegaPlus        ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
3090                 f2dHistAsMCResROmegaPlus         ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
3091                 f2dHistAsMCResPhiOmegaPlus       ->Fill( lmcPt, lDeltaPhiMcReco );
3092                 f2dHistAsMCptAntiprotonMCptOmegaPlus->Fill(lmcPt,lmcPtNegV0Dghter);
3093                 fHistV0CosineOfPointingAnglevsPtOmega->Fill(lmcPt,lV0CosineOfPointingAngle);
3094         }
3095         fHistV0toXiCosineOfPointingAngle->Fill(lV0CosineOfPointingAngleXi);
3096
3097         //------------------         
3098         // - Fill containers
3099        
3100         // - Filling the AliCFContainer (optimisation of topological selections + systematics)
3101         Double_t lContainerCutVars[19] = {0.0};
3102         lContainerCutVars[0]  = lDcaXiDaughters;
3103         lContainerCutVars[1]  = lDcaBachToPrimVertexXi;
3104         lContainerCutVars[2]  = lXiCosineOfPointingAngle;
3105         lContainerCutVars[3]  = lXiRadius;
3106         lContainerCutVars[4]  = lInvMassLambdaAsCascDghter;
3107         lContainerCutVars[5]  = lDcaV0DaughtersXi;
3108         lContainerCutVars[6]  = lV0CosineOfPointingAngleXi;
3109         lContainerCutVars[7]  = lV0RadiusXi;
3110         lContainerCutVars[8]  = lDcaV0ToPrimVertexXi;   
3111         lContainerCutVars[9]  = lDcaPosToPrimVertexXi;
3112         lContainerCutVars[10] = lDcaNegToPrimVertexXi;
3113         lContainerCutVars[13] = lmcPt;
3114         lContainerCutVars[16] = lctau;
3115         lContainerCutVars[17] = lctauV0;
3116         lContainerCutVars[18] = distTV0Xi;
3117         // All cases should be covered below
3118         if( lChargeXi < 0 && lAssoXiMinus    ) {
3119                      lContainerCutVars[11] = lInvMassXiMinus;
3120                      lContainerCutVars[12] = lInvMassOmegaMinus;//1.63;
3121                      lContainerCutVars[14] = lmcRapCasc;
3122                      lContainerCutVars[15] = -1.;
3123                      if ( lIsBachelorPionForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )    
3124                        fCFContAsCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
3125         }
3126         if( lChargeXi > 0 && lAssoXiPlus    ) {
3127                      lContainerCutVars[11] = lInvMassXiPlus;
3128                      lContainerCutVars[12] = lInvMassOmegaPlus;//1.26;
3129                      lContainerCutVars[14] = lmcRapCasc;
3130                      lContainerCutVars[15] = -1.; 
3131                      if ( lIsBachelorPionForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )    
3132                        fCFContAsCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
3133         }
3134         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3135                      lContainerCutVars[11] = lInvMassXiMinus;//1.63;
3136                      lContainerCutVars[12] = lInvMassOmegaMinus;
3137                      lContainerCutVars[14] = -1.;
3138                      lContainerCutVars[15] = lmcRapCasc;
3139                      if ( lIsBachelorKaonForTPC   && lIsPosProtonForTPC    && lIsNegPionForTPC )    
3140                        fCFContAsCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
3141         }
3142         if( lChargeXi > 0 && lAssoOmegaPlus  ) {
3143                      lContainerCutVars[11] = lInvMassXiPlus;//1.26;
3144                      lContainerCutVars[12] = lInvMassOmegaPlus;
3145                      lContainerCutVars[14] = -1.;
3146                      lContainerCutVars[15] = lmcRapCasc;
3147                      if ( lIsBachelorKaonForTPC   && lIsNegProtonForTPC    && lIsPosPionForTPC )    
3148                        fCFContAsCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
3149         }
3150         
3151         // - Filling the AliCFContainers related to PID
3152         Double_t lContainerPIDVars[3] = {0.0};
3153         
3154         // Xi Minus             
3155         if( lChargeXi < 0 && lAssoXiMinus ) {
3156                 lContainerPIDVars[0] = lmcPt;
3157                 lContainerPIDVars[1] = lInvMassXiMinus;
3158                 lContainerPIDVars[2] = lmcRapCasc;
3159                 // No PID
3160                 fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 0); // No PID
3161                 // TPC PID
3162                 if( lIsBachelorPionForTPC ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3163                 if( lIsBachelorPionForTPC && lIsPosProtonForTPC ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks   
3164                 if( lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3165                 // Combined PID
3166                 if( lIsBachelorPion ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3167                 if( lIsBachelorPion && lIsPosInXiProton ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3168                 if( lIsBachelorPion && lIsPosInXiProton && lIsNegInXiPion ) fCFContCascadePIDAsXiMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3169         }       
3170         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; 
3171         
3172         // Xi Plus              
3173         if( lChargeXi > 0 && lAssoXiPlus ) {
3174                lContainerPIDVars[0] = lmcPt;
3175                lContainerPIDVars[1] = lInvMassXiPlus;
3176                lContainerPIDVars[2] = lmcRapCasc;
3177                // No PID
3178                fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 0); // No PID
3179                // TPC PID
3180                if( lIsBachelorPionForTPC ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3181                if( lIsBachelorPionForTPC && lIsNegProtonForTPC ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3182                if( lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3183                // Combined PID
3184                if( lIsBachelorPion ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3185                if( lIsBachelorPion && lIsNegInXiProton ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3186                if( lIsBachelorPion && lIsNegInXiProton && lIsPosInXiPion ) fCFContCascadePIDAsXiPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3187         }       
3188         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; 
3189         
3190         // Omega Minus          
3191         if( lChargeXi < 0 && lAssoOmegaMinus ) {
3192                lContainerPIDVars[0] = lmcPt;
3193                lContainerPIDVars[1] = lInvMassOmegaMinus;
3194                lContainerPIDVars[2] = lmcRapCasc;               
3195                // No PID
3196                fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 0); // No PID
3197                // TPC PID
3198                if( lIsBachelorKaonForTPC ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3199                if( lIsBachelorKaonForTPC && lIsPosProtonForTPC ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3200                if( lIsBachelorKaonForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3201                // Combined PID
3202                if( lIsBachelorKaon ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3203                if( lIsBachelorKaon && lIsPosInOmegaProton ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3204                if( lIsBachelorKaon && lIsPosInOmegaProton && lIsNegInOmegaPion ) fCFContCascadePIDAsOmegaMinus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3205         }       
3206         lContainerPIDVars[0] = 0.; lContainerPIDVars[1] = 0.; lContainerPIDVars[2] = 0.; 
3207         
3208         // Omega Plus           
3209         if( lChargeXi > 0 && lAssoOmegaPlus) {
3210                lContainerPIDVars[0] = lmcPt;
3211                lContainerPIDVars[1] = lInvMassOmegaPlus;
3212                lContainerPIDVars[2] = lmcRapCasc;               
3213                // No PID
3214                fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 0); // No PID
3215                // TPC PID
3216                if( lIsBachelorKaonForTPC ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 1); // TPC PID / 4-#sigma cut on Bachelor track
3217                      if( lIsBachelorKaonForTPC && lIsNegProtonForTPC ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 2); // TPC PID / 4-#sigma cut on Bachelor+Baryon tracks
3218                      if( lIsBachelorKaonForTPC && lIsNegProtonForTPC && lIsPosPionForTPC ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 3); // TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks
3219                      // Combined PID
3220                      if( lIsBachelorKaon ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 4); // Comb. PID / Bachelor
3221                      if( lIsBachelorKaon && lIsNegInOmegaProton ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 5); // Comb. PID / Bachelor+Baryon
3222                      if( lIsBachelorKaon && lIsNegInOmegaProton && lIsPosInOmegaPion ) fCFContCascadePIDAsOmegaPlus->Fill(lContainerPIDVars, 6); // Comb. PID / Bachelor+Baryon+Meson
3223         }       
3224         
3225    }// End of loop over reconstructed cascades
3226  
3227    fHistnAssoXiMinus->Fill(nAssoXiMinus);
3228    fHistnAssoXiPlus->Fill(nAssoXiPlus);
3229    fHistnAssoOmegaMinus->Fill(nAssoOmegaMinus);
3230    fHistnAssoOmegaPlus->Fill(nAssoOmegaPlus);  
3231  
3232    // Post output data.
3233    PostData(1, fListHistCascade);
3234    PostData(2, fCFContCascadePIDAsXiMinus);
3235    PostData(3, fCFContCascadePIDAsXiPlus);
3236    PostData(4, fCFContCascadePIDAsOmegaMinus);
3237    PostData(5, fCFContCascadePIDAsOmegaPlus);
3238    PostData(6, fCFContAsCascadeCuts);
3239
3240 }      
3241
3242
3243 //________________________________________________________________________
3244 void AliAnalysisTaskCheckPerformanceCascadepp276::Terminate(Option_t *) {
3245   // Draw result to the screen
3246   // Called once at the end of the query
3247         
3248  /* TList *cRetrievedList = 0x0;
3249   cRetrievedList = (TList*)GetOutputData(1);
3250   if(!cRetrievedList) {
3251         Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadepp276 : ouput data container list not available\n");
3252         return;
3253   }     
3254         
3255   fHistTrackMultiplicityBeforeAnySel = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistTrackMultiplicityBeforeAnySel")  );
3256   if (!fHistTrackMultiplicityBeforeAnySel) {
3257     Printf("ERROR - AliAnalysisTaskCheckPerformanceCascadepp276 : fHistTrackMultiplicityBeforeAnySel not available");
3258     return;
3259   }
3260   
3261    
3262   TCanvas *canCheckPerformanceCascade = new TCanvas("AliAnalysisTaskCheckPerformanceCascadepp276","Multiplicity",10,10,510,510);
3263   canCheckPerformanceCascade->cd(1)->SetLogy();
3264
3265   fHistTrackMultiplicityBeforeAnySel->SetMarkerStyle(22);
3266   fHistTrackMultiplicityBeforeAnySel->DrawCopy("E");
3267  */
3268 }