Bug fix in AliExternalTrackParam::Invert()
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskCheckPerformanceCascade.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 //               AliAnalysisTaskCheckPerformanceCascade class
16 //            This task is for a performance study of cascade identification.
17 //            It works with MC info and ESD.
18 //              Use with AOD tree = under development
19 //            Origin   : A.Maire Mar2009, antonin.maire@ires.in2p3.fr
20 //            Modified : A.Maire Nov2010, antonin.maire@ires.in2p3.fr
21 //-----------------------------------------------------------------
22
23
24 #include <Riostream.h>
25
26 #include "TList.h"
27 #include "TFile.h"
28 #include "TH1F.h"
29 #include "TH2F.h"
30 #include "TVector3.h"
31 #include "TCanvas.h"
32 #include "TParticle.h"
33 #include "TMath.h"
34
35 #include "AliLog.h"
36 #include "AliHeader.h"
37 #include "AliMCEvent.h"
38 #include "AliStack.h"
39 #include "AliMultiplicity.h"
40 #include "AliInputEventHandler.h"
41 #include "AliAnalysisManager.h"
42
43 #include "AliCFContainer.h"
44 #include "AliESDpid.h"
45 #include "AliESDtrackCuts.h"
46 //   #include "AliV0vertexer.h"
47 //   #include "AliCascadeVertexer.h"
48 #include "AliESDEvent.h"
49 #include "AliESDcascade.h"
50
51 #include "AliAODEvent.h"
52
53 #include "AliAnalysisTaskCheckPerformanceCascade.h"
54
55 ClassImp(AliAnalysisTaskCheckPerformanceCascade)
56
57
58
59      //_____Dummy constructor________________________________________________________________
60 AliAnalysisTaskCheckPerformanceCascade::AliAnalysisTaskCheckPerformanceCascade() 
61 : AliAnalysisTaskSE(), // <- take care to AliAnalysisTask( empty )
62   fDebugCascade(0), fAnalysisType("ESD"), fTriggerMaskType("kMB"), fCollidingSystems(0), fESDpid(0), fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
63     fkRerunV0CascVertexers         (0),
64     fkQualityCutZprimVtxPos        (kTRUE),
65     fkRejectEventPileUp            (kTRUE),
66     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
67     fkQualityCutTPCrefit           (kTRUE),
68     fkQualityCut80TPCcls           (kTRUE),
69     fkIsDataRecoWith1PadTPCCluster (kTRUE),
70     fkExtraSelections              (0),
71     
72         // - Cascade part initialisation
73 fListHistCascade(0),
74     fHistMCTrackMultiplicity(0), 
75        // - Resolution of the multiplicity estimator
76     f2dHistRecoPrimTrckMultVsMCMult(0), f2dHistRecoEstimateMultVsMCMult(0),
77
78     fHistEtaGenProton(0),
79     fHistEtaGenAntiProton(0),
80     
81    // Xi-
82    fHistEtaGenCascXiMinus(0),
83    f2dHistGenPtVsGenYGenXiMinus(0),
84    
85     fHistThetaGenCascXiMinus(0),
86     f2dHistGenPtVsGenYFdblXiMinus(0),
87     
88     fHistThetaLambdaXiMinus(0), 
89     fHistThetaBachXiMinus(0),
90     
91     fHistThetaMesDghterXiMinus(0), 
92     fHistThetaBarDghterXiMinus(0),
93     
94     fHistPtBachXiMinus(0),
95     fHistPtMesDghterXiMinus(0),
96     fHistPtBarDghterXiMinus(0),
97    
98    
99    // Xi+
100    fHistEtaGenCascXiPlus(0),
101    f2dHistGenPtVsGenYGenXiPlus(0),
102    
103     fHistThetaGenCascXiPlus(0), 
104     f2dHistGenPtVsGenYFdblXiPlus(0),
105     
106     fHistThetaLambdaXiPlus(0), 
107     fHistThetaBachXiPlus(0),
108     
109     fHistThetaMesDghterXiPlus(0), 
110     fHistThetaBarDghterXiPlus(0),
111     
112     fHistPtBachXiPlus(0),
113     fHistPtMesDghterXiPlus(0),
114     fHistPtBarDghterXiPlus(0),
115    
116    // Omega-
117    fHistEtaGenCascOmegaMinus(0),
118    f2dHistGenPtVsGenYGenOmegaMinus(0),
119    
120     fHistThetaGenCascOmegaMinus(0),
121     f2dHistGenPtVsGenYFdblOmegaMinus(0),
122     
123     fHistThetaLambdaOmegaMinus(0), 
124     fHistThetaBachOmegaMinus(0),
125     
126     fHistThetaMesDghterOmegaMinus(0), 
127     fHistThetaBarDghterOmegaMinus(0),
128     
129     fHistPtBachOmegaMinus(0),
130     fHistPtMesDghterOmegaMinus(0),
131     fHistPtBarDghterOmegaMinus(0),
132    
133    
134    // Omega+
135    fHistEtaGenCascOmegaPlus(0),
136    f2dHistGenPtVsGenYGenOmegaPlus(0),
137    
138     fHistThetaGenCascOmegaPlus(0),
139     f2dHistGenPtVsGenYFdblOmegaPlus(0),
140     
141     fHistThetaLambdaOmegaPlus(0), 
142     fHistThetaBachOmegaPlus(0),
143     
144     fHistThetaMesDghterOmegaPlus(0), 
145     fHistThetaBarDghterOmegaPlus(0),
146     
147     fHistPtBachOmegaPlus(0),
148     fHistPtMesDghterOmegaPlus(0),
149     fHistPtBarDghterOmegaPlus(0),
150
151 // Part 2 - Association to MC
152         
153     fHistMassXiMinus(0),
154     fHistMassXiPlus(0),
155     fHistMassOmegaMinus(0),
156     fHistMassOmegaPlus(0),
157     
158         // - Effective mass histos with combined PID
159     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
160     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
161         
162         // - PID Probability versus MC Pt(bachelor track)
163     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
164     
165         // - Effective mass histos with perfect MC PID on the bachelor
166     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
167     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
168         
169     
170         // - Effective mass histos for the cascade candidates associated with MC
171     fHistAsMCMassXiMinus(0),            
172     fHistAsMCMassXiPlus(0),             
173     fHistAsMCMassOmegaMinus(0),
174     fHistAsMCMassOmegaPlus(0),
175     
176         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
177     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
178     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
179     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0),
180     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
181            
182         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
183     f2dHistAsMCGenPtVsGenYXiMinus(0),
184     f2dHistAsMCGenPtVsGenYXiPlus(0),
185     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
186     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
187     
188         // - Generated Eta of the the cascade candidates associated with MC
189     fHistAsMCGenEtaXiMinus(0),
190     fHistAsMCGenEtaXiPlus(0),
191     fHistAsMCGenEtaOmegaMinus(0),
192     fHistAsMCGenEtaOmegaPlus(0),
193         
194         // - Resolution in Pt as function of generated Pt
195     f2dHistAsMCResPtXiMinus(0),         
196     f2dHistAsMCResPtXiPlus(0),          
197     f2dHistAsMCResPtOmegaMinus(0),
198     f2dHistAsMCResPtOmegaPlus(0),       
199         
200         // - Resolution in R(2D) as function of generated R
201     f2dHistAsMCResRXiMinus(0),          
202     f2dHistAsMCResRXiPlus(0),           
203     f2dHistAsMCResROmegaMinus(0),
204     f2dHistAsMCResROmegaPlus(0),
205
206         // - Resolution in phi as function of generated Pt
207     f2dHistAsMCResPhiXiMinus(0),
208     f2dHistAsMCResPhiXiPlus(0),
209     f2dHistAsMCResPhiOmegaMinus(0),
210     f2dHistAsMCResPhiOmegaPlus(0),
211
212         // - Correlation in Pt between the cascade and its (anti)proton daughter
213     f2dHistAsMCPtProtonVsPtXiMinus(0),
214     f2dHistAsMCPtAntiProtonVsPtXiPlus(0),
215     f2dHistAsMCPtProtonVsPtOmegaMinus(0),
216     f2dHistAsMCPtAntiProtonVsPtOmegaPlus(0),
217
218     fCFContCascadePIDAsXiMinus(0),
219     fCFContCascadePIDAsXiPlus(0),
220     fCFContCascadePIDAsOmegaMinus(0),
221     fCFContCascadePIDAsOmegaPlus(0),
222
223     fCFContAsCascadeCuts(0)
224
225 {
226 // Dummy constructor
227         for(Int_t iAlephIdx   = 0; iAlephIdx   < 5; iAlephIdx++   ) { fAlephParameters [iAlephIdx]    = -1.; }
228         for(Int_t iV0selIdx   = 0; iV0selIdx   < 7; iV0selIdx++   ) { fV0Sels          [iV0selIdx   ] = -1.; }
229         for(Int_t iCascSelIdx = 0; iCascSelIdx < 8; iCascSelIdx++ ) { fCascSels        [iCascSelIdx ] = -1.; }
230 }
231      
232        
233      
234      
235 //_____Non-default Constructor________________________________________________________________
236 AliAnalysisTaskCheckPerformanceCascade::AliAnalysisTaskCheckPerformanceCascade(const char *name) 
237   : AliAnalysisTaskSE(name),
238     fDebugCascade(0), fAnalysisType("ESD"), fTriggerMaskType("kMB"), fCollidingSystems(0), fESDpid(0), fESDtrackCuts(0), /*fPaveTextBookKeeping(0),*/
239     fkRerunV0CascVertexers         (0),
240     fkQualityCutZprimVtxPos        (kTRUE),
241     fkRejectEventPileUp            (kTRUE),
242     fkQualityCutNoTPConlyPrimVtx   (kTRUE),
243     fkQualityCutTPCrefit           (kTRUE),
244     fkQualityCut80TPCcls           (kTRUE),
245     fkIsDataRecoWith1PadTPCCluster (kTRUE),
246     fkExtraSelections              (0),
247       
248         // - Cascade part initialisation
249 fListHistCascade(0),
250     fHistMCTrackMultiplicity(0), 
251        // - Resolution of the multiplicity estimator
252     f2dHistRecoPrimTrckMultVsMCMult(0), f2dHistRecoEstimateMultVsMCMult(0),
253
254     fHistEtaGenProton(0),
255     fHistEtaGenAntiProton(0),
256     
257 // Xi-
258    fHistEtaGenCascXiMinus(0),
259    f2dHistGenPtVsGenYGenXiMinus(0),
260    
261     fHistThetaGenCascXiMinus(0),
262     f2dHistGenPtVsGenYFdblXiMinus(0),
263     
264     fHistThetaLambdaXiMinus(0), 
265     fHistThetaBachXiMinus(0),
266     
267     fHistThetaMesDghterXiMinus(0), 
268     fHistThetaBarDghterXiMinus(0),
269     
270     fHistPtBachXiMinus(0),
271     fHistPtMesDghterXiMinus(0),
272     fHistPtBarDghterXiMinus(0),
273    
274    
275    // Xi+
276    fHistEtaGenCascXiPlus(0),
277    f2dHistGenPtVsGenYGenXiPlus(0),
278    
279     fHistThetaGenCascXiPlus(0), 
280     f2dHistGenPtVsGenYFdblXiPlus(0),
281     
282     fHistThetaLambdaXiPlus(0), 
283     fHistThetaBachXiPlus(0),
284     
285     fHistThetaMesDghterXiPlus(0), 
286     fHistThetaBarDghterXiPlus(0),
287     
288     fHistPtBachXiPlus(0),
289     fHistPtMesDghterXiPlus(0),
290     fHistPtBarDghterXiPlus(0),
291    
292    // Omega-
293    fHistEtaGenCascOmegaMinus(0),
294    f2dHistGenPtVsGenYGenOmegaMinus(0),
295    
296     fHistThetaGenCascOmegaMinus(0),
297     f2dHistGenPtVsGenYFdblOmegaMinus(0),
298     
299     fHistThetaLambdaOmegaMinus(0), 
300     fHistThetaBachOmegaMinus(0),
301     
302     fHistThetaMesDghterOmegaMinus(0), 
303     fHistThetaBarDghterOmegaMinus(0),
304     
305     fHistPtBachOmegaMinus(0),
306     fHistPtMesDghterOmegaMinus(0),
307     fHistPtBarDghterOmegaMinus(0),
308    
309    
310    // Omega+
311    fHistEtaGenCascOmegaPlus(0),
312    f2dHistGenPtVsGenYGenOmegaPlus(0),
313    
314     fHistThetaGenCascOmegaPlus(0),
315     f2dHistGenPtVsGenYFdblOmegaPlus(0),
316     
317     fHistThetaLambdaOmegaPlus(0), 
318     fHistThetaBachOmegaPlus(0),
319     
320     fHistThetaMesDghterOmegaPlus(0), 
321     fHistThetaBarDghterOmegaPlus(0),
322     
323     fHistPtBachOmegaPlus(0),
324     fHistPtMesDghterOmegaPlus(0),
325     fHistPtBarDghterOmegaPlus(0),
326
327 // Part 2 - Association to MC
328         
329     fHistMassXiMinus(0),
330     fHistMassXiPlus(0),
331     fHistMassOmegaMinus(0),
332     fHistMassOmegaPlus(0),
333     
334         // - Effective mass histos with combined PID
335     fHistMassWithCombPIDXiMinus(0), fHistMassWithCombPIDXiPlus(0),
336     fHistMassWithCombPIDOmegaMinus(0), fHistMassWithCombPIDOmegaPlus(0),
337
338         // - PID Probability versus MC Pt(bachelor track)
339     f2dHistPIDprobaKaonVsMCPtBach(0), f2dHistPIDprobaPionVsMCPtBach(0),
340     
341         // - Effective mass histos with perfect MC PID on the bachelor
342     fHistMassWithMcPIDXiMinus(0), fHistMassWithMcPIDXiPlus(0),
343     fHistMassWithMcPIDOmegaMinus(0), fHistMassWithMcPIDOmegaPlus(0),
344         
345     
346         // - Effective mass histos for the cascade candidates associated with MC
347     fHistAsMCMassXiMinus(0),            
348     fHistAsMCMassXiPlus(0),             
349     fHistAsMCMassOmegaMinus(0),
350     fHistAsMCMassOmegaPlus(0),
351     
352         // - Generated Pt Vs generated y, for the cascade candidates associated with MC + Info Comb. PID
353     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus(0),
354     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus(0),
355     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus(0), 
356     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus(0),
357     
358         // - Generated Pt Vs generated y, for the cascade candidates associated with MC
359     f2dHistAsMCGenPtVsGenYXiMinus(0),
360     f2dHistAsMCGenPtVsGenYXiPlus(0),
361     f2dHistAsMCGenPtVsGenYOmegaMinus(0),
362     f2dHistAsMCGenPtVsGenYOmegaPlus(0),
363     
364         // - Generated Eta of the the cascade candidates associated with MC
365     fHistAsMCGenEtaXiMinus(0),
366     fHistAsMCGenEtaXiPlus(0),
367     fHistAsMCGenEtaOmegaMinus(0),
368     fHistAsMCGenEtaOmegaPlus(0),
369         
370         // - Resolution in Pt as function of generated Pt
371     f2dHistAsMCResPtXiMinus(0),         
372     f2dHistAsMCResPtXiPlus(0),          
373     f2dHistAsMCResPtOmegaMinus(0),
374     f2dHistAsMCResPtOmegaPlus(0),       
375         
376         // - Resolution in R(2D) as function of generated R
377     f2dHistAsMCResRXiMinus(0),          
378     f2dHistAsMCResRXiPlus(0),           
379     f2dHistAsMCResROmegaMinus(0),
380     f2dHistAsMCResROmegaPlus(0),
381
382        // - Resolution in phi as function of generated Pt
383     f2dHistAsMCResPhiXiMinus(0),
384     f2dHistAsMCResPhiXiPlus(0),
385     f2dHistAsMCResPhiOmegaMinus(0),
386     f2dHistAsMCResPhiOmegaPlus(0),
387
388         // - Correlation in Pt between the cascade and its (anti)proton daughter
389     f2dHistAsMCPtProtonVsPtXiMinus(0),
390     f2dHistAsMCPtAntiProtonVsPtXiPlus(0),
391     f2dHistAsMCPtProtonVsPtOmegaMinus(0),
392     f2dHistAsMCPtAntiProtonVsPtOmegaPlus(0),
393
394     fCFContCascadePIDAsXiMinus(0),
395     fCFContCascadePIDAsXiPlus(0),
396     fCFContCascadePIDAsOmegaMinus(0),
397     fCFContCascadePIDAsOmegaPlus(0),
398
399     fCFContAsCascadeCuts(0)
400
401 {
402   // Constructor
403
404   // Define input and output slots here
405   // Input slot #0 works with a TChain
406   // Output slot #1 writes into a TList container (cascade)
407         
408         for(Int_t iAlephIdx   = 0; iAlephIdx   < 5; iAlephIdx++   ) { fAlephParameters [iAlephIdx]    = -1.; }
409         
410         // New Loose : 1st step for the 7 TeV pp analysis
411         
412         fV0Sels[0] =  33.  ;  // max allowed chi2
413         fV0Sels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
414         fV0Sels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
415         fV0Sels[3] =   2.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
416         fV0Sels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
417         fV0Sels[5] =   1.0 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
418         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
419         
420         fCascSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
421         fCascSels[1] =   0.05 ;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
422         fCascSels[2] =   0.010;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
423         fCascSels[3] =   0.03 ;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
424         fCascSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
425         fCascSels[5] =   0.95 ;  // min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
426         fCascSels[6] =   0.4  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
427         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
428         
429         
430         // Hyper Loose "à la 900 GeV 2009 data", with lower cosine of pointing angle for Xi (0.95 down to 0.82) = 900 GeV paper
431         /*
432         fV0Sels[0] =  33.  ;  // max allowed chi2
433         fV0Sels[1] =   0.001; // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
434         fV0Sels[2] =   0.001; // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
435         fV0Sels[3] =   5.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
436         fV0Sels[4] =   0.0 ;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
437         fV0Sels[5] =   0.1 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
438         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
439         
440         fCascSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
441         fCascSels[1] =   0.001;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
442         fCascSels[2] =   0.008;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
443         fCascSels[3] =   0.001;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
444         fCascSels[4] =   5.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
445         fCascSels[5] =   0.82 ;  //FIXME min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
446         fCascSels[6] =   0.1  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
447         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
448         */
449         
450         //New default vtxR (http://alisoft.cern.ch/viewvc?view=rev&root=AliRoot&revision=40955, 5 May 2010)        
451         /*
452         fV0Sels[0] =  33.  ;  // max allowed chi2
453         fV0Sels[1] =   0.05;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
454         fV0Sels[2] =   0.05;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
455         fV0Sels[3] =   1.5 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
456         fV0Sels[4] =   0.9 ;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
457         fV0Sels[5] =   0.2 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
458         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
459         
460         fCascSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
461         fCascSels[1] =   0.01 ;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
462         fCascSels[2] =   0.008;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
463         fCascSels[3] =   0.01 ;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
464         fCascSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
465         fCascSels[5] =   0.98 ;  // min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
466         fCascSels[6] =   0.2  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
467         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
468         */
469         
470   DefineOutput(1, TList::Class());
471  
472 }
473
474
475 AliAnalysisTaskCheckPerformanceCascade::~AliAnalysisTaskCheckPerformanceCascade()
476 {
477   //
478   // Destructor
479   //
480
481   // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
482   // They will be deleted when fListCascade is deleted by the TSelector dtor
483   // Because of TList::SetOwner()
484
485   if (fListHistCascade)      { delete fListHistCascade;     fListHistCascade = 0x0;  }  
486   if (fESDpid)               { delete fESDpid;              fESDpid = 0x0;} // fESDpid is not stored into the TList
487   if (fESDtrackCuts)         { delete fESDtrackCuts;        fESDtrackCuts = 0x0; }
488   /*if (fPaveTextBookKeeping)  { delete fPaveTextBookKeeping; fPaveTextBookKeeping = 0x0; } // fPaveTextBookKeeping is not stored into the TList*/
489 }
490
491
492 //________________________________________________________________________
493 void AliAnalysisTaskCheckPerformanceCascade::UserCreateOutputObjects()
494 {
495   // Create histograms
496   // Called once
497
498         
499    // Option for AliLog
500         AliLog::SetGlobalLogLevel(AliLog::kError); 
501         // to suppress the extensive info prompted by a run with MC                     
502
503    // Definition of the output datamembers      
504    fListHistCascade = new TList();
505    fListHistCascade->SetOwner(); // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
506         
507 if(! fESDpid){
508
509          if(fkIsDataRecoWith1PadTPCCluster){
510                 // Home made parameterization for LHC10f6a production = p+p 7 TeV
511                 fAlephParameters[0] = 0.04;
512                 fAlephParameters[1] = 17.5;
513                 fAlephParameters[2] = 3.4e-09;
514                 fAlephParameters[3] = 2.15;
515                 fAlephParameters[4] = 3.91720e+00;
516                 
517                 // Home made parameterization for LHC10e13 production = p+p 900 GeV/c
518          }
519          else {
520                 // Reasonable parameters extracted for p-p simulation (LHC09a4) - A.Kalweit
521                  // fAlephParameters[0] = 4.23232575531564326e+00/50;//50*0.76176e-1;
522                  // fAlephParameters[1] = 8.68482806165147636e+00;//10.632; 
523                  // fAlephParameters[2] = 1.34000000000000005e-05;//0.13279e-4;
524                  // fAlephParameters[3] = 2.30445734159456084e+00;//1.8631;
525                  // fAlephParameters[4] = 2.25624744086878559e+00;//1.9479;
526                 
527                 // Param for LHC09d10 prod - A.Kalweit
528                 fAlephParameters[0] = 2.15898e+00/50.;
529                 fAlephParameters[1] = 1.75295e+01;
530                 fAlephParameters[2] = 3.40030e-09;
531                 fAlephParameters[3] = 1.96178e+00;
532                 fAlephParameters[4] = 3.91720e+00; 
533          }
534          Printf("CheckPerfCascade - Check Aleph Param in case of MC Data   (fAlephParameters[3] = %f) (To be compared with : 2.15 for 1-pad-cluster prod. / 1.96178 otherwise)\n",  fAlephParameters[3]);
535         
536   fESDpid = new AliESDpid();
537   fESDpid->GetTPCResponse().SetBetheBlochParameters( fAlephParameters[0],
538                                                      fAlephParameters[1],
539                                                      fAlephParameters[2],
540                                                      fAlephParameters[3],
541                                                      fAlephParameters[4] );
542 }
543
544
545 if(! fESDtrackCuts ){
546       fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // Std definition of primary (see kTRUE argument) tracks for 2010
547       fESDtrackCuts->SetEtaRange(-0.8,+0.8);
548       fESDtrackCuts->SetPtRange(0.15, 1e10);
549       Printf("CheckCascade - ESDtrackCuts set up to 2010 std ITS-TPC cuts...");
550 }
551
552
553 /*
554 if( !fPaveTextBookKeeping){
555         fPaveTextBookKeeping = new TPaveText(0.1, 0.1, 0.9, 0.9,"NDC");
556         fPaveTextBookKeeping->SetName("fPaveTextBookKeeping");
557         fPaveTextBookKeeping->SetBorderSize(0);
558         fPaveTextBookKeeping->SetTextAlign(12);
559         fPaveTextBookKeeping->SetFillColor(kWhite);
560         fPaveTextBookKeeping->SetTextFont(42);        // regular Arial or Helvetica,
561         fPaveTextBookKeeping->SetTextColor(kGray+3);
562         
563         
564         fPaveTextBookKeeping->AddText( "Task CHECK PERFORMANCE CASCADE analysis" );
565         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
566         fPaveTextBookKeeping->AddText( Form("AnalysisType : %s ", fAnalysisType.Data() ));
567         if(!fCollidingSystems)  fPaveTextBookKeeping->AddText("Colliding system : p-p collisions ");
568         else                    fPaveTextBookKeeping->AddText("Colliding system : A-A collisions ");
569
570         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
571     
572         if(fkRerunV0CascVertexers){
573                 fPaveTextBookKeeping->AddText("A.1. With V0 vertexer : ");
574                 fPaveTextBookKeeping->AddText( Form("  - V0 #chi^{2} _________________ <  %.3f ",               fV0Sels[0] ));
575                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 1^{st} daughter) ___ >  %.3f     cm ",  fV0Sels[1] ));
576                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ 2^{nd} daughter) __  >  %.3f     cm",   fV0Sels[2] ));
577                 fPaveTextBookKeeping->AddText( Form("  - DCA between V0 daughters ___ <  %.3f      cm",         fV0Sels[3] ));
578                 fPaveTextBookKeeping->AddText( Form("  - cos(V0 pointing angle) _______ >  %.3f ",              fV0Sels[4] ));
579                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ >  %.3f             cm", fV0Sels[5] ));
580                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(V0 decay) ________ <  %.3f         cm",     fV0Sels[6] ));
581                 
582                 fPaveTextBookKeeping->AddText(" "); 
583                 
584                 fPaveTextBookKeeping->AddText("A.2. With Casc. vertexer : ");
585                 fPaveTextBookKeeping->AddText( Form("  - Casc. #chi^{2} ______________  <  %.3f ",                               fCascSels[0] ));
586                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ V0) _________ >  %.3f    cm",                            fCascSels[1] ));
587                 fPaveTextBookKeeping->AddText( Form("  - | M_{#Lambda}(reco) - M_{#Lambda}(pdg) | _______ <  %.3f    GeV/c^{2}", fCascSels[2] ));
588                 fPaveTextBookKeeping->AddText( Form("  - DCA(prim. Vtx/ Bach) _______ >  %.3f    cm",                            fCascSels[3] ));
589                 fPaveTextBookKeeping->AddText( Form("  - DCA between Bach/ #Lambda ______ <  %.3f    cm",                        fCascSels[4] ));
590                 fPaveTextBookKeeping->AddText( Form("  - cos(Casc. pointing angle) ____ >  %.3f ",                               fCascSels[5] ));
591                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ >  %.3f         cm",                     fCascSels[6] ));
592                 fPaveTextBookKeeping->AddText( Form("  - R_{transv}(Casc. decay) ______ <  %.3f     cm",                         fCascSels[7] ));
593         }
594         else{   fPaveTextBookKeeping->AddText("A. No rerunning of the V0/Casc. vertexers ... See std cuts in (AliRoot+Rec.C) used for this prod. cycle");}
595
596         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
597         
598         if(fkQualityCutZprimVtxPos)      fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = ON  ");
599         else                             fPaveTextBookKeeping->AddText("B. Quality Cut(prim. Vtx z-Pos)    = Off ");
600         if(fkQualityCutNoTPConlyPrimVtx) fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = ON  ");
601         else                             fPaveTextBookKeeping->AddText("C. Quality Cut(No TPC prim. vtx) = Off ");
602         if(fkQualityCutTPCrefit)         fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = ON  ");
603         else                             fPaveTextBookKeeping->AddText("D. Quality Cut(TPCrefit)               = Off ");
604         if(fkQualityCut80TPCcls)         fPaveTextBookKeeping->AddText("E. Quality Cut(80 TPC clusters)   = ON  ");
605         else                             fPaveTextBookKeeping->AddText("E. Quality Cut(80 TPC clusters)   = Off ");
606         if(fkExtraSelections)            fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = ON  ");
607         else                             fPaveTextBookKeeping->AddText("F. Extra Analysis Selections         = Off ");
608
609         fPaveTextBookKeeping->AddText("- - - - - - - - - - - ");
610
611         fPaveTextBookKeeping->AddText("G. TPC Aleph Param : ");
612         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [0] =  %.5g", fAlephParameters[0] ));
613         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [1] =  %.5g", fAlephParameters[1] ));
614         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [2] =  %.5g", fAlephParameters[2] ));
615         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [3] =  %.5g", fAlephParameters[3] ));
616         fPaveTextBookKeeping->AddText( Form("   - fAlephParam [4] =  %.5g", fAlephParameters[4] ));
617         
618         fListHistCascade->Add(fPaveTextBookKeeping);
619 }       
620 */
621                 
622   // - General
623   
624   if (!fHistMCTrackMultiplicity) {
625      fHistMCTrackMultiplicity = new TH1F("fHistMCTrackMultiplicity", "MC Track Multiplicity;Number of MC tracks;Events", 100, 0, 500);
626    //  fHistMCTrackMultiplicity = new TH1F("fHistMCTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 200, 0, 40000); //HERE
627   fListHistCascade->Add(fHistMCTrackMultiplicity);
628   }
629   
630     // - Resolution of the multiplicity estimator
631   if(! f2dHistRecoPrimTrckMultVsMCMult){
632        f2dHistRecoPrimTrckMultVsMCMult = new TH2F("f2dHistRecoPrimTrckMultVsMCMult", "Resolution of the multiplicity estimator (prim. tracks in |#eta| < 0.8); Reco Multiplicity (prim. tracks); MC multiplicity (gen. part. in |#eta| < 0.8)", 120, 0., 120., 300, 0., 300.);
633        fListHistCascade->Add(f2dHistRecoPrimTrckMultVsMCMult);
634   }
635   
636   if(! f2dHistRecoEstimateMultVsMCMult){
637        f2dHistRecoEstimateMultVsMCMult = new TH2F("f2dHistRecoEstimateMultVsMCMult", "Resolution of the multiplicity estimator (EstimateMult. in |#eta| < 1.0); Reco Multiplicity (tr(ITS-TPC)+ITSsa+tracklets); MC multiplicity (gen. part. in |#eta| < 1.0)", 160, 0., 160., 300, 0., 300.);
638        fListHistCascade->Add(f2dHistRecoEstimateMultVsMCMult);
639   }
640   
641   if (!fHistEtaGenProton) {
642      fHistEtaGenProton = new TH1F("fHistEtaGenProton", "#eta of any gen. p^{+};#eta;Number of prim. protons", 200, -10, 10);
643      fListHistCascade->Add(fHistEtaGenProton);
644   }
645   
646   if (!fHistEtaGenAntiProton) {
647      fHistEtaGenAntiProton = new TH1F("fHistEtaGenAntiProton", "#eta of any gen. #bar{p}^{-};#eta;Number of prim. #bar{p}", 200, -10, 10);
648      fListHistCascade->Add(fHistEtaGenAntiProton);
649   }
650   
651
652
653
654
655
656   
657   //--------
658   // I - Xi- 
659   // - Pseudo-Rapidity distribution
660   if (!fHistEtaGenCascXiMinus) {
661      fHistEtaGenCascXiMinus = new TH1F("fHistEtaGenCascXiMinus", "#eta of any gen. #Xi^{-};#eta;Number of Casc", 200, -10, 10);
662      fListHistCascade->Add(fHistEtaGenCascXiMinus);
663   }
664   
665   if (!f2dHistGenPtVsGenYGenXiMinus) {
666      f2dHistGenPtVsGenYGenXiMinus = new TH2F("f2dHistGenPtVsGenYGenXiMinus", "MC P_{t} Vs MC Y of Gen #Xi^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
667      fListHistCascade->Add(f2dHistGenPtVsGenYGenXiMinus);
668   }
669   
670   
671   // - Info at the generation level of multi-strange particle
672   
673   if (!fHistThetaGenCascXiMinus) {
674      fHistThetaGenCascXiMinus = new TH1F("fHistThetaGenCascXiMinus", "#theta of gen. #Xi^{-};#theta;Number of Casc.", 200, -10, 190);
675      fListHistCascade->Add(fHistThetaGenCascXiMinus);
676   }
677
678   if (!f2dHistGenPtVsGenYFdblXiMinus) {
679      f2dHistGenPtVsGenYFdblXiMinus = new TH2F("f2dHistGenPtVsGenYFdblXiMinus", "MC P_{t} Vs MC Y of findable Gen #Xi^{-}; Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
680      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiMinus);
681   }
682   
683                 // - Theta distribution the daughters (control plots)
684   
685   if (!fHistThetaLambdaXiMinus) {
686      fHistThetaLambdaXiMinus = new TH1F("fHistThetaLambdaXiMinus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda^0", 200, -10, 190);
687      fListHistCascade->Add(fHistThetaLambdaXiMinus);
688   }
689
690   if (!fHistThetaBachXiMinus) {
691      fHistThetaBachXiMinus = new TH1F("fHistThetaBachXiMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
692      fListHistCascade->Add(fHistThetaBachXiMinus);
693   }
694   
695   if (!fHistThetaMesDghterXiMinus) {
696      fHistThetaMesDghterXiMinus = new TH1F("fHistThetaMesDghterXiMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
697      fListHistCascade->Add(fHistThetaMesDghterXiMinus);
698   }
699   
700   if (!fHistThetaBarDghterXiMinus) {
701      fHistThetaBarDghterXiMinus = new TH1F("fHistThetaBarDghterXiMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
702      fListHistCascade->Add(fHistThetaBarDghterXiMinus);
703   }
704  
705                 // - Pt distribution (control plots)
706     
707   if (!fHistPtBachXiMinus) {
708      fHistPtBachXiMinus = new TH1F("fHistPtBachXiMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
709      fListHistCascade->Add(fHistPtBachXiMinus);
710   }
711   
712   if (!fHistPtMesDghterXiMinus) {
713      fHistPtMesDghterXiMinus = new TH1F("fHistPtMesDghterXiMinus", "p_{t} of gen. Meson #Lambda dghter;pt_{MesDght};Number of Mes.", 200, 0, 10);
714      fListHistCascade->Add(fHistPtMesDghterXiMinus);
715   }
716     
717   if (!fHistPtBarDghterXiMinus) {
718      fHistPtBarDghterXiMinus = new TH1F("fHistPtBarDghterXiMinus", "p_{t} of gen. Baryon #Lambda dghter;pt_{BarDght};Number of Bar.", 200, 0, 10);
719      fListHistCascade->Add(fHistPtBarDghterXiMinus);
720   }
721   
722   
723   
724   //--------
725   // II - Xi+ 
726   // - Pseudo-Rapidity distribution
727   if (!fHistEtaGenCascXiPlus) {
728      fHistEtaGenCascXiPlus = new TH1F("fHistEtaGenCascXiPlus", "#eta of any gen. #bar{#Xi}^{+};#eta;Number of Casc", 200, -10, 10);
729      fListHistCascade->Add(fHistEtaGenCascXiPlus);
730   }
731   
732   if (!f2dHistGenPtVsGenYGenXiPlus) {
733      f2dHistGenPtVsGenYGenXiPlus = new TH2F("f2dHistGenPtVsGenYGenXiPlus", "MC P_{t} Vs MC Y of Gen #bar{#Xi}^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
734      fListHistCascade->Add(f2dHistGenPtVsGenYGenXiPlus);
735   }
736   
737   
738   // - Info at the generation level of multi-strange particle
739   
740   if (!fHistThetaGenCascXiPlus) {
741      fHistThetaGenCascXiPlus = new TH1F("fHistThetaGenCascXiPlus", "#theta of gen. #bar{#Xi}^{+};#theta;Number of Casc.", 200, -10, 190);
742      fListHistCascade->Add(fHistThetaGenCascXiPlus);
743   }
744  
745   if (!f2dHistGenPtVsGenYFdblXiPlus) {
746      f2dHistGenPtVsGenYFdblXiPlus = new TH2F("f2dHistGenPtVsGenYFdblXiPlus", "MC P_{t} Vs MC Y of findable Gen #bar{#Xi}^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
747      fListHistCascade->Add(f2dHistGenPtVsGenYFdblXiPlus);
748   }
749   
750                 // - Theta distribution the daughters (control plots)
751   
752   if (!fHistThetaLambdaXiPlus) {
753      fHistThetaLambdaXiPlus = new TH1F("fHistThetaLambdaXiPlus", "#theta of gen. #Lambda (Xi dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
754      fListHistCascade->Add(fHistThetaLambdaXiPlus);
755   }
756
757   if (!fHistThetaBachXiPlus) {
758      fHistThetaBachXiPlus = new TH1F("fHistThetaBachXiPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
759      fListHistCascade->Add(fHistThetaBachXiPlus);
760   }
761   
762   if (!fHistThetaMesDghterXiPlus) {
763      fHistThetaMesDghterXiPlus = new TH1F("fHistThetaMesDghterXiPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
764      fListHistCascade->Add(fHistThetaMesDghterXiPlus);
765   }
766   
767   if (!fHistThetaBarDghterXiPlus) {
768      fHistThetaBarDghterXiPlus = new TH1F("fHistThetaBarDghterXiPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
769      fListHistCascade->Add(fHistThetaBarDghterXiPlus);
770   }
771  
772                 // - Pt distribution (control plots)
773     
774   if (!fHistPtBachXiPlus) {
775      fHistPtBachXiPlus = new TH1F("fHistPtBachXiPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
776      fListHistCascade->Add(fHistPtBachXiPlus);
777   }
778   
779   if (!fHistPtMesDghterXiPlus) {
780      fHistPtMesDghterXiPlus = new TH1F("fHistPtMesDghterXiPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
781      fListHistCascade->Add(fHistPtMesDghterXiPlus);
782   }
783     
784   if (!fHistPtBarDghterXiPlus) {
785      fHistPtBarDghterXiPlus = new TH1F("fHistPtBarDghterXiPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
786      fListHistCascade->Add(fHistPtBarDghterXiPlus);
787   }
788   
789   
790   //---------
791   // III - Omega- 
792   // - Pseudo-Rapidity distribution
793   if (!fHistEtaGenCascOmegaMinus) {
794      fHistEtaGenCascOmegaMinus = new TH1F("fHistEtaGenCascOmegaMinus", "#eta of any gen. #Omega^{-};#eta;Number of Casc", 200, -10, 10);
795      fListHistCascade->Add(fHistEtaGenCascOmegaMinus);
796   }
797   
798   if (!f2dHistGenPtVsGenYGenOmegaMinus) {
799      f2dHistGenPtVsGenYGenOmegaMinus = new TH2F("f2dHistGenPtVsGenYGenOmegaMinus", "MC P_{t} Vs MC Y of Gen #Omega^{-} ;Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
800      fListHistCascade->Add(f2dHistGenPtVsGenYGenOmegaMinus);
801   }
802   
803   
804   // - Info at the generation level of multi-strange particle
805   
806   if (!fHistThetaGenCascOmegaMinus) {
807      fHistThetaGenCascOmegaMinus = new TH1F("fHistThetaGenCascOmegaMinus", "#theta of gen. #Omega^{-};#theta;Number of Casc.", 200, -10, 190);
808      fListHistCascade->Add(fHistThetaGenCascOmegaMinus);
809   }
810  
811   if (!f2dHistGenPtVsGenYFdblOmegaMinus) {
812      f2dHistGenPtVsGenYFdblOmegaMinus = new TH2F("f2dHistGenPtVsGenYFdblOmegaMinus", "MC P_{t} Vs MC Y of findable Gen #Omega^{-}; Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
813      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaMinus);
814   }
815   
816                 // - Theta distribution the daughters (control plots)
817   
818   if (!fHistThetaLambdaOmegaMinus) {
819      fHistThetaLambdaOmegaMinus = new TH1F("fHistThetaLambdaOmegaMinus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
820      fListHistCascade->Add(fHistThetaLambdaOmegaMinus);
821   }
822
823   if (!fHistThetaBachOmegaMinus) {
824      fHistThetaBachOmegaMinus = new TH1F("fHistThetaBachOmegaMinus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
825      fListHistCascade->Add(fHistThetaBachOmegaMinus);
826   }
827   
828   if (!fHistThetaMesDghterOmegaMinus) {
829      fHistThetaMesDghterOmegaMinus = new TH1F("fHistThetaMesDghterOmegaMinus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
830      fListHistCascade->Add(fHistThetaMesDghterOmegaMinus);
831   }
832   
833   if (!fHistThetaBarDghterOmegaMinus) {
834      fHistThetaBarDghterOmegaMinus = new TH1F("fHistThetaBarDghterOmegaMinus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
835      fListHistCascade->Add(fHistThetaBarDghterOmegaMinus);
836   }
837  
838                 // - Pt distribution (control plots)
839     
840   if (!fHistPtBachOmegaMinus) {
841      fHistPtBachOmegaMinus = new TH1F("fHistPtBachOmegaMinus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
842      fListHistCascade->Add(fHistPtBachOmegaMinus);
843   }
844   
845   if (!fHistPtMesDghterOmegaMinus) {
846      fHistPtMesDghterOmegaMinus = new TH1F("fHistPtMesDghterOmegaMinus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
847      fListHistCascade->Add(fHistPtMesDghterOmegaMinus);
848   }
849     
850   if (!fHistPtBarDghterOmegaMinus) {
851      fHistPtBarDghterOmegaMinus = new TH1F("fHistPtBarDghterOmegaMinus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
852      fListHistCascade->Add(fHistPtBarDghterOmegaMinus);
853   }
854   
855   
856   //---------
857   // IV - Omega+ 
858   // - Pseudo-Rapidity distribution
859   if (!fHistEtaGenCascOmegaPlus) {
860      fHistEtaGenCascOmegaPlus = new TH1F("fHistEtaGenCascOmegaPlus", "#eta of any gen. #bar{#Omega}^{+};#eta;Number of Casc", 200, -10, 10);
861      fListHistCascade->Add(fHistEtaGenCascOmegaPlus);
862   }
863   
864   if (!f2dHistGenPtVsGenYGenOmegaPlus) {
865      f2dHistGenPtVsGenYGenOmegaPlus = new TH2F("f2dHistGenPtVsGenYGenOmegaPlus", "MC P_{t} Vs MC Y of Gen #bar{#Omega}^{+} ;Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
866      fListHistCascade->Add(f2dHistGenPtVsGenYGenOmegaPlus);
867   }
868   
869   
870   
871   // - Info at the generation level of multi-strange particle
872   
873   if (!fHistThetaGenCascOmegaPlus) {
874      fHistThetaGenCascOmegaPlus = new TH1F("fHistThetaGenCascOmegaPlus", "#theta of gen. #bar{#Omega}^{+};#theta;Number of Casc.", 200, -10, 190);
875      fListHistCascade->Add(fHistThetaGenCascOmegaPlus);
876   }
877  
878   if (!f2dHistGenPtVsGenYFdblOmegaPlus) {
879      f2dHistGenPtVsGenYFdblOmegaPlus = new TH2F("f2dHistGenPtVsGenYFdblOmegaPlus", "MC P_{t} Vs MC Y of findable Gen #bar{#Omega}^{+}; Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
880      fListHistCascade->Add(f2dHistGenPtVsGenYFdblOmegaPlus);
881   }
882
883   
884                 // - Theta distribution the daughters (control plots)
885   
886   if (!fHistThetaLambdaOmegaPlus) {
887      fHistThetaLambdaOmegaPlus = new TH1F("fHistThetaLambdaOmegaPlus", "#theta of gen. #Lambda (Omega dghter);#theta_{#Lambda};Number of #Lambda", 200, -10, 190);
888      fListHistCascade->Add(fHistThetaLambdaOmegaPlus);
889   }
890
891   if (!fHistThetaBachOmegaPlus) {
892      fHistThetaBachOmegaPlus = new TH1F("fHistThetaBachOmegaPlus", "#theta of gen. Bach.;#theta_{Bach};Number of Bach.", 200, -10, 190);
893      fListHistCascade->Add(fHistThetaBachOmegaPlus);
894   }
895   
896   if (!fHistThetaMesDghterOmegaPlus) {
897      fHistThetaMesDghterOmegaPlus = new TH1F("fHistThetaMesDghterOmegaPlus", "#theta of gen. Meson #Lambda dghter;#theta_{MesDght};Number of Mes.", 200, -10, 190);
898      fListHistCascade->Add(fHistThetaMesDghterOmegaPlus);
899   }
900   
901   if (!fHistThetaBarDghterOmegaPlus) {
902      fHistThetaBarDghterOmegaPlus = new TH1F("fHistThetaBarDghterOmegaPlus", "#theta of gen. Baryon #Lambda dghter;#theta_{BarDght};Number of Bar.", 200, -10, 190);
903      fListHistCascade->Add(fHistThetaBarDghterOmegaPlus);
904   }
905  
906                 // - Pt distribution (control plots)
907     
908   if (!fHistPtBachOmegaPlus) {
909      fHistPtBachOmegaPlus = new TH1F("fHistPtBachOmegaPlus", "p_{t} of gen. Bach.;pt_{Bach};Number of Bach.", 200, 0, 10);
910      fListHistCascade->Add(fHistPtBachOmegaPlus);
911   }
912   
913   if (!fHistPtMesDghterOmegaPlus) {
914      fHistPtMesDghterOmegaPlus = new TH1F("fHistPtMesDghterOmegaPlus", "p_{t} of gen. Meson #Lambda dghter);pt_{MesDght};Number of Mes.", 200, 0, 10);
915      fListHistCascade->Add(fHistPtMesDghterOmegaPlus);
916   }
917     
918   if (!fHistPtBarDghterOmegaPlus) {
919      fHistPtBarDghterOmegaPlus = new TH1F("fHistPtBarDghterOmegaPlus", "p_{t} of gen. Baryon #Lambda dghter);pt_{BarDght};Number of Bar.", 200, 0, 10);
920      fListHistCascade->Add(fHistPtBarDghterOmegaPlus);
921   }
922     
923   
924 //--------------------------------------------------------------------------------
925 // Part 2 - Any reconstructed cascades + reconstructed cascades associated with MC
926   
927                 // - Effective mass histos for cascades candidates.
928   
929   if (! fHistMassXiMinus) {
930           fHistMassXiMinus = new TH1F("fHistMassXiMinus","#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
931           fListHistCascade->Add(fHistMassXiMinus);
932   }
933   
934   if (! fHistMassXiPlus) {
935           fHistMassXiPlus = new TH1F("fHistMassXiPlus","#bar{#Xi}^{+} candidates;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
936           fListHistCascade->Add(fHistMassXiPlus);
937   }
938
939   if (! fHistMassOmegaMinus) {
940           fHistMassOmegaMinus = new TH1F("fHistMassOmegaMinus","#Omega^{-} candidates;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
941           fListHistCascade->Add(fHistMassOmegaMinus);
942   }
943  
944   if (! fHistMassOmegaPlus) {
945           fHistMassOmegaPlus = new TH1F("fHistMassOmegaPlus","#bar{#Omega}^{+} candidates;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
946           fListHistCascade->Add(fHistMassOmegaPlus);
947   }
948   
949   
950   
951                 // - Effective mass histos with combined PID
952   
953   if (! fHistMassWithCombPIDXiMinus) {
954     fHistMassWithCombPIDXiMinus = new TH1F("fHistMassWithCombPIDXiMinus","#Xi^{-} candidates, with Bach. comb. PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
955     fListHistCascade->Add(fHistMassWithCombPIDXiMinus);
956   }
957   
958   if (! fHistMassWithCombPIDXiPlus) {
959     fHistMassWithCombPIDXiPlus = new TH1F("fHistMassWithCombPIDXiPlus","#bar{#Xi}^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
960     fListHistCascade->Add(fHistMassWithCombPIDXiPlus);
961   }
962
963   if (! fHistMassWithCombPIDOmegaMinus) {
964         fHistMassWithCombPIDOmegaMinus = new TH1F("fHistMassWithCombPIDOmegaMinus","#Omega^{-} candidates, with Bach. comb. PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
965     fListHistCascade->Add(fHistMassWithCombPIDOmegaMinus);
966   }
967  
968   if (! fHistMassWithCombPIDOmegaPlus) {
969         fHistMassWithCombPIDOmegaPlus = new TH1F("fHistMassWithCombPIDOmegaPlus","#bar{#Omega}^{+} candidates, with Bach. comb. PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
970     fListHistCascade->Add(fHistMassWithCombPIDOmegaPlus);
971   }
972   
973                 // - PID Probability versus MC Pt(bachelor track)
974   if(! f2dHistPIDprobaKaonVsMCPtBach ){
975         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 );
976         fListHistCascade->Add(f2dHistPIDprobaKaonVsMCPtBach);
977   }
978   
979   if(! f2dHistPIDprobaPionVsMCPtBach ){
980         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 );
981         fListHistCascade->Add(f2dHistPIDprobaPionVsMCPtBach);
982   }
983   
984   
985                 // - Effective mass histos with perfect MC PID on the bachelor
986   
987   if (! fHistMassWithMcPIDXiMinus) {
988     fHistMassWithMcPIDXiMinus = new TH1F("fHistMassWithMcPIDXiMinus","#Xi^{-} candidates, with Bach. MC PID;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
989     fListHistCascade->Add(fHistMassWithMcPIDXiMinus);
990   }
991   
992   if (! fHistMassWithMcPIDXiPlus) {
993     fHistMassWithMcPIDXiPlus = new TH1F("fHistMassWithMcPIDXiPlus","#bar{#Xi}^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
994     fListHistCascade->Add(fHistMassWithMcPIDXiPlus);
995   }
996
997   if (! fHistMassWithMcPIDOmegaMinus) {
998         fHistMassWithMcPIDOmegaMinus = new TH1F("fHistMassWithMcPIDOmegaMinus","#Omega^{-} candidates, with Bach. MC PID;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
999     fListHistCascade->Add(fHistMassWithMcPIDOmegaMinus);
1000   }
1001  
1002   if (! fHistMassWithMcPIDOmegaPlus) {
1003         fHistMassWithMcPIDOmegaPlus = new TH1F("fHistMassWithMcPIDOmegaPlus","#bar{#Omega}^{+} candidates, with Bach. MC PID;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1004     fListHistCascade->Add(fHistMassWithMcPIDOmegaPlus);
1005   }
1006   
1007   
1008                 // - Effective mass histos for cascades candidates ASSOCIATED with MC.
1009   
1010   if (! fHistAsMCMassXiMinus) {
1011           fHistAsMCMassXiMinus = new TH1F("fHistAsMCMassXiMinus","#Xi^{-} candidates associated to MC;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts", 400,1.2,2.0);
1012           fListHistCascade->Add(fHistAsMCMassXiMinus);
1013   }
1014   
1015   if (! fHistAsMCMassXiPlus) {
1016           fHistAsMCMassXiPlus = new TH1F("fHistAsMCMassXiPlus","#bar{#Xi}^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , #pi^{+} ) (GeV/c^{2});Counts",400,1.2,2.0);
1017           fListHistCascade->Add(fHistAsMCMassXiPlus);
1018   }
1019
1020   if (! fHistAsMCMassOmegaMinus) {
1021           fHistAsMCMassOmegaMinus = new TH1F("fHistAsMCMassOmegaMinus","#Omega^{-} candidates associated to MC;M( #Lambda , K^{-} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1022           fListHistCascade->Add(fHistAsMCMassOmegaMinus);
1023   }
1024  
1025   if (! fHistAsMCMassOmegaPlus) {
1026           fHistAsMCMassOmegaPlus = new TH1F("fHistAsMCMassOmegaPlus","#bar{#Omega}^{+} candidates associated to MC;M( #bar{#Lambda}^{0} , K^{+} ) (GeV/c^{2});Counts", 500,1.5,2.5);
1027           fListHistCascade->Add(fHistAsMCMassOmegaPlus);
1028   }
1029   
1030                 
1031                 // -  Generated Pt Vs generated Y of the cascade candidates associated with MC 
1032                 //     + having the proper maximum proba of combined PID for the bachelor
1033   
1034   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiMinus) {
1035      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);
1036      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiMinus);
1037   }
1038   
1039   if (!f2dHistAsMCandCombPIDGenPtVsGenYXiPlus) {
1040      f2dHistAsMCandCombPIDGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of #bar{#Xi}^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1041      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYXiPlus);
1042   }
1043   
1044   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus) {
1045      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);
1046      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus);
1047   }
1048   
1049   if (!f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus) {
1050      f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of #bar{#Omega}^{+} (associated+Bach.PID); Pt_{MC} (GeV/c); Y_{MC}", 200, 0., 10., 220, -1.1, 1.1);
1051      fListHistCascade->Add(f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus);
1052   }
1053   
1054
1055                 // - Generated Pt Vs Generated Y, for the cascade candidates associated with MC
1056   
1057   if (!f2dHistAsMCGenPtVsGenYXiMinus) {
1058         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);
1059           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiMinus );
1060   }
1061   
1062   if (!f2dHistAsMCGenPtVsGenYXiPlus) {
1063           f2dHistAsMCGenPtVsGenYXiPlus = new TH2F("f2dHistAsMCGenPtVsGenYXiPlus", "MC P_{t} Vs MC Y of gen. #bar{#Xi}^{+} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1064           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYXiPlus );
1065   }
1066   
1067   if (!f2dHistAsMCGenPtVsGenYOmegaMinus) {
1068           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);
1069           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaMinus );
1070   }
1071   
1072   if (!f2dHistAsMCGenPtVsGenYOmegaPlus) {
1073           f2dHistAsMCGenPtVsGenYOmegaPlus = new TH2F("f2dHistAsMCGenPtVsGenYOmegaPlus", "MC P_{t} Vs MC Y of gen. #bar{#Omega}^{+} (associated);Pt_{MC} (GeV/c); Rapidity, Y_{MC}",200, 0., 10., 220, -1.1, 1.1);
1074           fListHistCascade->Add(f2dHistAsMCGenPtVsGenYOmegaPlus );
1075   } 
1076   
1077   
1078                 // - Generated Eta of the the cascade candidates associated with MC
1079   if (!fHistAsMCGenEtaXiMinus) {
1080           fHistAsMCGenEtaXiMinus = new TH1F("fHistAsMCGenEtaXiMinus", "#eta of gen. #Xi^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1081           fListHistCascade->Add( fHistAsMCGenEtaXiMinus );
1082   }
1083   
1084   if (!fHistAsMCGenEtaXiPlus) {
1085           fHistAsMCGenEtaXiPlus = new TH1F("fHistAsMCGenEtaXiPlus", "#eta of gen. #bar{#Xi}^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1086           fListHistCascade->Add( fHistAsMCGenEtaXiPlus );
1087   }
1088   
1089   if (!fHistAsMCGenEtaOmegaMinus) {
1090           fHistAsMCGenEtaOmegaMinus = new TH1F("fHistAsMCGenEtaOmegaMinus", "#eta of gen. #Omega^{-} (associated);#eta;Number of Casc", 100, -5, 5);
1091           fListHistCascade->Add( fHistAsMCGenEtaOmegaMinus );
1092   }
1093   
1094   if (!fHistAsMCGenEtaOmegaPlus) {
1095           fHistAsMCGenEtaOmegaPlus = new TH1F("fHistAsMCGenEtaOmegaPlus", "#eta of gen. #bar{#Omega}^{+} (associated);#eta;Number of Casc", 100, -5, 5);
1096           fListHistCascade->Add( fHistAsMCGenEtaOmegaPlus );
1097   }
1098   
1099   
1100   
1101                 // - Resolution in Pt as function of generated Pt
1102   
1103   if(! f2dHistAsMCResPtXiMinus) {
1104           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);
1105           fListHistCascade->Add(f2dHistAsMCResPtXiMinus);
1106   }
1107   
1108   if(! f2dHistAsMCResPtXiPlus) {
1109           f2dHistAsMCResPtXiPlus = new TH2F( "f2dHistAsMCResPtXiPlus", "Resolution in Pt reconstruction for #bar{#Xi}^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1110           fListHistCascade->Add(f2dHistAsMCResPtXiPlus);
1111   }
1112   
1113   if(! f2dHistAsMCResPtOmegaMinus) {
1114           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);
1115           fListHistCascade->Add(f2dHistAsMCResPtOmegaMinus);
1116   }
1117   
1118   if(! f2dHistAsMCResPtOmegaPlus) {
1119           f2dHistAsMCResPtOmegaPlus = new TH2F( "f2dHistAsMCResPtOmegaPlus", "Resolution in Pt reconstruction for #bar{#Omega}^{+}; Pt_{MC} (GeV/c); (Pt_{reco} - Pt_{MC}) / Pt_{MC}", 200, 0., 10., 200, -0.1, 0.1);
1120           fListHistCascade->Add(f2dHistAsMCResPtOmegaPlus);
1121   }
1122   
1123                 // - Resolution in R(2D) as function of generated R
1124   
1125   if(! f2dHistAsMCResRXiMinus) {
1126           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);
1127           fListHistCascade->Add(f2dHistAsMCResRXiMinus);
1128   }
1129   
1130   if(! f2dHistAsMCResRXiPlus) {
1131           f2dHistAsMCResRXiPlus = new TH2F( "f2dHistAsMCResRXiPlus", "Resolution in transv. position for #bar{#Xi}^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1132           fListHistCascade->Add(f2dHistAsMCResRXiPlus);
1133   }
1134   
1135   if(! f2dHistAsMCResROmegaMinus) {
1136           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);
1137           fListHistCascade->Add(f2dHistAsMCResROmegaMinus);
1138   }
1139   
1140   if(! f2dHistAsMCResROmegaPlus) {
1141           f2dHistAsMCResROmegaPlus = new TH2F( "f2dHistAsMCResROmegaPlus", "Resolution in transv. position for #bar{#Omega}^{+}; R_{MC} (cm); (R_{reco} - R_{MC}) / R_{MC}", 450, 0., 45.0, 240, -0.3, 0.3);
1142           fListHistCascade->Add(f2dHistAsMCResROmegaPlus);
1143   }
1144   
1145                 // - Resolution in phi as function of generated Pt
1146     
1147   if(! f2dHistAsMCResPhiXiMinus) {
1148           f2dHistAsMCResPhiXiMinus = new TH2F( "f2dHistAsMCResPhiXiMinus", "Resolution in #phi for #Xi^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1149           fListHistCascade->Add(f2dHistAsMCResPhiXiMinus);
1150   }
1151   
1152   if(! f2dHistAsMCResPhiXiPlus) {
1153           f2dHistAsMCResPhiXiPlus = new TH2F( "f2dHistAsMCResPhiXiPlus", "Resolution in #phi for #bar{#Xi}^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1154           fListHistCascade->Add(f2dHistAsMCResPhiXiPlus);
1155   }
1156   
1157   if(! f2dHistAsMCResPhiOmegaMinus) {
1158           f2dHistAsMCResPhiOmegaMinus = new TH2F( "f2dHistAsMCResPhiOmegaMinus", "Resolution in #phi for #Omega^{-}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);  
1159           fListHistCascade->Add(f2dHistAsMCResPhiOmegaMinus);
1160   }
1161   
1162   if(! f2dHistAsMCResPhiOmegaPlus) {
1163           f2dHistAsMCResPhiOmegaPlus = new TH2F( "f2dHistAsMCResPhiOmegaPlus", "Resolution in #phi for #bar{#Omega}^{+}; Pt_{MC} (GeV/c); #phi(MC) - #phi(reco)   (deg)", 200, 0., 10., 60, -30., 30.);
1164           fListHistCascade->Add(f2dHistAsMCResPhiOmegaPlus);
1165   }
1166   
1167   
1168                 // - Correlation in Pt between the cascade and its (anti)proton daughter
1169   if(! f2dHistAsMCPtProtonVsPtXiMinus) {
1170         f2dHistAsMCPtProtonVsPtXiMinus = new TH2F( "f2dHistAsMCPtProtonVsPtXiMinus", "Correlation Pt(p) Vs Pt(#Xi^{-}), associated to MC; Pt_{MC}(p) (GeV/c); Pt_{MC}(#Xi^{-}) (GeV/c)", 100, 0., 10., 100, 0., 10.);
1171         fListHistCascade->Add(f2dHistAsMCPtProtonVsPtXiMinus);
1172   }
1173   
1174   if(! f2dHistAsMCPtAntiProtonVsPtXiPlus) {
1175         f2dHistAsMCPtAntiProtonVsPtXiPlus = new TH2F( "f2dHistAsMCPtAntiProtonVsPtXiPlus", "Correlation Pt(#bar{p}) Vs Pt(#bar{#Xi}^{+}), associated to MC; Pt_{MC}(#bar{p}) (GeV/c); Pt_{MC}(#bar{#Xi}^{+}) (GeV/c)", 100, 0., 10., 100, 0., 10.);
1176         fListHistCascade->Add(f2dHistAsMCPtAntiProtonVsPtXiPlus);
1177   }
1178   
1179   if(! f2dHistAsMCPtProtonVsPtOmegaMinus) {
1180         f2dHistAsMCPtProtonVsPtOmegaMinus = new TH2F( "f2dHistAsMCPtProtonVsPtOmegaMinus", "Correlation Pt(p) Vs Pt(#Omega^{-}), associated to MC; Pt_{MC}(p) (GeV/c); Pt_{MC}(#Omega^{-}) (GeV/c)", 100, 0., 10., 100, 0., 10.);
1181         fListHistCascade->Add(f2dHistAsMCPtProtonVsPtOmegaMinus);
1182   }
1183   
1184   
1185   if(! f2dHistAsMCPtAntiProtonVsPtOmegaPlus) {
1186         f2dHistAsMCPtAntiProtonVsPtOmegaPlus = new TH2F( "f2dHistAsMCPtAntiProtonVsPtOmegaPlus", "Correlation Pt(#bar{p}) Vs Pt(#bar{#Omega}^{+}), associated to MC; Pt_{MC}(#bar{p}) (GeV/c); Pt_{MC}(#bar{#Omega}^{+}) (GeV/c)", 100, 0., 10., 100, 0., 10.);
1187         fListHistCascade->Add(f2dHistAsMCPtAntiProtonVsPtOmegaPlus);
1188   }
1189   
1190   
1191   
1192                 // - PID container
1193 if(! fCFContCascadePIDAsXiMinus)  {
1194   const Int_t  lNbSteps      =  7 ;
1195   const Int_t  lNbVariables  =  4 ;
1196
1197   //array for the number of bins in each dimension :
1198   Int_t lNbBinsPerVar[4] = {0};
1199   lNbBinsPerVar[0] = 200;
1200   lNbBinsPerVar[1] = 75;
1201   lNbBinsPerVar[2] = 44;
1202   lNbBinsPerVar[3] = 250;
1203    
1204   
1205   fCFContCascadePIDAsXiMinus = new AliCFContainer("fCFContCascadePIDAsXiMinus","Pt_{cascade} Vs M_{#Xi^{-} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1206   
1207   //setting the bin limits (valid  for v4-18-10-AN)
1208   fCFContCascadePIDAsXiMinus->SetBinLimits(0,   0.0  ,  10.0  );        // Pt(Cascade)
1209   fCFContCascadePIDAsXiMinus->SetBinLimits(1,   1.25 ,   1.40 );        // Xi Effective mass
1210   fCFContCascadePIDAsXiMinus->SetBinLimits(2,  -1.1  ,   1.1  );        // Rapidity
1211   if(fCollidingSystems) 
1212         fCFContCascadePIDAsXiMinus->SetBinLimits(3, 0.0, 20000.0  );    // nTrackPrimaryMultiplicity
1213   else
1214         fCFContCascadePIDAsXiMinus->SetBinLimits(3, 0.0, 250.0  );     // nTrackPrimaryMultiplicity
1215   
1216   // Setting the step title : one per PID case
1217   fCFContCascadePIDAsXiMinus->SetStepTitle(0, "No PID");
1218   fCFContCascadePIDAsXiMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1219   fCFContCascadePIDAsXiMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1220   fCFContCascadePIDAsXiMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1221   fCFContCascadePIDAsXiMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1222   fCFContCascadePIDAsXiMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1223   fCFContCascadePIDAsXiMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1224   
1225   // Setting the variable title, per axis
1226   fCFContCascadePIDAsXiMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1227   fCFContCascadePIDAsXiMinus->SetVarTitle(1, "M( #Lambda , #pi^{-} ) (GeV/c^{2})");
1228   fCFContCascadePIDAsXiMinus->SetVarTitle(2, "Y_{#Xi}");
1229   fCFContCascadePIDAsXiMinus->SetVarTitle(3, "Primary Track Multiplicity");
1230   
1231   fListHistCascade->Add(fCFContCascadePIDAsXiMinus);
1232   
1233 }
1234
1235 if(! fCFContCascadePIDAsXiPlus)  {
1236   const Int_t  lNbSteps      =  7 ;
1237   const Int_t  lNbVariables  =  4 ;
1238
1239   //array for the number of bins in each dimension :
1240   Int_t lNbBinsPerVar[4] = {0};
1241   lNbBinsPerVar[0] = 200;
1242   lNbBinsPerVar[1] = 75;
1243   lNbBinsPerVar[2] = 44;
1244   lNbBinsPerVar[3] = 250;
1245    
1246   
1247   fCFContCascadePIDAsXiPlus = new AliCFContainer("fCFContCascadePIDAsXiPlus","Pt_{cascade} Vs M_{#bar{#Xi}^{+} candidates} Vs Y_{#Xi}", lNbSteps, lNbVariables, lNbBinsPerVar );
1248   
1249   
1250   //setting the bin limits (valid  for v4-18-10-AN)
1251   fCFContCascadePIDAsXiPlus->SetBinLimits(0,   0.0  ,  10.0  );         // Pt(Cascade)
1252   fCFContCascadePIDAsXiPlus->SetBinLimits(1,   1.25 ,   1.40 );         // Xi Effective mass
1253   fCFContCascadePIDAsXiPlus->SetBinLimits(2,  -1.1  ,   1.1  );         // Rapidity
1254   if(fCollidingSystems) 
1255         fCFContCascadePIDAsXiPlus->SetBinLimits(3, 0.0, 20000.0  );    // nTrackPrimaryMultiplicity
1256   else
1257         fCFContCascadePIDAsXiPlus->SetBinLimits(3, 0.0, 250.0  );     // nTrackPrimaryMultiplicity
1258   
1259   // Setting the step title : one per PID case
1260   fCFContCascadePIDAsXiPlus->SetStepTitle(0, "No PID");
1261   fCFContCascadePIDAsXiPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1262   fCFContCascadePIDAsXiPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1263   fCFContCascadePIDAsXiPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1264   fCFContCascadePIDAsXiPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1265   fCFContCascadePIDAsXiPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1266   fCFContCascadePIDAsXiPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1267   
1268   // Setting the variable title, per axis
1269   fCFContCascadePIDAsXiPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1270   fCFContCascadePIDAsXiPlus->SetVarTitle(1, "M( #Lambda , #pi^{+} ) (GeV/c^{2})");
1271   fCFContCascadePIDAsXiPlus->SetVarTitle(2, "Y_{#Xi}");
1272   fCFContCascadePIDAsXiPlus->SetVarTitle(3, "Primary Track Multiplicity");
1273   
1274   fListHistCascade->Add(fCFContCascadePIDAsXiPlus);
1275   
1276 }
1277
1278
1279 if(! fCFContCascadePIDAsOmegaMinus)  {
1280   const Int_t  lNbSteps      =  7 ;
1281   const Int_t  lNbVariables  =  4 ;
1282
1283   //array for the number of bins in each dimension :
1284   Int_t lNbBinsPerVar[4] = {0};
1285   lNbBinsPerVar[0] = 200;
1286   lNbBinsPerVar[1] = 60;
1287   lNbBinsPerVar[2] = 44;
1288   lNbBinsPerVar[3] = 250;
1289    
1290   
1291   fCFContCascadePIDAsOmegaMinus = new AliCFContainer("fCFContCascadePIDAsOmegaMinus","Pt_{cascade} Vs M_{#Omega^{-} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1292   
1293   
1294   //setting the bin limits (valid  for v4-18-10-AN)
1295   fCFContCascadePIDAsOmegaMinus->SetBinLimits(0,   0.0  ,  10.0  );     // Pt(Cascade)
1296   fCFContCascadePIDAsOmegaMinus->SetBinLimits(1,   1.62 ,   1.74 );     // Omega Effective mass
1297   fCFContCascadePIDAsOmegaMinus->SetBinLimits(2,  -1.1  ,   1.1  );     // Rapidity
1298   if(fCollidingSystems) 
1299         fCFContCascadePIDAsOmegaMinus->SetBinLimits(3, 0.0, 20000.0  );    // nTrackPrimaryMultiplicity
1300   else
1301         fCFContCascadePIDAsOmegaMinus->SetBinLimits(3, 0.0, 250.0  );     // nTrackPrimaryMultiplicity
1302   
1303   // Setting the step title : one per PID case
1304   fCFContCascadePIDAsOmegaMinus->SetStepTitle(0, "No PID");
1305   fCFContCascadePIDAsOmegaMinus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1306   fCFContCascadePIDAsOmegaMinus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1307   fCFContCascadePIDAsOmegaMinus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1308   fCFContCascadePIDAsOmegaMinus->SetStepTitle(4, "Comb. PID / Bachelor");
1309   fCFContCascadePIDAsOmegaMinus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1310   fCFContCascadePIDAsOmegaMinus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1311   
1312   // Setting the variable title, per axis
1313   fCFContCascadePIDAsOmegaMinus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1314   fCFContCascadePIDAsOmegaMinus->SetVarTitle(1, "M( #Lambda , K^{-} ) (GeV/c^{2})");
1315   fCFContCascadePIDAsOmegaMinus->SetVarTitle(2, "Y_{#Omega}");
1316   fCFContCascadePIDAsOmegaMinus->SetVarTitle(3, "Primary Track Multiplicity");
1317   
1318   fListHistCascade->Add(fCFContCascadePIDAsOmegaMinus);
1319   
1320 }
1321
1322 if(! fCFContCascadePIDAsOmegaPlus)  {
1323   const Int_t  lNbSteps      =  7 ;
1324   const Int_t  lNbVariables  =  4 ;
1325
1326   //array for the number of bins in each dimension :
1327   Int_t lNbBinsPerVar[4]= {0};
1328   lNbBinsPerVar[0] = 200;
1329   lNbBinsPerVar[1] = 60;
1330   lNbBinsPerVar[2] = 44;
1331   lNbBinsPerVar[3] = 250;
1332    
1333   
1334   fCFContCascadePIDAsOmegaPlus = new AliCFContainer("fCFContCascadePIDAsOmegaPlus","Pt_{cascade} Vs M_{#bar{#Omega}^{+} candidates} Vs Y_{#Omega}", lNbSteps, lNbVariables, lNbBinsPerVar );
1335   
1336   
1337   //setting the bin limits (valid  for v4-18-10-AN)
1338   fCFContCascadePIDAsOmegaPlus->SetBinLimits(0,   0.0  ,  10.0  );      // Pt(Cascade)
1339   fCFContCascadePIDAsOmegaPlus->SetBinLimits(1,   1.62 ,   1.74 );      // Omega Effective mass
1340   fCFContCascadePIDAsOmegaPlus->SetBinLimits(2,  -1.1  ,   1.1  );      // Rapidity
1341   if(fCollidingSystems) 
1342         fCFContCascadePIDAsOmegaPlus->SetBinLimits(3, 0.0, 20000.0  );    // nTrackPrimaryMultiplicity
1343   else
1344         fCFContCascadePIDAsOmegaPlus->SetBinLimits(3, 0.0, 250.0  );     // nTrackPrimaryMultiplicity
1345   
1346   // Setting the step title : one per PID case
1347   fCFContCascadePIDAsOmegaPlus->SetStepTitle(0, "No PID");
1348   fCFContCascadePIDAsOmegaPlus->SetStepTitle(1, "TPC PID / 4-#sigma cut on Bachelor track");
1349   fCFContCascadePIDAsOmegaPlus->SetStepTitle(2, "TPC PID / 4-#sigma cut on Bachelor+Baryon tracks");
1350   fCFContCascadePIDAsOmegaPlus->SetStepTitle(3, "TPC PID / 4-#sigma cut on Bachelor+Baryon+Meson tracks");
1351   fCFContCascadePIDAsOmegaPlus->SetStepTitle(4, "Comb. PID / Bachelor");
1352   fCFContCascadePIDAsOmegaPlus->SetStepTitle(5, "Comb. PID / Bachelor+Baryon");
1353   fCFContCascadePIDAsOmegaPlus->SetStepTitle(6, "Comb. PID / Bachelor+Baryon+Meson");
1354   
1355   // Setting the variable title, per axis
1356   fCFContCascadePIDAsOmegaPlus->SetVarTitle(0, "Pt_{cascade} (GeV/c)");
1357   fCFContCascadePIDAsOmegaPlus->SetVarTitle(1, "M( #Lambda , K^{+} ) (GeV/c^{2})");
1358   fCFContCascadePIDAsOmegaPlus->SetVarTitle(2, "Y_{#Omega}");
1359   fCFContCascadePIDAsOmegaPlus->SetVarTitle(3, "Primary Track Multiplicity");
1360   
1361   fListHistCascade->Add(fCFContCascadePIDAsOmegaPlus);
1362   
1363 }
1364
1365 // Part 3 : Towards the optimisation of topological selections -------
1366 if(! fCFContAsCascadeCuts){
1367    
1368         // Container meant to store all the relevant distributions corresponding to the cut variables.
1369         // So far, 20 variables have been identified.
1370         // The following will be done in quite a brut force way ... 
1371         // FIXME Improvement expected later (before Pb-Pb data at least)
1372         //          - Define a user binning to have less bins in each dimension
1373         //          - boolean for enabling/disbaling this CFContainer
1374   const Int_t  lNbSteps      =  4 ;
1375   const Int_t  lNbVariables  =  20 ;
1376   
1377   //array for the number of bins in each dimension :
1378   Int_t lNbBinsPerVar[20] = {0};
1379   lNbBinsPerVar[0]  = 25;
1380   lNbBinsPerVar[1]  = 25;
1381   lNbBinsPerVar[2]  = 20;
1382   lNbBinsPerVar[3]  = 40;
1383   lNbBinsPerVar[4]  = 30;
1384   lNbBinsPerVar[5]  = 25;
1385   
1386   lNbBinsPerVar[6]  = 20;
1387   lNbBinsPerVar[7]  = 40;
1388   lNbBinsPerVar[8]  = 40;
1389   lNbBinsPerVar[9]  = 25;
1390   lNbBinsPerVar[10] = 25;
1391   
1392   lNbBinsPerVar[11] = 75; // 2-MeV/c2 bins
1393   lNbBinsPerVar[12] = 60; // 2-MeV/c2 bins
1394   
1395   lNbBinsPerVar[13] = 100;
1396   
1397   lNbBinsPerVar[14] = 44; // 0.05 in rapidity units
1398   lNbBinsPerVar[15] = 44; // 0.05 in rapidity units
1399   
1400   lNbBinsPerVar[16] = 20;
1401  
1402   lNbBinsPerVar[17] = 50;
1403   lNbBinsPerVar[18] = 100;
1404   lNbBinsPerVar[19] = 24;
1405    
1406    
1407   fCFContAsCascadeCuts = new AliCFContainer("fCFContAsCascadeCuts","Cut Container for Asso. Cascades", lNbSteps, lNbVariables, lNbBinsPerVar );
1408   
1409   //0
1410   Double_t *lBinLim0  = new Double_t[ lNbBinsPerVar[0]+1 ];
1411         for(Int_t i=0; i< lNbBinsPerVar[0];i++)  lBinLim0[i]  = (Double_t)0.0   + (4.8  - 0.0 )/(lNbBinsPerVar[0]-1)  * (Double_t)i ;
1412         lBinLim0[ lNbBinsPerVar[0]  ] = 20.0;
1413   fCFContAsCascadeCuts -> SetBinLimits(0,  lBinLim0 );            // DcaXiDaughters : 0.0 to 5.0  
1414   delete [] lBinLim0;  
1415   //1
1416   Double_t *lBinLim1  = new Double_t[ lNbBinsPerVar[1]+1 ];
1417         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 ;
1418         lBinLim1[ lNbBinsPerVar[1]  ] = 100.0;
1419   fCFContAsCascadeCuts -> SetBinLimits(1,  lBinLim1 );            // DcaBachToPrimVertexXi : 0.0 to 0.25  
1420   delete [] lBinLim1;  
1421   //2 
1422   Double_t *lBinLim2  = new Double_t[ lNbBinsPerVar[2]+1 ];
1423         for(Int_t i=1; i< lNbBinsPerVar[2]+1;i++)   lBinLim2[i]  = (Double_t)0.81   + (1.0  - 0.81 )/(lNbBinsPerVar[2]-1)  * (Double_t) (i-1) ;   
1424         lBinLim2[0] = 0.0;
1425   fCFContAsCascadeCuts -> SetBinLimits(2,  lBinLim2 );            // XiCosineOfPointingAngle : 0.80 to 1.0        
1426   delete [] lBinLim2;  
1427   //3
1428   Double_t *lBinLim3  = new Double_t[ lNbBinsPerVar[3]+1 ];
1429         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 ;
1430         lBinLim3[ lNbBinsPerVar[3]  ] = 110.0;
1431   fCFContAsCascadeCuts -> SetBinLimits(3,  lBinLim3 );            // XiRadius : 0.0 to 4.0        
1432   delete [] lBinLim3;  
1433   //4
1434   fCFContAsCascadeCuts->SetBinLimits(4,    1.1  ,  1.13 );        // InvMassLambdaAsCascDghter
1435   //5
1436   Double_t *lBinLim5  = new Double_t[ lNbBinsPerVar[5]+1 ];
1437         for(Int_t i=0; i< lNbBinsPerVar[5];i++)   lBinLim5[i]  = (Double_t)0.0   + (4.8  - 0.0 )/(lNbBinsPerVar[5]-1)  * (Double_t)i ;
1438         lBinLim5[ lNbBinsPerVar[5]  ] = 20.0;
1439   fCFContAsCascadeCuts -> SetBinLimits(5,  lBinLim5 );            // DcaV0DaughtersXi : 0.0 to 5.0        
1440   delete [] lBinLim5;  
1441   
1442   
1443   //6
1444   Double_t *lBinLim6  = new Double_t[ lNbBinsPerVar[6]+1 ];
1445         for(Int_t i=1; i< lNbBinsPerVar[6]+1 ;i++)   lBinLim6[i]  = (Double_t)0.81   + (1.0  - 0.81 )/(lNbBinsPerVar[6]-1)  * (Double_t) (i-1) ;   
1446         lBinLim6[0] = 0.0;
1447   fCFContAsCascadeCuts -> SetBinLimits(6,  lBinLim6 );            // V0CosineOfPointingAngleXi : 0.80 to 1.0      
1448   delete [] lBinLim6;  
1449   //7
1450   Double_t *lBinLim7  = new Double_t[ lNbBinsPerVar[7]+1 ];
1451         for(Int_t i=0; i< lNbBinsPerVar[7];i++)   lBinLim7[i]  = (Double_t)0.0   + (7.8  - 0.0 )/(lNbBinsPerVar[7]-1)  * (Double_t)i ;
1452         lBinLim7[ lNbBinsPerVar[7]  ] = 100.0;
1453   fCFContAsCascadeCuts -> SetBinLimits(7,  lBinLim7 );            // V0RadiusXi : 0.0 to 8.0      
1454   delete [] lBinLim7;  
1455   //8
1456   Double_t *lBinLim8  = new Double_t[ lNbBinsPerVar[8]+1 ];
1457         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 ;
1458         lBinLim8[ lNbBinsPerVar[8]  ] = 100.0;
1459   fCFContAsCascadeCuts -> SetBinLimits(8,  lBinLim8 );            // DcaV0ToPrimVertexXi : 0.0 to 0.4     
1460   delete [] lBinLim8;  
1461   //9
1462   Double_t *lBinLim9  = new Double_t[ lNbBinsPerVar[9]+1 ];
1463         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 ;
1464         lBinLim9[ lNbBinsPerVar[9]  ] = 100.0;
1465   fCFContAsCascadeCuts -> SetBinLimits(9,  lBinLim9 );            // DcaPosToPrimVertexXi : 0.0 to 0.25   
1466   delete [] lBinLim9; 
1467   //10
1468   Double_t *lBinLim10  = new Double_t[ lNbBinsPerVar[10]+1 ];
1469         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 ;
1470         lBinLim10[ lNbBinsPerVar[10]  ] = 100.0;
1471   fCFContAsCascadeCuts -> SetBinLimits(10,  lBinLim10 );            // DcaPosToPrimVertexXi : 0.0 to 0.25 
1472   delete [] lBinLim10; 
1473   
1474   
1475   //11
1476   fCFContAsCascadeCuts->SetBinLimits(11,   1.25 ,  1.40 );      // InvMassXi
1477   fCFContAsCascadeCuts->SetBinLimits(12,   1.62  , 1.74 );      // InvMassOmega
1478   fCFContAsCascadeCuts->SetBinLimits(13,   0.0  , 10.0  );      // XiTransvMom 
1479   fCFContAsCascadeCuts->SetBinLimits(14,  -1.1  ,  1.1  );      // Y(Xi)
1480   fCFContAsCascadeCuts->SetBinLimits(15,  -1.1  ,  1.1  );      // Y(Omega)
1481   fCFContAsCascadeCuts->SetBinLimits(16, -10.0  , 10.0  );      // BestPrimaryVtxPosZ
1482   if(fCollidingSystems){
1483         fCFContAsCascadeCuts->SetBinLimits(17,   0.0, 10000.0  );    // nTrackPrimaryMultiplicity
1484         fCFContAsCascadeCuts->SetBinLimits(18,   0.0, 10000.0  );    // nITSandTPCtracksAndSPDtracklets
1485   }
1486   else{  
1487         //fCFContAsCascadeCuts->SetBinLimits(17,   0.0, 250.0  );     // nTrackPrimaryMultiplicity
1488         Double_t *lBinLim17  = new Double_t[ lNbBinsPerVar[17]+1 ];
1489                 lBinLim17[0] = 0;       lBinLim17[10] = 10;      lBinLim17[20] = 24;     lBinLim17[30] = 45;    lBinLim17[40] = 95;    lBinLim17[50] = 250;
1490                 lBinLim17[1] = 1;       lBinLim17[11] = 11;      lBinLim17[21] = 25;     lBinLim17[31] = 50;    lBinLim17[41] = 100;
1491                 lBinLim17[2] = 2;       lBinLim17[12] = 13;      lBinLim17[22] = 27;     lBinLim17[32] = 55;    lBinLim17[42] = 105;
1492                 lBinLim17[3] = 3;       lBinLim17[13] = 14;      lBinLim17[23] = 30;     lBinLim17[33] = 60;    lBinLim17[43] = 110;
1493                 lBinLim17[4] = 4;       lBinLim17[14] = 15;      lBinLim17[24] = 31;     lBinLim17[34] = 65;    lBinLim17[44] = 115;
1494                 lBinLim17[5] = 5;       lBinLim17[15] = 16;      lBinLim17[25] = 32;     lBinLim17[35] = 70;    lBinLim17[45] = 120;
1495                 lBinLim17[6] = 6;       lBinLim17[16] = 20;      lBinLim17[26] = 33;     lBinLim17[36] = 75;    lBinLim17[46] = 125;
1496                 lBinLim17[7] = 7;       lBinLim17[17] = 21;      lBinLim17[27] = 34;     lBinLim17[37] = 80;    lBinLim17[47] = 130;
1497                 lBinLim17[8] = 8;       lBinLim17[18] = 22;      lBinLim17[28] = 35;     lBinLim17[38] = 85;    lBinLim17[48] = 135;
1498                 lBinLim17[9] = 9;       lBinLim17[19] = 23;      lBinLim17[29] = 40;     lBinLim17[39] = 90;    lBinLim17[49] = 140;
1499                 
1500                 fCFContAsCascadeCuts -> SetBinLimits(17,  lBinLim17 );            // nTrackPrimaryMultiplicity : 0 to 250       
1501         delete [] lBinLim17;     
1502           
1503
1504         fCFContAsCascadeCuts->SetBinLimits(18,   0.0, 200.0  );     // nITSandTPCtracksAndSPDtracklets
1505   }
1506   fCFContAsCascadeCuts->SetBinLimits(19,  68.0  ,164.0  );      // BachTPCClusters
1507   
1508   
1509   // Regular binning definition (valid for v4-18-10-AN on)
1510   /*
1511   //setting the bin limits
1512   fCFContAsCascadeCuts->SetBinLimits(0,    0.0  ,  2.5  );      // DcaXiDaughters
1513   fCFContAsCascadeCuts->SetBinLimits(1,    0.0  ,  0.25 );      // DcaBachToPrimVertexXi
1514   fCFContAsCascadeCuts->SetBinLimits(2,    0.99 ,  1.0  );      // XiCosineOfPointingAngle
1515   fCFContAsCascadeCuts->SetBinLimits(3,    0.0  ,  4.0  );      // XiRadius
1516   fCFContAsCascadeCuts->SetBinLimits(4,    1.1  ,  1.15 );      // InvMassLambdaAsCascDghter
1517   fCFContAsCascadeCuts->SetBinLimits(5,    0.0  ,  1.0  );      // DcaV0DaughtersXi
1518   fCFContAsCascadeCuts->SetBinLimits(6,    0.98 ,  1.0  );      // V0CosineOfPointingAngleXi
1519   fCFContAsCascadeCuts->SetBinLimits(7,    0.0  , 20.0  );      // V0RadiusXi
1520   fCFContAsCascadeCuts->SetBinLimits(8,    0.0  ,  1.0  );      // DcaV0ToPrimVertexXi
1521   fCFContAsCascadeCuts->SetBinLimits(9,    0.0  ,  0.25 );      // DcaPosToPrimVertexXi
1522   fCFContAsCascadeCuts->SetBinLimits(10,   0.0  ,  0.25 );      // DcaNegToPrimVertexXi
1523   fCFContAsCascadeCuts->SetBinLimits(11,   1.25 ,  1.40 );      // InvMassXi
1524   fCFContAsCascadeCuts->SetBinLimits(12,   1.62 ,  1.74 );      // InvMassOmega
1525   fCFContAsCascadeCuts->SetBinLimits(13,   0.0  , 10.0  );      // pt_MC(Xi)
1526   fCFContAsCascadeCuts->SetBinLimits(14,  -1.1  ,  1.1  );      // Y_MC(Xi)
1527   fCFContAsCascadeCuts->SetBinLimits(15,  -1.1  ,  1.1  );      // Y_MC(Omega)
1528   fCFContAsCascadeCuts->SetBinLimits(16, -10.0  , 10.0  );      // BestPrimaryVtxPosZ
1529   if(fCollidingSystems){
1530         fCFContAsCascadeCuts->SetBinLimits(17,   0.0, 10000.0  );    // nTrackPrimaryMultiplicity
1531         fCFContAsCascadeCuts->SetBinLimits(18,   0.0, 10000.0  );    // nITSandTPCtracksAndSPDtracklets
1532   }
1533   else{  
1534         fCFContAsCascadeCuts->SetBinLimits(17,   0.0, 250.0  );     // nTrackPrimaryMultiplicity
1535         fCFContAsCascadeCuts->SetBinLimits(18,   0.0, 200.0  );     // nITSandTPCtracksAndSPDtracklets
1536   }
1537   fCFContAsCascadeCuts->SetBinLimits(19,  25.0  ,165.0  );      // BachTPCClusters
1538   */
1539   
1540   // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
1541   fCFContAsCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates associated to MC");
1542   fCFContAsCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates associated to MC");
1543   fCFContAsCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates associated to MC");
1544   fCFContAsCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates associated to MC");
1545   
1546   // Setting the variable title, per axis
1547   // fCFContAsCascadeCuts->SetVarTitle(40,  "Chi2Xi");
1548   fCFContAsCascadeCuts->SetVarTitle(0,  "Dca(XiDaughters) (cm)");
1549   fCFContAsCascadeCuts->SetVarTitle(1,  "Dca(Bach/PrimVertex) (cm)");
1550   fCFContAsCascadeCuts->SetVarTitle(2,  "cos(Xi pointing angle)");
1551   fCFContAsCascadeCuts->SetVarTitle(3,  "R_{2d}(Xi decay) (cm)");
1552   fCFContAsCascadeCuts->SetVarTitle(4,  "M_{#Lambda}(As Casc Dghter) (GeV/c^{2})");
1553    // fCFContAsCascadeCuts->SetVarTitle(40,  "V0Chi2Xi");
1554   fCFContAsCascadeCuts->SetVarTitle(5,  "Dca(V0 Daughters Xi) (cm)");
1555   
1556   fCFContAsCascadeCuts->SetVarTitle(6,  "cos(V0 pointing Angle) in Casc");
1557   fCFContAsCascadeCuts->SetVarTitle(7,  "R_{2d}(V0 decay) (cm)");
1558   fCFContAsCascadeCuts->SetVarTitle(8,  "Dca(V0/PrimVertex) (cm)");
1559   fCFContAsCascadeCuts->SetVarTitle(9,  "Dca(Pos/PrimVertex) (cm)");
1560   fCFContAsCascadeCuts->SetVarTitle(10, "Dca(Neg/PrimVertex) (cm)");
1561   
1562   fCFContAsCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
1563   fCFContAsCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
1564   
1565   fCFContAsCascadeCuts->SetVarTitle(13, "Pt_{MC}(Casc.) (GeV/c)");
1566   //fCFContAsCascadeCuts->SetVarTitle(40, "V0toXiCosineOfPointingAngle");
1567   
1568   fCFContAsCascadeCuts->SetVarTitle(14, "Y_{MC}(Xi)");
1569   fCFContAsCascadeCuts->SetVarTitle(15, "Y_{MC}(Omega)");
1570   
1571   fCFContAsCascadeCuts->SetVarTitle(16, "Z-position(BestPrimVtx) (cm)");
1572   
1573   fCFContAsCascadeCuts->SetVarTitle(17, "Primary Track Multiplicity");
1574   fCFContAsCascadeCuts->SetVarTitle(18, "(ITS+TPC tracks + SPD tracklets) Multiplicity");
1575   fCFContAsCascadeCuts->SetVarTitle(19, "Bach.TPC Clusters");
1576   
1577   fListHistCascade->Add(fCFContAsCascadeCuts);
1578 }
1579
1580
1581 PostData(1, fListHistCascade); 
1582 }// end CreateOutputObjects
1583
1584
1585
1586
1587
1588
1589 //________________________________________________________________________
1590 void AliAnalysisTaskCheckPerformanceCascade::UserExec(Option_t *) 
1591 {
1592         
1593   // Main loop
1594   // Called for each event
1595         
1596         AliESDEvent *lESDevent = 0x0;
1597         AliAODEvent *lAODevent = 0x0;
1598         AliMCEvent  *lMCevent  = 0x0; 
1599         AliStack    *lMCstack  = 0x0; 
1600         Int_t ncascades = -1;
1601         
1602         
1603   // Connect to the InputEvent  
1604   // After these lines, we should have an ESD/AOD event + the number of cascades in it.
1605                 
1606         if(fAnalysisType == "ESD"){
1607                 lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
1608                 if (!lESDevent) {
1609                         Printf("ERROR: lESDevent not available \n");
1610                         cout << "Name of the file with pb :" <<  CurrentFileName() << endl;  // or AliAnalysisTaskSE::CurrentFileName()
1611                         return;
1612                 }
1613         }
1614   
1615         else if(fAnalysisType == "AOD"){  
1616                 lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() ); 
1617                 if (!lAODevent) {
1618                         Printf("ERROR: lAODevent not available \n");
1619                         cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1620                         return;
1621                 }
1622         }
1623         
1624
1625         lMCevent = MCEvent();
1626         if (!lMCevent) {
1627                 Printf("ERROR: Could not retrieve MC event \n");
1628                 cout << "Name of the file with pb :" <<  CurrentFileName() << endl;     
1629                 return;
1630         }
1631
1632         lMCstack = lMCevent->Stack();
1633         if (!lMCstack) {
1634                 Printf("ERROR: Could not retrieve MC stack \n");
1635                 cout << "Name of the file with pb :" <<  CurrentFileName() << endl;
1636                 return;
1637                 
1638         }
1639         
1640         
1641    // Temporary way : AOD awareness of the code to be developed  FIXME
1642    if(fAnalysisType == "AOD") return;
1643    
1644    
1645    
1646    
1647   //-------------------------------------------------
1648   // 0 - Trigger managment + global event selection
1649   // NOTE : Check the availability of the proper trigger 
1650
1651    // Note : Presuppose the presence of AliPhysicsSelectionTask
1652    
1653         UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1654         Bool_t isSelected = 0;
1655         if(     fTriggerMaskType == "kMB")           isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1656         else if(fTriggerMaskType == "kHighMult")     isSelected = (maskIsSelected & AliVEvent::kHighMult) == AliVEvent::kHighMult;
1657         else                                         isSelected = 1; // default = select anyway (use case = run without Phys Selection task)
1658         
1659         if ( ! isSelected ) { 
1660                 PostData(1, fListHistCascade); 
1661                 return;
1662         }
1663         //else Printf("Event selected ... \n");
1664    
1665
1666   //-------------------------------------------------
1667   // 1 - Cascade vertexer (ESD)
1668         if(fkRerunV0CascVertexers){ // FIXME : relaunch V0 and Cascade vertexers
1669                 if(fAnalysisType == "ESD" ){
1670 //                         lESDevent->ResetCascades();
1671 //                         lESDevent->ResetV0s();
1672 // 
1673 //                         AliV0vertexer lV0vtxer;
1674 //                         AliCascadeVertexer lCascVtxer;
1675 // 
1676 //                         lV0vtxer.SetDefaultCuts(fV0Sels);
1677 //                         lCascVtxer.SetDefaultCuts(fCascSels);
1678 // 
1679 //                         lV0vtxer.Tracks2V0vertices(lESDevent);
1680 //                         lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
1681                 }
1682         }
1683
1684
1685   //------------------------------------------------
1686   // 2 - Preparing the general info about of the event = prim. Vtx + magnetic field (ESD)
1687   
1688
1689 //      if(fAnalysisType == "ESD" ){
1690
1691         // Magnetic field
1692                 const Double_t lMagneticField = lESDevent->GetMagneticField( );
1693
1694         // Prim vertex
1695                 const AliESDVertex *lPrimaryTrackingVtx = lESDevent->GetPrimaryVertexTracks();  // get the vtx stored in ESD found with tracks
1696                 const AliESDVertex *lPrimarySPDVtx      = lESDevent->GetPrimaryVertexSPD();     // get the vtx stored in ESD found with SPD tracklets
1697
1698                 const AliESDVertex *lPrimaryBestVtx     = lESDevent->GetPrimaryVertex();
1699                         // get the best primary vertex available for the event
1700                         // As done in AliCascadeVertexer, we keep the one which is the best one available.
1701                         // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1702                         Double_t lBestPrimaryVtxPos[3]   = {-100.0, -100.0, -100.0};
1703                 lPrimaryBestVtx->GetXYZ( lBestPrimaryVtxPos );
1704
1705                 // FIXME : quality cut on the z-position of the prim vertex.
1706                 if(fkQualityCutZprimVtxPos) {
1707                         if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) { 
1708                                 AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !"); 
1709                                 PostData(1, fListHistCascade); 
1710                                 return;
1711                         }
1712                 }
1713                 // FIXME : quality selection regarding pile-up rejection 
1714                 if(fkRejectEventPileUp) {
1715                         if(lESDevent->IsPileupFromSPDInMultBins() ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1716                                 AliWarning("Pb / Event tagged as pile-up by SPD... return !"); 
1717                                 PostData(1, fListHistCascade); 
1718                                 return; 
1719                         }
1720                 }
1721                 // FIXME : remove TPC-only primary vertex : retain only events with tracking + SPD vertex
1722                 if(fkQualityCutNoTPConlyPrimVtx) {
1723                         if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingVtx->GetStatus() ){
1724                                 AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1725                                 PostData(1, fListHistCascade); 
1726                                 return;
1727                         }
1728                 }
1729 //      }// if ESD
1730         
1731         
1732   //    cout << "Name of the accessed file :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1733
1734   //    cout << "Tree characteristics ..." << endl;
1735   //    fInputHandler->GetTree()->Print("toponly");
1736   //    fInputHandler->GetTree()->GetBranch("PrimaryVertex")->Print();
1737   //    fInputHandler->GetTree()->GetBranch("SPDVertex")->Print();
1738
1739
1740
1741   // ---------------------------------------------------------------
1742   // - Initialisation of the part dedicated to cascade vertices
1743
1744   if(fAnalysisType == "ESD")            ncascades = lESDevent->GetNumberOfCascades();
1745   else if(fAnalysisType == "AOD")       ncascades = lAODevent->GetNumberOfCascades();
1746         
1747   
1748   Int_t   nNumberOfMCPrimaries       = -1;
1749   Int_t   nMCPrimariesInEtaBelow0p8  =  0;
1750   Int_t   nMCPrimariesInEtaBelow1p0  =  0;
1751   
1752   Int_t   nTrackPrimaryMultiplicity        = -1;
1753   Int_t   nSPDTracklets                    =  0; // AliESDEvent::EstimateMultiplicity will re-initialise the value to 0
1754   Int_t   nITSandTPCtracksAndSPDtracklets  =  0; // AliESDEvent::EstimateMultiplicity will re-initialise the value to 0
1755   Int_t   nTracksITSSApure                 =  0; // AliESDEvent::EstimateMultiplicity will re-initialise the value to 0
1756   
1757         nNumberOfMCPrimaries       = lMCstack->GetNprimary();
1758         if(nNumberOfMCPrimaries < 1) return;
1759         
1760         nTrackPrimaryMultiplicity  = fESDtrackCuts->CountAcceptedTracks(lESDevent);
1761         //EstimateMultiplicity(Int_t &tracklets, Int_t &trITSTPC, Int_t &trITSSApure, Double_t eta, Bool_t useDCAFlag,Bool_t useV0Flag)
1762         lESDevent->EstimateMultiplicity( nSPDTracklets, nITSandTPCtracksAndSPDtracklets, nTracksITSSApure, 1.0, kTRUE, kTRUE);
1763
1764
1765         fHistMCTrackMultiplicity->Fill( nNumberOfMCPrimaries );
1766     
1767 //_____________________________________________________________________________ 
1768 // Part 1 - Loop over the MC primaries  
1769         
1770     for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) 
1771     {// This is the begining of the loop on primaries
1772       
1773         TParticle* lCurrentParticle = 0x0; 
1774                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1775         if(!lCurrentParticle){
1776                 Printf("MC Primary loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1777                 continue;
1778         }
1779         
1780         Double_t lEtaCurrentParticle = TMath::Abs( lCurrentParticle->Eta() );
1781         if( lEtaCurrentParticle < 1.0 ){
1782                 nMCPrimariesInEtaBelow1p0++;
1783                 if( lEtaCurrentParticle < 0.8 ) nMCPrimariesInEtaBelow0p8++;
1784         }
1785     }
1786     
1787     f2dHistRecoPrimTrckMultVsMCMult->Fill( nTrackPrimaryMultiplicity,       nMCPrimariesInEtaBelow0p8 );
1788     f2dHistRecoEstimateMultVsMCMult->Fill( nITSandTPCtracksAndSPDtracklets, nMCPrimariesInEtaBelow1p0 );
1789
1790
1791    // For proton
1792    /*
1793    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) 
1794         {// This is the begining of the loop on primaries, for protons
1795           
1796         TParticle* lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1797         if(!lCurrentParticle){
1798                 Printf("Proton loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1799                 continue;
1800                 
1801         }
1802         
1803         if( lCurrentParticle->GetPdgCode() == 2212 )
1804                 fHistEtaGenProton->Fill( lCurrentParticle->Eta() );
1805
1806         if( lCurrentParticle->GetPdgCode() == -2212 )
1807                 fHistEtaGenAntiProton->Fill( lCurrentParticle->Eta() );
1808         }// end loop over primary proton
1809    */
1810
1811       
1812        
1813 //_____________________________________________________________________________ 
1814 // Part 2 - Loop over the different types of GENERATED cascades (Xi-+, Omega-+) 
1815
1816         // - Initialisation of useful local variables
1817                 
1818         Int_t lPdgCodeCasc            = 0;
1819         Int_t lPdgCodeBach            = 0;
1820         Int_t lPdgCodeLambda          = 0;
1821         Int_t lPdgCodeDghtMesV0       = 0;
1822         Int_t lPdgCodeDghtBarV0       = 0;
1823         
1824         
1825         TH1F *lHistEtaGenCasc         = 0;      
1826         TH2F *l2dHistGenPtVsGenYGen   = 0;
1827                 
1828         TH1F *lHistThetaGenCasc       = 0;
1829         TH2F *l2dHistGenPtVsGenYFdbl  = 0;
1830         TH1F *lHistThetaLambda        = 0;
1831         TH1F *lHistThetaBach          = 0;
1832         TH1F *lHistThetaBarDghter     = 0;
1833         TH1F *lHistThetaMesDghter     = 0;
1834         TH1F *lHistPtBach             = 0;
1835         TH1F *lHistPtBarDghter        = 0;
1836         TH1F *lHistPtMesDghter        = 0;
1837
1838
1839 for(Int_t iCascType = 1; iCascType < 5; iCascType++)
1840 {
1841        
1842 switch (iCascType)
1843   {
1844     case 1: // Xi-
1845          lPdgCodeCasc       =   3312;  //Xi-
1846          lPdgCodeBach       =   -211;  //Pi-
1847          lPdgCodeLambda     =   3122;  //Lambda0
1848          lPdgCodeDghtMesV0  =   -211;  //Pi-
1849          lPdgCodeDghtBarV0  =   2212;  //Proton 
1850                 
1851                 // any Xi-
1852          lHistEtaGenCasc        = fHistEtaGenCascXiMinus;
1853          l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenXiMinus;
1854         
1855                 // cascades generated within acceptance (cut in pt + theta)
1856          lHistThetaGenCasc      = fHistThetaGenCascXiMinus;
1857          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiMinus;
1858          lHistThetaLambda       = fHistThetaLambdaXiMinus;
1859          lHistThetaBach         = fHistThetaBachXiMinus;
1860          lHistThetaBarDghter    = fHistThetaBarDghterXiMinus;
1861          lHistThetaMesDghter    = fHistThetaMesDghterXiMinus;
1862          lHistPtBach            = fHistPtBachXiMinus;
1863          lHistPtBarDghter       = fHistPtBarDghterXiMinus;
1864          lHistPtMesDghter       = fHistPtMesDghterXiMinus;
1865         break; 
1866            
1867     case 2: // Xi+
1868          lPdgCodeCasc        =  -3312;  //Xi+
1869          lPdgCodeBach        =    211;  //Pi+
1870          lPdgCodeLambda      =  -3122;  //AntiLambda0
1871          lPdgCodeDghtMesV0   =    211;  //Pi+
1872          lPdgCodeDghtBarV0   =  -2212;  //AntiProton  
1873          
1874                 // any Xi+
1875          lHistEtaGenCasc        = fHistEtaGenCascXiPlus;
1876          l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenXiPlus;
1877         
1878                 // cascades generated within acceptance (cut in pt + theta)      
1879          lHistThetaGenCasc      = fHistThetaGenCascXiPlus;
1880          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblXiPlus;
1881          lHistThetaLambda       = fHistThetaLambdaXiPlus;
1882          lHistThetaBach         = fHistThetaBachXiPlus;
1883          lHistThetaBarDghter    = fHistThetaBarDghterXiPlus;
1884          lHistThetaMesDghter    = fHistThetaMesDghterXiPlus;
1885          lHistPtBach            = fHistPtBachXiPlus;
1886          lHistPtBarDghter       = fHistPtBarDghterXiPlus;
1887          lHistPtMesDghter       = fHistPtMesDghterXiPlus;  
1888         break;
1889    
1890     case 3: // Omega-
1891          lPdgCodeCasc       =   3334;  //Omega-
1892          lPdgCodeBach       =   -321;  //K-
1893          lPdgCodeLambda     =   3122;  //Lambda0
1894          lPdgCodeDghtMesV0  =   -211;  //Pi-
1895          lPdgCodeDghtBarV0  =   2212;  //Proton 
1896          
1897                 // any Omega-
1898          lHistEtaGenCasc        = fHistEtaGenCascOmegaMinus;            
1899          l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenOmegaMinus;      
1900         
1901                 // cascades generated within acceptance (cut in pt + theta)
1902          lHistThetaGenCasc      = fHistThetaGenCascOmegaMinus;
1903          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaMinus;
1904          lHistThetaLambda       = fHistThetaLambdaOmegaMinus;
1905          lHistThetaBach         = fHistThetaBachOmegaMinus;
1906          lHistThetaBarDghter    = fHistThetaBarDghterOmegaMinus;
1907          lHistThetaMesDghter    = fHistThetaMesDghterOmegaMinus;
1908          lHistPtBach            = fHistPtBachOmegaMinus;
1909          lHistPtBarDghter       = fHistPtBarDghterOmegaMinus;
1910          lHistPtMesDghter       = fHistPtMesDghterOmegaMinus;   
1911         break;
1912     
1913     case 4:  // Omega+
1914          lPdgCodeCasc       =  -3334;  //Omega+
1915          lPdgCodeBach       =    321;  //K+
1916          lPdgCodeLambda     =  -3122;  //AntiLambda0
1917          lPdgCodeDghtMesV0  =    211;  //Pi+
1918          lPdgCodeDghtBarV0  =  -2212;  //AntiProton 
1919          
1920                 // any Omega+
1921          lHistEtaGenCasc        = fHistEtaGenCascOmegaPlus;
1922          l2dHistGenPtVsGenYGen  = f2dHistGenPtVsGenYGenOmegaPlus;               
1923                 
1924                 // cascades generated within acceptance (cut in pt + theta)
1925          lHistThetaGenCasc      = fHistThetaGenCascOmegaPlus;
1926          l2dHistGenPtVsGenYFdbl = f2dHistGenPtVsGenYFdblOmegaPlus;
1927          lHistThetaLambda       = fHistThetaLambdaOmegaPlus;
1928          lHistThetaBach         = fHistThetaBachOmegaPlus;
1929          lHistThetaBarDghter    = fHistThetaBarDghterOmegaPlus;
1930          lHistThetaMesDghter    = fHistThetaMesDghterOmegaPlus;
1931          lHistPtBach            = fHistPtBachOmegaPlus;
1932          lHistPtBarDghter       = fHistPtBarDghterOmegaPlus;
1933          lHistPtMesDghter       = fHistPtMesDghterOmegaPlus;  
1934         break;
1935
1936   }// end switch cascade
1937
1938
1939    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < nNumberOfMCPrimaries; iCurrentLabelStack++) 
1940     {// This is the begining of the loop on primaries
1941       
1942         TParticle* lCurrentParticle = 0x0; 
1943                    lCurrentParticle = lMCstack->Particle( iCurrentLabelStack );
1944         if(!lCurrentParticle){
1945                 Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1946                 continue;
1947                 
1948         }
1949          
1950         if( lCurrentParticle->GetPdgCode() == lPdgCodeCasc ){  // Here !
1951                 //cout << "Xi- within loop " << iCurrentLabelStack << "/ " << nNumberOfMCPrimaries << endl;
1952                 
1953                 // -  Xi level ... _____________________________________________________________
1954                 TParticle* xiMC = 0x0;
1955                            xiMC = lCurrentParticle;
1956                 if(!xiMC){
1957                         Printf("MC TParticle pointer to Cascade = 0x0 ! Skip ...");
1958                         continue;
1959                 
1960                 }
1961                 
1962                 // Fill the first histos : = any generated Xi, not necessarily within the acceptance
1963                 Double_t lRapXiMC = 0.5*TMath::Log((xiMC->Energy() + xiMC->Pz()) / (xiMC->Energy() - xiMC->Pz() +1.e-13));
1964                 
1965                 lHistEtaGenCasc         ->Fill( xiMC->Eta() );   
1966                 l2dHistGenPtVsGenYGen   ->Fill( xiMC->Pt(), lRapXiMC  );        
1967                         
1968                 
1969                 
1970                 // Check the emission of particle stays within the acceptance of the detector (cut in theta)
1971                 if( xiMC->Theta() < TMath::Pi()/4.0  ||    xiMC->Theta() > 3.0*TMath::Pi()/4.0 ) continue;      
1972                 if( xiMC->GetNDaughters() != 2) continue;
1973                 if( xiMC->GetDaughter(0) < 0 )  continue;
1974                 if( xiMC->GetDaughter(1) < 0 )  continue;
1975                 
1976                         TParticle* lDght0ofXi = lMCstack->Particle(  xiMC->GetDaughter(0) );
1977                         TParticle* lDght1ofXi = lMCstack->Particle(  xiMC->GetDaughter(1) );
1978                         
1979                 TParticle* lLambda = 0;
1980                 TParticle* lBach   = 0;
1981                         
1982                 // Xi - Case 1
1983                         if(     lDght0ofXi->GetPdgCode() == lPdgCodeLambda   &&  // Here !
1984                                 lDght1ofXi->GetPdgCode() == lPdgCodeBach ){      // Here !
1985                                 
1986                                 lLambda = lDght0ofXi;
1987                                 lBach   = lDght1ofXi;
1988                         }// end if dghter 0 = Lambda and    dghter 1 = Pi-  
1989                         
1990                 // Xi - Case 2
1991                         else if( lDght0ofXi->GetPdgCode() == lPdgCodeBach  &&      // Here !
1992                                  lDght1ofXi->GetPdgCode() == lPdgCodeLambda ){     // Here !
1993                         
1994                                 lBach   = lDght0ofXi;
1995                                 lLambda = lDght1ofXi;
1996                         }//  end if dghter 0 = Pi-  and   dghter 1 = Lambda
1997                         
1998                 // V0 otherwise - Case 3        
1999                         else continue;
2000                 
2001                         // Check the emission of particle stays within the acceptance of the detector (cut in pt + theta)
2002                         if( lLambda->Theta() < TMath::Pi()/4.0  ||    lLambda->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2003                         if( lBach->Theta() < TMath::Pi()/4.0    ||    lBach->Theta() > 3.0*TMath::Pi()/4.0 )   continue;
2004                 
2005                         if( lBach->Pt() < 0.150 ) continue; //FIXME : maybe tuned for Xi but not for K- from Omega ...
2006                         
2007                 
2008                 
2009                 // -  V0 level ... _____________________________________________________________
2010                 TParticle* lDghtBarV0 = 0;
2011                 TParticle* lDghtMesV0 = 0;
2012                 
2013                 if( lLambda->GetNDaughters() != 2 )  continue;
2014                 if( lLambda->GetDaughter(0) < 0 )    continue;
2015                 if( lLambda->GetDaughter(1) < 0 )    continue;
2016                 
2017                 
2018                 TParticle* lDght0ofLambda = lMCstack->Particle(  lLambda->GetDaughter(0) );
2019                 TParticle* lDght1ofLambda = lMCstack->Particle(  lLambda->GetDaughter(1) );
2020                         
2021                 // V0 - Case 1
2022                         if(     lDght0ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 &&    // Here !
2023                                 lDght1ofLambda->GetPdgCode() == lPdgCodeDghtMesV0 ){    // Here !
2024                         
2025                                 lDghtBarV0 = lDght0ofLambda;
2026                                 lDghtMesV0 = lDght1ofLambda;
2027                         }// end if dghter 0 = Proton  and   dghter 1 = Pi-  
2028                         
2029                 // V0 - Case 2
2030                         else if( lDght0ofLambda->GetPdgCode() == lPdgCodeDghtMesV0  &&     // Here !
2031                                  lDght1ofLambda->GetPdgCode() == lPdgCodeDghtBarV0 ){      // Here !
2032                         
2033                                 lDghtMesV0 = lDght0ofLambda;
2034                                 lDghtBarV0 = lDght1ofLambda;
2035                         }//  end if dghter 0 = Pi-  and   dghter 1 = proton
2036                         
2037                 // V0 otherwise - Case 3
2038                         else continue;
2039         
2040                         
2041                         // Check the emission of particle stays within the acceptance of the detector
2042                         if( lDghtBarV0->Theta() < TMath::Pi()/4.0  ||  lDghtBarV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2043                         if( lDghtMesV0->Theta() < TMath::Pi()/4.0  ||  lDghtMesV0->Theta() > 3.0*TMath::Pi()/4.0 ) continue;
2044                         
2045                         if( lDghtBarV0->Pt() < 0.250 ) continue;
2046                         if( lDghtMesV0->Pt() < 0.150 ) continue;
2047                         
2048                         
2049                         
2050                 // - Just to know which file is currently open : locate the file containing Xi 
2051                 //cout << "Name of the file containing generated Xi :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() 
2052                 //                                                   <<  endl;  
2053                         
2054                 Double_t lRadToDeg = 180.0/TMath::Pi(); 
2055                         
2056                 // - Filling histos ... _________________________________________________________________       
2057                         lHistThetaGenCasc       ->Fill( lRadToDeg * xiMC->Theta()  );
2058                         l2dHistGenPtVsGenYFdbl  ->Fill( xiMC->Pt(), lRapXiMC );
2059                         
2060                         // - Fill theta histos for Lambda and Bach
2061                         lHistThetaLambda        ->Fill( lRadToDeg * lLambda->Theta() );
2062                         lHistThetaBach          ->Fill( lRadToDeg *   lBach->Theta() );
2063                         
2064                         // - Fill theta histos for V0 daughters
2065                         lHistThetaBarDghter     ->Fill( lRadToDeg * lDghtBarV0->Theta() );
2066                         lHistThetaMesDghter     ->Fill( lRadToDeg * lDghtMesV0->Theta() );
2067                         
2068                         // - Fill pt histos.
2069                         lHistPtBach             ->Fill(      lBach->Pt() );
2070                         lHistPtBarDghter        ->Fill( lDghtBarV0->Pt() );
2071                         lHistPtMesDghter        ->Fill( lDghtMesV0->Pt() );
2072                                                 
2073                 }// end if current particle = Xi-
2074              
2075      }// This is the end of the loop on primaries
2076      
2077 // - Re-initialisation of the local TH1F pointers
2078 lHistEtaGenCasc         = 0x0;
2079 l2dHistGenPtVsGenYGen   = 0x0;
2080
2081 lHistThetaGenCasc       = 0x0;
2082 l2dHistGenPtVsGenYFdbl  = 0x0;
2083 lHistThetaLambda        = 0x0;
2084 lHistThetaBach          = 0x0;
2085 lHistThetaBarDghter     = 0x0;
2086 lHistThetaMesDghter     = 0x0;
2087 lHistPtBach             = 0x0;
2088 lHistPtBarDghter        = 0x0;
2089 lHistPtMesDghter        = 0x0;  
2090
2091 } // end of loop over the different types of cascades (Xi-+, Omega-+)
2092         
2093  
2094  
2095 //__________________________________________________________________________    
2096 // Part 3 - Loop over the reconstructed candidates
2097   
2098
2099 for (Int_t iXi = 0; iXi < ncascades; iXi++) 
2100 {// This is the begining of the Cascade loop
2101                 
2102         AliESDcascade *xiESD = lESDevent->GetCascade(iXi);
2103         if (!xiESD) continue;
2104         
2105         // - Step II.1 : Connection to daughter tracks of the current cascade
2106         //-------------
2107                         
2108                 UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xiESD->GetPindex() );
2109                 UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xiESD->GetNindex() );
2110                 UInt_t lBachIdx         = (UInt_t) TMath::Abs( xiESD->GetBindex() );
2111                 // abs value not needed ; the index should always be positive (!= label ...)
2112                 
2113                 
2114         // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
2115         if(lBachIdx == lIdxNegXi) {
2116                 AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
2117         }
2118         if(lBachIdx == lIdxPosXi) {
2119                 AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
2120         }
2121       
2122         AliESDtrack *pTrackXi           = lESDevent->GetTrack( lIdxPosXi );
2123         AliESDtrack *nTrackXi           = lESDevent->GetTrack( lIdxNegXi );
2124         AliESDtrack *bachTrackXi        = lESDevent->GetTrack( lBachIdx  );
2125         if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
2126                 Printf("ERROR: Could not retrieve one of the 3 daughter tracks of the cascade ...");
2127                 continue;
2128         }
2129         
2130         Int_t lPosTPCClusters   = pTrackXi->GetTPCNcls();
2131         Int_t lNegTPCClusters   = nTrackXi->GetTPCNcls();
2132         Int_t lBachTPCClusters  = bachTrackXi->GetTPCNcls(); 
2133
2134                 // FIXME : rejection of a poor quality tracks
2135         if(fkQualityCutTPCrefit){
2136                 // 1 - Poor quality related to TPCrefit
2137                 ULong_t pStatus    = pTrackXi->GetStatus();
2138                 ULong_t nStatus    = nTrackXi->GetStatus();
2139                 ULong_t bachStatus = bachTrackXi->GetStatus();
2140                 if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
2141                 if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
2142                 if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
2143         }
2144         if(fkQualityCut80TPCcls){
2145                 // 2 - Poor quality related to TPC clusters
2146                 if(lPosTPCClusters  < 80) { AliWarning("Pb / V0 Pos. track has less than 80 TPC clusters ... continue!"); continue; }
2147                 if(lNegTPCClusters  < 80) { AliWarning("Pb / V0 Neg. track has less than 80 TPC clusters ... continue!"); continue; }
2148                 if(lBachTPCClusters < 80) { AliWarning("Pb / Bach.   track has less than 80 TPC clusters ... continue!"); continue; }
2149         }
2150         
2151         // - Step II.2 : Info over reconstructed cascades
2152         //------------- 
2153         
2154         Double_t lInvMassXiMinus    = 0.;
2155         Double_t lInvMassXiPlus     = 0.;
2156         Double_t lInvMassOmegaMinus = 0.;
2157         Double_t lInvMassOmegaPlus  = 0.;
2158         
2159         Double_t lV0quality = 0.;
2160         
2161         if( bachTrackXi->Charge() < 0 ) {
2162                 lV0quality = 0.;
2163                 xiESD->ChangeMassHypothesis(lV0quality , 3312);         
2164                         // Calculate the effective mass of the Xi- candidate. 
2165                         // pdg code 3312 = Xi-
2166                 lInvMassXiMinus = xiESD->GetEffMassXi();
2167                 
2168                 lV0quality = 0.;
2169                 xiESD->ChangeMassHypothesis(lV0quality , 3334);         
2170                         // Calculate the effective mass of the Xi- candidate. 
2171                         // pdg code 3334 = Omega-
2172                 lInvMassOmegaMinus = xiESD->GetEffMassXi();
2173                                         
2174                 lV0quality = 0.;
2175                 xiESD->ChangeMassHypothesis(lV0quality , 3312);         // Back to default hyp.
2176                 
2177         }
2178         
2179         if( bachTrackXi->Charge() >  0 ){
2180                 lV0quality = 0.;
2181                 xiESD->ChangeMassHypothesis(lV0quality , -3312);        
2182                         // Calculate the effective mass of the Xi+ candidate. 
2183                         // pdg code -3312 = Xi+
2184                 lInvMassXiPlus = xiESD->GetEffMassXi();
2185                 
2186                 lV0quality = 0.;
2187                 xiESD->ChangeMassHypothesis(lV0quality , -3334);        
2188                         // Calculate the effective mass of the Xi+ candidate. 
2189                         // pdg code -3334  = Omega+
2190                 lInvMassOmegaPlus = xiESD->GetEffMassXi();
2191                 
2192                 lV0quality = 0.;
2193                 xiESD->ChangeMassHypothesis(lV0quality , -3312);        // Back to "default" hyp.
2194         }
2195                 
2196         Double_t lChargeXi = xiESD->Charge();
2197         
2198         if( lChargeXi < 0 )     fHistMassXiMinus        ->Fill( lInvMassXiMinus );
2199         if( lChargeXi > 0 )     fHistMassXiPlus         ->Fill( lInvMassXiPlus );
2200         if( lChargeXi < 0 )     fHistMassOmegaMinus     ->Fill( lInvMassOmegaMinus );
2201         if( lChargeXi > 0 )     fHistMassOmegaPlus      ->Fill( lInvMassOmegaPlus );
2202         
2203         
2204         // - Step II.3 : PID info
2205         //-------------
2206         
2207         
2208         // 3.1 - PID Information
2209
2210         Bool_t   lIsPosInXiProton      = kFALSE;
2211         Bool_t   lIsPosInXiPion        = kFALSE;
2212         Bool_t   lIsPosInOmegaProton   = kFALSE;
2213         Bool_t   lIsPosInOmegaPion     = kFALSE;
2214         
2215         Bool_t   lIsNegInXiProton      = kFALSE;
2216         Bool_t   lIsNegInXiPion        = kFALSE;
2217         Bool_t   lIsNegInOmegaProton   = kFALSE;
2218         Bool_t   lIsNegInOmegaPion     = kFALSE;
2219         
2220         Bool_t   lIsBachelorKaon       = kFALSE;
2221         Bool_t   lIsBachelorPion       = kFALSE; 
2222         
2223         Bool_t   lIsBachelorKaonForTPC = kFALSE; // For ESD only ...//FIXME : wait for availability in AOD
2224         Bool_t   lIsBachelorPionForTPC = kFALSE; // For ESD only ...
2225         Bool_t   lIsNegPionForTPC      = kFALSE; // For ESD only ...
2226         Bool_t   lIsPosPionForTPC      = kFALSE; // For ESD only ...
2227         Bool_t   lIsNegProtonForTPC    = kFALSE; // For ESD only ...
2228         Bool_t   lIsPosProtonForTPC    = kFALSE; // For ESD only ...
2229
2230         // 3.1.A - Combined PID
2231         // Reasonable guess for the priors for the cascade track sample (e-, mu, pi, K, p)
2232         Double_t lPriorsGuessXi[10]    = {0, 0, 2, 0, 1,  0,0,0,0,0};
2233         Double_t lPriorsGuessOmega[10] = {0, 0, 1, 1, 1,  0,0,0,0,0};
2234         
2235         // Combined VO-positive-daughter PID
2236         AliPID pPidXi;         pPidXi.SetPriors(    lPriorsGuessXi    , kTRUE); // kTRUE = for charged particle PID
2237         AliPID pPidOmega;      pPidOmega.SetPriors( lPriorsGuessOmega , kTRUE); // kTRUE = for charged particle PID
2238                 
2239         if( pTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2240                 Double_t r[10] = {0.}; pTrackXi->GetESDpid(r);
2241                 pPidXi.SetProbabilities(r);
2242                 pPidOmega.SetProbabilities(r);
2243                 
2244                 // Check if the V0 positive track is a proton (case for Xi-)
2245                 Double_t pproton = pPidXi.GetProbability(AliPID::kProton);
2246                 if (pproton > pPidXi.GetProbability(AliPID::kElectron) &&
2247                     pproton > pPidXi.GetProbability(AliPID::kMuon)     &&
2248                     pproton > pPidXi.GetProbability(AliPID::kPion)     &&
2249                     pproton > pPidXi.GetProbability(AliPID::kKaon)     )     lIsPosInXiProton = kTRUE;
2250                 
2251                 // Check if the V0 positive track is a pi+ (case for Xi+)
2252                 Double_t ppion = pPidXi.GetProbability(AliPID::kPion);
2253                 if (ppion > pPidXi.GetProbability(AliPID::kElectron) &&
2254                     ppion > pPidXi.GetProbability(AliPID::kMuon)     &&
2255                     ppion > pPidXi.GetProbability(AliPID::kKaon)     &&
2256                     ppion > pPidXi.GetProbability(AliPID::kProton)   )     lIsPosInXiPion = kTRUE;
2257                 
2258                 
2259                 // Check if the V0 positive track is a proton (case for Omega-)
2260                 pproton = 0.;
2261                     pproton = pPidOmega.GetProbability(AliPID::kProton);
2262                 if (pproton > pPidOmega.GetProbability(AliPID::kElectron) &&
2263                     pproton > pPidOmega.GetProbability(AliPID::kMuon)     &&
2264                     pproton > pPidOmega.GetProbability(AliPID::kPion)     &&
2265                     pproton > pPidOmega.GetProbability(AliPID::kKaon)     )  lIsPosInOmegaProton = kTRUE;
2266                 
2267                 // Check if the V0 positive track is a pi+ (case for Omega+)
2268                 ppion = 0.;
2269                     ppion = pPidOmega.GetProbability(AliPID::kPion);
2270                 if (ppion > pPidOmega.GetProbability(AliPID::kElectron) &&
2271                     ppion > pPidOmega.GetProbability(AliPID::kMuon)     &&
2272                     ppion > pPidOmega.GetProbability(AliPID::kKaon)     &&
2273                     ppion > pPidOmega.GetProbability(AliPID::kProton)   )    lIsPosInOmegaPion = kTRUE;
2274                 
2275         }// end if V0 positive track with existing combined PID 
2276         
2277         
2278         // Combined VO-negative-daughter PID
2279         AliPID nPidXi;         nPidXi.SetPriors(    lPriorsGuessXi    , kTRUE); // kTRUE = for charged particle PID
2280         AliPID nPidOmega;      nPidOmega.SetPriors( lPriorsGuessOmega , kTRUE); // kTRUE = for charged particle PID
2281                 
2282         if( nTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2283                 Double_t r[10] = {0.}; nTrackXi->GetESDpid(r);
2284                 nPidXi.SetProbabilities(r);
2285                 nPidOmega.SetProbabilities(r);
2286                 
2287                 // Check if the V0 negative track is a pi- (case for Xi-)
2288                 Double_t ppion = nPidXi.GetProbability(AliPID::kPion);
2289                 if (ppion > nPidXi.GetProbability(AliPID::kElectron) &&
2290                     ppion > nPidXi.GetProbability(AliPID::kMuon)     &&
2291                     ppion > nPidXi.GetProbability(AliPID::kKaon)     &&
2292                     ppion > nPidXi.GetProbability(AliPID::kProton)   )     lIsNegInXiPion = kTRUE;
2293
2294                 // Check if the V0 negative track is an anti-proton (case for Xi+)
2295                 Double_t pproton = nPidXi.GetProbability(AliPID::kProton);
2296                 if (pproton > nPidXi.GetProbability(AliPID::kElectron) &&
2297                     pproton > nPidXi.GetProbability(AliPID::kMuon)     &&
2298                     pproton > nPidXi.GetProbability(AliPID::kPion)     &&
2299                     pproton > nPidXi.GetProbability(AliPID::kKaon)     )     lIsNegInXiProton = kTRUE;
2300                 
2301                 // Check if the V0 negative track is a pi- (case for Omega-)
2302                 ppion = 0.;
2303                     ppion = nPidOmega.GetProbability(AliPID::kPion);
2304                 if (ppion > nPidOmega.GetProbability(AliPID::kElectron) &&
2305                     ppion > nPidOmega.GetProbability(AliPID::kMuon)     &&
2306                     ppion > nPidOmega.GetProbability(AliPID::kKaon)     &&
2307                     ppion > nPidOmega.GetProbability(AliPID::kProton)   )    lIsNegInOmegaPion = kTRUE;
2308                 
2309                 // Check if the V0 negative track is an anti-proton (case for Omega+)
2310                 pproton = 0.;
2311                     pproton = nPidOmega.GetProbability(AliPID::kProton);
2312                 if (pproton > nPidOmega.GetProbability(AliPID::kElectron) &&
2313                     pproton > nPidOmega.GetProbability(AliPID::kMuon)     &&
2314                     pproton > nPidOmega.GetProbability(AliPID::kPion)     &&
2315                     pproton > nPidOmega.GetProbability(AliPID::kKaon)     )  lIsNegInOmegaProton = kTRUE;
2316                 
2317         }// end if V0 negative track with existing combined PID 
2318         
2319                 
2320         // Combined bachelor PID
2321         AliPID bachPidXi;      bachPidXi.SetPriors(    lPriorsGuessXi    , kTRUE); // kTRUE = for charged particle PID
2322         AliPID bachPidOmega;   bachPidOmega.SetPriors( lPriorsGuessOmega , kTRUE); // kTRUE = for charged particle PID
2323         
2324         Double_t ppionBach = 0.0, pkaonBach = 0.0;
2325         
2326         if( bachTrackXi->IsOn(AliESDtrack::kESDpid) ){  // Combined PID exists
2327                 Double_t r[10] = {0.}; bachTrackXi->GetESDpid(r);
2328                 bachPidXi.SetProbabilities(r);
2329                 bachPidOmega.SetProbabilities(r);
2330                 // Check if the bachelor track is a pion
2331                     ppionBach = bachPidXi.GetProbability(AliPID::kPion);
2332                 if (ppionBach > bachPidXi.GetProbability(AliPID::kElectron) &&
2333                     ppionBach > bachPidXi.GetProbability(AliPID::kMuon)     &&
2334                     ppionBach > bachPidXi.GetProbability(AliPID::kKaon)     &&
2335                     ppionBach > bachPidXi.GetProbability(AliPID::kProton)   )     lIsBachelorPion = kTRUE;
2336                 // Check if the bachelor track is a kaon
2337                     pkaonBach = bachPidOmega.GetProbability(AliPID::kKaon);
2338                 if (pkaonBach > bachPidOmega.GetProbability(AliPID::kElectron) &&
2339                     pkaonBach > bachPidOmega.GetProbability(AliPID::kMuon)     &&
2340                     pkaonBach > bachPidOmega.GetProbability(AliPID::kPion)     &&
2341                     pkaonBach > bachPidOmega.GetProbability(AliPID::kProton)   )  lIsBachelorKaon = kTRUE;      
2342         }// end if bachelor track with existing combined PID
2343         
2344         
2345         // 3.1.B - TPC PID : 4-sigma bands on Bethe-Bloch curve
2346         
2347         // Bachelor
2348         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
2349         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
2350         
2351         // Negative V0 daughter
2352         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
2353         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
2354         
2355         // Positive V0 daughter
2356         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
2357         if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
2358         
2359         /*
2360         const AliExternalTrackParam *pInnerWallTrackXi    = pTrackXi    ->GetInnerParam(); // Do not use GetTPCInnerWall
2361         const AliExternalTrackParam *nInnerWallTrackXi    = nTrackXi    ->GetInnerParam();
2362         const AliExternalTrackParam *bachInnerWallTrackXi = bachTrackXi ->GetInnerParam();
2363         if(pInnerWallTrackXi && nInnerWallTrackXi && bachInnerWallTrackXi ){
2364                 
2365                 Double_t pMomInnerWall    = pInnerWallTrackXi   ->GetP();
2366                 Double_t nMomInnerWall    = nInnerWallTrackXi   ->GetP();
2367                 Double_t bachMomInnerWall = bachInnerWallTrackXi->GetP();
2368                 
2369                 // Bachelor
2370                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 3)                              lIsBachelorPionForTPC = kTRUE;
2371                 if (bachMomInnerWall < 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 5) lIsBachelorKaonForTPC = kTRUE;
2372                 if (bachMomInnerWall > 0.350  && TMath::Abs(fESDpid->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 3) lIsBachelorKaonForTPC = kTRUE;
2373                 
2374                 // Negative V0 daughter
2375                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 3  )                           lIsNegPionForTPC   = kTRUE;
2376                 if (nMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 5 )   lIsNegProtonForTPC = kTRUE;
2377                 if (nMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( nTrackXi,AliPID::kProton ) ) < 3 )   lIsNegProtonForTPC = kTRUE;
2378                 
2379                 // Positive V0 daughter
2380                 if (TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 3 )                            lIsPosPionForTPC   = kTRUE;
2381                 if (pMomInnerWall < 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 5)     lIsPosProtonForTPC = kTRUE;
2382                 if (pMomInnerWall > 0.6  && TMath::Abs(fESDpid->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 3)     lIsPosProtonForTPC = kTRUE;
2383                 
2384         }
2385         */
2386                 
2387         // Combined PID TH1s
2388         if( lChargeXi < 0 && lIsBachelorPion )    fHistMassWithCombPIDXiMinus     ->Fill( lInvMassXiMinus    );
2389         if( lChargeXi > 0 && lIsBachelorPion )    fHistMassWithCombPIDXiPlus      ->Fill( lInvMassXiPlus     );
2390         if( lChargeXi < 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
2391         if( lChargeXi > 0 && lIsBachelorKaon )    fHistMassWithCombPIDOmegaPlus   ->Fill( lInvMassOmegaPlus  );
2392          
2393         
2394         // 3.2 - PID proba Vs Pt(Bach)
2395         Int_t      lblBachForPID  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2396         TParticle* mcBachForPID   = lMCstack->Particle( lblBachForPID );
2397         Double_t   lmcPtBach      = mcBachForPID->Pt();
2398         
2399         if(lIsBachelorPion)   f2dHistPIDprobaPionVsMCPtBach->Fill( lmcPtBach, ppionBach );
2400         if(lIsBachelorKaon)   f2dHistPIDprobaKaonVsMCPtBach->Fill( lmcPtBach, pkaonBach );
2401         
2402                         
2403         // 3.3 - MC perfect PID
2404         Bool_t   lIsBachelorMCPiMinus  = kFALSE;
2405         Bool_t   lIsBachelorMCPiPlus   = kFALSE;
2406         Bool_t   lIsBachelorMCKMinus   = kFALSE;
2407         Bool_t   lIsBachelorMCKPlus    = kFALSE;        
2408         
2409         if( mcBachForPID->GetPdgCode() == -211) lIsBachelorMCPiMinus = kTRUE;
2410         if( mcBachForPID->GetPdgCode() ==  211) lIsBachelorMCPiPlus  = kTRUE;
2411         if( mcBachForPID->GetPdgCode() == -321) lIsBachelorMCKMinus  = kTRUE;
2412         if( mcBachForPID->GetPdgCode() ==  321) lIsBachelorMCKPlus   = kTRUE;
2413         
2414         if( lChargeXi < 0 && lIsBachelorMCPiMinus )    fHistMassWithMcPIDXiMinus     ->Fill( lInvMassXiMinus );
2415         if( lChargeXi > 0 && lIsBachelorMCPiPlus  )    fHistMassWithMcPIDXiPlus      ->Fill( lInvMassXiPlus );
2416         if( lChargeXi < 0 && lIsBachelorMCKMinus  )    fHistMassWithMcPIDOmegaMinus  ->Fill( lInvMassOmegaMinus );
2417         if( lChargeXi > 0 && lIsBachelorMCKPlus   )    fHistMassWithMcPIDOmegaPlus   ->Fill( lInvMassOmegaPlus );
2418         
2419         
2420         
2421         
2422         // - Step II.4 : MC association (care : lots of "continue;" below this line)
2423         //------------- 
2424         
2425         Bool_t lAssoXiMinus    = kFALSE;
2426         Bool_t lAssoXiPlus     = kFALSE;
2427         Bool_t lAssoOmegaMinus = kFALSE;
2428         Bool_t lAssoOmegaPlus  = kFALSE;
2429         
2430         
2431         if(fDebug > 5)
2432                 cout    << "MC EventNumber : " << lMCevent->Header()->GetEvent() 
2433                         << " / MC event Number in Run : " << lMCevent->Header()->GetEventNrInRun() << endl;
2434         
2435         // - Step 4.1 : level of the V0 daughters
2436                 
2437         Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrackXi->GetLabel() );  
2438                 // Abs value = needed ! question of quality track association ...
2439         Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrackXi->GetLabel() );
2440                 
2441         TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
2442         TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
2443         
2444
2445         // - Step 4.2 : level of the Xi daughters
2446                 
2447         Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ; 
2448         Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
2449         
2450                 if( lblMotherPosV0Dghter != lblMotherNegV0Dghter) continue; // same mother
2451                 if( lblMotherPosV0Dghter < 0 ) continue; // mother != primary (!= -1)
2452                 if( lblMotherNegV0Dghter < 0 ) continue;
2453                                         
2454
2455                 // mothers = Lambda candidate ... a priori
2456         
2457         TParticle* mcMotherPosV0Dghter = lMCstack->Particle( lblMotherPosV0Dghter );
2458         TParticle* mcMotherNegV0Dghter = lMCstack->Particle( lblMotherNegV0Dghter );
2459
2460         Int_t      lblBach  = (Int_t) TMath::Abs( bachTrackXi->GetLabel() );
2461         TParticle* mcBach   = lMCstack->Particle( lblBach );    
2462                                 
2463
2464         // - Step 4.3 : level of Xi candidate
2465         
2466         Int_t lblGdMotherPosV0Dghter =   mcMotherPosV0Dghter->GetFirstMother() ;
2467         Int_t lblGdMotherNegV0Dghter =   mcMotherNegV0Dghter->GetFirstMother() ;
2468                         
2469                 if( lblGdMotherPosV0Dghter != lblGdMotherNegV0Dghter ) continue;
2470                 if( lblGdMotherPosV0Dghter < 0 ) continue; // primary lambda ...
2471                 if( lblGdMotherNegV0Dghter < 0 ) continue; // primary lambda ...
2472                         
2473                 
2474                 // Gd mothers = Xi candidate ... a priori
2475         
2476         TParticle* mcGdMotherPosV0Dghter = lMCstack->Particle( lblGdMotherPosV0Dghter );
2477         TParticle* mcGdMotherNegV0Dghter = lMCstack->Particle( lblGdMotherNegV0Dghter );
2478                                         
2479         Int_t lblMotherBach = (Int_t) TMath::Abs( mcBach->GetFirstMother()  );
2480         
2481                 if( lblMotherBach != lblGdMotherPosV0Dghter ) continue; //same mother for bach and V0 daughters
2482         
2483         TParticle* mcMotherBach = lMCstack->Particle( lblMotherBach );
2484         
2485                 
2486         // - Step 4.4 : Manage boolean for association
2487         
2488         if( mcMotherBach                ->GetPdgCode() ==   3312 &&
2489             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3312 &&
2490             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3312)       lAssoXiMinus = kTRUE;
2491         
2492         else if( mcMotherBach           ->GetPdgCode() ==  -3312 &&
2493             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3312 &&
2494             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3312)       lAssoXiPlus = kTRUE;
2495         
2496         else if( mcMotherBach           ->GetPdgCode() ==   3334 &&
2497             mcGdMotherPosV0Dghter       ->GetPdgCode() ==   3334 &&
2498             mcGdMotherNegV0Dghter       ->GetPdgCode() ==   3334)       lAssoOmegaMinus = kTRUE;
2499                 
2500         else if( mcMotherBach           ->GetPdgCode() ==  -3334 &&
2501             mcGdMotherPosV0Dghter       ->GetPdgCode() ==  -3334 &&
2502             mcGdMotherNegV0Dghter       ->GetPdgCode() ==  -3334)       lAssoOmegaPlus = kTRUE;
2503         
2504         
2505         
2506         if(!lAssoXiMinus && !lAssoXiPlus && !lAssoOmegaMinus && !lAssoOmegaPlus) continue; // no association, skip the rest of the code 
2507         // If a proper association  exists ...
2508                 
2509         if(fDebug > 4){
2510                 cout << "XiMinus    = " << lAssoXiMinus    << endl;
2511                 cout << "XiPlus     = " << lAssoXiPlus     << endl;
2512                 cout << "OmegaMinus = " << lAssoOmegaMinus << endl;
2513                 cout << "OmegaPlus  = " << lAssoOmegaPlus  << endl 
2514                      << "----"          << endl;        
2515         }
2516
2517
2518         if(fDebug > 5){
2519                 cout << endl;
2520                 cout << "- V0 daughters - " << endl;
2521                 cout << "     + V0 Pos. / Label : " << lblPosV0Dghter 
2522                 << " - Pdg Code : " << mcPosV0Dghter->GetTitle() << endl;
2523                 cout << "     - V0 Neg. / Label : " << lblNegV0Dghter 
2524                 << " - Pdg Code : " << mcNegV0Dghter->GetTitle() << endl;
2525                 
2526                 cout << "- Xi daughters - " << endl;
2527                 cout << "     + V0 Pos. mother / Label : " << lblMotherPosV0Dghter 
2528                 << " - Pdg Code : " << mcMotherPosV0Dghter->GetTitle() << endl;
2529                 cout << "     - V0 Neg. mother / Label : " << lblMotherNegV0Dghter 
2530                 << " - Pdg Code : " << mcMotherNegV0Dghter->GetTitle() << endl;
2531                 
2532                 cout << "     --  Bach. / Label :" << lblBach 
2533                 << " -  Pdg Code : " << mcBach->GetTitle() << endl;
2534                 
2535                 cout << "- Xi candidate -" << endl;
2536                 cout << "    +  V0 Pos. Gd Mother / Label : " << lblGdMotherPosV0Dghter 
2537                 << " - Pdg Code : " << mcGdMotherPosV0Dghter->GetTitle() << endl;
2538                 cout << "    -  V0 Neg. Gd Mother / Label : "  << lblGdMotherNegV0Dghter 
2539                 << " - Pdg Code : "<< mcGdMotherNegV0Dghter->GetTitle() << endl;
2540                 
2541                 cout << "    --  Mother Bach. / Label : " << lblMotherBach 
2542                 << " - Pdg Code    : " << mcMotherBach->GetTitle() << endl;
2543                 cout << endl;
2544         }
2545
2546         
2547         // - Step 5 : Plots around the cascade candidates associated with MC
2548         //------------- 
2549         
2550         Double_t lmcPt             = mcMotherBach->Pt();
2551         Double_t lmcRapCasc        = 0.5*TMath::Log( (mcMotherBach->Energy() + mcMotherBach->Pz()) / 
2552                                                      (mcMotherBach->Energy() - mcMotherBach->Pz() +1.e-13) );
2553         Double_t lmcPtBaryon       = (mcNegV0Dghter->GetPdgCode() == -2212 ) ? mcNegV0Dghter->Pt() : mcPosV0Dghter->Pt(); // spot the baryon daughter : pbar or p
2554         
2555         
2556         Double_t lmcEta            = mcMotherBach->Eta();
2557         Double_t lmcTransvRadius   = mcBach->R(); // to get the decay point of Xi, = the production vertex of Bachelor ...
2558         
2559         TVector3 lmcTVect3Mom( mcMotherBach->Px(), mcMotherBach->Py(), mcMotherBach->Pz() );
2560
2561         Double_t lrecoPt           = xiESD->Pt();
2562         Double_t lrecoTransvRadius = TMath::Sqrt( xiESD->Xv() * xiESD->Xv() + xiESD->Yv() * xiESD->Yv() );
2563         
2564         TVector3 lrecoTVect3Mom( xiESD->Px(), xiESD->Py(), xiESD->Pz()  );
2565         Double_t lDeltaPhiMcReco   = lmcTVect3Mom.DeltaPhi( lrecoTVect3Mom ) * 180.0/TMath::Pi();
2566
2567         
2568         // - Histos for the cascade candidates associated with MC
2569         
2570         if( lChargeXi < 0 && lAssoXiMinus){     
2571                 fHistAsMCMassXiMinus          ->Fill( lInvMassXiMinus  );
2572                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiMinus->Fill( lmcPt, lmcRapCasc );
2573                 f2dHistAsMCGenPtVsGenYXiMinus ->Fill( lmcPt, lmcRapCasc);
2574                 fHistAsMCGenEtaXiMinus        ->Fill( lmcEta           );
2575                 f2dHistAsMCResPtXiMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2576                 f2dHistAsMCResRXiMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2577                 f2dHistAsMCResPhiXiMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
2578                 f2dHistAsMCPtProtonVsPtXiMinus->Fill( lmcPtBaryon, lmcPt );
2579         }
2580         
2581         else if( lChargeXi > 0 && lAssoXiPlus){ 
2582                 fHistAsMCMassXiPlus                     ->Fill( lInvMassXiPlus   );
2583                 if(lIsBachelorPion)     f2dHistAsMCandCombPIDGenPtVsGenYXiPlus->Fill( lmcPt, lmcRapCasc );
2584                 f2dHistAsMCGenPtVsGenYXiPlus            ->Fill( lmcPt, lmcRapCasc);
2585                 fHistAsMCGenEtaXiPlus                   ->Fill( lmcEta           );
2586                 f2dHistAsMCResPtXiPlus                  ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2587                 f2dHistAsMCResRXiPlus                   ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2588                 f2dHistAsMCResPhiXiPlus                 ->Fill( lmcPt, lDeltaPhiMcReco );
2589                 f2dHistAsMCPtAntiProtonVsPtXiPlus       ->Fill( lmcPtBaryon, lmcPt );
2590     
2591         }
2592         
2593         else if( lChargeXi < 0 && lAssoOmegaMinus){     
2594                 fHistAsMCMassOmegaMinus          ->Fill( lInvMassOmegaMinus );
2595                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaMinus->Fill( lmcPt, lmcRapCasc );
2596                 f2dHistAsMCGenPtVsGenYOmegaMinus ->Fill( lmcPt, lmcRapCasc  );
2597                 fHistAsMCGenEtaOmegaMinus        ->Fill( lmcEta             );
2598                 f2dHistAsMCResPtOmegaMinus       ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2599                 f2dHistAsMCResROmegaMinus        ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2600                 f2dHistAsMCResPhiOmegaMinus      ->Fill( lmcPt, lDeltaPhiMcReco );
2601                 f2dHistAsMCPtProtonVsPtOmegaMinus->Fill( lmcPtBaryon, lmcPt );
2602         }
2603         
2604         else if( lChargeXi > 0 && lAssoOmegaPlus){      
2605                 fHistAsMCMassOmegaPlus                  ->Fill( lInvMassOmegaPlus );
2606                 if(lIsBachelorKaon)     f2dHistAsMCandCombPIDGenPtVsGenYOmegaPlus->Fill( lmcPt, lmcRapCasc );
2607                 f2dHistAsMCGenPtVsGenYOmegaPlus         ->Fill( lmcPt, lmcRapCasc );
2608                 fHistAsMCGenEtaOmegaPlus                ->Fill( lmcEta            );
2609                 f2dHistAsMCResPtOmegaPlus               ->Fill( lmcPt,           (lrecoPt - lmcPt)/ lmcPt );
2610                 f2dHistAsMCResROmegaPlus                ->Fill( lmcTransvRadius, (lrecoTransvRadius - lmcTransvRadius)/ lmcTransvRadius    );
2611                 f2dHistAsMCResPhiOmegaPlus              ->Fill( lmcPt, lDeltaPhiMcReco );
2612                 f2dHistAsMCPtAntiProtonVsPtOmegaPlus    ->Fill( lmcPtBaryon, lmcPt );
2613         }
2614         
2615         
2616         // - Step 6 : Containers = Cascade cuts + PID
2617         //------------- 
2618
2619         // Double_t lChi2Xi         = -1. ;
2620         Double_t lDcaXiDaughters            = -1. ;
2621         Double_t lDcaBachToPrimVertexXi     = -1. ;
2622         Double_t lXiCosineOfPointingAngle   = -1. ;
2623         Double_t lPosXi[3]                  = { -1000.0, -1000.0, -1000.0 };
2624         Double_t lXiRadius                  = -1000. ;
2625         
2626         Double_t lInvMassLambdaAsCascDghter = 0.;
2627         Double_t lDcaV0DaughtersXi          = -1.;
2628         // Double_t lV0Chi2Xi               = -1. ;
2629         Double_t lV0CosineOfPointingAngleXi = -1.;
2630         Double_t lPosV0Xi[3]                = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
2631         Double_t lV0RadiusXi                = -1000.;
2632         Double_t lDcaV0ToPrimVertexXi       = -1.;
2633         Double_t lDcaPosToPrimVertexXi      = -1.;
2634         Double_t lDcaNegToPrimVertexXi      = -1.;
2635
2636         //Int_t    nTrackWithTPCrefitMultiplicity  =  0;
2637
2638         
2639         // 6.2 - Definition of the needed variables
2640         
2641         //lChi2Xi                       = xiESD->GetChi2Xi();
2642         lDcaXiDaughters                 = xiESD->GetDcaXiDaughters();
2643         lDcaBachToPrimVertexXi          = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0], 
2644                                                                          lBestPrimaryVtxPos[1], 
2645                                                                          lMagneticField  ) ); 
2646                                         // NOTE : AliExternalTrackParam::GetD returns an algebraic value
2647         lXiCosineOfPointingAngle        = xiESD->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
2648                                                                                   lBestPrimaryVtxPos[1],
2649                                                                                   lBestPrimaryVtxPos[2] );
2650                                         // Take care : the best available vertex should be used (like in AliCascadeVertexer)
2651                 xiESD->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] ); 
2652         lXiRadius                       = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
2653         lInvMassLambdaAsCascDghter      = xiESD->GetEffMass(); 
2654                                         // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
2655         lDcaV0DaughtersXi               = xiESD->GetDcaV0Daughters();
2656         // lV0Chi2Xi                    = xiESD->GetChi2V0();
2657         lV0CosineOfPointingAngleXi      = xiESD->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
2658                                                                              lBestPrimaryVtxPos[1],
2659                                                                              lBestPrimaryVtxPos[2] );
2660                 xiESD->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] ); 
2661         lV0RadiusXi                     = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
2662         
2663         lDcaV0ToPrimVertexXi            = xiESD->GetD( lBestPrimaryVtxPos[0], 
2664                                                 lBestPrimaryVtxPos[1], 
2665                                                 lBestPrimaryVtxPos[2] );
2666         
2667         lDcaPosToPrimVertexXi           = TMath::Abs( pTrackXi   ->GetD( lBestPrimaryVtxPos[0],
2668                                                                          lBestPrimaryVtxPos[1],
2669                                                                          lMagneticField  )     );
2670         
2671         lDcaNegToPrimVertexXi           = TMath::Abs( nTrackXi   ->GetD( lBestPrimaryVtxPos[0], 
2672                                                                          lBestPrimaryVtxPos[1],