]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/STRANGENESS/Cascades/AliAnalysisTaskExtractCascade.cxx
Logic Fix
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / Cascades / AliAnalysisTaskExtractCascade.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
17 //
18 // Modified version of AliAnalysisTaskCheckCascade.cxx.
19 // This is a 'hybrid' output version, in that it uses a classic TTree
20 // ROOT object to store the candidates, plus a couple of histograms filled on
21 // a per-event basis for storing variables too numerous to put in a tree. 
22 //
23 //
24 //  --- Algorithm Description 
25 //   1. Loop over primaries in stack to acquire generated charged Xi
26 //   2. Loop over stack to find Cascades, fill TH3Fs "PrimRawPt"s for Efficiency
27 //   3. Perform Physics Selection
28 //   4. Perform Primary Vertex |z|<10cm selection
29 //   5. Perform Primary Vertex NoTPCOnly vertexing selection
30 //   6. Perform Pileup Rejection
31 //   7. Analysis Loops: 
32 //    7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
33 //
34 //  Please Report Any Bugs! 
35 //
36 //   --- David Dobrigkeit Chinellato
37 //        (david.chinellato@gmail.com)
38 //
39 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
40
41 class TTree;
42 class TParticle;
43 class TVector3;
44
45 //class AliMCEventHandler;
46 //class AliMCEvent;
47 //class AliStack;
48
49 class AliESDVertex;
50 class AliAODVertex;
51 class AliESDv0;
52 class AliAODv0;
53
54 #include <Riostream.h>
55 #include "TList.h"
56 #include "TH1.h"
57 #include "TH2.h"
58 #include "TH3.h"
59 #include "TFile.h"
60 #include "THnSparse.h"
61 #include "TVector3.h"
62 #include "TCanvas.h"
63 #include "TMath.h"
64 #include "TLegend.h"
65 #include "AliLog.h"
66
67 #include "AliESDEvent.h"
68 #include "AliAODEvent.h"
69 #include "AliV0vertexer.h"
70 #include "AliCascadeVertexer.h"
71 #include "AliESDpid.h"
72 #include "AliESDtrack.h"
73 #include "AliESDtrackCuts.h"
74 #include "AliInputEventHandler.h"
75 #include "AliAnalysisManager.h"
76 #include "AliMCEventHandler.h"
77 #include "AliMCEvent.h"
78 #include "AliStack.h"
79
80 #include "AliV0vertexer.h"
81 #include "AliCascadeVertexer.h"
82
83 #include "AliCFContainer.h"
84 #include "AliMultiplicity.h"
85 #include "AliAODMCParticle.h"
86 #include "AliESDcascade.h"
87 #include "AliAODcascade.h"
88 #include "AliESDUtils.h"
89 #include "AliGenEventHeader.h"
90 #include "AliAnalysisUtils.h"
91
92 #include "AliAnalysisTaskExtractCascade.h"
93
94 using std::cout;
95 using std::endl;
96
97 ClassImp(AliAnalysisTaskExtractCascade)
98
99 AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade() 
100   : AliAnalysisTaskSE(), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), fUtils(0),
101    fkIsNuclear   ( kFALSE ), 
102    fkSwitchINT7  ( kFALSE ),
103    fCentralityEstimator("V0M"),
104    fkpAVertexSelection( kFALSE ),
105    fEtaRefMult ( 0.5 ),
106    fkRunVertexers ( kFALSE ),
107    fkDebugMode (kTRUE),
108    fkSelectPeripheral (kFALSE),
109 //------------------------------------------------
110 // Tree Variables
111 //------------------------------------------------
112
113    fTreeCascVarCharge(0),
114    fTreeCascVarMassAsXi(0),
115    fTreeCascVarMassAsOmega(0),
116    fTreeCascVarPt(0),
117    fTreeCascVarPtMC(0),
118    fTreeCascVarRapMC(0),
119    fTreeCascVarRapXi(0),
120    fTreeCascVarRapOmega(0),
121    fTreeCascVarNegEta(0),
122    fTreeCascVarPosEta(0),
123    fTreeCascVarBachEta(0),
124    fTreeCascVarDCACascDaughters(0),
125    fTreeCascVarDCABachToPrimVtx(0),
126    fTreeCascVarDCAV0Daughters(0),
127    fTreeCascVarDCAV0ToPrimVtx(0),
128    fTreeCascVarDCAPosToPrimVtx(0),
129    fTreeCascVarDCANegToPrimVtx(0),
130    fTreeCascVarCascCosPointingAngle(0),
131    fTreeCascVarCascRadius(0),
132    fTreeCascVarV0Mass(0),
133    fTreeCascVarV0CosPointingAngle(0),
134    fTreeCascVarV0CosPointingAngleSpecial(0),
135    fTreeCascVarV0Radius(0),
136    fTreeCascVarLeastNbrClusters(0),
137    fTreeCascVarMultiplicity(0),
138    fTreeCascVarMultiplicityV0A(0),
139    fTreeCascVarMultiplicityZNA(0),
140    fTreeCascVarMultiplicityTRK(0),
141    fTreeCascVarMultiplicitySPD(0),
142    fTreeCascVarDistOverTotMom(0),
143    fTreeCascVarPID(0),
144    fTreeCascVarPIDBachelor(0),
145    fTreeCascVarPIDNegative(0),
146    fTreeCascVarPIDPositive(0),
147    fTreeCascVarBachTransMom(0),
148    fTreeCascVarPosTransMom(0),
149    fTreeCascVarNegTransMom(0),
150    fTreeCascVarPosTransMomMC(0),
151    fTreeCascVarNegTransMomMC(0),
152    fTreeCascVarNegNSigmaPion(0),
153    fTreeCascVarNegNSigmaProton(0),
154    fTreeCascVarPosNSigmaPion(0),
155    fTreeCascVarPosNSigmaProton(0),
156    fTreeCascVarBachNSigmaPion(0),
157    fTreeCascVarBachNSigmaKaon(0),
158
159   fTreeCascVarkITSRefitBachelor(0),
160   fTreeCascVarkITSRefitNegative(0),
161   fTreeCascVarkITSRefitPositive(0),
162
163 //Debug information
164 //Part A: EbyE info, Run number
165 fTreeCascVarRunNumber(0),
166 fTreeCascVarEventNumber(0),
167
168 //Part B: Shared Clusters
169 fTreeCascVarNegClusters(0),
170 fTreeCascVarPosClusters(0),
171 fTreeCascVarBachClusters(0),
172 fTreeCascVarNegSharedClusters(0),
173 fTreeCascVarPosSharedClusters(0),
174 fTreeCascVarBachSharedClusters(0),
175
176 //Part C: All momenta
177 fTreeCascVarNegPx(0),
178 fTreeCascVarNegPy(0),
179 fTreeCascVarNegPz(0),
180 fTreeCascVarPosPx(0),
181 fTreeCascVarPosPy(0),
182 fTreeCascVarPosPz(0),
183 fTreeCascVarBachPx(0),
184 fTreeCascVarBachPy(0),
185 fTreeCascVarBachPz(0),
186
187 fTreeCascVarV0DecayX(0),
188 fTreeCascVarV0DecayY(0),
189 fTreeCascVarV0DecayZ(0),
190
191 fTreeCascVarCascadeDecayX(0),
192 fTreeCascVarCascadeDecayY(0),
193 fTreeCascVarCascadeDecayZ(0),
194
195 fTreeCascVarBadCascadeJai(0),
196 fTreeCascVarDeltaDCA(0),
197
198 //------------------------------------------------
199 // HISTOGRAMS
200 // --- Filled on an Event-by-event basis
201 //------------------------------------------------
202    fHistV0MultiplicityBeforeTrigSel(0),
203    fHistV0MultiplicityForTrigEvt(0), 
204    fHistV0MultiplicityForSelEvt(0),
205    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
206    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
207    fHistMultiplicityBeforeTrigSel(0),
208    fHistMultiplicityForTrigEvt(0),
209    fHistMultiplicity(0),
210    fHistMultiplicityNoTPCOnly(0),
211    fHistMultiplicityNoTPCOnlyNoPileup(0),
212
213 //V0A Centrality
214 fHistMultiplicityV0ABeforeTrigSel(0),
215 fHistMultiplicityV0AForTrigEvt(0),
216 fHistMultiplicityV0A(0),
217 fHistMultiplicityV0ANoTPCOnly(0),
218 fHistMultiplicityV0ANoTPCOnlyNoPileup(0),
219
220 //ZNA Centrality
221 fHistMultiplicityZNABeforeTrigSel(0),
222 fHistMultiplicityZNAForTrigEvt(0),
223 fHistMultiplicityZNA(0),
224 fHistMultiplicityZNANoTPCOnly(0),
225 fHistMultiplicityZNANoTPCOnlyNoPileup(0),
226
227 //TRK Centrality
228 fHistMultiplicityTRKBeforeTrigSel(0),
229 fHistMultiplicityTRKForTrigEvt(0),
230 fHistMultiplicityTRK(0),
231 fHistMultiplicityTRKNoTPCOnly(0),
232 fHistMultiplicityTRKNoTPCOnlyNoPileup(0),
233
234 //SPD Centrality
235 fHistMultiplicitySPDBeforeTrigSel(0),
236 fHistMultiplicitySPDForTrigEvt(0),
237 fHistMultiplicitySPD(0),
238 fHistMultiplicitySPDNoTPCOnly(0),
239 fHistMultiplicitySPDNoTPCOnlyNoPileup(0),
240
241    fHistPVx(0),
242    fHistPVy(0),
243    fHistPVz(0),
244    fHistPVxAnalysis(0),
245    fHistPVyAnalysis(0),
246    fHistPVzAnalysis(0)
247 {
248   // Dummy Constructor
249 }
250
251 AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade(const char *name) 
252   : AliAnalysisTaskSE(name), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), fUtils(0),
253    fkIsNuclear   ( kFALSE ), 
254    fkSwitchINT7  ( kFALSE ),
255    fCentralityEstimator("V0M"),
256    fkpAVertexSelection( kFALSE ),
257    fEtaRefMult ( 0.5 ),
258    fkRunVertexers ( kFALSE ),
259    fkDebugMode (kTRUE),
260    fkSelectPeripheral (kFALSE),
261 //------------------------------------------------
262 // Tree Variables
263 //------------------------------------------------
264
265    fTreeCascVarCharge(0),
266    fTreeCascVarMassAsXi(0),
267    fTreeCascVarMassAsOmega(0),
268    fTreeCascVarPt(0),
269    fTreeCascVarPtMC(0),
270    fTreeCascVarRapMC(0),
271    fTreeCascVarRapXi(0),
272    fTreeCascVarRapOmega(0),
273    fTreeCascVarNegEta(0),
274    fTreeCascVarPosEta(0),
275    fTreeCascVarBachEta(0),
276    fTreeCascVarDCACascDaughters(0),
277    fTreeCascVarDCABachToPrimVtx(0),
278    fTreeCascVarDCAV0Daughters(0),
279    fTreeCascVarDCAV0ToPrimVtx(0),
280    fTreeCascVarDCAPosToPrimVtx(0),
281    fTreeCascVarDCANegToPrimVtx(0),
282    fTreeCascVarCascCosPointingAngle(0),
283    fTreeCascVarCascRadius(0),
284    fTreeCascVarV0Mass(0),
285    fTreeCascVarV0CosPointingAngle(0),
286    fTreeCascVarV0CosPointingAngleSpecial(0),
287    fTreeCascVarV0Radius(0),
288    fTreeCascVarLeastNbrClusters(0),
289    fTreeCascVarMultiplicity(0),
290    fTreeCascVarMultiplicityV0A(0),
291    fTreeCascVarMultiplicityZNA(0),
292    fTreeCascVarMultiplicityTRK(0),
293    fTreeCascVarMultiplicitySPD(0),
294    fTreeCascVarDistOverTotMom(0),
295    fTreeCascVarPID(0),
296    fTreeCascVarPIDBachelor(0),
297    fTreeCascVarPIDNegative(0),
298    fTreeCascVarPIDPositive(0),
299    fTreeCascVarBachTransMom(0),
300    fTreeCascVarPosTransMom(0),
301    fTreeCascVarNegTransMom(0),
302    fTreeCascVarPosTransMomMC(0),
303    fTreeCascVarNegTransMomMC(0),
304    fTreeCascVarNegNSigmaPion(0),
305    fTreeCascVarNegNSigmaProton(0),
306    fTreeCascVarPosNSigmaPion(0),
307    fTreeCascVarPosNSigmaProton(0),
308    fTreeCascVarBachNSigmaPion(0),
309    fTreeCascVarBachNSigmaKaon(0),
310
311    fTreeCascVarkITSRefitBachelor(0),
312    fTreeCascVarkITSRefitNegative(0),
313    fTreeCascVarkITSRefitPositive(0),
314
315     //Debug information
316     //Part A: EbyE info, Run number
317     fTreeCascVarRunNumber(0),
318     fTreeCascVarEventNumber(0),
319
320     //Part B: Shared Clusters
321     fTreeCascVarNegClusters(0),
322     fTreeCascVarPosClusters(0),
323     fTreeCascVarBachClusters(0),
324     fTreeCascVarNegSharedClusters(0),
325     fTreeCascVarPosSharedClusters(0),
326     fTreeCascVarBachSharedClusters(0),
327
328     //Part C: All momenta
329     fTreeCascVarNegPx(0),
330     fTreeCascVarNegPy(0),
331     fTreeCascVarNegPz(0),
332     fTreeCascVarPosPx(0),
333     fTreeCascVarPosPy(0),
334     fTreeCascVarPosPz(0),
335     fTreeCascVarBachPx(0),
336     fTreeCascVarBachPy(0),
337     fTreeCascVarBachPz(0),
338
339     fTreeCascVarV0DecayX(0),
340     fTreeCascVarV0DecayY(0),
341     fTreeCascVarV0DecayZ(0),
342
343     fTreeCascVarCascadeDecayX(0),
344     fTreeCascVarCascadeDecayY(0),
345     fTreeCascVarCascadeDecayZ(0),
346     
347     fTreeCascVarBadCascadeJai(0), 
348     fTreeCascVarDeltaDCA(0),
349
350
351 //------------------------------------------------
352 // HISTOGRAMS
353 // --- Filled on an Event-by-event basis
354 //------------------------------------------------
355    fHistV0MultiplicityBeforeTrigSel(0),
356    fHistV0MultiplicityForTrigEvt(0), 
357    fHistV0MultiplicityForSelEvt(0),
358    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
359    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
360    fHistMultiplicityBeforeTrigSel(0),
361    fHistMultiplicityForTrigEvt(0),
362    fHistMultiplicity(0),
363    fHistMultiplicityNoTPCOnly(0),
364    fHistMultiplicityNoTPCOnlyNoPileup(0),
365
366 //V0A Centrality
367 fHistMultiplicityV0ABeforeTrigSel(0),
368 fHistMultiplicityV0AForTrigEvt(0),
369 fHistMultiplicityV0A(0),
370 fHistMultiplicityV0ANoTPCOnly(0),
371 fHistMultiplicityV0ANoTPCOnlyNoPileup(0),
372
373 //ZNA Centrality
374 fHistMultiplicityZNABeforeTrigSel(0),
375 fHistMultiplicityZNAForTrigEvt(0),
376 fHistMultiplicityZNA(0),
377 fHistMultiplicityZNANoTPCOnly(0),
378 fHistMultiplicityZNANoTPCOnlyNoPileup(0),
379
380 //TRK Centrality
381 fHistMultiplicityTRKBeforeTrigSel(0),
382 fHistMultiplicityTRKForTrigEvt(0),
383 fHistMultiplicityTRK(0),
384 fHistMultiplicityTRKNoTPCOnly(0),
385 fHistMultiplicityTRKNoTPCOnlyNoPileup(0),
386
387 //SPD Centrality
388 fHistMultiplicitySPDBeforeTrigSel(0),
389 fHistMultiplicitySPDForTrigEvt(0),
390 fHistMultiplicitySPD(0),
391 fHistMultiplicitySPDNoTPCOnly(0),
392 fHistMultiplicitySPDNoTPCOnlyNoPileup(0),
393
394    fHistPVx(0),
395    fHistPVy(0),
396    fHistPVz(0),
397    fHistPVxAnalysis(0),
398    fHistPVyAnalysis(0),
399    fHistPVzAnalysis(0)
400 {
401    // Constructor
402
403    //Set Variables for re-running the cascade vertexers (as done for MS paper)
404         
405         // New Loose : 1st step for the 7 TeV pp analysis
406         
407         fV0VertexerSels[0] =  33.  ;  // max allowed chi2
408         fV0VertexerSels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
409         fV0VertexerSels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
410         fV0VertexerSels[3] =   2.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
411         fV0VertexerSels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
412         fV0VertexerSels[5] =   1.0 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
413         fV0VertexerSels[6] = 200.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
414         
415         fCascadeVertexerSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
416         fCascadeVertexerSels[1] =   0.05 ;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
417         fCascadeVertexerSels[2] =   0.010;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
418         fCascadeVertexerSels[3] =   0.03 ;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
419         fCascadeVertexerSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
420         fCascadeVertexerSels[5] =   0.95 ;  // min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
421         fCascadeVertexerSels[6] =   0.4  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
422         fCascadeVertexerSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
423         
424    // Output slot #0 writes into a TList container (Cascade)
425    DefineOutput(1, TList::Class());
426    DefineOutput(2, TTree::Class());
427 }
428
429
430 AliAnalysisTaskExtractCascade::~AliAnalysisTaskExtractCascade()
431 {
432 //------------------------------------------------
433 // DESTRUCTOR
434 //------------------------------------------------
435
436    if (fListHist){
437       delete fListHist;
438       fListHist = 0x0;
439    }
440    if (fTreeCascade){
441       delete fTreeCascade;
442       fTreeCascade = 0x0;
443    }
444     //cleanup esd track cuts object too...
445    if (fESDtrackCuts){
446     delete fESDtrackCuts;
447     fESDtrackCuts = 0x0; 
448   }
449   if (fUtils){
450     delete fUtils;
451     fUtils = 0x0;
452   }
453
454 }
455
456 //________________________________________________________________________
457 void AliAnalysisTaskExtractCascade::UserCreateOutputObjects()
458 {
459    OpenFile(2); 
460    // Called once
461
462 //------------------------------------------------
463
464    fTreeCascade = new TTree("fTreeCascade","CascadeCandidates");
465
466 //------------------------------------------------
467 // fTreeCascade Branch definitions - Cascade Tree
468 //------------------------------------------------
469
470 //------------------------------------------------
471 // fTreeCascade Branch definitions
472 //------------------------------------------------
473
474 //-----------BASIC-INFO---------------------------
475 /* 1*/          fTreeCascade->Branch("fTreeCascVarCharge",&fTreeCascVarCharge,"fTreeCascVarCharge/I");  
476 /* 2*/          fTreeCascade->Branch("fTreeCascVarMassAsXi",&fTreeCascVarMassAsXi,"fTreeCascVarMassAsXi/F");
477 /* 3*/          fTreeCascade->Branch("fTreeCascVarMassAsOmega",&fTreeCascVarMassAsOmega,"fTreeCascVarMassAsOmega/F");
478 /* 4*/          fTreeCascade->Branch("fTreeCascVarPt",&fTreeCascVarPt,"fTreeCascVarPt/F");
479 /* 5*/          fTreeCascade->Branch("fTreeCascVarRapXi",&fTreeCascVarRapXi,"fTreeCascVarRapXi/F");
480 /* 6*/          fTreeCascade->Branch("fTreeCascVarRapOmega",&fTreeCascVarRapOmega,"fTreeCascVarRapOmega/F");
481 /* 7*/          fTreeCascade->Branch("fTreeCascVarNegEta",&fTreeCascVarNegEta,"fTreeCascVarNegEta/F");
482 /* 8*/          fTreeCascade->Branch("fTreeCascVarPosEta",&fTreeCascVarPosEta,"fTreeCascVarPosEta/F");
483 /* 9*/          fTreeCascade->Branch("fTreeCascVarBachEta",&fTreeCascVarBachEta,"fTreeCascVarBachEta/F");
484 //-----------INFO-FOR-CUTS------------------------
485 /*10*/          fTreeCascade->Branch("fTreeCascVarDCACascDaughters",&fTreeCascVarDCACascDaughters,"fTreeCascVarDCACascDaughters/F");
486 /*11*/          fTreeCascade->Branch("fTreeCascVarDCABachToPrimVtx",&fTreeCascVarDCABachToPrimVtx,"fTreeCascVarDCABachToPrimVtx/F");
487 /*12*/          fTreeCascade->Branch("fTreeCascVarDCAV0Daughters",&fTreeCascVarDCAV0Daughters,"fTreeCascVarDCAV0Daughters/F");
488 /*13*/          fTreeCascade->Branch("fTreeCascVarDCAV0ToPrimVtx",&fTreeCascVarDCAV0ToPrimVtx,"fTreeCascVarDCAV0ToPrimVtx/F");
489 /*14*/          fTreeCascade->Branch("fTreeCascVarDCAPosToPrimVtx",&fTreeCascVarDCAPosToPrimVtx,"fTreeCascVarDCAPosToPrimVtx/F");
490 /*15*/          fTreeCascade->Branch("fTreeCascVarDCANegToPrimVtx",&fTreeCascVarDCANegToPrimVtx,"fTreeCascVarDCANegToPrimVtx/F");
491 /*16*/          fTreeCascade->Branch("fTreeCascVarCascCosPointingAngle",&fTreeCascVarCascCosPointingAngle,"fTreeCascVarCascCosPointingAngle/F");
492 /*17*/          fTreeCascade->Branch("fTreeCascVarCascRadius",&fTreeCascVarCascRadius,"fTreeCascVarCascRadius/F");
493 /*18*/          fTreeCascade->Branch("fTreeCascVarV0Mass",&fTreeCascVarV0Mass,"fTreeCascVarV0Mass/F");
494 /*19*/          fTreeCascade->Branch("fTreeCascVarV0CosPointingAngle",&fTreeCascVarV0CosPointingAngle,"fTreeCascVarV0CosPointingAngle/F");
495 /*19*/          fTreeCascade->Branch("fTreeCascVarV0CosPointingAngleSpecial",&fTreeCascVarV0CosPointingAngleSpecial,"fTreeCascVarV0CosPointingAngleSpecial/F");
496 /*20*/          fTreeCascade->Branch("fTreeCascVarV0Radius",&fTreeCascVarV0Radius,"fTreeCascVarV0Radius/F");
497 /*21*/          fTreeCascade->Branch("fTreeCascVarLeastNbrClusters",&fTreeCascVarLeastNbrClusters,"fTreeCascVarLeastNbrClusters/I");
498 //-----------MULTIPLICITY-INFO--------------------
499 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicity",&fTreeCascVarMultiplicity,"fTreeCascVarMultiplicity/I");
500 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityV0A",&fTreeCascVarMultiplicityV0A,"fTreeCascVarMultiplicityV0A/I");
501 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityZNA",&fTreeCascVarMultiplicityZNA,"fTreeCascVarMultiplicityZNA/I");
502 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityTRK",&fTreeCascVarMultiplicityTRK,"fTreeCascVarMultiplicityTRK/I");
503 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicitySPD",&fTreeCascVarMultiplicitySPD,"fTreeCascVarMultiplicitySPD/I");
504 //-----------DECAY-LENGTH-INFO--------------------
505 /*23*/          fTreeCascade->Branch("fTreeCascVarDistOverTotMom",&fTreeCascVarDistOverTotMom,"fTreeCascVarDistOverTotMom/F");
506 //------------------------------------------------
507 /*24*/          fTreeCascade->Branch("fTreeCascVarNegNSigmaPion",&fTreeCascVarNegNSigmaPion,"fTreeCascVarNegNSigmaPion/F");
508 /*25*/          fTreeCascade->Branch("fTreeCascVarNegNSigmaProton",&fTreeCascVarNegNSigmaProton,"fTreeCascVarNegNSigmaProton/F");
509 /*26*/          fTreeCascade->Branch("fTreeCascVarPosNSigmaPion",&fTreeCascVarPosNSigmaPion,"fTreeCascVarPosNSigmaPion/F");
510 /*27*/          fTreeCascade->Branch("fTreeCascVarPosNSigmaProton",&fTreeCascVarPosNSigmaProton,"fTreeCascVarPosNSigmaProton/F");
511 /*28*/          fTreeCascade->Branch("fTreeCascVarBachNSigmaPion",&fTreeCascVarBachNSigmaPion,"fTreeCascVarBachNSigmaPion/F");
512 /*29*/          fTreeCascade->Branch("fTreeCascVarBachNSigmaKaon",&fTreeCascVarBachNSigmaKaon,"fTreeCascVarBachNSigmaKaon/F");
513     
514         //Commented out: not needed since all momenta provided! (less info)
515 /*30*/          //fTreeCascade->Branch("fTreeCascVarBachTransMom",&fTreeCascVarBachTransMom,"fTreeCascVarBachTransMom/F");
516 /*30*/          //fTreeCascade->Branch("fTreeCascVarPosTransMom",&fTreeCascVarPosTransMom,"fTreeCascVarPosTransMom/F");
517 /*31*/          //fTreeCascade->Branch("fTreeCascVarNegTransMom",&fTreeCascVarNegTransMom,"fTreeCascVarNegTransMom/F");
518
519 /*29*/          fTreeCascade->Branch("fTreeCascVarkITSRefitBachelor",&fTreeCascVarkITSRefitBachelor,"fTreeCascVarkITSRefitBachelor/O");
520 /*29*/          fTreeCascade->Branch("fTreeCascVarkITSRefitNegative",&fTreeCascVarkITSRefitNegative,"fTreeCascVarkITSRefitNegative/O");
521 /*29*/          fTreeCascade->Branch("fTreeCascVarkITSRefitPositive",&fTreeCascVarkITSRefitPositive,"fTreeCascVarkITSRefitPositive/O");
522
523     //-----------Debugging information----------------
524
525   if(fkDebugMode){ 
526       //Only save this if requested - can be turned off
527
528       //Part A: Event-by-event, run-by-run debugging
529       fTreeCascade->Branch("fTreeCascVarRunNumber",&fTreeCascVarRunNumber,"fTreeCascVarRunNumber/I");
530       fTreeCascade->Branch("fTreeCascVarEventNumber",&fTreeCascVarEventNumber,"fTreeCascVarEventNumber/l");
531       
532       //Part B: Shared Clusters for all daughter tracks
533       fTreeCascade->Branch("fTreeCascVarNegClusters",&fTreeCascVarNegClusters,"fTreeCascVarNegClusters/I");
534       fTreeCascade->Branch("fTreeCascVarPosClusters",&fTreeCascVarPosClusters,"fTreeCascVarPosClusters/I");
535       fTreeCascade->Branch("fTreeCascVarBachClusters",&fTreeCascVarBachClusters,"fTreeCascVarBachClusters/I");
536       fTreeCascade->Branch("fTreeCascVarNegSharedClusters",&fTreeCascVarNegSharedClusters,"fTreeCascVarNegSharedClusters/I");
537       fTreeCascade->Branch("fTreeCascVarPosSharedClusters",&fTreeCascVarPosSharedClusters,"fTreeCascVarPosSharedClusters/I");
538       fTreeCascade->Branch("fTreeCascVarBachSharedClusters",&fTreeCascVarBachSharedClusters,"fTreeCascVarBachSharedClusters/I");
539       
540       //Part C: All Momenta of all daughters
541       fTreeCascade->Branch("fTreeCascVarNegPx",&fTreeCascVarNegPx,"fTreeCascVarNegPx/F");
542       fTreeCascade->Branch("fTreeCascVarNegPy",&fTreeCascVarNegPy,"fTreeCascVarNegPy/F");
543       fTreeCascade->Branch("fTreeCascVarNegPz",&fTreeCascVarNegPz,"fTreeCascVarNegPz/F");
544       fTreeCascade->Branch("fTreeCascVarPosPx",&fTreeCascVarPosPx,"fTreeCascVarPosPx/F");
545       fTreeCascade->Branch("fTreeCascVarPosPy",&fTreeCascVarPosPy,"fTreeCascVarPosPy/F");
546       fTreeCascade->Branch("fTreeCascVarPosPz",&fTreeCascVarPosPz,"fTreeCascVarPosPz/F");
547       fTreeCascade->Branch("fTreeCascVarBachPx",&fTreeCascVarBachPx,"fTreeCascVarBachPx/F");
548       fTreeCascade->Branch("fTreeCascVarBachPy",&fTreeCascVarBachPy,"fTreeCascVarBachPy/F");
549       fTreeCascade->Branch("fTreeCascVarBachPz",&fTreeCascVarBachPz,"fTreeCascVarBachPz/F");
550       
551       //Part D: Decay positions
552       fTreeCascade->Branch("fTreeCascVarV0DecayX",&fTreeCascVarV0DecayX,"fTreeCascVarV0DecayX/F");
553       fTreeCascade->Branch("fTreeCascVarV0DecayY",&fTreeCascVarV0DecayY,"fTreeCascVarV0DecayY/F");
554       fTreeCascade->Branch("fTreeCascVarV0DecayZ",&fTreeCascVarV0DecayZ,"fTreeCascVarV0DecayZ/F");
555       fTreeCascade->Branch("fTreeCascVarCascadeDecayX",&fTreeCascVarCascadeDecayX,"fTreeCascVarCascadeDecayX/F");
556       fTreeCascade->Branch("fTreeCascVarCascadeDecayY",&fTreeCascVarCascadeDecayY,"fTreeCascVarCascadeDecayY/F");
557       fTreeCascade->Branch("fTreeCascVarCascadeDecayZ",&fTreeCascVarCascadeDecayZ,"fTreeCascVarCascadeDecayZ/F");
558       
559       fTreeCascade->Branch("fTreeCascVarBadCascadeJai",&fTreeCascVarBadCascadeJai,"fTreeCascVarBadCascadeJai/O");
560       fTreeCascade->Branch("fTreeCascVarDeltaDCA",&fTreeCascVarDeltaDCA,"fTreeCascVarDeltaDCA/F");
561       //------------------------------------------------
562   }
563     
564 //------------------------------------------------
565 // Particle Identification Setup
566 //------------------------------------------------
567
568    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
569    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
570    fPIDResponse = inputHandler->GetPIDResponse();
571
572 // Multiplicity 
573
574   if(! fESDtrackCuts ){
575     fESDtrackCuts = new AliESDtrackCuts();
576   }
577   if(! fUtils ){
578     fUtils = new AliAnalysisUtils();
579   }
580
581 //------------------------------------------------
582 // V0 Multiplicity Histograms
583 //------------------------------------------------
584
585    // Create histograms
586    OpenFile(1);
587    fListHist = new TList();
588    fListHist->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
589
590
591    if(! fHistV0MultiplicityBeforeTrigSel) {
592       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", 
593          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", 
594          25, 0, 25);
595       fListHist->Add(fHistV0MultiplicityBeforeTrigSel);
596    }
597            
598    if(! fHistV0MultiplicityForTrigEvt) {
599       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", 
600          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", 
601          25, 0, 25);
602       fListHist->Add(fHistV0MultiplicityForTrigEvt);
603    }
604
605    if(! fHistV0MultiplicityForSelEvt) {
606       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", 
607          "V0s per event;Nbr of V0s/Evt;Events", 
608          25, 0, 25);
609       fListHist->Add(fHistV0MultiplicityForSelEvt);
610    }
611
612    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
613       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", 
614          "V0s per event;Nbr of V0s/Evt;Events", 
615          25, 0, 25);
616       fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
617    }
618    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
619       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", 
620          "V0s per event;Nbr of V0s/Evt;Events", 
621          25, 0, 25);
622       fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
623    }
624
625 //------------------------------------------------
626 // Track Multiplicity Histograms
627 //------------------------------------------------
628
629    if(! fHistMultiplicityBeforeTrigSel) {
630       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", 
631          "Tracks per event;Nbr of Tracks;Events", 
632          200, 0, 200);          
633       fListHist->Add(fHistMultiplicityBeforeTrigSel);
634    }
635    if(! fHistMultiplicityForTrigEvt) {
636       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", 
637          "Tracks per event;Nbr of Tracks;Events", 
638          200, 0, 200);          
639       fListHist->Add(fHistMultiplicityForTrigEvt);
640    }
641    if(! fHistMultiplicity) {
642       fHistMultiplicity = new TH1F("fHistMultiplicity", 
643          "Tracks per event;Nbr of Tracks;Events", 
644          200, 0, 200);          
645       fListHist->Add(fHistMultiplicity);
646    }
647    if(! fHistMultiplicityNoTPCOnly) {
648       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", 
649          "Tracks per event;Nbr of Tracks;Events", 
650          200, 0, 200);          
651       fListHist->Add(fHistMultiplicityNoTPCOnly);
652    }
653    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
654       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", 
655          "Tracks per event;Nbr of Tracks;Events", 
656          200, 0, 200);          
657       fListHist->Add(fHistMultiplicityNoTPCOnlyNoPileup);
658    }
659   
660   
661   //V0A Centrality (if PbPb / pPb)
662   if(! fHistMultiplicityV0ABeforeTrigSel) {
663     fHistMultiplicityV0ABeforeTrigSel = new TH1F("fHistMultiplicityV0ABeforeTrigSel",
664                                                  "Centrality Distribution: V0A;V0A Centrality;Events",
665                                                  200, 0, 200);
666     fListHist->Add(fHistMultiplicityV0ABeforeTrigSel);
667   }
668   if(! fHistMultiplicityV0AForTrigEvt) {
669     fHistMultiplicityV0AForTrigEvt = new TH1F("fHistMultiplicityV0AForTrigEvt",
670                                               "Centrality Distribution: V0A;V0A Centrality;Events",
671                                               200, 0, 200);
672     fListHist->Add(fHistMultiplicityV0AForTrigEvt);
673   }
674   if(! fHistMultiplicityV0A) {
675     fHistMultiplicityV0A = new TH1F("fHistMultiplicityV0A",
676                                     "Centrality Distribution: V0A;V0A Centrality;Events",
677                                     200, 0, 200);
678     fListHist->Add(fHistMultiplicityV0A);
679   }
680   if(! fHistMultiplicityV0ANoTPCOnly) {
681     fHistMultiplicityV0ANoTPCOnly = new TH1F("fHistMultiplicityV0ANoTPCOnly",
682                                              "Centrality Distribution: V0A;V0A Centrality;Events",
683                                              200, 0, 200);
684     fListHist->Add(fHistMultiplicityV0ANoTPCOnly);
685   }
686   if(! fHistMultiplicityV0ANoTPCOnlyNoPileup) {
687     fHistMultiplicityV0ANoTPCOnlyNoPileup = new TH1F("fHistMultiplicityV0ANoTPCOnlyNoPileup",
688                                                      "Centrality Distribution: V0A;V0A Centrality;Events",
689                                                      200, 0, 200);
690     fListHist->Add(fHistMultiplicityV0ANoTPCOnlyNoPileup);
691   }
692   
693   //ZNA Centrality (if PbPb / pPb)
694   if(! fHistMultiplicityZNABeforeTrigSel) {
695     fHistMultiplicityZNABeforeTrigSel = new TH1F("fHistMultiplicityZNABeforeTrigSel",
696                                                  "Centrality Distribution: ZNA;ZNA Centrality;Events",
697                                                  200, 0, 200);
698     fListHist->Add(fHistMultiplicityZNABeforeTrigSel);
699   }
700   if(! fHistMultiplicityZNAForTrigEvt) {
701     fHistMultiplicityZNAForTrigEvt = new TH1F("fHistMultiplicityZNAForTrigEvt",
702                                               "Centrality Distribution: ZNA;ZNA Centrality;Events",
703                                               200, 0, 200);
704     fListHist->Add(fHistMultiplicityZNAForTrigEvt);
705   }
706   if(! fHistMultiplicityZNA) {
707     fHistMultiplicityZNA = new TH1F("fHistMultiplicityZNA",
708                                     "Centrality Distribution: ZNA;ZNA Centrality;Events",
709                                     200, 0, 200);
710     fListHist->Add(fHistMultiplicityZNA);
711   }
712   if(! fHistMultiplicityZNANoTPCOnly) {
713     fHistMultiplicityZNANoTPCOnly = new TH1F("fHistMultiplicityZNANoTPCOnly",
714                                              "Centrality Distribution: ZNA;ZNA Centrality;Events",
715                                              200, 0, 200);
716     fListHist->Add(fHistMultiplicityZNANoTPCOnly);
717   }
718   if(! fHistMultiplicityZNANoTPCOnlyNoPileup) {
719     fHistMultiplicityZNANoTPCOnlyNoPileup = new TH1F("fHistMultiplicityZNANoTPCOnlyNoPileup",
720                                                      "Centrality Distribution: ZNA;ZNA Centrality;Events",
721                                                      200, 0, 200);
722     fListHist->Add(fHistMultiplicityZNANoTPCOnlyNoPileup);
723   }
724   
725   //TRK Centrality (if PbPb / pPb)
726   if(! fHistMultiplicityTRKBeforeTrigSel) {
727     fHistMultiplicityTRKBeforeTrigSel = new TH1F("fHistMultiplicityTRKBeforeTrigSel",
728                                                  "Centrality Distribution: TRK;TRK Centrality;Events",
729                                                  200, 0, 200);
730     fListHist->Add(fHistMultiplicityTRKBeforeTrigSel);
731   }
732   if(! fHistMultiplicityTRKForTrigEvt) {
733     fHistMultiplicityTRKForTrigEvt = new TH1F("fHistMultiplicityTRKForTrigEvt",
734                                               "Centrality Distribution: TRK;TRK Centrality;Events",
735                                               200, 0, 200);
736     fListHist->Add(fHistMultiplicityTRKForTrigEvt);
737   }
738   if(! fHistMultiplicityTRK) {
739     fHistMultiplicityTRK = new TH1F("fHistMultiplicityTRK",
740                                     "Centrality Distribution: TRK;TRK Centrality;Events",
741                                     200, 0, 200);
742     fListHist->Add(fHistMultiplicityTRK);
743   }
744   if(! fHistMultiplicityTRKNoTPCOnly) {
745     fHistMultiplicityTRKNoTPCOnly = new TH1F("fHistMultiplicityTRKNoTPCOnly",
746                                              "Centrality Distribution: TRK;TRK Centrality;Events",
747                                              200, 0, 200);
748     fListHist->Add(fHistMultiplicityTRKNoTPCOnly);
749   }
750   if(! fHistMultiplicityTRKNoTPCOnlyNoPileup) {
751     fHistMultiplicityTRKNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityTRKNoTPCOnlyNoPileup",
752                                                      "Centrality Distribution: TRK;TRK Centrality;Events",
753                                                      200, 0, 200);
754     fListHist->Add(fHistMultiplicityTRKNoTPCOnlyNoPileup);
755   }
756   
757   //SPD Centrality (if PbPb / pPb)
758   if(! fHistMultiplicitySPDBeforeTrigSel) {
759     fHistMultiplicitySPDBeforeTrigSel = new TH1F("fHistMultiplicitySPDBeforeTrigSel",
760                                                  "Centrality Distribution: SPD;SPD Centrality;Events",
761                                                  200, 0, 200);
762     fListHist->Add(fHistMultiplicitySPDBeforeTrigSel);
763   }
764   if(! fHistMultiplicitySPDForTrigEvt) {
765     fHistMultiplicitySPDForTrigEvt = new TH1F("fHistMultiplicitySPDForTrigEvt",
766                                               "Centrality Distribution: SPD;SPD Centrality;Events",
767                                               200, 0, 200);
768     fListHist->Add(fHistMultiplicitySPDForTrigEvt);
769   }
770   if(! fHistMultiplicitySPD) {
771     fHistMultiplicitySPD = new TH1F("fHistMultiplicitySPD",
772                                     "Centrality Distribution: SPD;SPD Centrality;Events",
773                                     200, 0, 200);
774     fListHist->Add(fHistMultiplicitySPD);
775   }
776   if(! fHistMultiplicitySPDNoTPCOnly) {
777     fHistMultiplicitySPDNoTPCOnly = new TH1F("fHistMultiplicitySPDNoTPCOnly",
778                                              "Centrality Distribution: SPD;SPD Centrality;Events",
779                                              200, 0, 200);
780     fListHist->Add(fHistMultiplicitySPDNoTPCOnly);
781   }
782   if(! fHistMultiplicitySPDNoTPCOnlyNoPileup) {
783     fHistMultiplicitySPDNoTPCOnlyNoPileup = new TH1F("fHistMultiplicitySPDNoTPCOnlyNoPileup",
784                                                      "Centrality Distribution: SPD;SPD Centrality;Events",
785                                                      200, 0, 200);
786     fListHist->Add(fHistMultiplicitySPDNoTPCOnlyNoPileup);
787   }
788
789
790 //----------------------------------
791 // Primary Vertex Position Histos
792 //----------------------------------
793
794    if(! fHistPVx) {
795          fHistPVx = new TH1F("fHistPVx", 
796             "PV x position;Nbr of Evts;x", 
797             2000, -0.5, 0.5);       
798       fListHist->Add(fHistPVx);
799    }
800    if(! fHistPVy) {
801          fHistPVy = new TH1F("fHistPVy", 
802             "PV y position;Nbr of Evts;y", 
803             2000, -0.5, 0.5);       
804       fListHist->Add(fHistPVy);
805    }
806    if(! fHistPVz) {
807          fHistPVz = new TH1F("fHistPVz", 
808             "PV z position;Nbr of Evts;z", 
809             400, -20, 20);       
810       fListHist->Add(fHistPVz);
811    }
812
813    if(! fHistPVxAnalysis) {
814          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis", 
815             "PV x position;Nbr of Evts;x", 
816             2000, -0.5, 0.5);       
817       fListHist->Add(fHistPVxAnalysis);
818    }
819    if(! fHistPVyAnalysis) {
820          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis", 
821             "PV y position;Nbr of Evts;y", 
822             2000, -0.5, 0.5);       
823       fListHist->Add(fHistPVyAnalysis);
824    }
825    if(! fHistPVzAnalysis) {
826          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis", 
827             "PV z position;Nbr of Evts;z", 
828             400, -20, 20);       
829       fListHist->Add(fHistPVzAnalysis);
830    }
831
832    //List of Histograms: Normal
833    PostData(1, fListHist);
834
835    //TTree Object: Saved to base directory. Should cache to disk while saving. 
836    //(Important to avoid excessive memory usage, particularly when merging)
837    PostData(2, fTreeCascade);
838
839 }// end UserCreateOutputObjects
840
841
842 //________________________________________________________________________
843 void AliAnalysisTaskExtractCascade::UserExec(Option_t *) 
844 {
845   // Main loop
846   // Called for each event
847
848    AliESDEvent *lESDevent = 0x0;
849
850    Int_t    lNumberOfV0s                      = -1;
851    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
852    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
853    Double_t lMagneticField                 = -10.;
854    
855   // Connect to the InputEvent   
856   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
857
858    // Appropriate for ESD analysis! 
859       
860    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
861    if (!lESDevent) {
862       AliWarning("ERROR: lESDevent not available \n");
863       return;
864    }
865
866     //--- Acquisition of exact event ID
867    fTreeCascVarRunNumber = lESDevent->GetRunNumber();
868    fTreeCascVarEventNumber =
869     ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
870       ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
871         ((ULong64_t)lESDevent->GetBunchCrossNumber() )  );
872
873         
874 //------------------------------------------------
875 // Multiplicity Information Acquistion
876 //------------------------------------------------
877   
878    //REVISED multiplicity estimator after 'multiplicity day' (2011)
879    Int_t lMultiplicity = -100;
880   Int_t lMultiplicityV0A = -100;
881   Int_t lMultiplicityZNA = -100;
882   Int_t lMultiplicityTRK = -100;
883   Int_t lMultiplicitySPD = -100;
884
885    //testing purposes
886    if(fkIsNuclear == kFALSE) lMultiplicity =  fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,  fEtaRefMult );
887
888    //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
889    //---> Warning: Experimental
890    if(fkIsNuclear == kTRUE){ 
891       AliCentrality* centrality;
892       centrality = lESDevent->GetCentrality();
893       lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( fCentralityEstimator.Data() ) ) );
894       lMultiplicityV0A = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "V0A" ) ) );
895       lMultiplicityZNA = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "ZNA" ) ) );
896       lMultiplicityTRK = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "TRK" ) ) );
897       lMultiplicitySPD = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "CL1" ) ) );
898       if (centrality->GetQuality()>1) {
899         PostData(1, fListHist);
900         PostData(2, fTreeCascade);
901         return;
902       }
903    }
904     
905     if( fkSelectPeripheral ){
906         if( lMultiplicity < 60 || lMultiplicity >= 80. ){
907             //Event is outside 60-80% centrality in V0M!
908             PostData(1, fListHist);
909             PostData(2, fTreeCascade);
910             return;
911         }
912     }
913   
914    //Set variable for filling tree afterwards!
915    //---> pp case......: GetReferenceMultiplicity
916    //---> Pb-Pb case...: Centrality by V0M
917
918   fTreeCascVarMultiplicity = lMultiplicity;
919   fTreeCascVarMultiplicityV0A = lMultiplicityV0A;
920   fTreeCascVarMultiplicityZNA = lMultiplicityZNA;
921   fTreeCascVarMultiplicityTRK = lMultiplicityTRK;
922   fTreeCascVarMultiplicitySPD = lMultiplicitySPD;
923
924   fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
925   fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
926   fHistMultiplicityV0ABeforeTrigSel->Fill ( lMultiplicityV0A );
927   fHistMultiplicityZNABeforeTrigSel->Fill ( lMultiplicityZNA );
928   fHistMultiplicityTRKBeforeTrigSel->Fill ( lMultiplicityTRK );
929   fHistMultiplicitySPDBeforeTrigSel->Fill ( lMultiplicitySPD );
930   
931 //------------------------------------------------
932 // Physics Selection
933 //------------------------------------------------
934
935    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
936    Bool_t isSelected = 0;
937    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
938
939    //pA triggering: CINT7
940    if( fkSwitchINT7 ) isSelected = (maskIsSelected & AliVEvent::kINT7) == AliVEvent::kINT7;
941
942    //Standard Min-Bias Selection
943    if ( ! isSelected ) { 
944         PostData(1, fListHist);
945         PostData(2, fTreeCascade);
946       return;
947    }
948
949 //------------------------------------------------
950 // Rerun cascade vertexer! 
951 //------------------------------------------------
952
953   if( fkRunVertexers ){ 
954     lESDevent->ResetCascades();
955     lESDevent->ResetV0s();
956
957     AliV0vertexer lV0vtxer;
958     AliCascadeVertexer lCascVtxer;
959                   
960     lV0vtxer.SetDefaultCuts(fV0VertexerSels);
961     lCascVtxer.SetDefaultCuts(fCascadeVertexerSels);
962
963     lV0vtxer.Tracks2V0vertices(lESDevent);
964     lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
965   }
966 //------------------------------------------------
967 // After Trigger Selection
968 //------------------------------------------------
969
970    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
971   
972    //Set variable for filling tree afterwards!
973    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
974    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
975    fHistMultiplicityV0AForTrigEvt       ->Fill( lMultiplicityV0A  );
976    fHistMultiplicityZNAForTrigEvt       ->Fill( lMultiplicityZNA  );
977    fHistMultiplicityTRKForTrigEvt       ->Fill( lMultiplicityTRK  );
978    fHistMultiplicitySPDForTrigEvt       ->Fill( lMultiplicitySPD  );
979
980 //------------------------------------------------
981 // Getting: Primary Vertex + MagField Info
982 //------------------------------------------------
983
984    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
985    // get the vtx stored in ESD found with tracks
986    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
987         
988    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();      
989    // get the best primary vertex available for the event
990    // As done in AliCascadeVertexer, we keep the one which is the best one available.
991    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
992    // This one will be used for next calculations (DCA essentially)
993    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
994
995    Double_t lPrimaryVtxPosition[3];
996    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
997    lPrimaryVtxPosition[0] = primaryVtx->GetX();
998    lPrimaryVtxPosition[1] = primaryVtx->GetY();
999    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
1000    fHistPVx->Fill( lPrimaryVtxPosition[0] );
1001    fHistPVy->Fill( lPrimaryVtxPosition[1] );
1002    fHistPVz->Fill( lPrimaryVtxPosition[2] );
1003
1004   //------------------------------------------------
1005   // Primary Vertex Requirements Section:
1006   //  ---> pp and PbPb: Only requires |z|<10cm
1007   //  ---> pPb: all requirements checked at this stage
1008   //------------------------------------------------
1009   
1010   //Roberto's PV selection criteria, implemented 17th April 2013
1011   
1012   /* vertex selection */
1013   Bool_t fHasVertex = kFALSE;
1014   
1015   const AliESDVertex *vertex = lESDevent->GetPrimaryVertexTracks();
1016   if (vertex->GetNContributors() < 1) {
1017     vertex = lESDevent->GetPrimaryVertexSPD();
1018     if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
1019     else fHasVertex = kTRUE;
1020     TString vtxTyp = vertex->GetTitle();
1021     Double_t cov[6]={0};
1022     vertex->GetCovarianceMatrix(cov);
1023     Double_t zRes = TMath::Sqrt(cov[5]);
1024     if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fHasVertex = kFALSE;
1025   }
1026   else fHasVertex = kTRUE;
1027   
1028   //Is First event in chunk rejection: Still present!
1029   if(fkpAVertexSelection==kTRUE && fHasVertex == kFALSE) {
1030     AliWarning("Pb / | PV does not satisfy selection criteria!");
1031     PostData(1, fListHist);
1032     PostData(2, fTreeCascade);
1033     return;
1034   }
1035   
1036   //Is First event in chunk rejection: Still present!
1037   if(fkpAVertexSelection==kTRUE && fUtils->IsFirstEventInChunk(lESDevent)) {
1038     AliWarning("Pb / | This is the first event in the chunk!");
1039     PostData(1, fListHist);
1040     PostData(2, fTreeCascade);
1041     return;
1042   }
1043   
1044   //17 April Fix: Always do primary vertex Z selection, after pA vertex selection from Roberto
1045   if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0) {
1046     AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1047     PostData(1, fListHist);
1048     PostData(2, fTreeCascade);
1049     return;
1050   }
1051   
1052   
1053   lMagneticField = lESDevent->GetMagneticField( );
1054   fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
1055   fHistMultiplicity->Fill(lMultiplicity);
1056   fHistMultiplicityV0A->Fill(lMultiplicityV0A);
1057   fHistMultiplicityZNA->Fill(lMultiplicityZNA);
1058   fHistMultiplicityTRK->Fill(lMultiplicityTRK);
1059   fHistMultiplicitySPD->Fill(lMultiplicitySPD);
1060
1061 //------------------------------------------------
1062 // SKIP: Events with well-established PVtx
1063 //------------------------------------------------
1064         
1065    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
1066    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
1067    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() && fkpAVertexSelection==kFALSE ){
1068       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1069         PostData(1, fListHist);
1070         PostData(2, fTreeCascade);
1071       return;
1072    }
1073    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
1074    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
1075    fHistMultiplicityV0ANoTPCOnly->Fill(lMultiplicityV0A);
1076    fHistMultiplicityZNANoTPCOnly->Fill(lMultiplicityZNA);
1077    fHistMultiplicityTRKNoTPCOnly->Fill(lMultiplicityTRK);
1078    fHistMultiplicitySPDNoTPCOnly->Fill(lMultiplicitySPD);
1079
1080 //------------------------------------------------
1081 // Pileup Rejection Studies
1082 //------------------------------------------------
1083
1084    // FIXME : quality selection regarding pile-up rejection 
1085    if(lESDevent->IsPileupFromSPD() && !fkIsNuclear ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1086       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
1087         PostData(1, fListHist);
1088         PostData(2, fTreeCascade);
1089       return;
1090    }
1091    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
1092    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
1093    fHistMultiplicityV0ANoTPCOnlyNoPileup->Fill(lMultiplicityV0A);
1094    fHistMultiplicityZNANoTPCOnlyNoPileup->Fill(lMultiplicityZNA);
1095    fHistMultiplicityTRKNoTPCOnlyNoPileup->Fill(lMultiplicityTRK);
1096    fHistMultiplicitySPDNoTPCOnlyNoPileup->Fill(lMultiplicitySPD);
1097
1098    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
1099    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){ 
1100       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
1101       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
1102       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
1103    }
1104
1105 //------------------------------------------------
1106 // MAIN CASCADE LOOP STARTS HERE
1107 //------------------------------------------------
1108 // Code Credit: Antonin Maire (thanks^100)
1109 // ---> This is an adaptation
1110
1111   Long_t ncascades = 0;
1112         ncascades = lESDevent->GetNumberOfCascades();
1113
1114   for (Int_t iXi = 0; iXi < ncascades; iXi++){
1115     //------------------------------------------------
1116     // Initializations
1117     //------------------------------------------------  
1118           //Double_t lTrkgPrimaryVtxRadius3D = -500.0;
1119           //Double_t lBestPrimaryVtxRadius3D = -500.0;
1120           fTreeCascVarBadCascadeJai = kFALSE ;
1121
1122           // - 1st part of initialisation : variables needed to store AliESDCascade data members
1123           Double_t lEffMassXi      = 0. ;
1124           //Double_t lChi2Xi         = -1. ;
1125           Double_t lDcaXiDaughters = -1. ;
1126           Double_t lXiCosineOfPointingAngle = -1. ;
1127           Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
1128           Double_t lXiRadius = -1000. ;
1129           
1130           // - 2nd part of initialisation : Nbr of clusters within TPC for the 3 daughter cascade tracks
1131           Int_t    lPosTPCClusters    = -1; // For ESD only ...//FIXME : wait for availability in AOD
1132           Int_t    lNegTPCClusters    = -1; // For ESD only ...
1133           Int_t    lBachTPCClusters   = -1; // For ESD only ...
1134                         
1135           // - 3rd part of initialisation : about V0 part in cascades
1136           Double_t lInvMassLambdaAsCascDghter = 0.;
1137           //Double_t lV0Chi2Xi         = -1. ;
1138           Double_t lDcaV0DaughtersXi = -1.;
1139                 
1140           Double_t lDcaBachToPrimVertexXi = -1., lDcaV0ToPrimVertexXi = -1.;
1141           Double_t lDcaPosToPrimVertexXi  = -1.;
1142           Double_t lDcaNegToPrimVertexXi  = -1.;
1143           Double_t lV0CosineOfPointingAngleXi = -1. ;
1144           Double_t lV0CosineOfPointingAngleXiSpecial = -1. ;
1145           Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
1146           Double_t lV0RadiusXi = -1000.0;
1147           Double_t lV0quality  = 0.;
1148         
1149           // - 4th part of initialisation : Effective masses
1150           Double_t lInvMassXiMinus    = 0.;
1151           Double_t lInvMassXiPlus     = 0.;
1152           Double_t lInvMassOmegaMinus = 0.;
1153           Double_t lInvMassOmegaPlus  = 0.;
1154     
1155           // - 6th part of initialisation : extra info for QA
1156           Double_t lXiMomX       = 0. , lXiMomY = 0., lXiMomZ = 0.;
1157           Double_t lXiTransvMom  = 0. ;
1158           Double_t lXiTransvMomMC= 0. ;
1159           Double_t lXiTotMom     = 0. ;
1160                 
1161           Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
1162           //Double_t lBachTransvMom  = 0.;
1163           //Double_t lBachTotMom     = 0.;
1164
1165     fTreeCascVarNegNSigmaPion   = -100;
1166     fTreeCascVarNegNSigmaProton = -100;
1167     fTreeCascVarPosNSigmaPion   = -100;
1168     fTreeCascVarPosNSigmaProton = -100;
1169     fTreeCascVarBachNSigmaPion  = -100;
1170     fTreeCascVarBachNSigmaKaon  = -100;
1171         
1172           Short_t  lChargeXi = -2;
1173           //Double_t lV0toXiCosineOfPointingAngle = 0. ;
1174         
1175           Double_t lRapXi   = -20.0, lRapOmega = -20.0, lRapMC = -20.0; //  lEta = -20.0, lTheta = 360., lPhi = 720. ;
1176           //Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
1177             
1178     // -------------------------------------
1179     // II.ESD - Calculation Part dedicated to Xi vertices (ESD)
1180     
1181           AliESDcascade *xi = lESDevent->GetCascade(iXi);
1182           if (!xi) continue;
1183         
1184           
1185                   // - II.Step 1 : around primary vertex
1186                   //-------------
1187           //lTrkgPrimaryVtxRadius3D = TMath::Sqrt(  lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
1188           //                                        lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
1189           //                                        lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
1190
1191           //lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
1192           //                                        lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1193           //                                        lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1194
1195                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)        
1196                 //-------------
1197           lV0quality = 0.;
1198           xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
1199
1200           lEffMassXi                    = xi->GetEffMassXi();
1201           //lChi2Xi                         = xi->GetChi2Xi();
1202           lDcaXiDaughters       = xi->GetDcaXiDaughters();
1203           lXiCosineOfPointingAngle                  = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
1204                                                                                  lBestPrimaryVtxPos[1],
1205                                                                                  lBestPrimaryVtxPos[2] );
1206                   // Take care : the best available vertex should be used (like in AliCascadeVertexer)
1207         
1208           xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] ); 
1209           lXiRadius                     = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );           
1210
1211     fTreeCascVarCascadeDecayX = lPosXi[0];
1212     fTreeCascVarCascadeDecayY = lPosXi[1];
1213     fTreeCascVarCascadeDecayZ = lPosXi[2];
1214
1215                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
1216                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
1217                 //-------------
1218                 
1219         UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
1220         UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
1221         UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
1222                 // Care track label can be negative in MC production (linked with the track quality)
1223                 // However = normally, not the case for track index ...
1224           
1225           // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
1226           if(lBachIdx == lIdxNegXi) {
1227                   AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
1228           }
1229     if(lBachIdx == lIdxPosXi) {
1230         AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
1231           }
1232           
1233           AliESDtrack *pTrackXi         = lESDevent->GetTrack( lIdxPosXi );
1234           AliESDtrack *nTrackXi         = lESDevent->GetTrack( lIdxNegXi );
1235           AliESDtrack *bachTrackXi      = lESDevent->GetTrack( lBachIdx );
1236
1237           if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1238                   AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
1239                   continue;
1240           }
1241
1242       fTreeCascVarPosEta = pTrackXi->Eta();
1243       fTreeCascVarNegEta = nTrackXi->Eta();
1244       fTreeCascVarBachEta = bachTrackXi->Eta();
1245       
1246       //Save shared clusters information
1247       fTreeCascVarNegSharedClusters = nTrackXi->GetTPCnclsS(0,159);
1248       fTreeCascVarPosSharedClusters = pTrackXi->GetTPCnclsS(0,159);
1249       fTreeCascVarBachSharedClusters = bachTrackXi->GetTPCnclsS(0,159);
1250       
1251       Double_t lBMom[3], lNMom[3], lPMom[3];
1252       xi->GetBPxPyPz( lBMom[0], lBMom[1], lBMom[2] );
1253       xi->GetPPxPyPz( lPMom[0], lPMom[1], lPMom[2] );
1254       xi->GetNPxPyPz( lNMom[0], lNMom[1], lNMom[2] );
1255       
1256       //Save all momentum information
1257       fTreeCascVarNegPx = lNMom[0];
1258       fTreeCascVarNegPy = lNMom[1];
1259       fTreeCascVarNegPz = lNMom[2];
1260       fTreeCascVarPosPx = lPMom[0];
1261       fTreeCascVarPosPy = lPMom[1];
1262       fTreeCascVarPosPz = lPMom[2];
1263       fTreeCascVarBachPx = lBMom[0];
1264       fTreeCascVarBachPy = lBMom[1];
1265       fTreeCascVarBachPz = lBMom[2];
1266       
1267       fTreeCascVarBachTransMom = TMath::Sqrt( lBMom[0]*lBMom[0] + lBMom[1]*lBMom[1] );
1268       fTreeCascVarPosTransMom  = TMath::Sqrt( lPMom[0]*lPMom[0] + lPMom[1]*lPMom[1] );
1269       fTreeCascVarNegTransMom  = TMath::Sqrt( lNMom[0]*lNMom[0] + lNMom[1]*lNMom[1] );
1270   
1271     //------------------------------------------------
1272     // TPC dEdx information 
1273     //------------------------------------------------
1274     fTreeCascVarNegNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kPion   );
1275     fTreeCascVarNegNSigmaProton = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kProton );
1276     fTreeCascVarPosNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kPion );
1277     fTreeCascVarPosNSigmaProton = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kProton );
1278     fTreeCascVarBachNSigmaPion  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kPion );
1279     fTreeCascVarBachNSigmaKaon  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kKaon );
1280
1281     //------------------------------------------------
1282     // TPC Number of clusters info
1283     // --- modified to save the smallest number 
1284     // --- of TPC clusters for the 3 tracks
1285     //------------------------------------------------
1286               
1287           lPosTPCClusters   = pTrackXi->GetTPCNcls();
1288           lNegTPCClusters   = nTrackXi->GetTPCNcls();
1289           lBachTPCClusters  = bachTrackXi->GetTPCNcls(); 
1290           
1291       fTreeCascVarNegClusters = lNegTPCClusters;
1292       fTreeCascVarPosClusters = lPosTPCClusters;
1293       fTreeCascVarBachClusters = lBachTPCClusters;
1294
1295     // 1 - Poor quality related to TPCrefit
1296           ULong_t pStatus    = pTrackXi->GetStatus();
1297           ULong_t nStatus    = nTrackXi->GetStatus();
1298           ULong_t bachStatus = bachTrackXi->GetStatus();
1299
1300     fTreeCascVarkITSRefitBachelor = kTRUE; 
1301     fTreeCascVarkITSRefitNegative = kTRUE; 
1302     fTreeCascVarkITSRefitPositive = kTRUE; 
1303
1304     if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
1305     if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
1306     if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
1307
1308     //Extra Debug Information: booleans for ITS refit
1309     if ((pStatus&AliESDtrack::kITSrefit)    == 0) { fTreeCascVarkITSRefitPositive = kFALSE; }
1310     if ((nStatus&AliESDtrack::kITSrefit)    == 0) { fTreeCascVarkITSRefitNegative = kFALSE; }
1311     if ((bachStatus&AliESDtrack::kITSrefit) == 0) { fTreeCascVarkITSRefitBachelor = kFALSE; }
1312
1313           // 2 - Poor quality related to TPC clusters: lowest cut of 70 clusters
1314     if(lPosTPCClusters  < 70) { AliWarning("Pb / V0 Pos. track has less than 70 TPC clusters ... continue!"); continue; }
1315           if(lNegTPCClusters  < 70) { AliWarning("Pb / V0 Neg. track has less than 70 TPC clusters ... continue!"); continue; }
1316           if(lBachTPCClusters < 70) { AliWarning("Pb / Bach.   track has less than 70 TPC clusters ... continue!"); continue; }
1317           Int_t leastnumberofclusters = 1000; 
1318           if( lPosTPCClusters < leastnumberofclusters ) leastnumberofclusters = lPosTPCClusters;
1319           if( lNegTPCClusters < leastnumberofclusters ) leastnumberofclusters = lNegTPCClusters;
1320           if( lBachTPCClusters < leastnumberofclusters ) leastnumberofclusters = lBachTPCClusters;
1321
1322           lInvMassLambdaAsCascDghter    = xi->GetEffMass();
1323           // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
1324           lDcaV0DaughtersXi             = xi->GetDcaV0Daughters(); 
1325           //lV0Chi2Xi                   = xi->GetChi2V0();
1326         
1327           lV0CosineOfPointingAngleXi    = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
1328                                                                             lBestPrimaryVtxPos[1],
1329                                                                             lBestPrimaryVtxPos[2] );
1330     //Modification: V0 CosPA wrt to Cascade decay vertex
1331           lV0CosineOfPointingAngleXiSpecial     = xi->GetV0CosineOfPointingAngle( lPosXi[0],
1332                                                                             lPosXi[1],
1333                                                                             lPosXi[2] );
1334
1335           lDcaV0ToPrimVertexXi          = xi->GetD( lBestPrimaryVtxPos[0], 
1336                                                       lBestPrimaryVtxPos[1], 
1337                                                       lBestPrimaryVtxPos[2] );
1338                 
1339           lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD(       lBestPrimaryVtxPos[0], 
1340                                                                 lBestPrimaryVtxPos[1], 
1341                                                                 lMagneticField  ) ); 
1342                                           // Note : AliExternalTrackParam::GetD returns an algebraic value ...
1343                 
1344           xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] ); 
1345           lV0RadiusXi           = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1346         
1347     fTreeCascVarV0DecayX = lPosV0Xi[0];
1348     fTreeCascVarV0DecayY = lPosV0Xi[1];
1349     fTreeCascVarV0DecayZ = lPosV0Xi[2];
1350
1351           lDcaPosToPrimVertexXi         = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0], 
1352                                                                 lBestPrimaryVtxPos[1], 
1353                                                                 lMagneticField  )     ); 
1354         
1355           lDcaNegToPrimVertexXi         = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0], 
1356                                                                 lBestPrimaryVtxPos[1], 
1357                                                                 lMagneticField  )     ); 
1358                 
1359           // - II.Step 4 : around effective masses (ESD)
1360           // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1361                 
1362           if( bachTrackXi->Charge() < 0 )       {
1363                   lV0quality = 0.;
1364                   xi->ChangeMassHypothesis(lV0quality , 3312);  
1365                           // Calculate the effective mass of the Xi- candidate. 
1366                           // pdg code 3312 = Xi-
1367                   lInvMassXiMinus = xi->GetEffMassXi();
1368                 
1369                   lV0quality = 0.;
1370                   xi->ChangeMassHypothesis(lV0quality , 3334);  
1371                           // Calculate the effective mass of the Xi- candidate. 
1372                           // pdg code 3334 = Omega-
1373                   lInvMassOmegaMinus = xi->GetEffMassXi();
1374                                         
1375                   lV0quality = 0.;
1376                   xi->ChangeMassHypothesis(lV0quality , 3312);  // Back to default hyp.
1377           }// end if negative bachelor
1378         
1379         
1380           if( bachTrackXi->Charge() >  0 ){
1381                   lV0quality = 0.;
1382                   xi->ChangeMassHypothesis(lV0quality , -3312);         
1383                           // Calculate the effective mass of the Xi+ candidate. 
1384                           // pdg code -3312 = Xi+
1385                   lInvMassXiPlus = xi->GetEffMassXi();
1386                 
1387                   lV0quality = 0.;
1388                   xi->ChangeMassHypothesis(lV0quality , -3334);         
1389                           // Calculate the effective mass of the Xi+ candidate. 
1390                           // pdg code -3334  = Omega+
1391                   lInvMassOmegaPlus = xi->GetEffMassXi();
1392                 
1393                   lV0quality = 0.;
1394                   xi->ChangeMassHypothesis(lV0quality , -3312);         // Back to "default" hyp.
1395           }// end if positive bachelor
1396                   // - II.Step 6 : extra info for QA (ESD)
1397                   // miscellaneous pieces of info that may help regarding data quality assessment.
1398                   //-------------
1399
1400           xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1401                   lXiTransvMom          = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1402                   lXiTotMom     = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1403                 
1404           xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
1405                   //lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1406                   //lBachTotMom         = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1407
1408           lChargeXi = xi->Charge();
1409
1410           //lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1411         
1412           lRapXi    = xi->RapXi();
1413           lRapOmega = xi->RapOmega();
1414           //lEta      = xi->Eta();
1415           //lTheta    = xi->Theta() *180.0/TMath::Pi();
1416           //lPhi      = xi->Phi()   *180.0/TMath::Pi();
1417           //lAlphaXi  = xi->AlphaXi();
1418           //lPtArmXi  = xi->PtArmXi();
1419           
1420   //------------------------------------------------
1421   // Jai Salzwedel's femto-cut: better V0 exists
1422   //------------------------------------------------                      
1423
1424   fTreeCascVarDeltaDCA = -100;
1425   Float_t DCAV0DaughtersDiff = -100; 
1426   for (Int_t iv0=0; iv0<lESDevent->GetNumberOfV0s(); iv0++) {
1427     AliESDv0 *v0 = lESDevent->GetV0(iv0);
1428     UInt_t posV0TrackIdx = (UInt_t) v0->GetPindex();
1429     UInt_t negV0TrackIdx = (UInt_t) v0->GetNindex();
1430     if ((posV0TrackIdx == lIdxPosXi) && (negV0TrackIdx == lIdxNegXi)) continue;
1431       // if both tracks are the same ones as the cascades V0 daughter tracks, then the V0 belongs to the cascade being analysed; so avoid it
1432     if ((posV0TrackIdx == lIdxPosXi) || (negV0TrackIdx == lIdxNegXi)) {
1433     DCAV0DaughtersDiff = lDcaV0DaughtersXi - v0->GetDcaV0Daughters();
1434     if( fTreeCascVarDeltaDCA < DCAV0DaughtersDiff ) fTreeCascVarDeltaDCA = DCAV0DaughtersDiff;
1435       if ( lDcaV0DaughtersXi > v0->GetDcaV0Daughters() )  {    // DCA comparison criterion
1436         fTreeCascVarBadCascadeJai = kTRUE;
1437       } //end DCA comparison 
1438     } // end shares a daughter check 
1439   } //end V0 loop 
1440
1441   //------------------------------------------------
1442   // Set Variables for adding to tree
1443   //------------------------------------------------            
1444         
1445 /* 1*/          fTreeCascVarCharge      = lChargeXi;
1446 /* 2*/          if(lInvMassXiMinus!=0)    fTreeCascVarMassAsXi = lInvMassXiMinus;
1447 /* 2*/          if(lInvMassXiPlus!=0)     fTreeCascVarMassAsXi = lInvMassXiPlus;
1448 /* 3*/          if(lInvMassOmegaMinus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaMinus;
1449 /* 3*/          if(lInvMassOmegaPlus!=0)  fTreeCascVarMassAsOmega = lInvMassOmegaPlus;
1450 /* 4*/          fTreeCascVarPt = lXiTransvMom;
1451 /* 4*/          fTreeCascVarPtMC = lXiTransvMomMC;
1452 /* 5*/          fTreeCascVarRapXi = lRapXi ;
1453 /* 5*/          fTreeCascVarRapMC = lRapMC ;
1454 /* 6*/          fTreeCascVarRapOmega = lRapOmega ;
1455 /* 7*/          fTreeCascVarDCACascDaughters = lDcaXiDaughters;
1456 /* 8*/          fTreeCascVarDCABachToPrimVtx = lDcaBachToPrimVertexXi;
1457 /* 9*/          fTreeCascVarDCAV0Daughters = lDcaV0DaughtersXi;
1458 /*10*/          fTreeCascVarDCAV0ToPrimVtx = lDcaV0ToPrimVertexXi;
1459 /*11*/          fTreeCascVarDCAPosToPrimVtx = lDcaPosToPrimVertexXi;
1460 /*12*/          fTreeCascVarDCANegToPrimVtx = lDcaNegToPrimVertexXi;
1461 /*13*/          fTreeCascVarCascCosPointingAngle = lXiCosineOfPointingAngle;
1462 /*14*/          fTreeCascVarCascRadius = lXiRadius;
1463 /*15*/          fTreeCascVarV0Mass = lInvMassLambdaAsCascDghter;
1464 /*16*/          fTreeCascVarV0CosPointingAngle = lV0CosineOfPointingAngleXi;
1465 /*16*/          fTreeCascVarV0CosPointingAngleSpecial = lV0CosineOfPointingAngleXiSpecial;
1466 /*17*/          fTreeCascVarV0Radius = lV0RadiusXi;
1467 /*20*/          fTreeCascVarLeastNbrClusters = leastnumberofclusters;
1468 /*21*/          fTreeCascVarMultiplicity = lMultiplicity; //multiplicity, whatever that may be
1469
1470 /*23*/          fTreeCascVarDistOverTotMom = TMath::Sqrt(
1471                                                 TMath::Power( lPosXi[0] - lBestPrimaryVtxPos[0] , 2) +
1472                                                 TMath::Power( lPosXi[1] - lBestPrimaryVtxPos[1] , 2) +
1473                                                 TMath::Power( lPosXi[2] - lBestPrimaryVtxPos[2] , 2)
1474                                         );
1475 /*23*/          fTreeCascVarDistOverTotMom /= (lXiTotMom+1e-13);
1476
1477 //All vars not specified here: specified elsewhere!
1478
1479 //------------------------------------------------
1480 // Fill Tree! 
1481 //------------------------------------------------
1482
1483 // The conditional is meant to decrease excessive
1484 // memory usage! Be careful when loosening the 
1485 // cut!
1486
1487   //Xi    Mass window: 150MeV wide
1488   //Omega mass window: 150MeV wide
1489
1490   if( (fTreeCascVarMassAsXi<1.32+0.075&&fTreeCascVarMassAsXi>1.32-0.075) ||
1491       (fTreeCascVarMassAsOmega<1.68+0.075&&fTreeCascVarMassAsOmega>1.68-0.075) ){
1492       fTreeCascade->Fill();
1493   }
1494
1495 //------------------------------------------------
1496 // Fill tree over.
1497 //------------------------------------------------
1498
1499         }// end of the Cascade loop (ESD or AOD)
1500
1501    // Post output data.
1502    PostData(1, fListHist);
1503    PostData(2, fTreeCascade);
1504 }
1505
1506 //________________________________________________________________________
1507 void AliAnalysisTaskExtractCascade::Terminate(Option_t *)
1508 {
1509    // Draw result to the screen
1510    // Called once at the end of the query
1511
1512    TList *cRetrievedList = 0x0;
1513    cRetrievedList = (TList*)GetOutputData(1);
1514    if(!cRetrievedList){
1515       Printf("ERROR - AliAnalysisTaskExtractCascade : ouput data container list not available\n");
1516       return;
1517    }    
1518         
1519    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1520    if (!fHistV0MultiplicityForTrigEvt) {
1521       Printf("ERROR - AliAnalysisTaskExtractCascade : fHistV0MultiplicityForTrigEvt not available");
1522       return;
1523    }
1524   
1525    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractCascade","V0 Multiplicity",10,10,510,510);
1526    canCheck->cd(1)->SetLogy();
1527
1528    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1529    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1530 }
1531
1532 //----------------------------------------------------------------------------
1533
1534 Double_t AliAnalysisTaskExtractCascade::MyRapidity(Double_t rE, Double_t rPz) const
1535 {
1536    // Local calculation for rapidity
1537    Double_t ReturnValue = -100;
1538    if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){ 
1539       ReturnValue =  0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1540    }
1541    return ReturnValue;
1542