AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskFlowEPCascade.cxx
1 #include <iostream>
2
3 #include "TFile.h"
4 #include "TChain.h"
5 #include "TList.h"
6 #include "TH1F.h"
7 #include "TH1I.h"
8 #include "TH2F.h"
9 #include "TH3F.h"
10 #include "TProfile.h"
11 #include "TProfile2D.h"
12 #include "TF1.h"
13
14 #include "AliAnalysisTaskSE.h"
15 #include "AliAnalysisManager.h"
16
17 #include "AliESDEvent.h"
18 #include "AliESDv0.h"
19 #include "AliESDcascade.h"
20 #include "AliVVertex.h"
21
22 #include "AliESDtrack.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliESDVZERO.h"
25 #include "AliESDUtils.h"
26
27 #include "AliAODEvent.h"
28
29 #include "AliESDInputHandler.h"
30 #include "AliCentrality.h"
31 #include "AliMultiplicity.h"
32
33 #include "AliFlowTrackSimple.h"
34 //#include "AliFlowEventCuts.h"
35 #include "AliFlowTrackCuts.h"
36
37 #include "AliTPCPIDResponse.h"
38 #include "AliTOFPIDResponse.h"
39 #include "AliPIDResponse.h"
40
41 #include "AliOADBContainer.h"
42
43 #include "AliAnalysisTaskFlowEPCascade.h"
44
45 #include "TMath.h"
46
47 #include "AliLog.h"
48
49 using namespace std;
50
51 ClassImp(AliAnalysisTaskFlowEPCascade)
52
53 //_______________________________________________________________
54   AliAnalysisTaskFlowEPCascade::AliAnalysisTaskFlowEPCascade():
55     AliAnalysisTaskSE()
56     ,fCutsDau (0x0)
57     ,fPIDResponse (0x0)
58     ,fMinCent(0)
59     ,fMaxCent(0)
60     ,fVtxCut(10.)
61     ,fOADB(0x0)
62     ,fRun(-1)
63     ,fICent(-1)
64     ,fMultV0(0x0)
65     ,fV0Cpol(100)
66     ,fV0Apol(100)
67     ,fHistList       (0x0)
68     ,fhEvent (0x0)
69     ,fhEPangleVZero (0x0)
70     ,fhEPangleV0A(0x0)
71     ,fhEPangleV0C(0x0)
72     ,fhEPangleTPC(0x0)
73     ,fh1Chi2Xi(0x0)
74     ,fh1DCAXiDaughters(0x0)
75     ,fh1DCABachToPrimVertex(0x0)
76     ,fh1XiCosOfPointingAngle(0x0)
77     ,fh1XiRadius(0x0)
78     ,fh1MassLambda(0x0)
79     ,fh1V0Chi2(0x0)
80     ,fh1V0CosOfPointingAngle(0x0)
81     ,fh1V0Radius(0x0)
82     ,fh1DcaV0DaughtersXi(0x0)
83     ,fh1DcaV0ToPrimVertex(0x0)
84     ,fh1DCAPosToPrimVertex(0x0)
85     ,fh1DCANegToPrimVertex(0x0)
86     ,fh1MassXiMinus(0x0)
87     ,fh1MassXiPlus(0x0)
88     ,fh1MassOmegaMinus(0x0)
89     ,fh1MassOmegaPlus(0x0)
90     ,fh1MassXi(0x0)
91     ,fh1MassOmega(0x0)
92     ,fh1XiPt(0x0)
93     ,fh1XiP(0x0)
94     ,fh1XiBachPt(0x0)
95     ,fh1XiBachP(0x0)
96     ,fh1ChargeXi(0x0)
97     ,fh1V0toXiCosOfPointingAngle(0x0)
98     ,fh1PhiXi(0x0)
99     ,fh2Armenteros(0x0)
100     ,fh2MassLambdaVsMassXiMinus(0x0)
101     ,fh2MassXiVsMassOmegaMinus(0x0)
102     ,fh2MassLambdaVsMassXiPlus(0x0)
103     ,fh2MassXiVsMassOmegaPlus(0x0)
104     ,fh2XiRadiusVsMassXiMinus(0x0)
105     ,fh2XiRadiusVsMassXiPlus(0x0)
106     ,fh2XiRadiusVsMassOmegaMinus(0x0)
107     ,fh2XiRadiusVsMassOmegaPlus(0x0)
108     ,fh2TPCdEdxOfCascDghters(0x0)
109     ,fh2MassVsPtXiMinus(0x0)
110     ,fh2MassVsPtXiPlus(0x0)
111     ,fh2MassVsPtXiAll(0x0)
112     ,fh2MassVsPtOmegaMinus(0x0)
113     ,fh2MassVsPtOmegaPlus(0x0)
114     ,fh2MassVsPtOmegaAll(0x0)
115     ,fhXiRapidity (0x0)
116     ,fhOmegaRapidity (0x0)
117     ,fProfResolution (0x0)
118     ,fh1DistToVtxZAfter(0x0)
119     ,fh1DistToVtxXYAfter(0x0)
120     ,fh2DistToVtxZBeforeVsAfter(0x0)
121     ,fh2DistToVtxXYBeforeVsAfter(0x0)
122     ,fh2PxBeforeVsAfter(0x0)
123     ,fh2PyBeforeVsAfter(0x0)
124     ,fh2PhiPosBeforeVsAfter(0x0)
125     ,fh2PhiNegBeforeVsAfter(0x0)
126 {
127
128   for(Int_t i = 0; i < 3; i++){
129     fProfXiV2PtV0A[i] = NULL;
130     fProfOmegaV2PtV0A[i] = NULL;
131     fProfXiSinePtV0A[i] = NULL;
132     fProfOmegaSinePtV0A[i] = NULL;
133
134     fProfXiV2PtV0C[i] = NULL;
135     fProfOmegaV2PtV0C[i] = NULL;
136     fProfXiSinePtV0C[i] = NULL;
137     fProfOmegaSinePtV0C[i] = NULL;
138
139     fProfXiV2Pt[i] = NULL;
140     fProfOmegaV2Pt[i] = NULL;
141     fProfXiSinePt[i] = NULL;
142     fProfOmegaSinePt[i] = NULL;
143     
144     fProf2dXiV2PtV0A[i] = NULL;
145     fProf2dOmegaV2PtV0A[i] = NULL;
146     fProf2dXiV2PtV0C[i] = NULL;
147     fProf2dOmegaV2PtV0C[i] = NULL;
148     fProf2dXiV2Pt[i] = NULL;
149     fProf2dOmegaV2Pt[i] = NULL;
150   }
151
152   for(int i=0; i!=3; ++i)
153     for(int j=0; j!=2; ++j) {
154       fXiBands[i][j] = 0;
155       fOmegaBands[i][j] = 0;
156     }
157
158   for(Int_t i = 0; i != 2; ++i)
159     for(Int_t j = 0; j != 2; ++j)
160       for(Int_t iC = 0; iC < 9; iC++){
161         fMeanQ[iC][i][j] = 0.;
162         fWidthQ[iC][i][j] = 0.;
163     }
164
165 }
166
167 //_________________________________________________________________________
168 AliAnalysisTaskFlowEPCascade::
169 AliAnalysisTaskFlowEPCascade(const char *name, double centMin, double centMax,
170                              double xis[3][2],
171                              double omegas[3][2]) 
172   : AliAnalysisTaskSE(name)
173   ,fCutsDau (0x0)
174   ,fPIDResponse (0x0)
175   ,fMinCent(centMin)
176   ,fMaxCent(centMax)
177   ,fVtxCut(10.)
178   ,fOADB(0x0)
179   ,fRun(-1)
180   ,fICent(-1)
181   ,fMultV0(0x0)
182   ,fV0Cpol(100)
183   ,fV0Apol(100)
184   ,fHistList       (0x0)
185   ,fhEvent (0x0)
186   ,fhEPangleVZero (0x0)
187   ,fhEPangleV0A(0x0)
188   ,fhEPangleV0C(0x0)
189   ,fhEPangleTPC(0x0)
190   ,fh1Chi2Xi(0x0)
191   ,fh1DCAXiDaughters(0x0)
192   ,fh1DCABachToPrimVertex(0x0)
193   ,fh1XiCosOfPointingAngle(0x0)
194   ,fh1XiRadius(0x0)
195   ,fh1MassLambda(0x0)
196   ,fh1V0Chi2(0x0)
197   ,fh1V0CosOfPointingAngle(0x0)
198   ,fh1V0Radius(0x0)
199   ,fh1DcaV0DaughtersXi(0x0)
200   ,fh1DcaV0ToPrimVertex(0x0)
201   ,fh1DCAPosToPrimVertex(0x0)
202   ,fh1DCANegToPrimVertex(0x0)
203   ,fh1MassXiMinus(0x0)
204   ,fh1MassXiPlus(0x0)
205   ,fh1MassOmegaMinus(0x0)
206   ,fh1MassOmegaPlus(0x0)
207   ,fh1MassXi(0x0)
208   ,fh1MassOmega(0x0)
209   ,fh1XiPt(0x0)
210   ,fh1XiP(0x0)
211   ,fh1XiBachPt(0x0)
212   ,fh1XiBachP(0x0)
213   ,fh1ChargeXi(0x0)
214   ,fh1V0toXiCosOfPointingAngle(0x0)
215   ,fh1PhiXi(0x0)
216   ,fh2Armenteros(0x0)
217   ,fh2MassLambdaVsMassXiMinus(0x0)
218   ,fh2MassXiVsMassOmegaMinus(0x0)
219   ,fh2MassLambdaVsMassXiPlus(0x0)
220   ,fh2MassXiVsMassOmegaPlus(0x0)
221   ,fh2XiRadiusVsMassXiMinus(0x0)
222   ,fh2XiRadiusVsMassXiPlus(0x0)
223   ,fh2XiRadiusVsMassOmegaMinus(0x0)
224   ,fh2XiRadiusVsMassOmegaPlus(0x0)
225   ,fh2TPCdEdxOfCascDghters(0x0)
226   ,fh2MassVsPtXiMinus(0x0)
227   ,fh2MassVsPtXiPlus(0x0)
228   ,fh2MassVsPtXiAll(0x0)
229   ,fh2MassVsPtOmegaMinus(0x0)
230   ,fh2MassVsPtOmegaPlus(0x0)
231   ,fh2MassVsPtOmegaAll(0x0)
232   ,fhXiRapidity (0x0)
233   ,fhOmegaRapidity (0x0)
234   ,fProfResolution (0x0)
235   ,fh1DistToVtxZAfter(0x0)
236   ,fh1DistToVtxXYAfter(0x0)
237   ,fh2DistToVtxZBeforeVsAfter(0x0)
238   ,fh2DistToVtxXYBeforeVsAfter(0x0)
239   ,fh2PxBeforeVsAfter(0x0)
240   ,fh2PyBeforeVsAfter(0x0)
241   ,fh2PhiPosBeforeVsAfter(0x0)
242   ,fh2PhiNegBeforeVsAfter(0x0)
243 {
244
245   for(Int_t i = 0; i < 3; i++){
246     fProfXiV2PtV0A[i] = NULL;
247     fProfOmegaV2PtV0A[i] = NULL;
248     fProfXiSinePtV0A[i] = NULL;
249     fProfOmegaSinePtV0A[i] = NULL;
250
251     fProfXiV2PtV0C[i] = NULL;
252     fProfOmegaV2PtV0C[i] = NULL;
253     fProfXiSinePtV0C[i] = NULL;
254     fProfOmegaSinePtV0C[i] = NULL;
255
256     fProfXiV2Pt[i] = NULL;
257     fProfOmegaV2Pt[i] = NULL;
258     fProfXiSinePt[i] = NULL;
259     fProfOmegaSinePt[i] = NULL;
260
261     fProf2dXiV2PtV0A[i] = NULL;
262     fProf2dOmegaV2PtV0A[i] = NULL;
263     fProf2dXiV2PtV0C[i] = NULL;
264     fProf2dOmegaV2PtV0C[i] = NULL;
265     fProf2dXiV2Pt[i] = NULL;
266     fProf2dOmegaV2Pt[i] = NULL;
267
268   }
269
270   for(int i=0; i!=3; ++i)
271     for(int j=0; j!=2; ++j) {
272       fXiBands[i][j] = xis[i][j];
273       fOmegaBands[i][j] = omegas[i][j];
274     }
275
276   for(Int_t i = 0; i != 2; ++i)
277     for(Int_t j = 0; j != 2; ++j)
278       for(Int_t iC = 0; iC <9; iC++){
279         fMeanQ[iC][i][j] = 0.;
280         fWidthQ[iC][i][j] = 0.;
281       }
282   
283   DefineInput( 0,TChain::Class());
284   //DefineInput (1, TList::Class());
285   DefineOutput(1, TList::Class());
286 }
287
288 void AliAnalysisTaskFlowEPCascade::UserCreateOutputObjects()
289 {
290   cout<<"AliAnalysisTaskFlowEPCascade::UserCreateOutputObjects()"<<endl;
291   
292   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
293   AliInputEventHandler* inputHandler
294     = (AliInputEventHandler*) (man->GetInputEventHandler());
295   fPIDResponse = inputHandler->GetPIDResponse();
296
297   TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
298   fOADB = TFile::Open(oadbfilename.Data());
299
300   if(!fOADB){
301     printf("OADB file %s cannot be opened\n",oadbfilename.Data());
302     return;
303   }
304
305   fHistList = new TList();
306   fHistList ->SetOwner();
307   
308   //Add histograms to the List
309   fhEvent = new TH1I("Event", "Number of Events", 3, 0, 3);
310   fHistList->Add(fhEvent);
311   
312   fhEPangleVZero = new TH1F("hEPangleVZero", 
313                        "EP from VZERO; #Psi; Number of Events",
314                        15, 0., TMath::Pi());
315   fHistList->Add(fhEPangleVZero);
316
317   fhEPangleV0A = new TH1F("hEPangleV0A",
318                             "EP from V0A; #Psi; Number of Events",
319                             15, 0., TMath::Pi());
320   fHistList->Add(fhEPangleV0A);
321
322   fhEPangleV0C = new TH1F("hEPangleV0C",
323                             "EP from V0C; #Psi; Number of Events",
324                             15, 0., TMath::Pi());
325   fHistList->Add(fhEPangleV0C);
326   
327
328   fhEPangleTPC = new TH1F("hEPangleTPC",
329                           "EP from TPC; #Psi; Number of Events",
330                           15, 0., TMath::Pi());
331   fHistList->Add(fhEPangleTPC);
332
333   fh1Chi2Xi = new TH1F("Chi2Xi", 
334                        "Cascade #chi^{2}; #chi^{2}; Number of Cascades", 
335                        160, 0, 160);
336   fHistList->Add(fh1Chi2Xi);
337   
338   fh1DCAXiDaughters 
339     = new TH1F( "DcaXiDaughters",  
340                 "DCA between Xi Daughters; DCA (cm); Number of Cascades", 
341                 100, 0., 0.5);
342   fHistList->Add(fh1DCAXiDaughters);
343
344   fh1DCABachToPrimVertex
345     = new TH1F("DcaBachToPrimVertex", 
346                "DCA of Bach. to Prim. Vertex; DCA (cm);Number of Cascades", 
347                250, 0., 2.5);
348   fHistList->Add(fh1DCABachToPrimVertex);
349
350   fh1XiCosOfPointingAngle
351     = new TH1F("XiCosineOfPointingAngle",
352                "Cos of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis", 
353                200, 0.99, 1.0);
354   fHistList->Add(fh1XiCosOfPointingAngle);
355   
356   fh1XiRadius = new TH1F("XiRadius",  
357                          "Casc. decay transv. radius; r (cm); Counts" , 
358                          1050, 0., 105.0 );
359   fHistList->Add(fh1XiRadius);
360   
361   fh1MassLambda 
362     = new TH1F("MassLambdaAsCascDghter",
363                "#Lambda assoc. to Casc. candidates; Eff. Mass (GeV/c^{2}); Counts", 300,1.00,1.3);
364   fHistList->Add(fh1MassLambda);
365   
366   fh1V0Chi2 = new TH1F("V0Chi2Xi", 
367                        "V0 #chi^{2}, in cascade; #chi^{2};Counts", 
368                        160, 0, 40);
369   fHistList->Add(fh1V0Chi2);
370   
371   fh1V0CosOfPointingAngle 
372     = new TH1F("V0CosOfPointingAngleXi", 
373                "Cos of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts", 
374                200, 0.98, 1.0);
375   fHistList->Add(fh1V0CosOfPointingAngle);
376   
377   fh1V0Radius  = new TH1F("V0RadiusXi", 
378                           "V0 decay radius, in cascade; radius (cm); Counts", 
379                           1050, 0., 105.0);
380   fHistList->Add(fh1V0Radius);
381   
382   fh1DcaV0DaughtersXi = new TH1F("DcaV0DaughtersXi", 
383                                  "DCA between V0 daughters, in cascade;DCA (cm);Number of V0s", 120, 0., 0.6);
384   fHistList->Add(fh1DcaV0DaughtersXi);
385   
386   fh1DcaV0ToPrimVertex = new TH1F("DcaV0ToPrimVertexXi", 
387                                   "DCA of V0 to Prim. Vertex, in cascade;DCA (cm);Number of Cascades", 200, 0., 1.);
388   fHistList->Add(fh1DcaV0ToPrimVertex);
389     
390    fh1DCAPosToPrimVertex =
391      new TH1F("DcaPosToPrimVertexXi", 
392               "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts", 
393               300, 0, 3);
394    fHistList->Add(fh1DCAPosToPrimVertex);
395
396   fh1DCANegToPrimVertex 
397     =  new TH1F("DcaNegToPrimVertexXi", 
398                 "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts", 
399                 300, 0, 3);
400   fHistList->Add(fh1DCANegToPrimVertex);
401
402   fh1MassXiMinus 
403     = new TH1F("MassXiMinus", 
404                "#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts",
405                1600, 1.2, 2.0);
406   fHistList->Add(fh1MassXiMinus);
407
408
409   fh1MassXiPlus 
410     = new TH1F("MassXiPlus", 
411                "#Xi^{+} candidates;M(#bar{#Lambda}^{0}, #pi^{+}) (GeV/c^{2});Counts",
412                1600, 1.2, 2.0);
413   fHistList->Add(fh1MassXiPlus);
414
415   fh1MassXi 
416     = new TH1F("MassXi",
417                "#Xi candidates;M(#bar{#Lambda}, #pi) (GeV/c^{2});Counts",
418                1600, 1.2, 2.0);
419   fHistList->Add(fh1MassXi);
420
421   fh1MassOmegaMinus 
422     = new TH1F("MassOmegaMinus", 
423                "#Omega^{-} candidates; M(#Lambda, K^{-})(GeV/c^{2});Counts",
424                2000, 1.5, 2.5);
425   fHistList->Add(fh1MassOmegaMinus);
426
427   fh1MassOmega 
428     = new TH1F("MassOmega",
429                "#Omega candidates; M(#Lambda, K)(GeV/c^{2});Counts",
430                2000, 1.5, 2.5);
431   fHistList->Add(fh1MassOmega);
432   
433   fh1MassOmegaPlus 
434     =  new TH1F("MassOmegaPlus", 
435                 "#Omega^{+} candidates;M(#bar{#Lambda}^{0}, K^{+})(GeV/c^{2});Counts", 2000, 1.5, 2.5);
436   fHistList->Add(fh1MassOmegaPlus);
437   
438   fh1XiPt 
439     = new TH1F("XiPt" , 
440                "#Xi Pt (cand. around the mass peak);p_{t}(#Xi)(GeV/c);Counts", 
441                100, 0.0, 10.0);
442   fHistList->Add(fh1XiPt);
443
444   fh1XiP 
445     =  new TH1F("XiTotMom", 
446                 "#Xi momentum (cand. around the mass peak); p_{tot}(#Xi)(GeV/c); Counts", 
447                 150, 0.0, 15.0);
448   fHistList->Add(fh1XiP);
449
450   fh1XiBachPt = new TH1F("XiBachPt", 
451                          "#Xi Bach. transverse momentum (cand. around the mass peak) ; p_{t}(Bach.) (GeV/c); Counts", 
452                          100, 0.0, 5.0);
453   fHistList->Add(fh1XiBachPt);
454
455   fh1XiBachP = new TH1F("BachTotMomXi", 
456                         "#Xi Bach. momentum (cand. around the mass peak); p_{tot}(Bach.) (GeV/c); Counts", 100, 0.0, 5.0);
457   fHistList->Add(fh1XiBachP);
458   
459   fh1ChargeXi = new TH1F("ChargeXi", 
460                          "Charge of casc. candidates; Sign; Counts", 
461                          5, -2.0, 3.0);
462   fHistList->Add(fh1ChargeXi);
463   
464   fh1V0toXiCosOfPointingAngle 
465     = new TH1F("V0toXiCosineOfPointingAngle", 
466                "Cos. of V0 Ptng Angl Xi vtx; Cos(V0 Point. Angl / Xi vtx); Counts", 
467                100, 0.99, 1.0);
468   fHistList->Add(fh1V0toXiCosOfPointingAngle);
469   
470   fh1PhiXi 
471     = new TH1F("PhiXi", 
472                "#phi of #Xi candidates (around the mass peak); #phi (deg); Counts", 64, 0., 6.4);
473   fHistList->Add(fh1PhiXi);
474
475   fh2Armenteros 
476     = new TH2F("Armenteros", 
477                "#alpha_{Arm}(casc. cand.) Vs Pt_{Arm}(casc. cand.); #alpha_{Arm} ; Pt_{Arm} (GeV/c)", 140, -1.2, 1.2, 300, 0., 0.3);
478   fHistList->Add(fh2Armenteros);
479
480    fh2XiRadiusVsMassXiMinus 
481      = new TH2F( "XiRadiusVsEffMassXiMinus", 
482                  "Transv. R_{Xi Decay} Vs M_{#Xi^{-} candidates}; r_{cascade} (cm); M( #Lambda , #pi^{-} ) (GeV/c^{2}) ", 
483                  450, 0., 45.0, 1600, 1.2, 2.0);
484   fHistList->Add(fh2XiRadiusVsMassXiMinus);
485
486   fh2XiRadiusVsMassXiPlus 
487     = new TH2F( "XiRadiusVsEffMassXiPlus",
488                 "Transv. R_{Xi Decay} Vs M_{#Xi^{+} candidates}; r_{cascade} (cm); M( #Lambda , #pi^{+} ) (GeV/c^{2}) ",
489                 450, 0., 45.0, 1600, 1.2, 2.0);
490   fHistList->Add(fh2XiRadiusVsMassXiPlus);
491
492   fh2XiRadiusVsMassOmegaMinus
493     = new TH2F( "XiRadiusVsEffMassOmegaMinus", 
494                 "Transv. R_{Xi Decay} Vs M_{#Omega^{-} candidates}; r_{cascade} (cm); M( #Lambda , K^{-} ) (GeV/c^{2})",  
495                 450, 0., 45.0, 2000, 1.5, 2.5);
496   fHistList->Add(fh2XiRadiusVsMassOmegaMinus);
497
498   fh2XiRadiusVsMassOmegaPlus 
499     = new TH2F( "XiRadiusVsEffMassOmegaPlus",
500                 "Transv. R_{Xi Decay} Vs M_{#Omega^{+} candidates}; r_{cascade} (cm); M( #Lambda , K^{+} ) (GeV/c^{2}) ",
501                 450, 0., 45.0, 2000, 1.5, 2.5);
502   fHistList->Add(fh2XiRadiusVsMassOmegaPlus);
503
504   fh2MassLambdaVsMassXiMinus 
505     = new TH2F( "f2dHistEffMassLambdaVsEffMassXiMinus", 
506                 "M_{#Lambda} Vs M_{#Xi^{-} candidates} ; Inv. M_{#Lambda^{0}}\
507                  (GeV/c^{2}); M(#Lambda, #pi^{-}) (GeV/c^{2})", 
508                 300, 1.1, 1.13, 1600, 1.2, 2.0);
509   fHistList->Add(fh2MassLambdaVsMassXiMinus);
510
511   fh2MassXiVsMassOmegaMinus 
512     = new TH2F( "EffMassXiVsEffMassOmegaMinus", 
513                 "M_{#Xi^{-} candidates} Vs M_{#Omega^{-} candidates} ; M( #Lambda , #pi^{-} ) (GeV/c^{2}) ; M( #Lambda , K^{-} ) (GeV/c^{2})", 
514                 1600, 1.2, 2.0, 2000, 1.5, 2.5);
515   fHistList->Add(fh2MassXiVsMassOmegaMinus);
516
517   fh2MassLambdaVsMassXiPlus
518     = new TH2F("EffMassLambdaVsEffMassXiPlus", 
519                "M_{#Lambda} Vs M_{#Xi^{+} candidates}; Inv. M_{#Lambda^{0}}(GeV/c^{2}); M( #Lambda , #pi^{+} ) (GeV/c^{2})", 
520                300, 1.1,1.13, 1600, 1.2, 2.0);
521   fHistList->Add(fh2MassLambdaVsMassXiPlus);
522
523   fh2MassXiVsMassOmegaPlus 
524     = new TH2F("EffMassXiVsEffMassOmegaPlus", 
525                "M_{#Xi^{+} candidates} Vs M_{#Omega^{+} candidates} ; M( #Lambda, #pi^{+} ) (GeV/c^{2}) ; M( #Lambda , K^{+} ) (GeV/c^{2})", 
526                1600, 1.2, 2.0, 2000, 1.5, 2.5);
527   fHistList->Add(fh2MassXiVsMassOmegaPlus);
528
529   fh2TPCdEdxOfCascDghters
530     = new TH2F( "TPCdEdxOfCascDghters",
531                 "TPC dE/dx of the cascade daughters; charge x || #vec{p}_{TPC inner wall}(Casc. daughter) || (GeV/c); TPC signal (ADC) ", 
532                 2000, -10.0, 10.0, 450, 0., 900.);
533   fHistList->Add(fh2TPCdEdxOfCascDghters);
534
535   fh2MassVsPtXiMinus 
536     = new TH2F("MassVsPtXiMinus", 
537                "M_{#Xi^{-} candidates} vs Pt; Pt (GeV/c); M(#Lambda, #pi^{-}) (GeV/c^{2})",   
538                100, 0., 10., 1600, 1.2, 2.0);
539   fHistList->Add(fh2MassVsPtXiMinus);
540
541   fh2MassVsPtXiPlus
542     = new TH2F("MassVsPtXiPlus",
543                "M_{#Xi^{+} candidates} vs Pt; Pt (GeV/c); M(#Lambda, #pi^{+})(GeV/c^{2})",
544                100, 0., 10., 1600, 1.2, 2.0);
545   fHistList->Add(fh2MassVsPtXiPlus);
546
547   fh2MassVsPtXiAll 
548     = new TH2F("MassVsPtXiAll", 
549                "M_{#Xi candidates} vs Pt; Pt (GeV/c); M(#Lambda, #pi) (GeV/c^{2})", 
550                100, 0., 10., 1600, 1.2, 2.0);
551   fHistList->Add(fh2MassVsPtXiAll);
552
553   fh2MassVsPtOmegaMinus 
554     = new TH2F("MassVsPtOmegaMinus", 
555                "M_{#Omega^{-} candidates} vs Pt; Pt (GeV/c); M(#Lambda, K^{-}) (GeV/c^{2})", 
556                100, 0., 10., 2000, 1.5, 2.5);
557   fHistList->Add(fh2MassVsPtOmegaMinus);
558
559   fh2MassVsPtOmegaPlus 
560     = new TH2F("MassVsPtOmegaPlus",
561                "M_{#Omega^{+} candidates} vs Pt; Pt (GeV/c); M(#Lambda, K^{+}) (GeV/c^{2})", 
562                100, 0., 10., 2000, 1.5, 2.5);
563   fHistList->Add(fh2MassVsPtOmegaPlus);
564
565   fh2MassVsPtOmegaAll
566     = new TH2F("MassVsPtOmegaAll",
567                "M_{#Omega candidates} vs Pt; Pt (GeV/c); M(#Lambda, K^{+}) (GeV/c^{2})",
568                100, 0., 10., 2000, 1.5, 2.5);
569   fHistList->Add(fh2MassVsPtOmegaAll);
570
571   fhXiRapidity = new TH1F("hXiRapidity", 
572                           "#Xi rapidity distribution before rap. cut; y; Number of counts",
573                           20, -1., 1.);
574   fHistList->Add(fhXiRapidity);
575   
576   fhOmegaRapidity = new TH1F("hOmegaRapidity",
577                              "#Omega rapidity distr. before rap. cut; y; Number of counts", 
578                              20, -1., 1.);
579   fHistList->Add(fhOmegaRapidity);
580
581   for(Int_t i = 0; i < 3; i++){
582     fProfXiV2PtV0A[i] = new TProfile(Form("hProfXiV2PtV0A%d", i), 
583                                      "; p_{T}[GeV/c]; v_{2}", 
584                                      100, 0., 10.);
585     fHistList->Add(fProfXiV2PtV0A[i]);
586
587     fProfOmegaV2PtV0A[i] = new TProfile(Form("hProfOmegaV2PtV0A%d", i),
588                                         "; p_{T}[GeV/c]; v_{2}",
589                                         100, 0., 10.);
590     fHistList->Add(fProfOmegaV2PtV0A[i]);
591
592     fProfXiSinePtV0A[i] = new TProfile(Form("hProfXiSinePtV0A%d", i),
593                                        ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
594                                        100, 0., 10.);
595     fHistList->Add(fProfXiSinePtV0A[i]);
596
597     fProfOmegaSinePtV0A[i] 
598       = new TProfile(Form("hProfOmegaSinePtV0A%d", i),
599                      "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
600                      100, 0., 10.);
601     fHistList->Add(fProfOmegaSinePtV0A[i]);
602     
603
604     fProfXiV2PtV0C[i] = new TProfile(Form("hProfXiV2PtV0C%d", i), 
605                                      "; p_{T}[GeV/c]; v_{2}", 
606                                      100, 0., 10.);
607     fHistList->Add(fProfXiV2PtV0C[i]);
608
609     fProfOmegaV2PtV0C[i] = new TProfile(Form("hProfOmegaV2PtV0C%d", i),
610                                         "; p_{T}[GeV/c]; v_{2}",
611                                         100, 0., 10.);
612     fHistList->Add(fProfOmegaV2PtV0C[i]);
613
614     fProfXiSinePtV0C[i] = new TProfile(Form("hProfXiSinePtV0C%d", i),
615                                        ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
616                                        100, 0., 10.);
617     fHistList->Add(fProfXiSinePtV0C[i]);
618
619     fProfOmegaSinePtV0C[i] 
620       = new TProfile(Form("hProfOmegaSinePtV0C%d", i),
621                      "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
622                      100, 0., 10.);
623     fHistList->Add(fProfOmegaSinePtV0C[i]);
624
625
626     fProfXiV2Pt[i] = new TProfile(Form("hProfXiV2Pt%d", i), 
627                                   "; p_{T}[GeV/c]; v_{2}", 
628                                   100, 0., 10.);
629     fHistList->Add(fProfXiV2Pt[i]);
630
631     fProfOmegaV2Pt[i] = new TProfile(Form("hProfOmegaV2Pt%d", i),
632                                   "; p_{T}[GeV/c]; v_{2}",
633                                   100, 0., 10.);
634     fHistList->Add(fProfOmegaV2Pt[i]);
635
636     fProfXiSinePt[i] = new TProfile(Form("hProfXiSinePt%d", i),
637                                   ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
638                                   100, 0., 10.);
639     fHistList->Add(fProfXiSinePt[i]);
640
641     fProfOmegaSinePt[i] = new TProfile(Form("hProfOmegaSinePt%d", i),
642                                      "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
643                                      100, 0., 10.);
644     fHistList->Add(fProfOmegaSinePt[i]);
645
646     fProf2dXiV2PtV0A[i] = new TProfile2D(Form("hProf2dXiV2PtV0A%d", i), 
647                                          "; p_{T}[GeV/c]; #Psi; v_{2}",
648                                          20, 0., 10., 
649                                          15, 0., TMath::Pi());
650     fHistList->Add(fProf2dXiV2PtV0A[i]);
651
652     fProf2dOmegaV2PtV0A[i] = new TProfile2D(Form("hProf2dOmegaV2PtV0A%d", i),
653                                             "; p_{T}[GeV/c]; #Psi; v_{2}",
654                                             20, 0., 10.,
655                                             15, 0., 
656                                             TMath::Pi());
657     fHistList->Add(fProf2dOmegaV2PtV0A[i]);
658
659     fProf2dXiV2PtV0C[i] = new TProfile2D(Form("hProf2dXiV2PtV0C%d", i),
660                                          "; p_{T}[GeV/c]; #Psi; v_{2}",
661                                          20, 0., 10.,
662                                          15, 0., TMath::Pi());
663     fHistList->Add(fProf2dXiV2PtV0C[i]);
664   
665     fProf2dOmegaV2PtV0C[i] = new TProfile2D(Form("hProf2dOmegaV2PtV0C%d", i),
666                                             "; p_{T}[GeV/c]; #Psi; v_{2}",
667                                             20, 0., 10.,
668                                             15, 0.,
669                                             TMath::Pi());
670     fHistList->Add(fProf2dOmegaV2PtV0C[i]);
671     
672     fProf2dXiV2Pt[i] = new TProfile2D(Form("hProf2dXiV2Pt%d", i),
673                                       "; p_{T}[GeV/c]; #Psi; v_{2}",
674                                       20, 0., 10.,
675                                       15, 0., TMath::Pi());
676     fHistList->Add(fProf2dXiV2Pt[i]);
677
678     fProf2dOmegaV2Pt[i] = new TProfile2D(Form("hProf2dOmegaV2Pt%d", i),
679                                          "; p_{T}[GeV/c]; #Psi; v_{2}",
680                                          20, 0., 10.,
681                                          15, 0.,
682                                          TMath::Pi());
683     fHistList->Add(fProf2dOmegaV2Pt[i]);
684   }
685   
686   fh1DistToVtxZAfter = new TH1F("DistToVtxZAfter",
687                                 "Distance to vtx z after propagation to vtx; z [cm]",
688                                 100, -5., 5.);
689   fHistList->Add(fh1DistToVtxZAfter);
690   
691   fh1DistToVtxXYAfter = new TH1F("DistToVtxXYAfter",
692                                  "Distance to vtx xy after propagation to vtx",
693                                  500, 0., 50.);
694   fHistList->Add(fh1DistToVtxXYAfter);
695   
696   fh2DistToVtxZBeforeVsAfter
697     = new TH2F("DistToVtxZBeforeVsAfter",
698                "Distance to vtx z before vs after propagation; Distance before [cm]; Distance after [cm]",
699                500, -50., 50., 100, -5., 5.);
700   fHistList->Add(fh2DistToVtxZBeforeVsAfter);
701   
702   fh2DistToVtxXYBeforeVsAfter
703     = new TH2F("DistToVtxXYBeforeVsAfter",
704                "Distance to vtx xy before vs after propagation; Distance before [cm]; Distance after [cm]",
705                500, 0., 50, 500, 0., 50);
706   fHistList->Add(fh2DistToVtxXYBeforeVsAfter);
707   
708   fh2PxBeforeVsAfter
709     = new TH2F("PxBeforeVsAfter",
710                "Px before vs after propagation; Px [GeV/c]; Px' [GeV/c]",
711                200, -10., 10, 200, -10., 10.);
712   fHistList->Add(fh2PxBeforeVsAfter);
713   
714   fh2PyBeforeVsAfter
715     = new TH2F("PyBeforeVsAfter",
716                "Py before vs after propagation; Py [GeV/c]; Py' [GeV/c]",
717                200, -10., 10, 200, -10., 10.);
718   fHistList->Add(fh2PyBeforeVsAfter);
719   
720   fh2PhiPosBeforeVsAfter
721     = new TH2F("PhiPosBeforeVsAfter",
722                "Phi for positively charged candidates before vs after propagation; #phi; #phi'",
723                360, 0., 2.0*TMath::Pi(), 360, 0., 2.0*TMath::Pi());
724   fHistList->Add(fh2PhiPosBeforeVsAfter);
725   
726   fh2PhiNegBeforeVsAfter
727     = new TH2F("PhiNegBeforeVsAfter",
728                "Phi for negatively charged candidates before vs after propagation; #phi; #phi'",
729                360, 0., 2.0*TMath::Pi(), 360, 0., 2.0*TMath::Pi());
730   fHistList->Add(fh2PhiNegBeforeVsAfter);
731   
732   fProfResolution = new TProfile("hProfResolution", 
733                                  "correlations between event planes", 
734                                  4, 0.5, 4.5);
735   (fProfResolution->GetXaxis())
736     ->SetBinLabel(1, "<cos[2(#Psi^{V0A})-#Psi^{V0C}]>");
737   (fProfResolution->GetXaxis())
738     ->SetBinLabel(2, "<cos[2(#Psi^{V0A})-#Psi^{TPC}]>");
739   (fProfResolution->GetXaxis())
740     ->SetBinLabel(3, "<cos[2(#Psi^{V0C})-#Psi^{TPC}]>");
741   (fProfResolution->GetXaxis())
742     ->SetBinLabel(4, "<cos[2(#Psi^{ZDCA})-#Psi^{ZDCC}]>");
743   fHistList->Add(fProfResolution);
744   
745   PostData(1, fHistList);
746 }
747
748 void AliAnalysisTaskFlowEPCascade::UserExec(Option_t *) 
749 {
750
751   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
752   AliAODEvent *fAOD=dynamic_cast<AliAODEvent*>(InputEvent());
753
754   if(fESD){
755     
756     Info("UserExec", "This task doesn't work with ESD!");
757     //fhEvent->Fill(0);
758     
759     // if (!fFlowEventCuts->IsSelected(InputEvent())) return;
760     
761     // fhEvent->Fill(2);
762     //    ReadFromESDv0(fESD);
763   }
764   else if(fAOD){
765     //routine for AOD info
766
767     fhEvent->Fill(0);
768
769     //At the momment the cutting class does not handle AOD event properly
770     //so we are doing the cuts explicitly here
771     AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
772     if(!aodHeader) AliFatal("Not a standard AOD");
773     if(!aodHeader) return;
774     AliCentrality *centrality = aodHeader->GetCentralityP();
775     if(!centrality) return;
776     Double_t cent = centrality->GetCentralityPercentile("V0M" );
777     if(cent<fMinCent||cent>=fMaxCent) return; //centrality cut
778
779     if(cent >= 0 && cent < 5) fICent = 0;
780     else if(cent >= 5 && cent < 10) fICent = 1;
781     else if(cent >= 10 && cent < 20) fICent = 2;
782     else if(cent >= 20 && cent < 30) fICent = 3;
783     else if(cent >= 30 && cent < 40) fICent = 4;
784     else if(cent >= 40 && cent < 50) fICent = 5;
785     else if(cent >= 50 && cent < 60) fICent = 6;
786     else if(cent >= 60 && cent < 70) fICent = 7;
787     else if(cent >= 70 && cent < 80) fICent = 8;
788     else 
789       return;
790
791     Double_t zvtx = fAOD->GetPrimaryVertex()->GetZ();
792     if(TMath::Abs(zvtx) > fVtxCut) return; //vertex cut
793     
794     fhEvent->Fill(2);
795     ReadFromAODv0(fAOD);
796   }
797   
798   return;
799 }
800
801 /*
802 void AliAnalysisTaskFlowEPCascade::ReadFromESDv0(AliESDEvent *fESD){
803
804   //  fEPcalib->CalculateEventPlanes(fESD);
805
806   //Get EP informations
807   //Double_t psiVZero = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0AC);
808   //fhEPangleVZero->Fill(psiVZero);
809
810   //Double_t psiV0A = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0A);
811   //Double_t psiV0C = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0C);
812   //Double_t psiTPC = fEPcalib->GetPsi(2, AliAnaEPcalib::kTPC);
813
814   //  Double_t psiZDC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCAC);
815   //fhEPangleZDC->Fill(psiZDC);
816   
817   //Double_t psiZDCA = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCA);
818   //Double_t psiZDCC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCC);
819
820   AliEventplane * ep = fESD->GetEventplane();
821   Double_t psiTPC = ep->GetEventplane("Q", fESD, 2);
822   
823   Int_t run = fESD->GetRunNumber();  
824   if(run != fRun){
825     // Load the calibrations run dependent
826     OpenInfoCalbration(run);
827     fRun=run;
828   }
829
830   //fill profile for resolution estimation
831   fProfResolution->Fill(1, TMath::Cos(2.*(psiV0A - psiV0C)));
832   fProfResolution->Fill(2, TMath::Cos(2.*(psiV0A - psiTPC)));
833   fProfResolution->Fill(3, TMath::Cos(2.*(psiV0C - psiTPC)));
834   //fProfResolution->Fill(4, TMath::Cos(2.*(psiZDCA - psiZDCC)));
835
836   //check Xi candidates
837   Double_t trkPrimaryVtxPos[3] = {-100., -100., -100.};
838   Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
839   int nCascades=fESD->GetNumberOfCascades();
840   
841   const AliESDVertex *primaryTrackingESDVtx = fESD->GetPrimaryVertexTracks();
842   primaryTrackingESDVtx->GetXYZ(trkPrimaryVtxPos);
843   
844   const AliESDVertex *primaryBestESDVtx = fESD->GetPrimaryVertex();
845   primaryBestESDVtx->GetXYZ(bestPrimaryVtxPos);
846
847   Double_t b = fESD->GetMagneticField();
848   
849   for(int i = 0; i != nCascades; ++i) {
850     Double_t effMassXi = 0.;
851     Double_t chi2Xi = -1.;
852     Double_t dcaXiDaughters = -1.;
853     Double_t XiCosOfPointingAngle = -1.;
854     Double_t posXi[3] = {-1000., -1000., -1000.};
855     Double_t XiRadius = -1000.;
856     
857     Double_t invMassLambdaAsCascDghter = 0.;
858     Double_t V0Chi2Xi = -1.;
859     Double_t dcaV0DaughtersXi = -1.;
860     
861     Double_t dcaBachToPrimaryVtxXi = -1.;
862     Double_t dcaV0ToPrimaryVtxXi = -1.;
863     Double_t dcaPosToPrimaryVtxXi = -1.;
864     Double_t dcaNegToPrimaryVtxXi = -1.;
865     Double_t V0CosOfPointingAngleXi = -1.;
866     Double_t posV0Xi[3] = {-1000., -1000., -1000.};
867     Double_t V0RadiusXi = -1000.;
868     Double_t V0quality = 0.;
869
870     Double_t invMassXiMinus = 0.;
871     Double_t invMassXiPlus = 0.;
872     Double_t invMassOmegaMinus = 0.;
873     Double_t invMassOmegaPlus = 0.;
874     
875     Bool_t isPosInXiProton = kFALSE;
876     Bool_t isPosInXiPion = kFALSE;
877     Bool_t isPosInOmegaProton = kFALSE;
878     Bool_t isPosInOmegaPion = kFALSE;
879     
880     Bool_t isNegInXiProton = kFALSE;
881     Bool_t isNegInXiPion = kFALSE;
882     Bool_t isNegInOmegaProton = kFALSE;
883     Bool_t isNegInOmegaPion = kFALSE;
884
885     Bool_t isBachelorKaon = kFALSE;
886     Bool_t isBachelorPion = kFALSE;
887     
888     Bool_t isBachelorKaonForTPC = kFALSE;
889     Bool_t isBachelorPionForTPC = kFALSE;
890     Bool_t isNegPionForTPC = kFALSE;
891     Bool_t isPosPionForTPC = kFALSE;
892     Bool_t isNegProtonForTPC = kFALSE;
893     Bool_t isPosProtonForTPC = kFALSE;
894
895     Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
896     Double_t XiPt = 0.;
897     Double_t XiPtot = 0.;
898
899     Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
900     Double_t bachPt = 0.;
901     Double_t bachPtot = 0.;
902     
903     //Short_t chargeXi = -2;
904     Double_t V0toXiCosOfPointingAngle = 0.;
905     
906     Double_t rapXi = -20.;
907     Double_t rapOmega = -20.;
908     Double_t phi = 6.3;
909     Double_t alphaXi = -200.;
910     Double_t ptArmXi = -200.;
911
912     Double_t distToVtxZBefore = -999.;
913     Double_t distToVtxZAfter = -999.;
914     Double_t distToVtxXYBefore = -999.;
915     Double_t distToVtxXYAfter = -999.;
916     Double_t XiPAfter[3] = {-999., -999., -999.};
917     Double_t phiAfter = -999.;
918
919     AliESDcascade *xi = fESD->GetCascade(i);
920     if(!xi) continue;
921     
922     if(xi->Charge()<0)
923       xi->ChangeMassHypothesis(V0quality, 3312); // Xi- hypothesis
924     else if(xi->Charge() > 0)
925       xi->ChangeMassHypothesis(V0quality, -3312);
926     else continue;
927
928     effMassXi = xi->GetEffMassXi();
929     chi2Xi = xi->GetChi2Xi();
930     dcaXiDaughters = xi->GetDcaXiDaughters();
931     XiCosOfPointingAngle 
932       = xi->GetCascadeCosineOfPointingAngle(bestPrimaryVtxPos[0],
933                                             bestPrimaryVtxPos[1],
934                                             bestPrimaryVtxPos[2]);
935     xi->GetXYZcascade(posXi[0], posXi[1], posXi[2]);
936     XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
937                            +posXi[1]*posXi[1]
938                            +posXi[2]*posXi[2]);
939     
940     UInt_t idxPosXi = (UInt_t)TMath::Abs(xi->GetPindex());
941     UInt_t idxNegXi = (UInt_t)TMath::Abs(xi->GetNindex());
942     UInt_t idxBach = (UInt_t)TMath::Abs(xi->GetBindex());
943
944     if(idxBach == idxPosXi || idxBach == idxNegXi) continue;
945
946     AliESDtrack *pTrkXi = fESD->GetTrack(idxPosXi);
947     AliESDtrack *nTrkXi = fESD->GetTrack(idxNegXi);
948     AliESDtrack *bTrkXi = fESD->GetTrack(idxBach);
949
950     if( !pTrkXi || !nTrkXi || !bTrkXi ) continue;
951
952     if( !fCutsDau->IsSelected(pTrkXi) 
953         || !fCutsDau->IsSelected(nTrkXi)
954         || !fCutsDau->IsSelected(bTrkXi) ) continue;
955
956     const AliExternalTrackParam *pExtTrk = pTrkXi->GetInnerParam();
957     const AliExternalTrackParam *nExtTrk = nTrkXi->GetInnerParam();
958     const AliExternalTrackParam *bExtTrk = bTrkXi->GetInnerParam();
959     
960     if(pExtTrk && pTrkXi->IsOn(AliESDtrack::kTPCin)){
961       fh2TPCdEdxOfCascDghters->Fill(pExtTrk->GetP()*pExtTrk->Charge(), 
962                                     pTrkXi->GetTPCsignal());
963     }
964     if(nExtTrk && nTrkXi->IsOn(AliESDtrack::kTPCin)){
965       fh2TPCdEdxOfCascDghters->Fill(nExtTrk->GetP()*nExtTrk->Charge(),
966                                     nTrkXi->GetTPCsignal());
967     }
968     if(bExtTrk && bTrkXi->IsOn(AliESDtrack::kTPCin)){
969       fh2TPCdEdxOfCascDghters->Fill(bExtTrk->GetP()*bExtTrk->Charge(),
970                                     bTrkXi->GetTPCsignal());
971     }
972     
973
974     invMassLambdaAsCascDghter = xi->GetEffMass(); // from V0
975     dcaV0DaughtersXi = xi->GetDcaV0Daughters();
976     V0Chi2Xi = xi->GetChi2V0();
977     
978     V0CosOfPointingAngleXi 
979       = xi->GetV0CosineOfPointingAngle(bestPrimaryVtxPos[0],
980                                        bestPrimaryVtxPos[1],
981                                        bestPrimaryVtxPos[2]);
982     dcaV0ToPrimaryVtxXi = xi->GetD(bestPrimaryVtxPos[0], 
983                                    bestPrimaryVtxPos[1],
984                                    bestPrimaryVtxPos[2]);
985     dcaBachToPrimaryVtxXi = TMath::Abs(bTrkXi->GetD(bestPrimaryVtxPos[0],
986                                                     bestPrimaryVtxPos[1],
987                                                     bestPrimaryVtxPos[2]));
988     
989     //V0
990     xi->GetXYZ(posV0Xi[0], posV0Xi[1], posV0Xi[2]);
991     V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
992                              +posV0Xi[1]*posV0Xi[1]
993                              +posV0Xi[2]*posV0Xi[2]);
994     dcaPosToPrimaryVtxXi = TMath::Abs(pTrkXi->GetD(bestPrimaryVtxPos[0],
995                                                    bestPrimaryVtxPos[1],
996                                                    bestPrimaryVtxPos[2]));
997     dcaNegToPrimaryVtxXi = TMath::Abs(nTrkXi->GetD(bestPrimaryVtxPos[0],
998                                                    bestPrimaryVtxPos[1],
999                                                    bestPrimaryVtxPos[2]));
1000     
1001     //apply cuts
1002     //    if(XiRadius < 1.5 || XiRadius > 100.) continue;
1003     //if(dcaXiDaughters > 0.3) continue;
1004     //if(XiCosOfPointingAngle < 0.999) continue;
1005     //if(dcaV0ToPrimaryVtxXi < 0.03) continue;
1006     //if(dcaBachToPrimaryVtxXi < 0.03) continue;
1007     
1008     //V0 mass cut?
1009     //if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.006) continue;
1010
1011     //if(dcaV0DaughtersXi > .6) continue;
1012     //if(V0CosOfPointingAngleXi > 0.9999) continue;
1013     //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1014     //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1015
1016     //if(V0RadiusXi < 6.0 || V0RadiusXi > 100) continue;
1017     
1018     //other cuts?
1019
1020     // change mass hypothesis to cover all the possibilities
1021     if(bTrkXi->Charge()<0){
1022       V0quality = 0.;
1023       xi->ChangeMassHypothesis(V0quality, 3312); //Xi- hyp.
1024       invMassXiMinus = xi->GetEffMassXi();
1025
1026       V0quality = 0.;
1027       xi->ChangeMassHypothesis(V0quality, 3334); //Omega- hyp.
1028       invMassOmegaMinus = xi->GetEffMassXi();
1029
1030       V0quality = 0.;
1031       xi->ChangeMassHypothesis(V0quality, 3312); //back to default hyp.
1032     }
1033
1034     if(bTrkXi->Charge() > 0){
1035       V0quality = 0.;
1036       xi->ChangeMassHypothesis(V0quality, -3312); //anti-Xi- hyp.
1037       invMassXiPlus = xi->GetEffMassXi();
1038
1039       V0quality = 0.;
1040       xi->ChangeMassHypothesis(V0quality, -3334); //anti-Omega- hyp.
1041       invMassOmegaPlus = xi->GetEffMassXi();
1042
1043       V0quality = 0.;
1044       xi->ChangeMassHypothesis(V0quality, -3312); //back to default hyp.
1045     }
1046
1047     //PID on the daughter tracks
1048     //A - Combined PID
1049     //Resonable guess the priors for the cascade track sample
1050     //(e, mu, pi, K, p)
1051     Double_t priorsGuessXi[5] = {0, 0, 2, 0, 1};
1052     Double_t priorsGuessOmega[5] = {0, 0, 1, 1, 1};
1053
1054     //Combined bachelor-daughter PID
1055     AliPID pidXi;
1056     pidXi.SetPriors(priorsGuessXi);
1057     AliPID pidOmega;
1058     pidOmega.SetPriors(priorsGuessOmega);
1059     
1060     if(pTrkXi->IsOn(AliESDtrack::kESDpid)){// combined PID exists
1061       Double_t r[10] = {0.};
1062       pTrkXi->GetESDpid(r);
1063       pidXi.SetProbabilities(r);
1064       pidOmega.SetProbabilities(r);
1065
1066       //Check if the V0 postive track is proton (case for Xi-)
1067       Double_t pProton = pidXi.GetProbability(AliPID::kProton);
1068       if(pProton > pidXi.GetProbability(AliPID::kElectron)
1069          && pProton > pidXi.GetProbability(AliPID::kMuon)
1070          && pProton > pidXi.GetProbability(AliPID::kPion)
1071          && pProton > pidXi.GetProbability(AliPID::kKaon))
1072         isPosInXiProton = kTRUE;
1073       
1074       //Check if the V0 postive track is a pi+ (case for Xi+)
1075       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1076       if(pPion > pidXi.GetProbability(AliPID::kElectron)
1077          && pPion > pidXi.GetProbability(AliPID::kMuon)
1078          && pPion > pidXi.GetProbability(AliPID::kKaon)
1079          && pPion > pidXi.GetProbability(AliPID::kProton))
1080         isPosInXiPion = kTRUE;
1081       // Check if the V0 positive track is a proton (case for Omega-)
1082       pProton = pidOmega.GetProbability(AliPID::kProton);
1083       if(pProton > pidOmega.GetProbability(AliPID::kElectron)
1084          && pProton > pidOmega.GetProbability(AliPID::kMuon)
1085          && pProton > pidOmega.GetProbability(AliPID::kPion)
1086          && pProton > pidOmega.GetProbability(AliPID::kKaon))
1087         isPosInOmegaProton = kTRUE;
1088     
1089       // Check if the V0 positive track is a pi+ (case for Omega+)
1090       pPion =  pidOmega.GetProbability(AliPID::kPion);
1091       if(pPion > pidOmega.GetProbability(AliPID::kElectron)
1092          && pPion > pidOmega.GetProbability(AliPID::kMuon)
1093          && pPion > pidOmega.GetProbability(AliPID::kKaon)
1094          && pPion > pidOmega.GetProbability(AliPID::kProton))
1095         isPosInOmegaPion = kTRUE;
1096     }
1097     
1098     //Combined V0-negative-daughter PID
1099     pidXi.SetPriors(priorsGuessXi);
1100     pidOmega.SetPriors(priorsGuessOmega);
1101     if(nTrkXi->IsOn(AliESDtrack::kESDpid)){
1102       Double_t r[10] = {0.};
1103       nTrkXi->GetESDpid(r);
1104       pidXi.SetProbabilities(r);
1105       pidOmega.SetProbabilities(r);
1106       
1107       // Check if the V0 negative track is a pi- (case for Xi-)
1108       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1109       if(pPion > pidXi.GetProbability(AliPID::kElectron)
1110          && pPion > pidXi.GetProbability(AliPID::kMuon)
1111          && pPion >  pidXi.GetProbability(AliPID::kKaon)
1112          && pPion > pidXi.GetProbability(AliPID::kProton))
1113         isNegInXiPion = kTRUE;
1114
1115       // Check if the V0 negative track is an anti-proton (case for Xi+)
1116       Double_t pProton = pidXi.GetProbability(AliPID::kProton);
1117       if(pProton > pidXi.GetProbability(AliPID::kElectron)
1118          && pProton > pidXi.GetProbability(AliPID::kMuon)
1119          && pProton > pidXi.GetProbability(AliPID::kPion) 
1120          && pProton > pidXi.GetProbability(AliPID::kKaon))
1121         isNegInXiProton = kTRUE;
1122       
1123       // Check if the V0 negative track is a pi- (case for Omega-)
1124       pPion = pidOmega.GetProbability(AliPID::kPion);
1125       if(pPion > pidOmega.GetProbability(AliPID::kElectron)
1126          && pPion > pidOmega.GetProbability(AliPID::kMuon)
1127          && pPion > pidOmega.GetProbability(AliPID::kKaon)
1128          && pPion > pidOmega.GetProbability(AliPID::kProton))
1129         isNegInOmegaPion = kTRUE;
1130       
1131       // Check if the V0 negative track is an anti-proton (case for Omega+)   
1132       pProton =  pidOmega.GetProbability(AliPID::kProton);
1133       if(pProton > pidOmega.GetProbability(AliPID::kElectron)
1134          && pProton > pidOmega.GetProbability(AliPID::kMuon)
1135          && pProton > pidOmega.GetProbability(AliPID::kPion)
1136          && pProton > pidOmega.GetProbability(AliPID::kKaon))
1137         isNegInOmegaProton = kTRUE;
1138     }
1139
1140     // Combined bachelor PID
1141     pidXi.SetPriors(priorsGuessXi);
1142     pidOmega.SetPriors(priorsGuessOmega);
1143     if(bTrkXi->IsOn(AliESDtrack::kESDpid)){//Combined PID exists
1144       Double_t r[10] = {0.};
1145       bTrkXi->GetESDpid(r);
1146       pidXi.SetProbabilities(r);
1147       pidOmega.SetProbabilities(r);
1148       
1149       //Check if the bachelor track is a pion
1150       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1151       if(pPion >  pidXi.GetProbability(AliPID::kElectron)
1152          && pPion >  pidXi.GetProbability(AliPID::kMuon)
1153          && pPion > pidXi.GetProbability(AliPID::kKaon)
1154          && pPion > pidXi.GetProbability(AliPID::kProton))
1155         isBachelorPion = kTRUE;
1156      
1157       // Check if the bachelor track is a kaon
1158       Double_t pKaon = pidOmega.GetProbability(AliPID::kKaon);
1159       if(pKaon > pidOmega.GetProbability(AliPID::kElectron)
1160          && pKaon > pidOmega.GetProbability(AliPID::kMuon)
1161          && pKaon > pidOmega.GetProbability(AliPID::kPion)
1162          && pKaon > pidOmega.GetProbability(AliPID::kProton))
1163         isBachelorKaon = kTRUE;
1164     }
1165
1166     //B - TPC PID: 3-sigma bands on Bethe-Bloch curve
1167     //Bachelor
1168     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
1169       isBachelorKaonForTPC = kTRUE;
1170     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
1171       isBachelorPionForTPC = kTRUE;
1172
1173     //Negative V0 daughter
1174     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
1175       isNegPionForTPC = kTRUE;
1176     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
1177       isNegProtonForTPC = kTRUE;
1178     
1179     //Positive V0 daughter
1180     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
1181       isPosPionForTPC = kTRUE;
1182     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
1183       isPosProtonForTPC = kTRUE;
1184
1185     //    //Extra QA information
1186     xi->GetPxPyPz(XiPx, XiPy, XiPz);
1187     XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
1188     XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
1189     
1190     xi->GetBPxPyPz(bachPx, bachPy, bachPz);
1191     bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1192     bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1193
1194      V0toXiCosOfPointingAngle 
1195        = xi->GetV0CosineOfPointingAngle(posXi[0], posXi[1], posXi[2]);
1196     rapXi = xi->RapXi();
1197     rapOmega = xi->RapOmega();
1198     phi = xi->Phi();
1199     alphaXi = xi->AlphaXi();
1200     ptArmXi = xi->PtArmXi();
1201
1202     XiPAfter[0] = XiPx;
1203     XiPAfter[1] = XiPy;
1204     XiPAfter[2] = XiPz;
1205
1206     distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1207     distToVtxXYBefore
1208       = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1209                     *(posXi[0] - bestPrimaryVtxPos[0])
1210                     +(posXi[1] - bestPrimaryVtxPos[1])
1211                     *(posXi[1] - bestPrimaryVtxPos[1]));
1212
1213     //propagation to the best primary vertex to determine the momentum   
1214     Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->Charge());
1215     distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
1216     distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1217                                    *(posXi[0] - bestPrimaryVtxPos[0])
1218                                    +(posXi[1] - bestPrimaryVtxPos[1])
1219                                    *(posXi[1] - bestPrimaryVtxPos[1]));
1220     phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
1221
1222     fh1DistToVtxZAfter->Fill(distToVtxZAfter);
1223     fh1DistToVtxXYAfter->Fill(distToVtxXYAfter);
1224     fh2DistToVtxZBeforeVsAfter->Fill(distToVtxZBefore, distToVtxZAfter);
1225     fh2DistToVtxXYBeforeVsAfter->Fill(distToVtxXYBefore, distToVtxXYAfter);
1226     fh2PxBeforeVsAfter->Fill(XiPx, XiPAfter[0]);
1227     fh2PyBeforeVsAfter->Fill(XiPy, XiPAfter[1]);
1228     if(xi->Charge()>0)
1229       fh2PhiPosBeforeVsAfter->Fill(phi, phiAfter);
1230     else if(xi->Charge()<0)
1231       fh2PhiNegBeforeVsAfter->Fill(phi, phiAfter);
1232
1233
1234     if( isBachelorPion && 
1235         ( (xi->Charge()<0 && isPosInXiProton && isNegInXiPion)
1236           || (xi->Charge()>0 && isNegInXiProton && isPosInXiPion ) ) )
1237       {//Xi candidate
1238         //for default hypothesis
1239         fh1Chi2Xi->Fill(chi2Xi);
1240         fh1DCAXiDaughters->Fill(dcaXiDaughters);
1241         fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1242         fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1243         fh1XiRadius->Fill(XiRadius);
1244         
1245         //V0
1246         fh1MassLambda->Fill(invMassLambdaAsCascDghter);
1247         fh1V0Chi2->Fill(V0Chi2Xi);
1248         fh1DcaV0DaughtersXi->Fill(dcaV0DaughtersXi);
1249         fh1V0CosOfPointingAngle->Fill(V0CosOfPointingAngleXi);
1250         fh1V0Radius->Fill(V0RadiusXi);
1251         fh1DcaV0ToPrimVertex->Fill(dcaV0ToPrimaryVtxXi);
1252         fh1DCAPosToPrimVertex->Fill(dcaPosToPrimaryVtxXi);
1253         fh1DCANegToPrimVertex->Fill(dcaNegToPrimaryVtxXi);
1254         fh1ChargeXi->Fill(xi->Charge());
1255         fh1V0toXiCosOfPointingAngle->Fill(V0toXiCosOfPointingAngle);
1256        
1257         if ( TMath::Abs( invMassXiMinus-1.3217 ) < 0.012 
1258              || TMath::Abs( invMassXiPlus-1.3217 ) < 0.012) 
1259           {// One InvMass should be different from 0             
1260             fh1XiPt->Fill(XiPt);
1261             fh1XiP->Fill(XiPtot);
1262             
1263             fh1XiBachPt->Fill(bachPt);
1264             fh1XiBachP->Fill(bachPtot);
1265             fh1PhiXi->Fill( xi->Phi() );
1266           }
1267         fh2Armenteros->Fill(alphaXi, ptArmXi);
1268       }
1269
1270     if(xi->Charge()<0){
1271       fh1MassXiMinus->Fill(invMassXiMinus);
1272       fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1273       fh1MassXi->Fill(invMassXiMinus);
1274       fh1MassOmega->Fill(invMassOmegaMinus);
1275       
1276       fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1277                                        invMassXiMinus);
1278       fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus, 
1279                                       invMassOmegaMinus);
1280       fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1281       fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1282     }
1283
1284     if(xi->Charge() > 0){
1285       fh1MassXiPlus->Fill(invMassXiPlus);
1286       fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1287       fh1MassXi->Fill(invMassXiPlus);
1288       fh1MassOmega->Fill(invMassOmegaPlus);
1289       
1290       fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter, 
1291                                       invMassXiPlus);
1292       fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1293                                      invMassOmegaPlus);
1294       fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1295       fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1296     }
1297     
1298
1299     Double_t phiNew 
1300       = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. ) ? 
1301       0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1302     Double_t phiV0 = phiNew;
1303     phiV0 -= psiVZero;
1304     //    if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1305     //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1306
1307     Double_t phiV0A = phiNew;
1308     phiV0A -= psiV0A;
1309     //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1310     //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1311
1312     Double_t phiV0C = phiNew;
1313     phiV0C -= psiV0C;
1314     //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1315     //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1316
1317     //PID cuts with TPC cuts
1318     if(xi->Charge() < 0){
1319       if(isPosProtonForTPC
1320          && isNegPionForTPC){
1321         if(isBachelorPionForTPC){
1322           //Xi
1323           fhXiRapidity->Fill(rapXi);
1324           if(TMath::Abs(rapXi) < 0.8){
1325             fh2MassVsPtXiMinus->Fill(XiPt, invMassXiMinus);
1326             fh2MassVsPtXiAll->Fill(XiPt, invMassXiMinus);
1327             
1328             for(int r=0; r!=3; ++r) {
1329               if(invMassXiMinus > fXiBands[r][0] 
1330                  && invMassXiMinus < fXiBands[r][1]){
1331
1332                 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1333                 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1334
1335                 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1336                 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1337
1338                 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1339                 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1340               }
1341             }
1342               
1343           }
1344         }
1345     
1346         if(isBachelorKaonForTPC){
1347           //Omega
1348           fhOmegaRapidity->Fill(rapOmega);
1349           if(TMath::Abs(rapOmega) < 0.8){
1350             fh2MassVsPtOmegaMinus->Fill(XiPt, invMassOmegaMinus);
1351             fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaMinus);
1352
1353             for(int r=0; r!=3; ++r) {
1354               if(invMassOmegaMinus > fOmegaBands[r][0]
1355                  && invMassOmegaMinus < fOmegaBands[r][1]){
1356
1357                 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1358                 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1359                 
1360                 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1361                 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1362
1363                 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1364                 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1365               }
1366             }
1367             
1368
1369           }
1370         }
1371       }
1372     }
1373   
1374     if(xi->Charge() > 0){
1375       if(isNegProtonForTPC
1376          && isPosPionForTPC){
1377         if(isBachelorPionForTPC){
1378           //Xi
1379           fhXiRapidity->Fill(rapXi);
1380           if(TMath::Abs(rapXi) < 0.8){
1381             fh2MassVsPtXiPlus->Fill(XiPt, invMassXiPlus);
1382             fh2MassVsPtXiAll->Fill(XiPt, invMassXiPlus);
1383             for(int r=0; r!=3; ++r) {
1384               if(invMassXiPlus > fXiBands[r][0]
1385                  && invMassXiPlus < fXiBands[r][1]){
1386                 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1387                 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1388
1389                 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1390                 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1391
1392                 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1393                 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1394               }
1395             }
1396
1397           }
1398         }
1399       
1400         if(isBachelorKaonForTPC){
1401           //Omega
1402           fhOmegaRapidity->Fill(rapOmega);
1403           if(TMath::Abs(rapOmega) < 0.8){
1404             fh2MassVsPtOmegaPlus->Fill(XiPt, invMassOmegaPlus);
1405             fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaPlus);
1406             
1407             for(int r=0; r!=3; ++r) {
1408               if(invMassOmegaPlus > fOmegaBands[r][0]
1409                  && invMassOmegaPlus < fOmegaBands[r][1]){
1410                 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1411                 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1412
1413                 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1414                 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1415                 
1416                 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1417                 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1418               }
1419             }
1420
1421
1422           }
1423         }
1424       }
1425     }
1426   
1427   }//end of cascade loop
1428
1429   //  fNEvent++;
1430
1431   PostData(1, fHistList);
1432 }
1433 */
1434
1435 void AliAnalysisTaskFlowEPCascade::ReadFromAODv0(AliAODEvent *fAOD){
1436
1437   AliEventplane * ep = ((AliVAODHeader*)fAOD->GetHeader())->GetEventplaneP();
1438   Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
1439   //  if(psiTPC > TMath::PiOver2()) 
1440   //  psiTPC -= TMath::Pi();
1441   fhEPangleTPC->Fill(psiTPC);
1442
1443   Int_t run = fAOD->GetRunNumber();
1444   if(run != fRun){
1445     // Load the calibrations run dependent
1446     OpenInfoCalbration(run);                    
1447     fRun=run;          
1448   }
1449
1450   //reset Q vector info       
1451   Double_t Qxa2 = 0, Qya2 = 0;
1452   Double_t Qxc2 = 0, Qyc2 = 0;
1453   //Double_t Qx2 = 0, Qy2 = 0;
1454
1455   //V0 info   
1456   AliAODVZERO* aodV0 = fAOD->GetVZEROData();
1457
1458   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1459     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1460     Float_t multv0 = aodV0->GetMultiplicity(iv0);
1461     if (iv0 < 32){ // V0C 
1462       Qxc2 += TMath::Cos(2*phiV0)*multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1463       Qyc2 += TMath::Sin(2*phiV0)*multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1464     } else { // V0A
1465       Qxa2 += TMath::Cos(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1466       Qya2 += TMath::Sin(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1467     }
1468   }
1469   
1470   // Qx2 = Qxa2 + Qxc2;
1471   //Qy2 = Qya2 + Qya2;
1472
1473   //grab for each centrality the proper histo with the Qx and Qy 
1474   //to do the recentering                    
1475   Double_t Qxamean2 = fMeanQ[fICent][1][0];
1476   Double_t Qxarms2  = fWidthQ[fICent][1][0];
1477   Double_t Qyamean2 = fMeanQ[fICent][1][1];
1478   Double_t Qyarms2  = fWidthQ[fICent][1][1];
1479
1480   Double_t Qxcmean2 = fMeanQ[fICent][0][0];
1481   Double_t Qxcrms2  = fWidthQ[fICent][0][0];
1482   Double_t Qycmean2 = fMeanQ[fICent][0][1];
1483   Double_t Qycrms2  = fWidthQ[fICent][0][1];
1484
1485   Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1486   Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1487   Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1488   Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1489
1490   Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
1491     /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
1492   Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
1493     /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
1494
1495   Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
1496   Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
1497   Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
1498
1499   fhEPangleVZero->Fill(psiVZero);
1500   fhEPangleV0A->Fill(psiV0A);
1501   fhEPangleV0C->Fill(psiV0C);
1502
1503   //fill profile for resolution estimation
1504   fProfResolution->Fill(1, TMath::Cos(2.*(psiV0A - psiV0C)));
1505   fProfResolution->Fill(2, TMath::Cos(2.*(psiV0A - psiTPC)));
1506   fProfResolution->Fill(3, TMath::Cos(2.*(psiV0C - psiTPC)));
1507   
1508   Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
1509
1510   Double_t b = fAOD->GetMagneticField();
1511
1512   int nCascades=fAOD->GetNumberOfCascades();
1513   //Info("ReadFromAODv0", "# cascades = %d", nCascades);
1514
1515   const AliAODVertex *primaryBestAODVtx = fAOD->GetPrimaryVertex();
1516   primaryBestAODVtx->GetXYZ(bestPrimaryVtxPos);
1517   
1518   for(Int_t iXi = 0; iXi < nCascades; iXi++){
1519     // Double_t effMassXi = 0.;
1520     Double_t chi2Xi = -1.;
1521     Double_t dcaXiDaughters = -1.;
1522     Double_t XiCosOfPointingAngle = -1.;
1523     Double_t posXi[3] = {-1000., -1000., -1000.};
1524     Double_t XiRadius = -1000.;
1525     
1526     Double_t invMassLambdaAsCascDghter = 0.;
1527     Double_t V0Chi2Xi = -1.;
1528     Double_t dcaV0DaughtersXi = -1.;
1529     
1530     Double_t dcaBachToPrimaryVtxXi = -1.;
1531     Double_t dcaV0ToPrimaryVtxXi = -1.;
1532     Double_t dcaPosToPrimaryVtxXi = -1.;
1533     Double_t dcaNegToPrimaryVtxXi = -1.;
1534     Double_t V0CosOfPointingAngleXi = -1.;
1535     Double_t posV0Xi[3] = {-1000., -1000., -1000.};
1536     Double_t V0RadiusXi = -1000.;
1537     
1538     Double_t invMassXiMinus = 0.;
1539     Double_t invMassXiPlus = 0.;
1540     Double_t invMassOmegaMinus = 0.;
1541     Double_t invMassOmegaPlus = 0.;
1542
1543     /*
1544     Bool_t isPosInXiProton = kFALSE;
1545     Bool_t isPosInXiPion = kFALSE;
1546     Bool_t isPosInOmegaProton = kFALSE;
1547     Bool_t isPosInOmegaPion = kFALSE;
1548     
1549     Bool_t isNegInXiProton = kFALSE;
1550     Bool_t isNegInXiPion = kFALSE;
1551     Bool_t isNegInOmegaProton = kFALSE;
1552     Bool_t isNegInOmegaPion = kFALSE;
1553
1554     Bool_t isBachelorKaon = kFALSE;
1555     Bool_t isBachelorPion = kFALSE;
1556     */
1557
1558     Bool_t isBachelorKaonForTPC = kFALSE;
1559     Bool_t isBachelorPionForTPC = kFALSE;
1560     Bool_t isNegPionForTPC = kFALSE;
1561     Bool_t isPosPionForTPC = kFALSE;
1562     Bool_t isNegProtonForTPC = kFALSE;
1563     Bool_t isPosProtonForTPC = kFALSE;
1564
1565     Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
1566     Double_t XiPt = 0.;
1567     Double_t XiPtot = 0.;
1568     
1569     Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
1570     Double_t bachPt = 0.;
1571     Double_t bachPtot = 0.;
1572     
1573     //Short_t chargeXi = -2;
1574     Double_t V0toXiCosOfPointingAngle = 0.;
1575     
1576     Double_t rapXi = -20.;
1577     Double_t rapOmega = -20.;
1578     Double_t phi = 6.3;
1579     Double_t alphaXi = -200.;
1580     Double_t ptArmXi = -200.;
1581
1582     Double_t distToVtxZBefore = -999.;
1583     Double_t distToVtxZAfter = -999.;
1584     Double_t distToVtxXYBefore = -999.;
1585     Double_t distToVtxXYAfter = -999.;
1586     Double_t XiPAfter[3] = {-999., -999., -999.};
1587     Double_t phiAfter = -999.;
1588
1589     const AliAODcascade *xi = fAOD->GetCascade(iXi);
1590     if (!xi) continue;
1591     
1592     //effMassXi = xi->MassXi(); //default working hypothesis: Xi- decay
1593     chi2Xi = xi->Chi2Xi();
1594     dcaXiDaughters = xi->DcaXiDaughters();
1595     XiCosOfPointingAngle = xi->CosPointingAngleXi(bestPrimaryVtxPos[0],
1596                                                   bestPrimaryVtxPos[1],
1597                                                   bestPrimaryVtxPos[2]);
1598     posXi[0] = xi->DecayVertexXiX();
1599     posXi[1] = xi->DecayVertexXiY();
1600     posXi[2] = xi->DecayVertexXiZ();
1601     XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
1602                            +posXi[1]*posXi[1]
1603                            +posXi[2]*posXi[2]);
1604     
1605     AliAODTrack *pTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1606     AliAODTrack *nTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1607     AliAODTrack *bTrkXi 
1608       = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1609
1610     if(!pTrkXi || !nTrkXi || !bTrkXi) continue;
1611     
1612     UInt_t idxPosXi  = (UInt_t) TMath::Abs( pTrkXi->GetID() );
1613     UInt_t idxNegXi  = (UInt_t) TMath::Abs( nTrkXi->GetID() );
1614     UInt_t idxBach   = (UInt_t) TMath::Abs( bTrkXi->GetID() );
1615
1616     if(idxBach == idxNegXi || idxBach == idxPosXi) continue;
1617
1618     if( !fCutsDau->IsSelected(pTrkXi) 
1619         || !fCutsDau->IsSelected(nTrkXi)
1620         || !fCutsDau->IsSelected(bTrkXi) ) continue;
1621     
1622     if(pTrkXi->IsOn(AliESDtrack::kTPCin)){
1623       fh2TPCdEdxOfCascDghters->Fill(pTrkXi->P()*pTrkXi->Charge(), 
1624                                     pTrkXi->GetTPCsignal());
1625     }
1626     if( nTrkXi->IsOn(AliESDtrack::kTPCin) ){
1627       fh2TPCdEdxOfCascDghters->Fill(nTrkXi->P()*nTrkXi->Charge(),
1628                                     nTrkXi->GetTPCsignal());
1629     }
1630     if(bTrkXi->IsOn(AliESDtrack::kTPCin)){
1631       fh2TPCdEdxOfCascDghters->Fill(bTrkXi->P()*bTrkXi->Charge(),
1632                                     bTrkXi->GetTPCsignal());
1633     }
1634
1635
1636     if(xi->ChargeXi() < 0)
1637       invMassLambdaAsCascDghter = xi->MassLambda();
1638     else
1639       invMassLambdaAsCascDghter = xi->MassAntiLambda();
1640
1641     dcaV0DaughtersXi = xi->DcaV0Daughters();
1642     V0Chi2Xi = xi->Chi2V0();
1643     
1644      V0CosOfPointingAngleXi 
1645        = xi->CosPointingAngle(bestPrimaryVtxPos);
1646      dcaV0ToPrimaryVtxXi = xi->DcaV0ToPrimVertex();
1647      dcaBachToPrimaryVtxXi = xi->DcaBachToPrimVertex();
1648
1649      //V0
1650      posV0Xi[0] = xi->DecayVertexV0X();
1651      posV0Xi[1] = xi->DecayVertexV0Y();
1652      posV0Xi[2] = xi->DecayVertexV0Z();
1653      V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
1654                              +posV0Xi[1]*posV0Xi[1]
1655                               +posV0Xi[2]*posV0Xi[2]);
1656      dcaPosToPrimaryVtxXi = xi->DcaPosToPrimVertex();
1657      dcaNegToPrimaryVtxXi = xi->DcaNegToPrimVertex();
1658
1659
1660      //apply cuts ?
1661      // if(XiRadius < 1. || XiRadius > 100.) continue;
1662      //if(dcaXiDaughters > 0.1) continue;
1663      //if(XiCosOfPointingAngle < 0.999) continue;
1664      //if(dcaV0ToPrimaryVtxXi < 0.05) continue;
1665      //if(dcaBachToPrimaryVtxXi < 0.03) continue;
1666     
1667      //V0 mass cut?
1668      if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.01) continue;
1669
1670      //    if(dcaV0DaughtersXi > 1.) continue;
1671      //if(V0CosOfPointingAngleXi > 0.9999) continue;
1672      //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1673      //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1674
1675      //if(V0RadiusXi < 1.0 || V0RadiusXi > 100) continue;
1676     
1677      //other cuts?
1678
1679      if(xi->ChargeXi()<0){
1680        invMassXiMinus = xi->MassXi();
1681        invMassOmegaMinus = xi->MassOmega();
1682      }else{
1683        invMassXiPlus = xi->MassXi();
1684        invMassOmegaPlus = xi->MassOmega();
1685      }
1686
1687      /*
1688      if(pTrkXi->GetMostProbablePID() == AliAODTrack::kProton) {
1689        isPosInXiProton = kTRUE;
1690        isPosInOmegaProton = kTRUE;
1691      }
1692      if(pTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1693        isPosInXiPion = kTRUE;
1694        isPosInOmegaPion = kTRUE;
1695      }
1696     
1697      if(nTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1698        isNegInXiPion = kTRUE;
1699        isNegInOmegaPion = kTRUE;
1700      }
1701      if(nTrkXi->GetMostProbablePID() == AliAODTrack::kProton){
1702        isNegInXiProton = kTRUE;
1703        isNegInOmegaProton = kTRUE;
1704      }
1705
1706      if(bTrkXi->GetMostProbablePID() == AliAODTrack::kPion)
1707        isBachelorPion = kTRUE;
1708      if(bTrkXi->GetMostProbablePID() == AliAODTrack::kKaon)
1709        isBachelorKaon = kTRUE;
1710      */
1711
1712      //PID with TPC only: ??? Fix me
1713      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
1714        isBachelorKaonForTPC = kTRUE;
1715      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
1716        isBachelorPionForTPC = kTRUE;
1717
1718      //Negative V0 daughter
1719      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
1720        isNegPionForTPC = kTRUE;
1721      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
1722        isNegProtonForTPC = kTRUE;
1723     
1724      //Positive V0 daughter
1725      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
1726        isPosPionForTPC = kTRUE;
1727      if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
1728        isPosProtonForTPC = kTRUE;
1729
1730      //Extra QA information
1731      XiPx = xi->MomXiX();
1732      XiPy = xi->MomXiY();
1733      XiPz = xi->MomXiZ();
1734      XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
1735      XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
1736     
1737      bachPx = xi->MomBachX();
1738      bachPy = xi->MomBachY();
1739      bachPz = xi->MomBachZ();
1740
1741      bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1742      bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1743   
1744      V0toXiCosOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1745     
1746      rapXi = xi->RapXi();
1747      rapOmega = xi->RapOmega();
1748      phi = xi->Phi();
1749      alphaXi = xi->AlphaXi();
1750      ptArmXi = xi->PtArmXi();
1751
1752      distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1753      distToVtxXYBefore
1754        = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1755                      *(posXi[0] - bestPrimaryVtxPos[0])
1756                      +(posXi[1] - bestPrimaryVtxPos[1])
1757                      *(posXi[1] - bestPrimaryVtxPos[1]));
1758
1759      XiPAfter[0] = XiPx;
1760      XiPAfter[1] = XiPy;
1761      XiPAfter[2] = XiPz;
1762      
1763      //propagation to the best primary vertex to determine the momentum
1764      Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->ChargeXi());
1765      distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
1766      distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1767                                     *(posXi[0] - bestPrimaryVtxPos[0])
1768                                     +(posXi[1] - bestPrimaryVtxPos[1])
1769                                     *(posXi[1] - bestPrimaryVtxPos[1]));
1770      phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
1771      
1772      fh1DistToVtxZAfter->Fill(distToVtxZAfter);
1773      fh1DistToVtxXYAfter->Fill(distToVtxXYAfter);
1774      fh2DistToVtxZBeforeVsAfter->Fill(distToVtxZBefore, distToVtxZAfter);
1775      fh2DistToVtxXYBeforeVsAfter->Fill(distToVtxXYBefore, distToVtxXYAfter);
1776      fh2PxBeforeVsAfter->Fill(XiPx, XiPAfter[0]);
1777      fh2PyBeforeVsAfter->Fill(XiPy, XiPAfter[1]);
1778      if(xi->ChargeXi()>0)
1779        fh2PhiPosBeforeVsAfter->Fill(phi, phiAfter);
1780      else if(xi->ChargeXi()<0)
1781        fh2PhiNegBeforeVsAfter->Fill(phi, phiAfter);
1782
1783      if( (xi->ChargeXi() < 0 && isBachelorPionForTPC 
1784           && isPosProtonForTPC && isNegPionForTPC)
1785          || (xi->ChargeXi() > 0 && isBachelorPionForTPC 
1786              && isNegProtonForTPC && isPosPionForTPC))
1787        {//Xi candidate
1788         //for default hypothesis
1789          fh1Chi2Xi->Fill(chi2Xi);
1790          fh1DCAXiDaughters->Fill(dcaXiDaughters);
1791          fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1792          fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1793          fh1XiRadius->Fill(XiRadius);
1794         
1795          //V0
1796          fh1MassLambda->Fill(invMassLambdaAsCascDghter);
1797          fh1V0Chi2->Fill(V0Chi2Xi);
1798          fh1DcaV0DaughtersXi->Fill(dcaV0DaughtersXi);
1799          fh1V0CosOfPointingAngle->Fill(V0CosOfPointingAngleXi);
1800          fh1V0Radius->Fill(V0RadiusXi);
1801          fh1DcaV0ToPrimVertex->Fill(dcaV0ToPrimaryVtxXi);
1802          fh1DCAPosToPrimVertex->Fill(dcaPosToPrimaryVtxXi);
1803          fh1DCANegToPrimVertex->Fill(dcaNegToPrimaryVtxXi);
1804          fh1ChargeXi->Fill(xi->ChargeXi());
1805          fh1V0toXiCosOfPointingAngle->Fill(V0toXiCosOfPointingAngle);
1806        
1807          if ( TMath::Abs( invMassXiMinus-1.3217 ) < 0.012 
1808               || TMath::Abs( invMassXiPlus-1.3217 ) < 0.012) 
1809            {// One InvMass should be different from 0             
1810              fh1XiPt->Fill(XiPt);
1811              fh1XiP->Fill(XiPtot);
1812             
1813              fh1XiBachPt->Fill(bachPt);
1814              fh1XiBachP->Fill(bachPtot);
1815              fh1PhiXi->Fill( xi->Phi() );
1816            }
1817          fh2Armenteros->Fill(alphaXi, ptArmXi);
1818        }
1819   
1820      if(xi->ChargeXi()<0){
1821        fh1MassXiMinus->Fill(invMassXiMinus);
1822        fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1823        fh1MassXi->Fill(invMassXiMinus);
1824        fh1MassOmega->Fill(invMassOmegaMinus);
1825       
1826        fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1827                                         invMassXiMinus);
1828        fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus, 
1829                                        invMassOmegaMinus);
1830        fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1831        fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1832      }
1833      
1834      if(xi->ChargeXi() > 0){
1835        fh1MassXiPlus->Fill(invMassXiPlus);
1836        fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1837        fh1MassXi->Fill(invMassXiPlus);
1838        fh1MassOmega->Fill(invMassOmegaPlus);
1839       
1840        fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter, 
1841                                        invMassXiPlus);
1842        fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1843                                       invMassOmegaPlus);
1844        fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1845        fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1846      }
1847
1848
1849      //      Double_t phiNew 
1850      // = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. ) 
1851      // ? 0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1852       Double_t phiV0 = phiAfter;
1853       phiV0 -= psiVZero;
1854       //if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1855       //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1856
1857
1858       Double_t phiV0A = phiAfter;
1859       phiV0A -= psiV0A;
1860       //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1861       //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1862
1863       Double_t phiV0C = phiAfter;
1864       phiV0C -= psiV0C;
1865       //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1866       //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1867       
1868       //PID cuts with TPC cuts
1869       if(xi->ChargeXi() < 0){
1870         if(isPosProtonForTPC
1871            && isNegPionForTPC){
1872           if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1873             //Xi
1874             fh2MassVsPtXiMinus->Fill(XiPt, invMassXiMinus);
1875             fh2MassVsPtXiAll->Fill(XiPt, invMassXiMinus);
1876             fhXiRapidity->Fill(rapXi);
1877             for(int r=0; r!=3; ++r) {
1878               if(invMassXiMinus > fXiBands[r][0] 
1879                  && invMassXiMinus < fXiBands[r][1]){
1880                 
1881                 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1882                 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1883                 fProf2dXiV2PtV0A[r]->Fill(XiPt, psiV0A, TMath::Cos(2.*phiV0A));
1884
1885                 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1886                 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1887                 fProf2dXiV2PtV0C[r]->Fill(XiPt, psiV0C, TMath::Cos(2.*phiV0C));
1888
1889                 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1890                 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1891                 fProf2dXiV2Pt[r]->Fill(XiPt, psiVZero, TMath::Cos(2.*phiV0));
1892               }
1893             }
1894           }
1895           
1896           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1897             fh2MassVsPtOmegaMinus->Fill(XiPt, invMassOmegaMinus);
1898             fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaMinus); 
1899             fhOmegaRapidity->Fill(rapOmega);
1900             for(int r=0; r!=3; ++r) {
1901               if(invMassOmegaMinus > fOmegaBands[r][0]
1902                  && invMassOmegaMinus < fOmegaBands[r][1]){
1903                 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1904                 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1905                 fProf2dOmegaV2PtV0A[r]->Fill(XiPt, psiV0A, 
1906                                              TMath::Cos(2.*phiV0A));
1907
1908                 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1909                 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1910                 fProf2dOmegaV2PtV0C[r]->Fill(XiPt, psiV0C, 
1911                                              TMath::Cos(2.*phiV0C));
1912
1913                 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1914                 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1915                 fProf2dOmegaV2Pt[r]->Fill(XiPt, psiVZero, 
1916                                           TMath::Cos(2.*phiV0));
1917               }
1918             }
1919           }
1920         }
1921       }//if charge < 0
1922
1923       if(xi->ChargeXi() > 0){
1924         if(isNegProtonForTPC
1925            && isPosPionForTPC){
1926           if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1927             //Xi
1928             fh2MassVsPtXiPlus->Fill(XiPt, invMassXiPlus);
1929             fh2MassVsPtXiAll->Fill(XiPt, invMassXiPlus);
1930             fhXiRapidity->Fill(rapXi);
1931             for(int r=0; r!=3; ++r) {
1932               if(invMassXiPlus > fXiBands[r][0]
1933                  && invMassXiPlus < fXiBands[r][1]){
1934                 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1935                 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1936                 fProf2dXiV2PtV0A[r]->Fill(XiPt, psiV0A, TMath::Cos(2.*phiV0A));
1937
1938                 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1939                 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1940                 fProf2dXiV2PtV0C[r]->Fill(XiPt, psiV0C, TMath::Cos(2.*phiV0C));
1941
1942                 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1943                 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1944                 fProf2dXiV2Pt[r]->Fill(XiPt, psiVZero, TMath::Cos(2.*phiV0));
1945               }
1946             }
1947           }
1948
1949           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1950             //Omega
1951             fh2MassVsPtOmegaPlus->Fill(XiPt, invMassOmegaPlus);
1952             fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaPlus);
1953             fhOmegaRapidity->Fill(rapOmega);
1954             for(int r=0; r!=3; ++r) {
1955               if(invMassOmegaPlus > fOmegaBands[r][0]
1956                  && invMassOmegaPlus < fOmegaBands[r][1]){
1957                 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1958                 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1959                 fProf2dOmegaV2PtV0A[r]->Fill(XiPt, psiV0A,
1960                                              TMath::Cos(2.*phiV0A));
1961
1962
1963                 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1964                 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1965                 fProf2dOmegaV2PtV0C[r]->Fill(XiPt, psiV0C,
1966                                              TMath::Cos(2.*phiV0C));
1967
1968                 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1969                 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1970                 fProf2dOmegaV2Pt[r]->Fill(XiPt, psiVZero,
1971                                           TMath::Cos(2.*phiV0));
1972               }
1973             }
1974           }
1975         }
1976       }// if charge > 0
1977
1978   }//end of cascade loop
1979
1980   PostData(1, fHistList);
1981       
1982 }
1983
1984
1985 AliAnalysisTaskFlowEPCascade::~AliAnalysisTaskFlowEPCascade(){
1986   if(fHistList) delete fHistList;
1987 }
1988
1989 void AliAnalysisTaskFlowEPCascade::Terminate(Option_t *) 
1990 {
1991   /*
1992     fHistList = dynamic_cast<TList*> (GetOutputData(1));
1993     
1994     if (!fHistList) {
1995     printf("ERROR: Output tree not available\n");
1996     return;
1997     }
1998   */
1999 }      
2000
2001
2002 void AliAnalysisTaskFlowEPCascade::Propagate(Double_t vv[3],
2003                                              Double_t x[3],
2004                                              Double_t p[3],
2005                                              Double_t bz,
2006                                              Short_t sign){
2007   //Propagation to the primary vertex to determine the px and py 
2008   //x, p are the position and momentum as input and output
2009   //bz is the magnetic field along z direction
2010   //sign is the charge of particle for propagation
2011   
2012   Double_t pp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
2013   Double_t len = (vv[2]-x[2])*pp/p[2];
2014   Double_t a = -kB2C*bz*sign;
2015
2016   Double_t rho = a/pp;
2017   x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
2018   x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
2019   x[2] += p[2]*len/pp;
2020
2021   Double_t p0=p[0];
2022   p[0] = p0  *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
2023   p[1] = p[1]*TMath::Cos(rho*len) + p0  *TMath::Sin(rho*len);
2024 }
2025
2026
2027 void 
2028 AliAnalysisTaskFlowEPCascade::OpenInfoCalbration(Int_t run){
2029
2030   AliOADBContainer *cont = (AliOADBContainer*) fOADB->Get("hMultV0BefCorr");
2031   if(!cont){
2032     printf("OADB object hMultV0BefCorr is not available in the file\n");
2033     return;
2034   }
2035
2036   if(!(cont->GetObject(run))){
2037     printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
2038     run = 137366;
2039   }
2040   fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
2041
2042   TF1 *fpol0 = new TF1("fpol0","pol0");
2043   fMultV0->Fit(fpol0,"","",0,32);
2044   fV0Cpol = fpol0->GetParameter(0);
2045   fMultV0->Fit(fpol0,"","",32,64);
2046   fV0Apol = fpol0->GetParameter(0);
2047
2048   for(Int_t iside=0;iside<2;iside++){
2049     for(Int_t icoord=0;icoord<2;icoord++){
2050       for(Int_t i=0;i  < 9;i++){
2051         char namecont[100];
2052         if(iside==0 && icoord==0)
2053           snprintf(namecont,100,"hQxc2_%i", i);
2054         else if(iside==1 && icoord==0)
2055           snprintf(namecont,100,"hQxa2_%i", i);
2056         else if(iside==0 && icoord==1)
2057           snprintf(namecont,100,"hQyc2_%i", i);
2058         else if(iside==1 && icoord==1)
2059           snprintf(namecont,100,"hQya2_%i", i);
2060       
2061         cont = (AliOADBContainer*) fOADB->Get(namecont);
2062         if(!cont){ 
2063           printf("OADB object %s is not available in the file\n",namecont);
2064           return;
2065         }
2066         
2067         if(!(cont->GetObject(run))){
2068           printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
2069           run = 137366;
2070         }
2071         fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2072         fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2073       }
2074     }
2075   }
2076 }