AliPhysics  fceccc5 (fceccc5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFlowEPCascade.cxx
Go to the documentation of this file.
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_PHYSICS/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 }
AliAnalysisTaskFlowEPCascade::fh2MassXiVsMassOmegaMinus
TH2F * fh2MassXiVsMassOmegaMinus
Xi- effective mass vs V0 eff. mass.
Definition: AliAnalysisTaskFlowEPCascade.h:109
Double_t
double Double_t
Definition: External.C:58
TH2F
Definition: External.C:236
AliAnalysisTaskFlowEPCascade::fh2MassLambdaVsMassXiMinus
TH2F * fh2MassLambdaVsMassXiMinus
alpha vs ptArm for casc. candidate
Definition: AliAnalysisTaskFlowEPCascade.h:108
AliFlowTrackCuts.h
AliAnalysisTaskFlowEPCascade::fh1MassOmegaPlus
TH1F * fh1MassOmegaPlus
effective mass under Omega- hyp.
Definition: AliAnalysisTaskFlowEPCascade.h:91
centrality
centrality
Definition: AverageDmesonRaa.C:71
AliAnalysisTaskFlowEPCascade::Propagate
void Propagate(Double_t vv[3], Double_t x[3], Double_t p[3], Double_t bz, Short_t sign)
Definition: AliAnalysisTaskFlowEPCascade.cxx:2002
AliAnalysisTaskFlowEPCascade::fh2XiRadiusVsMassXiPlus
TH2F * fh2XiRadiusVsMassXiPlus
decay radius under Xi- hyp.
Definition: AliAnalysisTaskFlowEPCascade.h:114
AliFlowTrackCuts::IsSelected
virtual Bool_t IsSelected(TObject *obj, Int_t id=-666)
Definition: AliFlowTrackCuts.cxx:823
AliAnalysisTaskFlowEPCascade::fMultV0
TProfile * fMultV0
centrality bin number
Definition: AliAnalysisTaskFlowEPCascade.h:60
AliAnalysisTaskFlowEPCascade::fh1DCANegToPrimVertex
TH1F * fh1DCANegToPrimVertex
V0 positive daughter to prim. vertex.
Definition: AliAnalysisTaskFlowEPCascade.h:86
AliAnalysisTaskFlowEPCascade::fh1PhiXi
TH1F * fh1PhiXi
cos of pointing angle btw the V0 mom and the Xi-V0 vtx line
Definition: AliAnalysisTaskFlowEPCascade.h:104
ClassImp
ClassImp(AliAnalysisTaskFlowEPCascade) AliAnalysisTaskFlowEPCascade
Definition: AliAnalysisTaskFlowEPCascade.cxx:51
Int_t
int Int_t
Definition: External.C:63
TH1I
Definition: External.C:204
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliAnalysisTaskFlowEPCascade::fProf2dOmegaV2Pt
TProfile2D * fProf2dOmegaV2Pt[3]
three mass bands V0A&V0C
Definition: AliAnalysisTaskFlowEPCascade.h:151
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskFlowEPCascade::fProf2dOmegaV2PtV0A
TProfile2D * fProf2dOmegaV2PtV0A[3]
three mass bands V0A
Definition: AliAnalysisTaskFlowEPCascade.h:147
AliAnalysisTaskFlowEPCascade::fMeanQ
Float_t fMeanQ[9][2][2]
loaded by OADB
Definition: AliAnalysisTaskFlowEPCascade.h:63
AliAnalysisTaskFlowEPCascade::UserExec
virtual void UserExec(Option_t *option)
Definition: AliAnalysisTaskFlowEPCascade.cxx:748
AliAnalysisTaskFlowEPCascade.h
AliAnalysisTaskFlowEPCascade::fICent
Int_t fICent
current run checked to load VZERO calibration
Definition: AliAnalysisTaskFlowEPCascade.h:58
kProton
Definition: FastCentEstimators.C:25
AliAnalysisTaskFlowEPCascade::fh2MassVsPtXiMinus
TH2F * fh2MassVsPtXiMinus
dEdx with the cascade daughters
Definition: AliAnalysisTaskFlowEPCascade.h:120
AliFlowTrackSimple.h
AliAnalysisTaskFlowEPCascade::fh1MassXiMinus
TH1F * fh1MassXiMinus
V0 neg. daughter to prim. vertex.
Definition: AliAnalysisTaskFlowEPCascade.h:88
Short_t
short Short_t
Definition: External.C:23
centMax
Double_t centMax
Definition: MuonTrackingEfficiency.C:51
AliAnalysisTaskFlowEPCascade::fWidthQ
Float_t fWidthQ[9][2][2]
and recentering
Definition: AliAnalysisTaskFlowEPCascade.h:64
AliAnalysisTaskFlowEPCascade::fV0Cpol
Float_t fV0Cpol
object containing VZERO calibration info
Definition: AliAnalysisTaskFlowEPCascade.h:61
AliAnalysisTaskFlowEPCascade::fh1V0CosOfPointingAngle
TH1F * fh1V0CosOfPointingAngle
for V0 associated to a cascade
Definition: AliAnalysisTaskFlowEPCascade.h:81
AliAnalysisTaskFlowEPCascade::fProf2dOmegaV2PtV0C
TProfile2D * fProf2dOmegaV2PtV0C[3]
three mass bands V0C
Definition: AliAnalysisTaskFlowEPCascade.h:149
AliAnalysisTaskFlowEPCascade::fh1DistToVtxXYAfter
TH1F * fh1DistToVtxXYAfter
After propagation to vertex.
Definition: AliAnalysisTaskFlowEPCascade.h:156
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskFlowEPCascade::fh1MassXiPlus
TH1F * fh1MassXiPlus
effective mass under Xi- hyp.
Definition: AliAnalysisTaskFlowEPCascade.h:89
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskFlowEPCascade::fh1DCAPosToPrimVertex
TH1F * fh1DCAPosToPrimVertex
DCA of V0 to prim. vtx.
Definition: AliAnalysisTaskFlowEPCascade.h:85
AliAnalysisTaskFlowEPCascade::fProfOmegaV2PtV0A
TProfile * fProfOmegaV2PtV0A[3]
three mass bands V0A
Definition: AliAnalysisTaskFlowEPCascade.h:132
AliAnalysisTaskFlowEPCascade::fProfOmegaV2Pt
TProfile * fProfOmegaV2Pt[3]
three mass bands V0A&V0C
Definition: AliAnalysisTaskFlowEPCascade.h:142
AliAnalysisTaskFlowEPCascade::fProfOmegaV2PtV0C
TProfile * fProfOmegaV2PtV0C[3]
three mass bands V0C
Definition: AliAnalysisTaskFlowEPCascade.h:137
AliAnalysisTaskFlowEPCascade::fh2XiRadiusVsMassOmegaPlus
TH2F * fh2XiRadiusVsMassOmegaPlus
decay radius under Omega- hyp.
Definition: AliAnalysisTaskFlowEPCascade.h:116
AliAnalysisTaskFlowEPCascade::fh1V0Chi2
TH1F * fh1V0Chi2
mass of lambda as the cascade daughter
Definition: AliAnalysisTaskFlowEPCascade.h:80
centMin
Double_t centMin
Definition: MuonTrackingEfficiency.C:50