11 #include "TProfile2D.h"
14 #include "AliAnalysisTaskSE.h"
15 #include "AliAnalysisManager.h"
17 #include "AliESDEvent.h"
19 #include "AliESDcascade.h"
20 #include "AliVVertex.h"
22 #include "AliESDtrack.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliESDVZERO.h"
25 #include "AliESDUtils.h"
27 #include "AliAODEvent.h"
29 #include "AliESDInputHandler.h"
30 #include "AliCentrality.h"
31 #include "AliMultiplicity.h"
33 #include "AliFlowTrackSimple.h"
34 //#include "AliFlowEventCuts.h"
35 #include "AliFlowTrackCuts.h"
37 #include "AliTPCPIDResponse.h"
38 #include "AliTOFPIDResponse.h"
39 #include "AliPIDResponse.h"
41 #include "AliOADBContainer.h"
43 #include "AliAnalysisTaskFlowEPCascade.h"
51 ClassImp(AliAnalysisTaskFlowEPCascade)
53 //_______________________________________________________________
54 AliAnalysisTaskFlowEPCascade::AliAnalysisTaskFlowEPCascade():
74 ,fh1DCAXiDaughters(0x0)
75 ,fh1DCABachToPrimVertex(0x0)
76 ,fh1XiCosOfPointingAngle(0x0)
80 ,fh1V0CosOfPointingAngle(0x0)
82 ,fh1DcaV0DaughtersXi(0x0)
83 ,fh1DcaV0ToPrimVertex(0x0)
84 ,fh1DCAPosToPrimVertex(0x0)
85 ,fh1DCANegToPrimVertex(0x0)
88 ,fh1MassOmegaMinus(0x0)
89 ,fh1MassOmegaPlus(0x0)
97 ,fh1V0toXiCosOfPointingAngle(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)
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)
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;
134 fProfXiV2PtV0C[i] = NULL;
135 fProfOmegaV2PtV0C[i] = NULL;
136 fProfXiSinePtV0C[i] = NULL;
137 fProfOmegaSinePtV0C[i] = NULL;
139 fProfXiV2Pt[i] = NULL;
140 fProfOmegaV2Pt[i] = NULL;
141 fProfXiSinePt[i] = NULL;
142 fProfOmegaSinePt[i] = NULL;
144 fProf2dXiV2PtV0A[i] = NULL;
145 fProf2dOmegaV2PtV0A[i] = NULL;
146 fProf2dXiV2PtV0C[i] = NULL;
147 fProf2dOmegaV2PtV0C[i] = NULL;
148 fProf2dXiV2Pt[i] = NULL;
149 fProf2dOmegaV2Pt[i] = NULL;
152 for(int i=0; i!=3; ++i)
153 for(int j=0; j!=2; ++j) {
155 fOmegaBands[i][j] = 0;
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.;
167 //_________________________________________________________________________
168 AliAnalysisTaskFlowEPCascade::
169 AliAnalysisTaskFlowEPCascade(const char *name, double centMin, double centMax,
172 : AliAnalysisTaskSE(name)
186 ,fhEPangleVZero (0x0)
191 ,fh1DCAXiDaughters(0x0)
192 ,fh1DCABachToPrimVertex(0x0)
193 ,fh1XiCosOfPointingAngle(0x0)
197 ,fh1V0CosOfPointingAngle(0x0)
199 ,fh1DcaV0DaughtersXi(0x0)
200 ,fh1DcaV0ToPrimVertex(0x0)
201 ,fh1DCAPosToPrimVertex(0x0)
202 ,fh1DCANegToPrimVertex(0x0)
205 ,fh1MassOmegaMinus(0x0)
206 ,fh1MassOmegaPlus(0x0)
214 ,fh1V0toXiCosOfPointingAngle(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)
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)
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;
251 fProfXiV2PtV0C[i] = NULL;
252 fProfOmegaV2PtV0C[i] = NULL;
253 fProfXiSinePtV0C[i] = NULL;
254 fProfOmegaSinePtV0C[i] = NULL;
256 fProfXiV2Pt[i] = NULL;
257 fProfOmegaV2Pt[i] = NULL;
258 fProfXiSinePt[i] = NULL;
259 fProfOmegaSinePt[i] = NULL;
261 fProf2dXiV2PtV0A[i] = NULL;
262 fProf2dOmegaV2PtV0A[i] = NULL;
263 fProf2dXiV2PtV0C[i] = NULL;
264 fProf2dOmegaV2PtV0C[i] = NULL;
265 fProf2dXiV2Pt[i] = NULL;
266 fProf2dOmegaV2Pt[i] = NULL;
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];
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.;
283 DefineInput( 0,TChain::Class());
284 //DefineInput (1, TList::Class());
285 DefineOutput(1, TList::Class());
288 void AliAnalysisTaskFlowEPCascade::UserCreateOutputObjects()
290 cout<<"AliAnalysisTaskFlowEPCascade::UserCreateOutputObjects()"<<endl;
292 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
293 AliInputEventHandler* inputHandler
294 = (AliInputEventHandler*) (man->GetInputEventHandler());
295 fPIDResponse = inputHandler->GetPIDResponse();
297 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
298 fOADB = TFile::Open(oadbfilename.Data());
301 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
305 fHistList = new TList();
306 fHistList ->SetOwner();
308 //Add histograms to the List
309 fhEvent = new TH1I("Event", "Number of Events", 3, 0, 3);
310 fHistList->Add(fhEvent);
312 fhEPangleVZero = new TH1F("hEPangleVZero",
313 "EP from VZERO; #Psi; Number of Events",
314 15, 0., TMath::Pi());
315 fHistList->Add(fhEPangleVZero);
317 fhEPangleV0A = new TH1F("hEPangleV0A",
318 "EP from V0A; #Psi; Number of Events",
319 15, 0., TMath::Pi());
320 fHistList->Add(fhEPangleV0A);
322 fhEPangleV0C = new TH1F("hEPangleV0C",
323 "EP from V0C; #Psi; Number of Events",
324 15, 0., TMath::Pi());
325 fHistList->Add(fhEPangleV0C);
328 fhEPangleTPC = new TH1F("hEPangleTPC",
329 "EP from TPC; #Psi; Number of Events",
330 15, 0., TMath::Pi());
331 fHistList->Add(fhEPangleTPC);
333 fh1Chi2Xi = new TH1F("Chi2Xi",
334 "Cascade #chi^{2}; #chi^{2}; Number of Cascades",
336 fHistList->Add(fh1Chi2Xi);
339 = new TH1F( "DcaXiDaughters",
340 "DCA between Xi Daughters; DCA (cm); Number of Cascades",
342 fHistList->Add(fh1DCAXiDaughters);
344 fh1DCABachToPrimVertex
345 = new TH1F("DcaBachToPrimVertex",
346 "DCA of Bach. to Prim. Vertex; DCA (cm);Number of Cascades",
348 fHistList->Add(fh1DCABachToPrimVertex);
350 fh1XiCosOfPointingAngle
351 = new TH1F("XiCosineOfPointingAngle",
352 "Cos of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis",
354 fHistList->Add(fh1XiCosOfPointingAngle);
356 fh1XiRadius = new TH1F("XiRadius",
357 "Casc. decay transv. radius; r (cm); Counts" ,
359 fHistList->Add(fh1XiRadius);
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);
366 fh1V0Chi2 = new TH1F("V0Chi2Xi",
367 "V0 #chi^{2}, in cascade; #chi^{2};Counts",
369 fHistList->Add(fh1V0Chi2);
371 fh1V0CosOfPointingAngle
372 = new TH1F("V0CosOfPointingAngleXi",
373 "Cos of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts",
375 fHistList->Add(fh1V0CosOfPointingAngle);
377 fh1V0Radius = new TH1F("V0RadiusXi",
378 "V0 decay radius, in cascade; radius (cm); Counts",
380 fHistList->Add(fh1V0Radius);
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);
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);
390 fh1DCAPosToPrimVertex =
391 new TH1F("DcaPosToPrimVertexXi",
392 "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts",
394 fHistList->Add(fh1DCAPosToPrimVertex);
396 fh1DCANegToPrimVertex
397 = new TH1F("DcaNegToPrimVertexXi",
398 "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts",
400 fHistList->Add(fh1DCANegToPrimVertex);
403 = new TH1F("MassXiMinus",
404 "#Xi^{-} candidates;M( #Lambda , #pi^{-} ) (GeV/c^{2});Counts",
406 fHistList->Add(fh1MassXiMinus);
410 = new TH1F("MassXiPlus",
411 "#Xi^{+} candidates;M(#bar{#Lambda}^{0}, #pi^{+}) (GeV/c^{2});Counts",
413 fHistList->Add(fh1MassXiPlus);
417 "#Xi candidates;M(#bar{#Lambda}, #pi) (GeV/c^{2});Counts",
419 fHistList->Add(fh1MassXi);
422 = new TH1F("MassOmegaMinus",
423 "#Omega^{-} candidates; M(#Lambda, K^{-})(GeV/c^{2});Counts",
425 fHistList->Add(fh1MassOmegaMinus);
428 = new TH1F("MassOmega",
429 "#Omega candidates; M(#Lambda, K)(GeV/c^{2});Counts",
431 fHistList->Add(fh1MassOmega);
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);
440 "#Xi Pt (cand. around the mass peak);p_{t}(#Xi)(GeV/c);Counts",
442 fHistList->Add(fh1XiPt);
445 = new TH1F("XiTotMom",
446 "#Xi momentum (cand. around the mass peak); p_{tot}(#Xi)(GeV/c); Counts",
448 fHistList->Add(fh1XiP);
450 fh1XiBachPt = new TH1F("XiBachPt",
451 "#Xi Bach. transverse momentum (cand. around the mass peak) ; p_{t}(Bach.) (GeV/c); Counts",
453 fHistList->Add(fh1XiBachPt);
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);
459 fh1ChargeXi = new TH1F("ChargeXi",
460 "Charge of casc. candidates; Sign; Counts",
462 fHistList->Add(fh1ChargeXi);
464 fh1V0toXiCosOfPointingAngle
465 = new TH1F("V0toXiCosineOfPointingAngle",
466 "Cos. of V0 Ptng Angl Xi vtx; Cos(V0 Point. Angl / Xi vtx); Counts",
468 fHistList->Add(fh1V0toXiCosOfPointingAngle);
472 "#phi of #Xi candidates (around the mass peak); #phi (deg); Counts", 64, 0., 6.4);
473 fHistList->Add(fh1PhiXi);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
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);
571 fhXiRapidity = new TH1F("hXiRapidity",
572 "#Xi rapidity distribution before rap. cut; y; Number of counts",
574 fHistList->Add(fhXiRapidity);
576 fhOmegaRapidity = new TH1F("hOmegaRapidity",
577 "#Omega rapidity distr. before rap. cut; y; Number of counts",
579 fHistList->Add(fhOmegaRapidity);
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}",
585 fHistList->Add(fProfXiV2PtV0A[i]);
587 fProfOmegaV2PtV0A[i] = new TProfile(Form("hProfOmegaV2PtV0A%d", i),
588 "; p_{T}[GeV/c]; v_{2}",
590 fHistList->Add(fProfOmegaV2PtV0A[i]);
592 fProfXiSinePtV0A[i] = new TProfile(Form("hProfXiSinePtV0A%d", i),
593 ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
595 fHistList->Add(fProfXiSinePtV0A[i]);
597 fProfOmegaSinePtV0A[i]
598 = new TProfile(Form("hProfOmegaSinePtV0A%d", i),
599 "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
601 fHistList->Add(fProfOmegaSinePtV0A[i]);
604 fProfXiV2PtV0C[i] = new TProfile(Form("hProfXiV2PtV0C%d", i),
605 "; p_{T}[GeV/c]; v_{2}",
607 fHistList->Add(fProfXiV2PtV0C[i]);
609 fProfOmegaV2PtV0C[i] = new TProfile(Form("hProfOmegaV2PtV0C%d", i),
610 "; p_{T}[GeV/c]; v_{2}",
612 fHistList->Add(fProfOmegaV2PtV0C[i]);
614 fProfXiSinePtV0C[i] = new TProfile(Form("hProfXiSinePtV0C%d", i),
615 ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
617 fHistList->Add(fProfXiSinePtV0C[i]);
619 fProfOmegaSinePtV0C[i]
620 = new TProfile(Form("hProfOmegaSinePtV0C%d", i),
621 "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
623 fHistList->Add(fProfOmegaSinePtV0C[i]);
626 fProfXiV2Pt[i] = new TProfile(Form("hProfXiV2Pt%d", i),
627 "; p_{T}[GeV/c]; v_{2}",
629 fHistList->Add(fProfXiV2Pt[i]);
631 fProfOmegaV2Pt[i] = new TProfile(Form("hProfOmegaV2Pt%d", i),
632 "; p_{T}[GeV/c]; v_{2}",
634 fHistList->Add(fProfOmegaV2Pt[i]);
636 fProfXiSinePt[i] = new TProfile(Form("hProfXiSinePt%d", i),
637 ";p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
639 fHistList->Add(fProfXiSinePt[i]);
641 fProfOmegaSinePt[i] = new TProfile(Form("hProfOmegaSinePt%d", i),
642 "; p_{T}[GeV/c]; <sin[2(#phi-#Psi)]>",
644 fHistList->Add(fProfOmegaSinePt[i]);
646 fProf2dXiV2PtV0A[i] = new TProfile2D(Form("hProf2dXiV2PtV0A%d", i),
647 "; p_{T}[GeV/c]; #Psi; v_{2}",
649 15, 0., TMath::Pi());
650 fHistList->Add(fProf2dXiV2PtV0A[i]);
652 fProf2dOmegaV2PtV0A[i] = new TProfile2D(Form("hProf2dOmegaV2PtV0A%d", i),
653 "; p_{T}[GeV/c]; #Psi; v_{2}",
657 fHistList->Add(fProf2dOmegaV2PtV0A[i]);
659 fProf2dXiV2PtV0C[i] = new TProfile2D(Form("hProf2dXiV2PtV0C%d", i),
660 "; p_{T}[GeV/c]; #Psi; v_{2}",
662 15, 0., TMath::Pi());
663 fHistList->Add(fProf2dXiV2PtV0C[i]);
665 fProf2dOmegaV2PtV0C[i] = new TProfile2D(Form("hProf2dOmegaV2PtV0C%d", i),
666 "; p_{T}[GeV/c]; #Psi; v_{2}",
670 fHistList->Add(fProf2dOmegaV2PtV0C[i]);
672 fProf2dXiV2Pt[i] = new TProfile2D(Form("hProf2dXiV2Pt%d", i),
673 "; p_{T}[GeV/c]; #Psi; v_{2}",
675 15, 0., TMath::Pi());
676 fHistList->Add(fProf2dXiV2Pt[i]);
678 fProf2dOmegaV2Pt[i] = new TProfile2D(Form("hProf2dOmegaV2Pt%d", i),
679 "; p_{T}[GeV/c]; #Psi; v_{2}",
683 fHistList->Add(fProf2dOmegaV2Pt[i]);
686 fh1DistToVtxZAfter = new TH1F("DistToVtxZAfter",
687 "Distance to vtx z after propagation to vtx; z [cm]",
689 fHistList->Add(fh1DistToVtxZAfter);
691 fh1DistToVtxXYAfter = new TH1F("DistToVtxXYAfter",
692 "Distance to vtx xy after propagation to vtx",
694 fHistList->Add(fh1DistToVtxXYAfter);
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);
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);
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);
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);
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);
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);
732 fProfResolution = new TProfile("hProfResolution",
733 "correlations between event planes",
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);
745 PostData(1, fHistList);
748 void AliAnalysisTaskFlowEPCascade::UserExec(Option_t *)
751 AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
752 AliAODEvent *fAOD=dynamic_cast<AliAODEvent*>(InputEvent());
756 Info("UserExec", "This task doesn't work with ESD!");
759 // if (!fFlowEventCuts->IsSelected(InputEvent())) return;
762 // ReadFromESDv0(fESD);
765 //routine for AOD info
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
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;
790 Double_t zvtx = fAOD->GetPrimaryVertex()->GetZ();
791 if(TMath::Abs(zvtx) > fVtxCut) return; //vertex cut
801 void AliAnalysisTaskFlowEPCascade::ReadFromESDv0(AliESDEvent *fESD){
803 // fEPcalib->CalculateEventPlanes(fESD);
805 //Get EP informations
806 //Double_t psiVZero = fEPcalib->GetPsi(2, AliAnaEPcalib::kV0AC);
807 //fhEPangleVZero->Fill(psiVZero);
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);
813 // Double_t psiZDC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCAC);
814 //fhEPangleZDC->Fill(psiZDC);
816 //Double_t psiZDCA = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCA);
817 //Double_t psiZDCC = fEPcalib->GetPsi(1, AliAnaEPcalib::kZDCC);
819 AliEventplane * ep = fESD->GetEventplane();
820 Double_t psiTPC = ep->GetEventplane("Q", fESD, 2);
822 Int_t run = fESD->GetRunNumber();
824 // Load the calibrations run dependent
825 OpenInfoCalbration(run);
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)));
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();
840 const AliESDVertex *primaryTrackingESDVtx = fESD->GetPrimaryVertexTracks();
841 primaryTrackingESDVtx->GetXYZ(trkPrimaryVtxPos);
843 const AliESDVertex *primaryBestESDVtx = fESD->GetPrimaryVertex();
844 primaryBestESDVtx->GetXYZ(bestPrimaryVtxPos);
846 Double_t b = fESD->GetMagneticField();
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.;
856 Double_t invMassLambdaAsCascDghter = 0.;
857 Double_t V0Chi2Xi = -1.;
858 Double_t dcaV0DaughtersXi = -1.;
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.;
869 Double_t invMassXiMinus = 0.;
870 Double_t invMassXiPlus = 0.;
871 Double_t invMassOmegaMinus = 0.;
872 Double_t invMassOmegaPlus = 0.;
874 Bool_t isPosInXiProton = kFALSE;
875 Bool_t isPosInXiPion = kFALSE;
876 Bool_t isPosInOmegaProton = kFALSE;
877 Bool_t isPosInOmegaPion = kFALSE;
879 Bool_t isNegInXiProton = kFALSE;
880 Bool_t isNegInXiPion = kFALSE;
881 Bool_t isNegInOmegaProton = kFALSE;
882 Bool_t isNegInOmegaPion = kFALSE;
884 Bool_t isBachelorKaon = kFALSE;
885 Bool_t isBachelorPion = kFALSE;
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;
894 Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
896 Double_t XiPtot = 0.;
898 Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
899 Double_t bachPt = 0.;
900 Double_t bachPtot = 0.;
902 //Short_t chargeXi = -2;
903 Double_t V0toXiCosOfPointingAngle = 0.;
905 Double_t rapXi = -20.;
906 Double_t rapOmega = -20.;
908 Double_t alphaXi = -200.;
909 Double_t ptArmXi = -200.;
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.;
918 AliESDcascade *xi = fESD->GetCascade(i);
922 xi->ChangeMassHypothesis(V0quality, 3312); // Xi- hypothesis
923 else if(xi->Charge() > 0)
924 xi->ChangeMassHypothesis(V0quality, -3312);
927 effMassXi = xi->GetEffMassXi();
928 chi2Xi = xi->GetChi2Xi();
929 dcaXiDaughters = xi->GetDcaXiDaughters();
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]
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());
943 if(idxBach == idxPosXi || idxBach == idxNegXi) continue;
945 AliESDtrack *pTrkXi = fESD->GetTrack(idxPosXi);
946 AliESDtrack *nTrkXi = fESD->GetTrack(idxNegXi);
947 AliESDtrack *bTrkXi = fESD->GetTrack(idxBach);
949 if( !pTrkXi || !nTrkXi || !bTrkXi ) continue;
951 if( !fCutsDau->IsSelected(pTrkXi)
952 || !fCutsDau->IsSelected(nTrkXi)
953 || !fCutsDau->IsSelected(bTrkXi) ) continue;
955 const AliExternalTrackParam *pExtTrk = pTrkXi->GetInnerParam();
956 const AliExternalTrackParam *nExtTrk = nTrkXi->GetInnerParam();
957 const AliExternalTrackParam *bExtTrk = bTrkXi->GetInnerParam();
959 if(pExtTrk && pTrkXi->IsOn(AliESDtrack::kTPCin)){
960 fh2TPCdEdxOfCascDghters->Fill(pExtTrk->GetP()*pExtTrk->Charge(),
961 pTrkXi->GetTPCsignal());
963 if(nExtTrk && nTrkXi->IsOn(AliESDtrack::kTPCin)){
964 fh2TPCdEdxOfCascDghters->Fill(nExtTrk->GetP()*nExtTrk->Charge(),
965 nTrkXi->GetTPCsignal());
967 if(bExtTrk && bTrkXi->IsOn(AliESDtrack::kTPCin)){
968 fh2TPCdEdxOfCascDghters->Fill(bExtTrk->GetP()*bExtTrk->Charge(),
969 bTrkXi->GetTPCsignal());
973 invMassLambdaAsCascDghter = xi->GetEffMass(); // from V0
974 dcaV0DaughtersXi = xi->GetDcaV0Daughters();
975 V0Chi2Xi = xi->GetChi2V0();
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]));
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]));
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;
1008 //if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.006) continue;
1010 //if(dcaV0DaughtersXi > .6) continue;
1011 //if(V0CosOfPointingAngleXi > 0.9999) continue;
1012 //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1013 //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1015 //if(V0RadiusXi < 6.0 || V0RadiusXi > 100) continue;
1019 // change mass hypothesis to cover all the possibilities
1020 if(bTrkXi->Charge()<0){
1022 xi->ChangeMassHypothesis(V0quality, 3312); //Xi- hyp.
1023 invMassXiMinus = xi->GetEffMassXi();
1026 xi->ChangeMassHypothesis(V0quality, 3334); //Omega- hyp.
1027 invMassOmegaMinus = xi->GetEffMassXi();
1030 xi->ChangeMassHypothesis(V0quality, 3312); //back to default hyp.
1033 if(bTrkXi->Charge() > 0){
1035 xi->ChangeMassHypothesis(V0quality, -3312); //anti-Xi- hyp.
1036 invMassXiPlus = xi->GetEffMassXi();
1039 xi->ChangeMassHypothesis(V0quality, -3334); //anti-Omega- hyp.
1040 invMassOmegaPlus = xi->GetEffMassXi();
1043 xi->ChangeMassHypothesis(V0quality, -3312); //back to default hyp.
1046 //PID on the daughter tracks
1048 //Resonable guess the priors for the cascade track sample
1050 Double_t priorsGuessXi[5] = {0, 0, 2, 0, 1};
1051 Double_t priorsGuessOmega[5] = {0, 0, 1, 1, 1};
1053 //Combined bachelor-daughter PID
1055 pidXi.SetPriors(priorsGuessXi);
1057 pidOmega.SetPriors(priorsGuessOmega);
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);
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;
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;
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;
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);
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;
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;
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;
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;
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);
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;
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;
1165 //B - TPC PID: 3-sigma bands on Bethe-Bloch curve
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;
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;
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;
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);
1189 xi->GetBPxPyPz(bachPx, bachPy, bachPz);
1190 bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1191 bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1193 V0toXiCosOfPointingAngle
1194 = xi->GetV0CosineOfPointingAngle(posXi[0], posXi[1], posXi[2]);
1195 rapXi = xi->RapXi();
1196 rapOmega = xi->RapOmega();
1198 alphaXi = xi->AlphaXi();
1199 ptArmXi = xi->PtArmXi();
1205 distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1207 = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1208 *(posXi[0] - bestPrimaryVtxPos[0])
1209 +(posXi[1] - bestPrimaryVtxPos[1])
1210 *(posXi[1] - bestPrimaryVtxPos[1]));
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]);
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]);
1228 fh2PhiPosBeforeVsAfter->Fill(phi, phiAfter);
1229 else if(xi->Charge()<0)
1230 fh2PhiNegBeforeVsAfter->Fill(phi, phiAfter);
1233 if( isBachelorPion &&
1234 ( (xi->Charge()<0 && isPosInXiProton && isNegInXiPion)
1235 || (xi->Charge()>0 && isNegInXiProton && isPosInXiPion ) ) )
1237 //for default hypothesis
1238 fh1Chi2Xi->Fill(chi2Xi);
1239 fh1DCAXiDaughters->Fill(dcaXiDaughters);
1240 fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1241 fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1242 fh1XiRadius->Fill(XiRadius);
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);
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);
1262 fh1XiBachPt->Fill(bachPt);
1263 fh1XiBachP->Fill(bachPtot);
1264 fh1PhiXi->Fill( xi->Phi() );
1266 fh2Armenteros->Fill(alphaXi, ptArmXi);
1270 fh1MassXiMinus->Fill(invMassXiMinus);
1271 fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1272 fh1MassXi->Fill(invMassXiMinus);
1273 fh1MassOmega->Fill(invMassOmegaMinus);
1275 fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1277 fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus,
1279 fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1280 fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1283 if(xi->Charge() > 0){
1284 fh1MassXiPlus->Fill(invMassXiPlus);
1285 fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1286 fh1MassXi->Fill(invMassXiPlus);
1287 fh1MassOmega->Fill(invMassOmegaPlus);
1289 fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter,
1291 fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1293 fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1294 fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1299 = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. ) ?
1300 0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1301 Double_t phiV0 = phiNew;
1303 // if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1304 //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1306 Double_t phiV0A = phiNew;
1308 //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1309 //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1311 Double_t phiV0C = phiNew;
1313 //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1314 //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1316 //PID cuts with TPC cuts
1317 if(xi->Charge() < 0){
1318 if(isPosProtonForTPC
1319 && isNegPionForTPC){
1320 if(isBachelorPionForTPC){
1322 fhXiRapidity->Fill(rapXi);
1323 if(TMath::Abs(rapXi) < 0.8){
1324 fh2MassVsPtXiMinus->Fill(XiPt, invMassXiMinus);
1325 fh2MassVsPtXiAll->Fill(XiPt, invMassXiMinus);
1327 for(int r=0; r!=3; ++r) {
1328 if(invMassXiMinus > fXiBands[r][0]
1329 && invMassXiMinus < fXiBands[r][1]){
1331 fProfXiV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1332 fProfXiSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1334 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1335 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1337 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1338 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1345 if(isBachelorKaonForTPC){
1347 fhOmegaRapidity->Fill(rapOmega);
1348 if(TMath::Abs(rapOmega) < 0.8){
1349 fh2MassVsPtOmegaMinus->Fill(XiPt, invMassOmegaMinus);
1350 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaMinus);
1352 for(int r=0; r!=3; ++r) {
1353 if(invMassOmegaMinus > fOmegaBands[r][0]
1354 && invMassOmegaMinus < fOmegaBands[r][1]){
1356 fProfOmegaV2PtV0A[r]->Fill(XiPt, TMath::Cos(2.*phiV0A));
1357 fProfOmegaSinePtV0A[r]->Fill(XiPt, TMath::Sin(2.*phiV0A));
1359 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1360 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1362 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1363 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1373 if(xi->Charge() > 0){
1374 if(isNegProtonForTPC
1375 && isPosPionForTPC){
1376 if(isBachelorPionForTPC){
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));
1388 fProfXiV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1389 fProfXiSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1391 fProfXiV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1392 fProfXiSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1399 if(isBachelorKaonForTPC){
1401 fhOmegaRapidity->Fill(rapOmega);
1402 if(TMath::Abs(rapOmega) < 0.8){
1403 fh2MassVsPtOmegaPlus->Fill(XiPt, invMassOmegaPlus);
1404 fh2MassVsPtOmegaAll->Fill(XiPt, invMassOmegaPlus);
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));
1412 fProfOmegaV2PtV0C[r]->Fill(XiPt, TMath::Cos(2.*phiV0C));
1413 fProfOmegaSinePtV0C[r]->Fill(XiPt, TMath::Sin(2.*phiV0C));
1415 fProfOmegaV2Pt[r]->Fill(XiPt, TMath::Cos(2.*phiV0));
1416 fProfOmegaSinePt[r]->Fill(XiPt, TMath::Sin(2.*phiV0));
1426 }//end of cascade loop
1430 PostData(1, fHistList);
1434 void AliAnalysisTaskFlowEPCascade::ReadFromAODv0(AliAODEvent *fAOD){
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);
1442 Int_t run = fAOD->GetRunNumber();
1444 // Load the calibrations run dependent
1445 OpenInfoCalbration(run);
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;
1455 AliAODVZERO* aodV0 = fAOD->GetVZEROData();
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);
1464 Qxa2 += TMath::Cos(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1465 Qya2 += TMath::Sin(2*phiV0)*multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1469 // Qx2 = Qxa2 + Qxc2;
1470 //Qy2 = Qya2 + Qya2;
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];
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];
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;
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);
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.;
1498 fhEPangleVZero->Fill(psiVZero);
1499 fhEPangleV0A->Fill(psiV0A);
1500 fhEPangleV0C->Fill(psiV0C);
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)));
1507 Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
1509 Double_t b = fAOD->GetMagneticField();
1511 int nCascades=fAOD->GetNumberOfCascades();
1512 //Info("ReadFromAODv0", "# cascades = %d", nCascades);
1514 const AliAODVertex *primaryBestAODVtx = fAOD->GetPrimaryVertex();
1515 primaryBestAODVtx->GetXYZ(bestPrimaryVtxPos);
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.;
1525 Double_t invMassLambdaAsCascDghter = 0.;
1526 Double_t V0Chi2Xi = -1.;
1527 Double_t dcaV0DaughtersXi = -1.;
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.;
1537 Double_t invMassXiMinus = 0.;
1538 Double_t invMassXiPlus = 0.;
1539 Double_t invMassOmegaMinus = 0.;
1540 Double_t invMassOmegaPlus = 0.;
1543 Bool_t isPosInXiProton = kFALSE;
1544 Bool_t isPosInXiPion = kFALSE;
1545 Bool_t isPosInOmegaProton = kFALSE;
1546 Bool_t isPosInOmegaPion = kFALSE;
1548 Bool_t isNegInXiProton = kFALSE;
1549 Bool_t isNegInXiPion = kFALSE;
1550 Bool_t isNegInOmegaProton = kFALSE;
1551 Bool_t isNegInOmegaPion = kFALSE;
1553 Bool_t isBachelorKaon = kFALSE;
1554 Bool_t isBachelorPion = kFALSE;
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;
1564 Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
1566 Double_t XiPtot = 0.;
1568 Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
1569 Double_t bachPt = 0.;
1570 Double_t bachPtot = 0.;
1572 //Short_t chargeXi = -2;
1573 Double_t V0toXiCosOfPointingAngle = 0.;
1575 Double_t rapXi = -20.;
1576 Double_t rapOmega = -20.;
1578 Double_t alphaXi = -200.;
1579 Double_t ptArmXi = -200.;
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.;
1588 const AliAODcascade *xi = fAOD->GetCascade(iXi);
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]
1602 +posXi[2]*posXi[2]);
1604 AliAODTrack *pTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1605 AliAODTrack *nTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1607 = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1609 if(!pTrkXi || !nTrkXi || !bTrkXi) continue;
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() );
1615 if(idxBach == idxNegXi || idxBach == idxPosXi) continue;
1617 if( !fCutsDau->IsSelected(pTrkXi)
1618 || !fCutsDau->IsSelected(nTrkXi)
1619 || !fCutsDau->IsSelected(bTrkXi) ) continue;
1621 if(pTrkXi->IsOn(AliESDtrack::kTPCin)){
1622 fh2TPCdEdxOfCascDghters->Fill(pTrkXi->P()*pTrkXi->Charge(),
1623 pTrkXi->GetTPCsignal());
1625 if( nTrkXi->IsOn(AliESDtrack::kTPCin) ){
1626 fh2TPCdEdxOfCascDghters->Fill(nTrkXi->P()*nTrkXi->Charge(),
1627 nTrkXi->GetTPCsignal());
1629 if(bTrkXi->IsOn(AliESDtrack::kTPCin)){
1630 fh2TPCdEdxOfCascDghters->Fill(bTrkXi->P()*bTrkXi->Charge(),
1631 bTrkXi->GetTPCsignal());
1635 if(xi->ChargeXi() < 0)
1636 invMassLambdaAsCascDghter = xi->MassLambda();
1638 invMassLambdaAsCascDghter = xi->MassAntiLambda();
1640 dcaV0DaughtersXi = xi->DcaV0Daughters();
1641 V0Chi2Xi = xi->Chi2V0();
1643 V0CosOfPointingAngleXi
1644 = xi->CosPointingAngle(bestPrimaryVtxPos);
1645 dcaV0ToPrimaryVtxXi = xi->DcaV0ToPrimVertex();
1646 dcaBachToPrimaryVtxXi = xi->DcaBachToPrimVertex();
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();
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;
1667 if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.01) continue;
1669 // if(dcaV0DaughtersXi > 1.) continue;
1670 //if(V0CosOfPointingAngleXi > 0.9999) continue;
1671 //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1672 //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1674 //if(V0RadiusXi < 1.0 || V0RadiusXi > 100) continue;
1678 if(xi->ChargeXi()<0){
1679 invMassXiMinus = xi->MassXi();
1680 invMassOmegaMinus = xi->MassOmega();
1682 invMassXiPlus = xi->MassXi();
1683 invMassOmegaPlus = xi->MassOmega();
1687 if(pTrkXi->GetMostProbablePID() == AliAODTrack::kProton) {
1688 isPosInXiProton = kTRUE;
1689 isPosInOmegaProton = kTRUE;
1691 if(pTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1692 isPosInXiPion = kTRUE;
1693 isPosInOmegaPion = kTRUE;
1696 if(nTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1697 isNegInXiPion = kTRUE;
1698 isNegInOmegaPion = kTRUE;
1700 if(nTrkXi->GetMostProbablePID() == AliAODTrack::kProton){
1701 isNegInXiProton = kTRUE;
1702 isNegInOmegaProton = kTRUE;
1705 if(bTrkXi->GetMostProbablePID() == AliAODTrack::kPion)
1706 isBachelorPion = kTRUE;
1707 if(bTrkXi->GetMostProbablePID() == AliAODTrack::kKaon)
1708 isBachelorKaon = kTRUE;
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;
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;
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;
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);
1736 bachPx = xi->MomBachX();
1737 bachPy = xi->MomBachY();
1738 bachPz = xi->MomBachZ();
1740 bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1741 bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1743 V0toXiCosOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1745 rapXi = xi->RapXi();
1746 rapOmega = xi->RapOmega();
1748 alphaXi = xi->AlphaXi();
1749 ptArmXi = xi->PtArmXi();
1751 distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1753 = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1754 *(posXi[0] - bestPrimaryVtxPos[0])
1755 +(posXi[1] - bestPrimaryVtxPos[1])
1756 *(posXi[1] - bestPrimaryVtxPos[1]));
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]);
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);
1782 if( (xi->ChargeXi() < 0 && isBachelorPionForTPC
1783 && isPosProtonForTPC && isNegPionForTPC)
1784 || (xi->ChargeXi() > 0 && isBachelorPionForTPC
1785 && isNegProtonForTPC && isPosPionForTPC))
1787 //for default hypothesis
1788 fh1Chi2Xi->Fill(chi2Xi);
1789 fh1DCAXiDaughters->Fill(dcaXiDaughters);
1790 fh1DCABachToPrimVertex->Fill(dcaBachToPrimaryVtxXi);
1791 fh1XiCosOfPointingAngle->Fill(XiCosOfPointingAngle);
1792 fh1XiRadius->Fill(XiRadius);
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);
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);
1812 fh1XiBachPt->Fill(bachPt);
1813 fh1XiBachP->Fill(bachPtot);
1814 fh1PhiXi->Fill( xi->Phi() );
1816 fh2Armenteros->Fill(alphaXi, ptArmXi);
1819 if(xi->ChargeXi()<0){
1820 fh1MassXiMinus->Fill(invMassXiMinus);
1821 fh1MassOmegaMinus->Fill(invMassOmegaMinus);
1822 fh1MassXi->Fill(invMassXiMinus);
1823 fh1MassOmega->Fill(invMassOmegaMinus);
1825 fh2MassLambdaVsMassXiMinus->Fill(invMassLambdaAsCascDghter,
1827 fh2MassXiVsMassOmegaMinus->Fill(invMassXiMinus,
1829 fh2XiRadiusVsMassXiMinus->Fill(XiRadius, invMassXiMinus);
1830 fh2XiRadiusVsMassOmegaMinus->Fill(XiRadius, invMassOmegaMinus);
1833 if(xi->ChargeXi() > 0){
1834 fh1MassXiPlus->Fill(invMassXiPlus);
1835 fh1MassOmegaPlus->Fill(invMassOmegaPlus);
1836 fh1MassXi->Fill(invMassXiPlus);
1837 fh1MassOmega->Fill(invMassOmegaPlus);
1839 fh2MassLambdaVsMassXiPlus->Fill(invMassLambdaAsCascDghter,
1841 fh2MassXiVsMassOmegaPlus->Fill(invMassXiPlus,
1843 fh2XiRadiusVsMassXiPlus->Fill(XiRadius, invMassXiPlus);
1844 fh2XiRadiusVsMassOmegaPlus->Fill(XiRadius, invMassOmegaPlus);
1849 // = ( XiPAfter[0] == 0. && XiPAfter[1] == 0. )
1850 // ? 0.0 : TMath::ATan2(XiPAfter[1], XiPAfter[0]);
1851 Double_t phiV0 = phiAfter;
1853 //if(phiV0 < 0) phiV0 += 2.*TMath::Pi();
1854 //if(phiV0 > TMath::Pi()) phiV0 -= TMath::Pi();
1857 Double_t phiV0A = phiAfter;
1859 //if(phiV0A < 0) phiV0A += 2.*TMath::Pi();
1860 //if(phiV0A > TMath::Pi()) phiV0A -= TMath::Pi();
1862 Double_t phiV0C = phiAfter;
1864 //if(phiV0C < 0) phiV0C += 2.*TMath::Pi();
1865 //if(phiV0C > TMath::Pi()) phiV0C -= TMath::Pi();
1867 //PID cuts with TPC cuts
1868 if(xi->ChargeXi() < 0){
1869 if(isPosProtonForTPC
1870 && isNegPionForTPC){
1871 if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
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]){
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));
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));
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));
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));
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));
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));
1922 if(xi->ChargeXi() > 0){
1923 if(isNegProtonForTPC
1924 && isPosPionForTPC){
1925 if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
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));
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));
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));
1948 if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
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));
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));
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));
1977 }//end of cascade loop
1979 PostData(1, fHistList);
1984 AliAnalysisTaskFlowEPCascade::~AliAnalysisTaskFlowEPCascade(){
1985 if(fHistList) delete fHistList;
1988 void AliAnalysisTaskFlowEPCascade::Terminate(Option_t *)
1991 fHistList = dynamic_cast<TList*> (GetOutputData(1));
1994 printf("ERROR: Output tree not available\n");
2001 void AliAnalysisTaskFlowEPCascade::Propagate(Double_t vv[3],
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
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;
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;
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);
2027 AliAnalysisTaskFlowEPCascade::OpenInfoCalbration(Int_t run){
2029 AliOADBContainer *cont = (AliOADBContainer*) fOADB->Get("hMultV0BefCorr");
2031 printf("OADB object hMultV0BefCorr is not available in the file\n");
2035 if(!(cont->GetObject(run))){
2036 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
2039 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
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);
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++){
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);
2060 cont = (AliOADBContainer*) fOADB->Get(namecont);
2062 printf("OADB object %s is not available in the file\n",namecont);
2066 if(!(cont->GetObject(run))){
2067 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
2070 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2071 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();