]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliAnalysisTaskFlowEPCascade.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
refresh random seed on each node
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskFlowEPCascade.cxx
CommitLineData
47976489 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
49using namespace std;
50
51ClassImp(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//_________________________________________________________________________
168AliAnalysisTaskFlowEPCascade::
169AliAnalysisTaskFlowEPCascade(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
288void 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
748void 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 = fAOD->GetHeader();
772 if(!aodHeader) return;
773 AliCentrality *centrality = aodHeader->GetCentralityP();
774 if(!centrality) return;
775 Double_t cent = centrality->GetCentralityPercentile("V0M" );
776 if(cent<fMinCent||cent>=fMaxCent) return; //centrality cut
777
778 if(cent >= 0 && cent < 5) fICent = 0;
779 else if(cent >= 5 && cent < 10) fICent = 1;
780 else if(cent >= 10 && cent < 20) fICent = 2;
781 else if(cent >= 20 && cent < 30) fICent = 3;
782 else if(cent >= 30 && cent < 40) fICent = 4;
783 else if(cent >= 40 && cent < 50) fICent = 5;
784 else if(cent >= 50 && cent < 60) fICent = 6;
785 else if(cent >= 60 && cent < 70) fICent = 7;
786 else if(cent >= 70 && cent < 80) fICent = 8;
787 else
788 return;
789
790 Double_t zvtx = fAOD->GetPrimaryVertex()->GetZ();
791 if(TMath::Abs(zvtx) > fVtxCut) return; //vertex cut
792
793 fhEvent->Fill(2);
794 ReadFromAODv0(fAOD);
795 }
796
797 return;
798}
799
800/*
801void AliAnalysisTaskFlowEPCascade::ReadFromESDv0(AliESDEvent *fESD){
802
803 // fEPcalib->CalculateEventPlanes(fESD);
804
805 //Get EP informations
806 //Double_t psiVZero = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0AC);
807 //fhEPangleVZero->Fill(psiVZero);
808
809 //Double_t psiV0A = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0A);
810 //Double_t psiV0C = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0C);
811 //Double_t psiTPC = fEPcalib->GetPsi(2, AliAnaEPcalib::kTPC);
812
813 // Double_t psiZDC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCAC);
814 //fhEPangleZDC->Fill(psiZDC);
815
816 //Double_t psiZDCA = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCA);
817 //Double_t psiZDCC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCC);
818
819 AliEventplane * ep = fESD->GetEventplane();
820 Double_t psiTPC = ep->GetEventplane("Q", fESD, 2);
821
822 Int_t run = fESD->GetRunNumber();
823 if(run != fRun){
824 // Load the calibrations run dependent
825 OpenInfoCalbration(run);
826 fRun=run;
827 }
828
829 //fill profile for resolution estimation
830 fProfResolution->Fill(1, TMath::Cos(2.*(psiV0A - psiV0C)));
831 fProfResolution->Fill(2, TMath::Cos(2.*(psiV0A - psiTPC)));
832 fProfResolution->Fill(3, TMath::Cos(2.*(psiV0C - psiTPC)));
833 //fProfResolution->Fill(4, TMath::Cos(2.*(psiZDCA - psiZDCC)));
834
835 //check Xi candidates
836 Double_t trkPrimaryVtxPos[3] = {-100., -100., -100.};
837 Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
838 int nCascades=fESD->GetNumberOfCascades();
839
840 const AliESDVertex *primaryTrackingESDVtx = fESD->GetPrimaryVertexTracks();
841 primaryTrackingESDVtx->GetXYZ(trkPrimaryVtxPos);
842
843 const AliESDVertex *primaryBestESDVtx = fESD->GetPrimaryVertex();
844 primaryBestESDVtx->GetXYZ(bestPrimaryVtxPos);
845
846 Double_t b = fESD->GetMagneticField();
847
848 for(int i = 0; i != nCascades; ++i) {
849 Double_t effMassXi = 0.;
850 Double_t chi2Xi = -1.;
851 Double_t dcaXiDaughters = -1.;
852 Double_t XiCosOfPointingAngle = -1.;
853 Double_t posXi[3] = {-1000., -1000., -1000.};
854 Double_t XiRadius = -1000.;
855
856 Double_t invMassLambdaAsCascDghter = 0.;
857 Double_t V0Chi2Xi = -1.;
858 Double_t dcaV0DaughtersXi = -1.;
859
860 Double_t dcaBachToPrimaryVtxXi = -1.;
861 Double_t dcaV0ToPrimaryVtxXi = -1.;
862 Double_t dcaPosToPrimaryVtxXi = -1.;
863 Double_t dcaNegToPrimaryVtxXi = -1.;
864 Double_t V0CosOfPointingAngleXi = -1.;
865 Double_t posV0Xi[3] = {-1000., -1000., -1000.};
866 Double_t V0RadiusXi = -1000.;
867 Double_t V0quality = 0.;
868
869 Double_t invMassXiMinus = 0.;
870 Double_t invMassXiPlus = 0.;
871 Double_t invMassOmegaMinus = 0.;
872 Double_t invMassOmegaPlus = 0.;
873
874 Bool_t isPosInXiProton = kFALSE;
875 Bool_t isPosInXiPion = kFALSE;
876 Bool_t isPosInOmegaProton = kFALSE;
877 Bool_t isPosInOmegaPion = kFALSE;
878
879 Bool_t isNegInXiProton = kFALSE;
880 Bool_t isNegInXiPion = kFALSE;
881 Bool_t isNegInOmegaProton = kFALSE;
882 Bool_t isNegInOmegaPion = kFALSE;
883
884 Bool_t isBachelorKaon = kFALSE;
885 Bool_t isBachelorPion = kFALSE;
886
887 Bool_t isBachelorKaonForTPC = kFALSE;
888 Bool_t isBachelorPionForTPC = kFALSE;
889 Bool_t isNegPionForTPC = kFALSE;
890 Bool_t isPosPionForTPC = kFALSE;
891 Bool_t isNegProtonForTPC = kFALSE;
892 Bool_t isPosProtonForTPC = kFALSE;
893
894 Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
895 Double_t XiPt = 0.;
896 Double_t XiPtot = 0.;
897
898 Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
899 Double_t bachPt = 0.;
900 Double_t bachPtot = 0.;
901
902 //Short_t chargeXi = -2;
903 Double_t V0toXiCosOfPointingAngle = 0.;
904
905 Double_t rapXi = -20.;
906 Double_t rapOmega = -20.;
907 Double_t phi = 6.3;
908 Double_t alphaXi = -200.;
909 Double_t ptArmXi = -200.;
910
911 Double_t distToVtxZBefore = -999.;
912 Double_t distToVtxZAfter = -999.;
913 Double_t distToVtxXYBefore = -999.;
914 Double_t distToVtxXYAfter = -999.;
915 Double_t XiPAfter[3] = {-999., -999., -999.};
916 Double_t phiAfter = -999.;
917
918 AliESDcascade *xi = fESD->GetCascade(i);
919 if(!xi) continue;
920
921 if(xi->Charge()<0)
922 xi->ChangeMassHypothesis(V0quality, 3312); // Xi- hypothesis
923 else if(xi->Charge() > 0)
924 xi->ChangeMassHypothesis(V0quality, -3312);
925 else continue;
926
927 effMassXi = xi->GetEffMassXi();
928 chi2Xi = xi->GetChi2Xi();
929 dcaXiDaughters = xi->GetDcaXiDaughters();
930 XiCosOfPointingAngle
931 = xi->GetCascadeCosineOfPointingAngle(bestPrimaryVtxPos[0],
932 bestPrimaryVtxPos[1],
933 bestPrimaryVtxPos[2]);
934 xi->GetXYZcascade(posXi[0], posXi[1], posXi[2]);
935 XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
936 +posXi[1]*posXi[1]
937 +posXi[2]*posXi[2]);
938
939 UInt_t idxPosXi = (UInt_t)TMath::Abs(xi->GetPindex());
940 UInt_t idxNegXi = (UInt_t)TMath::Abs(xi->GetNindex());
941 UInt_t idxBach = (UInt_t)TMath::Abs(xi->GetBindex());
942
943 if(idxBach == idxPosXi || idxBach == idxNegXi) continue;
944
945 AliESDtrack *pTrkXi = fESD->GetTrack(idxPosXi);
946 AliESDtrack *nTrkXi = fESD->GetTrack(idxNegXi);
947 AliESDtrack *bTrkXi = fESD->GetTrack(idxBach);
948
949 if( !pTrkXi || !nTrkXi || !bTrkXi ) continue;
950
951 if( !fCutsDau->IsSelected(pTrkXi)
952 || !fCutsDau->IsSelected(nTrkXi)
953 || !fCutsDau->IsSelected(bTrkXi) ) continue;
954
955 const AliExternalTrackParam *pExtTrk = pTrkXi->GetInnerParam();
956 const AliExternalTrackParam *nExtTrk = nTrkXi->GetInnerParam();
957 const AliExternalTrackParam *bExtTrk = bTrkXi->GetInnerParam();
958
959 if(pExtTrk && pTrkXi->IsOn(AliESDtrack::kTPCin)){
960 fh2TPCdEdxOfCascDghters->Fill(pExtTrk->GetP()*pExtTrk->Charge(),
961 pTrkXi->GetTPCsignal());
962 }
963 if(nExtTrk && nTrkXi->IsOn(AliESDtrack::kTPCin)){
964 fh2TPCdEdxOfCascDghters->Fill(nExtTrk->GetP()*nExtTrk->Charge(),
965 nTrkXi->GetTPCsignal());
966 }
967 if(bExtTrk && bTrkXi->IsOn(AliESDtrack::kTPCin)){
968 fh2TPCdEdxOfCascDghters->Fill(bExtTrk->GetP()*bExtTrk->Charge(),
969 bTrkXi->GetTPCsignal());
970 }
971
972
973 invMassLambdaAsCascDghter = xi->GetEffMass(); // from V0
974 dcaV0DaughtersXi = xi->GetDcaV0Daughters();
975 V0Chi2Xi = xi->GetChi2V0();
976
977 V0CosOfPointingAngleXi
978 = xi->GetV0CosineOfPointingAngle(bestPrimaryVtxPos[0],
979 bestPrimaryVtxPos[1],
980 bestPrimaryVtxPos[2]);
981 dcaV0ToPrimaryVtxXi = xi->GetD(bestPrimaryVtxPos[0],
982 bestPrimaryVtxPos[1],
983 bestPrimaryVtxPos[2]);
984 dcaBachToPrimaryVtxXi = TMath::Abs(bTrkXi->GetD(bestPrimaryVtxPos[0],
985 bestPrimaryVtxPos[1],
986 bestPrimaryVtxPos[2]));
987
988 //V0
989 xi->GetXYZ(posV0Xi[0], posV0Xi[1], posV0Xi[2]);
990 V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
991 +posV0Xi[1]*posV0Xi[1]
992 +posV0Xi[2]*posV0Xi[2]);
993 dcaPosToPrimaryVtxXi = TMath::Abs(pTrkXi->GetD(bestPrimaryVtxPos[0],
994 bestPrimaryVtxPos[1],
995 bestPrimaryVtxPos[2]));
996 dcaNegToPrimaryVtxXi = TMath::Abs(nTrkXi->GetD(bestPrimaryVtxPos[0],
997 bestPrimaryVtxPos[1],
998 bestPrimaryVtxPos[2]));
999
1000 //apply cuts
1001 // if(XiRadius < 1.5 || XiRadius > 100.) continue;
1002 //if(dcaXiDaughters > 0.3) continue;
1003 //if(XiCosOfPointingAngle < 0.999) continue;
1004 //if(dcaV0ToPrimaryVtxXi < 0.03) continue;
1005 //if(dcaBachToPrimaryVtxXi < 0.03) continue;
1006
1007 //V0 mass cut?
1008 //if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.006) continue;
1009
1010 //if(dcaV0DaughtersXi > .6) continue;
1011 //if(V0CosOfPointingAngleXi > 0.9999) continue;
1012 //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1013 //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1014
1015 //if(V0RadiusXi < 6.0 || V0RadiusXi > 100) continue;
1016
1017 //other cuts?
1018
1019 // change mass hypothesis to cover all the possibilities
1020 if(bTrkXi->Charge()<0){
1021 V0quality = 0.;
1022 xi->ChangeMassHypothesis(V0quality, 3312); //Xi- hyp.
1023 invMassXiMinus = xi->GetEffMassXi();
1024
1025 V0quality = 0.;
1026 xi->ChangeMassHypothesis(V0quality, 3334); //Omega- hyp.
1027 invMassOmegaMinus = xi->GetEffMassXi();
1028
1029 V0quality = 0.;
1030 xi->ChangeMassHypothesis(V0quality, 3312); //back to default hyp.
1031 }
1032
1033 if(bTrkXi->Charge() > 0){
1034 V0quality = 0.;
1035 xi->ChangeMassHypothesis(V0quality, -3312); //anti-Xi- hyp.
1036 invMassXiPlus = xi->GetEffMassXi();
1037
1038 V0quality = 0.;
1039 xi->ChangeMassHypothesis(V0quality, -3334); //anti-Omega- hyp.
1040 invMassOmegaPlus = xi->GetEffMassXi();
1041
1042 V0quality = 0.;
1043 xi->ChangeMassHypothesis(V0quality, -3312); //back to default hyp.
1044 }
1045
1046 //PID on the daughter tracks
1047 //A - Combined PID
1048 //Resonable guess the priors for the cascade track sample
1049 //(e, mu, pi, K, p)
1050 Double_t priorsGuessXi[5] = {0, 0, 2, 0, 1};
1051 Double_t priorsGuessOmega[5] = {0, 0, 1, 1, 1};
1052
1053 //Combined bachelor-daughter PID
1054 AliPID pidXi;
1055 pidXi.SetPriors(priorsGuessXi);
1056 AliPID pidOmega;
1057 pidOmega.SetPriors(priorsGuessOmega);
1058
1059 if(pTrkXi->IsOn(AliESDtrack::kESDpid)){// combined PID exists
1060 Double_t r[10] = {0.};
1061 pTrkXi->GetESDpid(r);
1062 pidXi.SetProbabilities(r);
1063 pidOmega.SetProbabilities(r);
1064
1065 //Check if the V0 postive track is proton (case for Xi-)
1066 Double_t pProton = pidXi.GetProbability(AliPID::kProton);
1067 if(pProton > pidXi.GetProbability(AliPID::kElectron)
1068 && pProton > pidXi.GetProbability(AliPID::kMuon)
1069 && pProton > pidXi.GetProbability(AliPID::kPion)
1070 && pProton > pidXi.GetProbability(AliPID::kKaon))
1071 isPosInXiProton = kTRUE;
1072
1073 //Check if the V0 postive track is a pi+ (case for Xi+)
1074 Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1075 if(pPion > pidXi.GetProbability(AliPID::kElectron)
1076 && pPion > pidXi.GetProbability(AliPID::kMuon)
1077 && pPion > pidXi.GetProbability(AliPID::kKaon)
1078 && pPion > pidXi.GetProbability(AliPID::kProton))
1079 isPosInXiPion = kTRUE;
1080 // Check if the V0 positive track is a proton (case for Omega-)
1081 pProton = pidOmega.GetProbability(AliPID::kProton);
1082 if(pProton > pidOmega.GetProbability(AliPID::kElectron)
1083 && pProton > pidOmega.GetProbability(AliPID::kMuon)
1084 && pProton > pidOmega.GetProbability(AliPID::kPion)
1085 && pProton > pidOmega.GetProbability(AliPID::kKaon))
1086 isPosInOmegaProton = kTRUE;
1087
1088 // Check if the V0 positive track is a pi+ (case for Omega+)
1089 pPion = pidOmega.GetProbability(AliPID::kPion);
1090 if(pPion > pidOmega.GetProbability(AliPID::kElectron)
1091 && pPion > pidOmega.GetProbability(AliPID::kMuon)
1092 && pPion > pidOmega.GetProbability(AliPID::kKaon)
1093 && pPion > pidOmega.GetProbability(AliPID::kProton))
1094 isPosInOmegaPion = kTRUE;
1095 }
1096
1097 //Combined V0-negative-daughter PID
1098 pidXi.SetPriors(priorsGuessXi);
1099 pidOmega.SetPriors(priorsGuessOmega);
1100 if(nTrkXi->IsOn(AliESDtrack::kESDpid)){
1101 Double_t r[10] = {0.};
1102 nTrkXi->GetESDpid(r);
1103 pidXi.SetProbabilities(r);
1104 pidOmega.SetProbabilities(r);
1105
1106 // Check if the V0 negative track is a pi- (case for Xi-)
1107 Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1108 if(pPion > pidXi.GetProbability(AliPID::kElectron)
1109 && pPion > pidXi.GetProbability(AliPID::kMuon)
1110 && pPion > pidXi.GetProbability(AliPID::kKaon)
1111 && pPion > pidXi.GetProbability(AliPID::kProton))
1112 isNegInXiPion = kTRUE;
1113
1114 // Check if the V0 negative track is an anti-proton (case for Xi+)
1115 Double_t pProton = pidXi.GetProbability(AliPID::kProton);
1116 if(pProton > pidXi.GetProbability(AliPID::kElectron)
1117 && pProton > pidXi.GetProbability(AliPID::kMuon)
1118 && pProton > pidXi.GetProbability(AliPID::kPion)
1119 && pProton > pidXi.GetProbability(AliPID::kKaon))
1120 isNegInXiProton = kTRUE;
1121
1122 // Check if the V0 negative track is a pi- (case for Omega-)
1123 pPion = pidOmega.GetProbability(AliPID::kPion);
1124 if(pPion > pidOmega.GetProbability(AliPID::kElectron)
1125 && pPion > pidOmega.GetProbability(AliPID::kMuon)
1126 && pPion > pidOmega.GetProbability(AliPID::kKaon)
1127 && pPion > pidOmega.GetProbability(AliPID::kProton))
1128 isNegInOmegaPion = kTRUE;
1129
1130 // Check if the V0 negative track is an anti-proton (case for Omega+)
1131 pProton = pidOmega.GetProbability(AliPID::kProton);
1132 if(pProton > pidOmega.GetProbability(AliPID::kElectron)
1133 && pProton > pidOmega.GetProbability(AliPID::kMuon)
1134 && pProton > pidOmega.GetProbability(AliPID::kPion)
1135 && pProton > pidOmega.GetProbability(AliPID::kKaon))
1136 isNegInOmegaProton = kTRUE;
1137 }
1138
1139 // Combined bachelor PID
1140 pidXi.SetPriors(priorsGuessXi);
1141 pidOmega.SetPriors(priorsGuessOmega);
1142 if(bTrkXi->IsOn(AliESDtrack::kESDpid)){//Combined PID exists
1143 Double_t r[10] = {0.};
1144 bTrkXi->GetESDpid(r);
1145 pidXi.SetProbabilities(r);
1146 pidOmega.SetProbabilities(r);
1147
1148 //Check if the bachelor track is a pion
1149 Double_t pPion = pidXi.GetProbability(AliPID::kPion);
1150 if(pPion > pidXi.GetProbability(AliPID::kElectron)
1151 && pPion > pidXi.GetProbability(AliPID::kMuon)
1152 && pPion > pidXi.GetProbability(AliPID::kKaon)
1153 && pPion > pidXi.GetProbability(AliPID::kProton))
1154 isBachelorPion = kTRUE;
1155
1156 // Check if the bachelor track is a kaon
1157 Double_t pKaon = pidOmega.GetProbability(AliPID::kKaon);
1158 if(pKaon > pidOmega.GetProbability(AliPID::kElectron)
1159 && pKaon > pidOmega.GetProbability(AliPID::kMuon)
1160 && pKaon > pidOmega.GetProbability(AliPID::kPion)
1161 && pKaon > pidOmega.GetProbability(AliPID::kProton))
1162 isBachelorKaon = kTRUE;
1163 }
1164
1165 //B - TPC PID: 3-sigma bands on Bethe-Bloch curve
1166 //Bachelor
1167 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
1168 isBachelorKaonForTPC = kTRUE;
1169 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
1170 isBachelorPionForTPC = kTRUE;
1171
1172 //Negative V0 daughter
1173 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
1174 isNegPionForTPC = kTRUE;
1175 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
1176 isNegProtonForTPC = kTRUE;
1177
1178 //Positive V0 daughter
1179 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
1180 isPosPionForTPC = kTRUE;
1181 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
1182 isPosProtonForTPC = kTRUE;
1183
1184 // //Extra QA information
1185 xi->GetPxPyPz(XiPx, XiPy, XiPz);
1186 XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
1187 XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
1188
1189 xi->GetBPxPyPz(bachPx, bachPy, bachPz);
1190 bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1191 bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1192
1193 V0toXiCosOfPointingAngle
1194 = xi->GetV0CosineOfPointingAngle(posXi[0], posXi[1], posXi[2]);
1195 rapXi = xi->RapXi();
1196 rapOmega = xi->RapOmega();
1197 phi = xi->Phi();
1198 alphaXi = xi->AlphaXi();
1199 ptArmXi = xi->PtArmXi();
1200
1201 XiPAfter[0] = XiPx;
1202 XiPAfter[1] = XiPy;
1203 XiPAfter[2] = XiPz;
1204
1205 distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1206 distToVtxXYBefore
1207 = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1208 *(posXi[0] - bestPrimaryVtxPos[0])
1209 +(posXi[1] - bestPrimaryVtxPos[1])
1210 *(posXi[1] - bestPrimaryVtxPos[1]));
1211
1212 //propagation to the best primary vertex to determine the momentum
1213 Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->Charge());
1214 distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
1215 distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1216 *(posXi[0] - bestPrimaryVtxPos[0])
1217 +(posXi[1] - bestPrimaryVtxPos[1])
1218 *(posXi[1] - bestPrimaryVtxPos[1]));
1219 phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
1220
1221 fh1DistToVtxZAfter->Fill(distToVtxZAfter);
1222 fh1DistToVtxXYAfter->Fill(distToVtxXYAfter);
1223 fh2DistToVtxZBeforeVsAfter->Fill(distToVtxZBefore, distToVtxZAfter);
1224 fh2DistToVtxXYBeforeVsAfter->Fill(distToVtxXYBefore, distToVtxXYAfter);
1225 fh2PxBeforeVsAfter->Fill(XiPx, XiPAfter[0]);
1226 fh2PyBeforeVsAfter->Fill(XiPy, XiPAfter[1]);
1227 if(xi->Charge()>0)
1228 fh2PhiPosBeforeVsAfter->Fill(phi, phiAfter);
1229 else if(xi->Charge()<0)
1230 fh2PhiNegBeforeVsAfter->Fill(phi, phiAfter);
1231
1232
1233 if( isBachelorPion &&
1234 ( (xi->Charge()<0 && isPosInXiProton && isNegInXiPion)
1235 || (xi->Charge()>0 && isNegInXiProton && isPosInXiPion ) ) )
1236 {//Xi candidate
1237 //for default hypothesis
1238 fh1Chi2Xi->Fill(chi2Xi);
1239 fh1DCAXiDaughters->Fill(dcaXiDaughters);
1240 fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1241 fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1242 fh1XiRadius->Fill(XiRadius);
1243
1244 //V0
1245 fh1MassLambda->Fill(invMassLambdaAsCascDghter);
1246 fh1V0Chi2->Fill(V0Chi2Xi);
1247 fh1DcaV0DaughtersXi->Fill(dcaV0DaughtersXi);
1248 fh1V0CosOfPointingAngle->Fill(V0CosOfPointingAngleXi);
1249 fh1V0Radius->Fill(V0RadiusXi);
1250 fh1DcaV0ToPrimVertex->Fill(dcaV0ToPrimaryVtxXi);
1251 fh1DCAPosToPrimVertex->Fill(dcaPosToPrimaryVtxXi);
1252 fh1DCANegToPrimVertex->Fill(dcaNegToPrimaryVtxXi);
1253 fh1ChargeXi->Fill(xi->Charge());
1254 fh1V0toXiCosOfPointingAngle->Fill(V0toXiCosOfPointingAngle);
1255
1256 if ( TMath::Abs( invMassXiMinus-1.3217 ) < 0.012
1257 || TMath::Abs( invMassXiPlus-1.3217 ) < 0.012)
1258 {// One InvMass should be different from 0
1259 fh1XiPt->Fill(XiPt);
1260 fh1XiP->Fill(XiPtot);
1261
1262 fh1XiBachPt->Fill(bachPt);
1263 fh1XiBachP->Fill(bachPtot);
1264 fh1PhiXi->Fill( xi->Phi() );
1265 }
1266 fh2Armenteros->Fill(alphaXi, ptArmXi);
1267 }
1268
1269 if(xi->Charge()<0){
1270 fh1MassXiMinus->Fill(invMassXiMinus);
1271 fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1272 fh1MassXi->Fill(invMassXiMinus);
1273 fh1MassOmega->Fill(invMassOmegaMinus);
1274
1275 fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1276 invMassXiMinus);
1277 fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus,
1278 invMassOmegaMinus);
1279 fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1280 fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1281 }
1282
1283 if(xi->Charge() > 0){
1284 fh1MassXiPlus->Fill(invMassXiPlus);
1285 fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1286 fh1MassXi->Fill(invMassXiPlus);
1287 fh1MassOmega->Fill(invMassOmegaPlus);
1288
1289 fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter,
1290 invMassXiPlus);
1291 fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1292 invMassOmegaPlus);
1293 fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1294 fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1295 }
1296
1297
1298 Double_t phiNew
1299 = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. ) ?
1300 0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1301 Double_t phiV0 = phiNew;
1302 phiV0 -= psiVZero;
1303 // if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1304 //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1305
1306 Double_t phiV0A = phiNew;
1307 phiV0A -= psiV0A;
1308 //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1309 //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1310
1311 Double_t phiV0C = phiNew;
1312 phiV0C -= psiV0C;
1313 //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1314 //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1315
1316 //PID cuts with TPC cuts
1317 if(xi->Charge() < 0){
1318 if(isPosProtonForTPC
1319 && isNegPionForTPC){
1320 if(isBachelorPionForTPC){
1321 //Xi
1322 fhXiRapidity->Fill(rapXi);
1323 if(TMath::Abs(rapXi) < 0.8){
1324 fh2MassVsPtXiMinus->Fill(XiPt, invMassXiMinus);
1325 fh2MassVsPtXiAll->Fill(XiPt, invMassXiMinus);
1326
1327 for(int r=0; r!=3; ++r) {
1328 if(invMassXiMinus > fXiBands[r][0]
1329 && invMassXiMinus < fXiBands[r][1]){
1330
1331 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1332 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1333
1334 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1335 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1336
1337 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1338 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1339 }
1340 }
1341
1342 }
1343 }
1344
1345 if(isBachelorKaonForTPC){
1346 //Omega
1347 fhOmegaRapidity->Fill(rapOmega);
1348 if(TMath::Abs(rapOmega) < 0.8){
1349 fh2MassVsPtOmegaMinus->Fill(XiPt, invMassOmegaMinus);
1350 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaMinus);
1351
1352 for(int r=0; r!=3; ++r) {
1353 if(invMassOmegaMinus > fOmegaBands[r][0]
1354 && invMassOmegaMinus < fOmegaBands[r][1]){
1355
1356 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1357 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1358
1359 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1360 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1361
1362 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1363 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1364 }
1365 }
1366
1367
1368 }
1369 }
1370 }
1371 }
1372
1373 if(xi->Charge() > 0){
1374 if(isNegProtonForTPC
1375 && isPosPionForTPC){
1376 if(isBachelorPionForTPC){
1377 //Xi
1378 fhXiRapidity->Fill(rapXi);
1379 if(TMath::Abs(rapXi) < 0.8){
1380 fh2MassVsPtXiPlus->Fill(XiPt, invMassXiPlus);
1381 fh2MassVsPtXiAll->Fill(XiPt, invMassXiPlus);
1382 for(int r=0; r!=3; ++r) {
1383 if(invMassXiPlus > fXiBands[r][0]
1384 && invMassXiPlus < fXiBands[r][1]){
1385 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1386 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1387
1388 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1389 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1390
1391 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1392 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1393 }
1394 }
1395
1396 }
1397 }
1398
1399 if(isBachelorKaonForTPC){
1400 //Omega
1401 fhOmegaRapidity->Fill(rapOmega);
1402 if(TMath::Abs(rapOmega) < 0.8){
1403 fh2MassVsPtOmegaPlus->Fill(XiPt, invMassOmegaPlus);
1404 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaPlus);
1405
1406 for(int r=0; r!=3; ++r) {
1407 if(invMassOmegaPlus > fOmegaBands[r][0]
1408 && invMassOmegaPlus < fOmegaBands[r][1]){
1409 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1410 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1411
1412 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1413 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1414
1415 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1416 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1417 }
1418 }
1419
1420
1421 }
1422 }
1423 }
1424 }
1425
1426 }//end of cascade loop
1427
1428 // fNEvent++;
1429
1430 PostData(1, fHistList);
1431}
1432*/
1433
1434void AliAnalysisTaskFlowEPCascade::ReadFromAODv0(AliAODEvent *fAOD){
1435
1436 AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
1437 Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
1438 // if(psiTPC > TMath::PiOver2())
1439 // psiTPC -= TMath::Pi();
1440 fhEPangleTPC->Fill(psiTPC);
1441
1442 Int_t run = fAOD->GetRunNumber();
1443 if(run != fRun){
1444 // Load the calibrations run dependent
1445 OpenInfoCalbration(run);
1446 fRun=run;
1447 }
1448
1449 //reset Q vector info
1450 Double_t Qxa2 = 0, Qya2 = 0;
1451 Double_t Qxc2 = 0, Qyc2 = 0;
1452 //Double_t Qx2 = 0, Qy2 = 0;
1453
1454 //V0 info
1455 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
1456
1457 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1458 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1459 Float_t multv0 = aodV0->GetMultiplicity(iv0);
1460 if (iv0 < 32){ // V0C
1461 Qxc2 += TMath::Cos(2*phiV0)*multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1462 Qyc2 += TMath::Sin(2*phiV0)*multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1463 } else { // V0A
1464 Qxa2 += TMath::Cos(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1465 Qya2 += TMath::Sin(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1466 }
1467 }
1468
1469 // Qx2 = Qxa2 + Qxc2;
1470 //Qy2 = Qya2 + Qya2;
1471
1472 //grab for each centrality the proper histo with the Qx and Qy
1473 //to do the recentering
1474 Double_t Qxamean2 = fMeanQ[fICent][1][0];
1475 Double_t Qxarms2 = fWidthQ[fICent][1][0];
1476 Double_t Qyamean2 = fMeanQ[fICent][1][1];
1477 Double_t Qyarms2 = fWidthQ[fICent][1][1];
1478
1479 Double_t Qxcmean2 = fMeanQ[fICent][0][0];
1480 Double_t Qxcrms2 = fWidthQ[fICent][0][0];
1481 Double_t Qycmean2 = fMeanQ[fICent][0][1];
1482 Double_t Qycrms2 = fWidthQ[fICent][0][1];
1483
1484 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1485 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1486 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1487 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1488
1489 Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
1490 /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
1491 Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
1492 /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
1493
1494 Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
1495 Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
1496 Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
1497
1498 fhEPangleVZero->Fill(psiVZero);
1499 fhEPangleV0A->Fill(psiV0A);
1500 fhEPangleV0C->Fill(psiV0C);
1501
1502 //fill profile for resolution estimation
1503 fProfResolution->Fill(1, TMath::Cos(2.*(psiV0A - psiV0C)));
1504 fProfResolution->Fill(2, TMath::Cos(2.*(psiV0A - psiTPC)));
1505 fProfResolution->Fill(3, TMath::Cos(2.*(psiV0C - psiTPC)));
1506
1507 Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
1508
1509 Double_t b = fAOD->GetMagneticField();
1510
1511 int nCascades=fAOD->GetNumberOfCascades();
1512 //Info("ReadFromAODv0", "# cascades = %d", nCascades);
1513
1514 const AliAODVertex *primaryBestAODVtx = fAOD->GetPrimaryVertex();
1515 primaryBestAODVtx->GetXYZ(bestPrimaryVtxPos);
1516
1517 for(Int_t iXi = 0; iXi < nCascades; iXi++){
1518 // Double_t effMassXi = 0.;
1519 Double_t chi2Xi = -1.;
1520 Double_t dcaXiDaughters = -1.;
1521 Double_t XiCosOfPointingAngle = -1.;
1522 Double_t posXi[3] = {-1000., -1000., -1000.};
1523 Double_t XiRadius = -1000.;
1524
1525 Double_t invMassLambdaAsCascDghter = 0.;
1526 Double_t V0Chi2Xi = -1.;
1527 Double_t dcaV0DaughtersXi = -1.;
1528
1529 Double_t dcaBachToPrimaryVtxXi = -1.;
1530 Double_t dcaV0ToPrimaryVtxXi = -1.;
1531 Double_t dcaPosToPrimaryVtxXi = -1.;
1532 Double_t dcaNegToPrimaryVtxXi = -1.;
1533 Double_t V0CosOfPointingAngleXi = -1.;
1534 Double_t posV0Xi[3] = {-1000., -1000., -1000.};
1535 Double_t V0RadiusXi = -1000.;
1536
1537 Double_t invMassXiMinus = 0.;
1538 Double_t invMassXiPlus = 0.;
1539 Double_t invMassOmegaMinus = 0.;
1540 Double_t invMassOmegaPlus = 0.;
1541
1542 /*
1543 Bool_t isPosInXiProton = kFALSE;
1544 Bool_t isPosInXiPion = kFALSE;
1545 Bool_t isPosInOmegaProton = kFALSE;
1546 Bool_t isPosInOmegaPion = kFALSE;
1547
1548 Bool_t isNegInXiProton = kFALSE;
1549 Bool_t isNegInXiPion = kFALSE;
1550 Bool_t isNegInOmegaProton = kFALSE;
1551 Bool_t isNegInOmegaPion = kFALSE;
1552
1553 Bool_t isBachelorKaon = kFALSE;
1554 Bool_t isBachelorPion = kFALSE;
1555 */
1556
1557 Bool_t isBachelorKaonForTPC = kFALSE;
1558 Bool_t isBachelorPionForTPC = kFALSE;
1559 Bool_t isNegPionForTPC = kFALSE;
1560 Bool_t isPosPionForTPC = kFALSE;
1561 Bool_t isNegProtonForTPC = kFALSE;
1562 Bool_t isPosProtonForTPC = kFALSE;
1563
1564 Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
1565 Double_t XiPt = 0.;
1566 Double_t XiPtot = 0.;
1567
1568 Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
1569 Double_t bachPt = 0.;
1570 Double_t bachPtot = 0.;
1571
1572 //Short_t chargeXi = -2;
1573 Double_t V0toXiCosOfPointingAngle = 0.;
1574
1575 Double_t rapXi = -20.;
1576 Double_t rapOmega = -20.;
1577 Double_t phi = 6.3;
1578 Double_t alphaXi = -200.;
1579 Double_t ptArmXi = -200.;
1580
1581 Double_t distToVtxZBefore = -999.;
1582 Double_t distToVtxZAfter = -999.;
1583 Double_t distToVtxXYBefore = -999.;
1584 Double_t distToVtxXYAfter = -999.;
1585 Double_t XiPAfter[3] = {-999., -999., -999.};
1586 Double_t phiAfter = -999.;
1587
1588 const AliAODcascade *xi = fAOD->GetCascade(iXi);
1589 if (!xi) continue;
1590
1591 //effMassXi = xi->MassXi(); //default working hypothesis: Xi- decay
1592 chi2Xi = xi->Chi2Xi();
1593 dcaXiDaughters = xi->DcaXiDaughters();
1594 XiCosOfPointingAngle = xi->CosPointingAngleXi(bestPrimaryVtxPos[0],
1595 bestPrimaryVtxPos[1],
1596 bestPrimaryVtxPos[2]);
1597 posXi[0] = xi->DecayVertexXiX();
1598 posXi[1] = xi->DecayVertexXiY();
1599 posXi[2] = xi->DecayVertexXiZ();
1600 XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
1601 +posXi[1]*posXi[1]
1602 +posXi[2]*posXi[2]);
1603
1604 AliAODTrack *pTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1605 AliAODTrack *nTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1606 AliAODTrack *bTrkXi
1607 = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1608
1609 if(!pTrkXi || !nTrkXi || !bTrkXi) continue;
1610
1611 UInt_t idxPosXi = (UInt_t) TMath::Abs( pTrkXi->GetID() );
1612 UInt_t idxNegXi = (UInt_t) TMath::Abs( nTrkXi->GetID() );
1613 UInt_t idxBach = (UInt_t) TMath::Abs( bTrkXi->GetID() );
1614
1615 if(idxBach == idxNegXi || idxBach == idxPosXi) continue;
1616
1617 if( !fCutsDau->IsSelected(pTrkXi)
1618 || !fCutsDau->IsSelected(nTrkXi)
1619 || !fCutsDau->IsSelected(bTrkXi) ) continue;
1620
1621 if(pTrkXi->IsOn(AliESDtrack::kTPCin)){
1622 fh2TPCdEdxOfCascDghters->Fill(pTrkXi->P()*pTrkXi->Charge(),
1623 pTrkXi->GetTPCsignal());
1624 }
1625 if( nTrkXi->IsOn(AliESDtrack::kTPCin) ){
1626 fh2TPCdEdxOfCascDghters->Fill(nTrkXi->P()*nTrkXi->Charge(),
1627 nTrkXi->GetTPCsignal());
1628 }
1629 if(bTrkXi->IsOn(AliESDtrack::kTPCin)){
1630 fh2TPCdEdxOfCascDghters->Fill(bTrkXi->P()*bTrkXi->Charge(),
1631 bTrkXi->GetTPCsignal());
1632 }
1633
1634
1635 if(xi->ChargeXi() < 0)
1636 invMassLambdaAsCascDghter = xi->MassLambda();
1637 else
1638 invMassLambdaAsCascDghter = xi->MassAntiLambda();
1639
1640 dcaV0DaughtersXi = xi->DcaV0Daughters();
1641 V0Chi2Xi = xi->Chi2V0();
1642
1643 V0CosOfPointingAngleXi
1644 = xi->CosPointingAngle(bestPrimaryVtxPos);
1645 dcaV0ToPrimaryVtxXi = xi->DcaV0ToPrimVertex();
1646 dcaBachToPrimaryVtxXi = xi->DcaBachToPrimVertex();
1647
1648 //V0
1649 posV0Xi[0] = xi->DecayVertexV0X();
1650 posV0Xi[1] = xi->DecayVertexV0Y();
1651 posV0Xi[2] = xi->DecayVertexV0Z();
1652 V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
1653 +posV0Xi[1]*posV0Xi[1]
1654 +posV0Xi[2]*posV0Xi[2]);
1655 dcaPosToPrimaryVtxXi = xi->DcaPosToPrimVertex();
1656 dcaNegToPrimaryVtxXi = xi->DcaNegToPrimVertex();
1657
1658
1659 //apply cuts ?
1660 // if(XiRadius < 1. || XiRadius > 100.) continue;
1661 //if(dcaXiDaughters > 0.1) continue;
1662 //if(XiCosOfPointingAngle < 0.999) continue;
1663 //if(dcaV0ToPrimaryVtxXi < 0.05) continue;
1664 //if(dcaBachToPrimaryVtxXi < 0.03) continue;
1665
1666 //V0 mass cut?
1667 if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.01) continue;
1668
1669 // if(dcaV0DaughtersXi > 1.) continue;
1670 //if(V0CosOfPointingAngleXi > 0.9999) continue;
1671 //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1672 //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1673
1674 //if(V0RadiusXi < 1.0 || V0RadiusXi > 100) continue;
1675
1676 //other cuts?
1677
1678 if(xi->ChargeXi()<0){
1679 invMassXiMinus = xi->MassXi();
1680 invMassOmegaMinus = xi->MassOmega();
1681 }else{
1682 invMassXiPlus = xi->MassXi();
1683 invMassOmegaPlus = xi->MassOmega();
1684 }
1685
1686 /*
1687 if(pTrkXi->GetMostProbablePID() == AliAODTrack::kProton) {
1688 isPosInXiProton = kTRUE;
1689 isPosInOmegaProton = kTRUE;
1690 }
1691 if(pTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1692 isPosInXiPion = kTRUE;
1693 isPosInOmegaPion = kTRUE;
1694 }
1695
1696 if(nTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1697 isNegInXiPion = kTRUE;
1698 isNegInOmegaPion = kTRUE;
1699 }
1700 if(nTrkXi->GetMostProbablePID() == AliAODTrack::kProton){
1701 isNegInXiProton = kTRUE;
1702 isNegInOmegaProton = kTRUE;
1703 }
1704
1705 if(bTrkXi->GetMostProbablePID() == AliAODTrack::kPion)
1706 isBachelorPion = kTRUE;
1707 if(bTrkXi->GetMostProbablePID() == AliAODTrack::kKaon)
1708 isBachelorKaon = kTRUE;
1709 */
1710
1711 //PID with TPC only: ??? Fix me
1712 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
1713 isBachelorKaonForTPC = kTRUE;
1714 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
1715 isBachelorPionForTPC = kTRUE;
1716
1717 //Negative V0 daughter
1718 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
1719 isNegPionForTPC = kTRUE;
1720 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
1721 isNegProtonForTPC = kTRUE;
1722
1723 //Positive V0 daughter
1724 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
1725 isPosPionForTPC = kTRUE;
1726 if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
1727 isPosProtonForTPC = kTRUE;
1728
1729 //Extra QA information
1730 XiPx = xi->MomXiX();
1731 XiPy = xi->MomXiY();
1732 XiPz = xi->MomXiZ();
1733 XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
1734 XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
1735
1736 bachPx = xi->MomBachX();
1737 bachPy = xi->MomBachY();
1738 bachPz = xi->MomBachZ();
1739
1740 bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1741 bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1742
1743 V0toXiCosOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1744
1745 rapXi = xi->RapXi();
1746 rapOmega = xi->RapOmega();
1747 phi = xi->Phi();
1748 alphaXi = xi->AlphaXi();
1749 ptArmXi = xi->PtArmXi();
1750
1751 distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1752 distToVtxXYBefore
1753 = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1754 *(posXi[0] - bestPrimaryVtxPos[0])
1755 +(posXi[1] - bestPrimaryVtxPos[1])
1756 *(posXi[1] - bestPrimaryVtxPos[1]));
1757
1758 XiPAfter[0] = XiPx;
1759 XiPAfter[1] = XiPy;
1760 XiPAfter[2] = XiPz;
1761
1762 //propagation to the best primary vertex to determine the momentum
1763 Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->ChargeXi());
1764 distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
1765 distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1766 *(posXi[0] - bestPrimaryVtxPos[0])
1767 +(posXi[1] - bestPrimaryVtxPos[1])
1768 *(posXi[1] - bestPrimaryVtxPos[1]));
1769 phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
1770
1771 fh1DistToVtxZAfter->Fill(distToVtxZAfter);
1772 fh1DistToVtxXYAfter->Fill(distToVtxXYAfter);
1773 fh2DistToVtxZBeforeVsAfter->Fill(distToVtxZBefore, distToVtxZAfter);
1774 fh2DistToVtxXYBeforeVsAfter->Fill(distToVtxXYBefore, distToVtxXYAfter);
1775 fh2PxBeforeVsAfter->Fill(XiPx, XiPAfter[0]);
1776 fh2PyBeforeVsAfter->Fill(XiPy, XiPAfter[1]);
1777 if(xi->ChargeXi()>0)
1778 fh2PhiPosBeforeVsAfter->Fill(phi, phiAfter);
1779 else if(xi->ChargeXi()<0)
1780 fh2PhiNegBeforeVsAfter->Fill(phi, phiAfter);
1781
1782 if( (xi->ChargeXi() < 0 && isBachelorPionForTPC
1783 && isPosProtonForTPC && isNegPionForTPC)
1784 || (xi->ChargeXi() > 0 && isBachelorPionForTPC
1785 && isNegProtonForTPC && isPosPionForTPC))
1786 {//Xi candidate
1787 //for default hypothesis
1788 fh1Chi2Xi->Fill(chi2Xi);
1789 fh1DCAXiDaughters->Fill(dcaXiDaughters);
1790 fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1791 fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1792 fh1XiRadius->Fill(XiRadius);
1793
1794 //V0
1795 fh1MassLambda->Fill(invMassLambdaAsCascDghter);
1796 fh1V0Chi2->Fill(V0Chi2Xi);
1797 fh1DcaV0DaughtersXi->Fill(dcaV0DaughtersXi);
1798 fh1V0CosOfPointingAngle->Fill(V0CosOfPointingAngleXi);
1799 fh1V0Radius->Fill(V0RadiusXi);
1800 fh1DcaV0ToPrimVertex->Fill(dcaV0ToPrimaryVtxXi);
1801 fh1DCAPosToPrimVertex->Fill(dcaPosToPrimaryVtxXi);
1802 fh1DCANegToPrimVertex->Fill(dcaNegToPrimaryVtxXi);
1803 fh1ChargeXi->Fill(xi->ChargeXi());
1804 fh1V0toXiCosOfPointingAngle->Fill(V0toXiCosOfPointingAngle);
1805
1806 if ( TMath::Abs( invMassXiMinus-1.3217 ) < 0.012
1807 || TMath::Abs( invMassXiPlus-1.3217 ) < 0.012)
1808 {// One InvMass should be different from 0
1809 fh1XiPt->Fill(XiPt);
1810 fh1XiP->Fill(XiPtot);
1811
1812 fh1XiBachPt->Fill(bachPt);
1813 fh1XiBachP->Fill(bachPtot);
1814 fh1PhiXi->Fill( xi->Phi() );
1815 }
1816 fh2Armenteros->Fill(alphaXi, ptArmXi);
1817 }
1818
1819 if(xi->ChargeXi()<0){
1820 fh1MassXiMinus->Fill(invMassXiMinus);
1821 fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1822 fh1MassXi->Fill(invMassXiMinus);
1823 fh1MassOmega->Fill(invMassOmegaMinus);
1824
1825 fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1826 invMassXiMinus);
1827 fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus,
1828 invMassOmegaMinus);
1829 fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1830 fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1831 }
1832
1833 if(xi->ChargeXi() > 0){
1834 fh1MassXiPlus->Fill(invMassXiPlus);
1835 fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1836 fh1MassXi->Fill(invMassXiPlus);
1837 fh1MassOmega->Fill(invMassOmegaPlus);
1838
1839 fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter,
1840 invMassXiPlus);
1841 fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1842 invMassOmegaPlus);
1843 fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1844 fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1845 }
1846
1847
1848 // Double_t phiNew
1849 // = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. )
1850 // ? 0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1851 Double_t phiV0 = phiAfter;
1852 phiV0 -= psiVZero;
1853 //if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1854 //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1855
1856
1857 Double_t phiV0A = phiAfter;
1858 phiV0A -= psiV0A;
1859 //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1860 //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1861
1862 Double_t phiV0C = phiAfter;
1863 phiV0C -= psiV0C;
1864 //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1865 //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1866
1867 //PID cuts with TPC cuts
1868 if(xi->ChargeXi() < 0){
1869 if(isPosProtonForTPC
1870 && isNegPionForTPC){
1871 if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1872 //Xi
1873 fh2MassVsPtXiMinus->Fill(XiPt, invMassXiMinus);
1874 fh2MassVsPtXiAll->Fill(XiPt, invMassXiMinus);
1875 fhXiRapidity->Fill(rapXi);
1876 for(int r=0; r!=3; ++r) {
1877 if(invMassXiMinus > fXiBands[r][0]
1878 && invMassXiMinus < fXiBands[r][1]){
1879
1880 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1881 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1882 fProf2dXiV2PtV0A[r]->Fill(XiPt, psiV0A, TMath::Cos(2.*phiV0A));
1883
1884 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1885 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1886 fProf2dXiV2PtV0C[r]->Fill(XiPt, psiV0C, TMath::Cos(2.*phiV0C));
1887
1888 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1889 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1890 fProf2dXiV2Pt[r]->Fill(XiPt, psiVZero, TMath::Cos(2.*phiV0));
1891 }
1892 }
1893 }
1894
1895 if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1896 fh2MassVsPtOmegaMinus->Fill(XiPt, invMassOmegaMinus);
1897 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaMinus);
1898 fhOmegaRapidity->Fill(rapOmega);
1899 for(int r=0; r!=3; ++r) {
1900 if(invMassOmegaMinus > fOmegaBands[r][0]
1901 && invMassOmegaMinus < fOmegaBands[r][1]){
1902 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1903 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1904 fProf2dOmegaV2PtV0A[r]->Fill(XiPt, psiV0A,
1905 TMath::Cos(2.*phiV0A));
1906
1907 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1908 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1909 fProf2dOmegaV2PtV0C[r]->Fill(XiPt, psiV0C,
1910 TMath::Cos(2.*phiV0C));
1911
1912 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1913 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1914 fProf2dOmegaV2Pt[r]->Fill(XiPt, psiVZero,
1915 TMath::Cos(2.*phiV0));
1916 }
1917 }
1918 }
1919 }
1920 }//if charge < 0
1921
1922 if(xi->ChargeXi() > 0){
1923 if(isNegProtonForTPC
1924 && isPosPionForTPC){
1925 if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1926 //Xi
1927 fh2MassVsPtXiPlus->Fill(XiPt, invMassXiPlus);
1928 fh2MassVsPtXiAll->Fill(XiPt, invMassXiPlus);
1929 fhXiRapidity->Fill(rapXi);
1930 for(int r=0; r!=3; ++r) {
1931 if(invMassXiPlus > fXiBands[r][0]
1932 && invMassXiPlus < fXiBands[r][1]){
1933 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1934 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1935 fProf2dXiV2PtV0A[r]->Fill(XiPt, psiV0A, TMath::Cos(2.*phiV0A));
1936
1937 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1938 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1939 fProf2dXiV2PtV0C[r]->Fill(XiPt, psiV0C, TMath::Cos(2.*phiV0C));
1940
1941 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1942 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1943 fProf2dXiV2Pt[r]->Fill(XiPt, psiVZero, TMath::Cos(2.*phiV0));
1944 }
1945 }
1946 }
1947
1948 if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1949 //Omega
1950 fh2MassVsPtOmegaPlus->Fill(XiPt, invMassOmegaPlus);
1951 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaPlus);
1952 fhOmegaRapidity->Fill(rapOmega);
1953 for(int r=0; r!=3; ++r) {
1954 if(invMassOmegaPlus > fOmegaBands[r][0]
1955 && invMassOmegaPlus < fOmegaBands[r][1]){
1956 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1957 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1958 fProf2dOmegaV2PtV0A[r]->Fill(XiPt, psiV0A,
1959 TMath::Cos(2.*phiV0A));
1960
1961
1962 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1963 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1964 fProf2dOmegaV2PtV0C[r]->Fill(XiPt, psiV0C,
1965 TMath::Cos(2.*phiV0C));
1966
1967 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1968 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1969 fProf2dOmegaV2Pt[r]->Fill(XiPt, psiVZero,
1970 TMath::Cos(2.*phiV0));
1971 }
1972 }
1973 }
1974 }
1975 }// if charge > 0
1976
1977 }//end of cascade loop
1978
1979 PostData(1, fHistList);
1980
1981}
1982
1983
1984AliAnalysisTaskFlowEPCascade::~AliAnalysisTaskFlowEPCascade(){
1985 if(fHistList) delete fHistList;
1986}
1987
1988void AliAnalysisTaskFlowEPCascade::Terminate(Option_t *)
1989{
1990 /*
1991 fHistList = dynamic_cast<TList*> (GetOutputData(1));
1992
1993 if (!fHistList) {
1994 printf("ERROR: Output tree not available\n");
1995 return;
1996 }
1997 */
1998}
1999
2000
2001void AliAnalysisTaskFlowEPCascade::Propagate(Double_t vv[3],
2002 Double_t x[3],
2003 Double_t p[3],
2004 Double_t bz,
2005 Short_t sign){
2006 //Propagation to the primary vertex to determine the px and py
2007 //x, p are the position and momentum as input and output
2008 //bz is the magnetic field along z direction
2009 //sign is the charge of particle for propagation
2010
2011 Double_t pp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
2012 Double_t len = (vv[2]-x[2])*pp/p[2];
2013 Double_t a = -kB2C*bz*sign;
2014
2015 Double_t rho = a/pp;
2016 x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
2017 x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
2018 x[2] += p[2]*len/pp;
2019
2020 Double_t p0=p[0];
2021 p[0] = p0 *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
2022 p[1] = p[1]*TMath::Cos(rho*len) + p0 *TMath::Sin(rho*len);
2023}
2024
2025
2026void
2027AliAnalysisTaskFlowEPCascade::OpenInfoCalbration(Int_t run){
2028
2029 AliOADBContainer *cont = (AliOADBContainer*) fOADB->Get("hMultV0BefCorr");
2030 if(!cont){
2031 printf("OADB object hMultV0BefCorr is not available in the file\n");
2032 return;
2033 }
2034
2035 if(!(cont->GetObject(run))){
2036 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
2037 run = 137366;
2038 }
2039 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
2040
2041 TF1 *fpol0 = new TF1("fpol0","pol0");
2042 fMultV0->Fit(fpol0,"","",0,32);
2043 fV0Cpol = fpol0->GetParameter(0);
2044 fMultV0->Fit(fpol0,"","",32,64);
2045 fV0Apol = fpol0->GetParameter(0);
2046
2047 for(Int_t iside=0;iside<2;iside++){
2048 for(Int_t icoord=0;icoord<2;icoord++){
2049 for(Int_t i=0;i < 9;i++){
2050 char namecont[100];
2051 if(iside==0 && icoord==0)
2052 snprintf(namecont,100,"hQxc2_%i", i);
2053 else if(iside==1 && icoord==0)
2054 snprintf(namecont,100,"hQxa2_%i", i);
2055 else if(iside==0 && icoord==1)
2056 snprintf(namecont,100,"hQyc2_%i", i);
2057 else if(iside==1 && icoord==1)
2058 snprintf(namecont,100,"hQya2_%i", i);
2059
2060 cont = (AliOADBContainer*) fOADB->Get(namecont);
2061 if(!cont){
2062 printf("OADB object %s is not available in the file\n",namecont);
2063 return;
2064 }
2065
2066 if(!(cont->GetObject(run))){
2067 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
2068 run = 137366;
2069 }
2070 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2071 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2072 }
2073 }
2074 }
2075}