]> 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
Enlarging window for DCS DPs retrieval for short runs for GRP + Keeping connection...
[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
0a918d8d 771 AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
772 if(!aodHeader) AliFatal("Not a standard AOD");
47976489 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/*
802void 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
1435void AliAnalysisTaskFlowEPCascade::ReadFromAODv0(AliAODEvent *fAOD){
1436
0a918d8d 1437 AliEventplane * ep = ((AliVAODHeader*)fAOD->GetHeader())->GetEventplaneP();
47976489 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
1985AliAnalysisTaskFlowEPCascade::~AliAnalysisTaskFlowEPCascade(){
1986 if(fHistList) delete fHistList;
1987}
1988
1989void 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
2002void 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
2027void
2028AliAnalysisTaskFlowEPCascade::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}