PWG/FLOW/Tasks/AliAnalysisTaskFlowEPCascade.cxx

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