PWG/FLOW/Tasks/AliAnalysisTaskFlowCascade.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-2008,ALICE Experiment at CERN,All rights reserved.    *
   3 *                                                                         *
   4 * Author: The ALICE Off-line Project.                                     *
   5 * Contributors are mentioned in the code where appropriate.               *
   6 *                                                                         *
   7 * Permission to use,copy,modify and distribute this software and its      *
   8 * documentation strictly for non-commercial purposes is hereby granted    *
   9 * without fee,provided that the above copyright notice appears in all     *
  10 * copies and that both the copyright notice and this permission notice    *
  11 * appear in the supporting documentation. The authors make no claims      *
  12 * about the suitability of this software for any purpose. It is           *
  13 * provided "as is" without express or implied warranty.                   *
  14 **************************************************************************/
  15
  16 /////////////////////////////////////////////////////
  17 // AliAnalysisTaskFlowCascade:
  18 // Analysis task to select Xi and Omega candidates for flow analysis.
  19 //
  20 // Author: Zhong-Bao.Yin@cern.ch
  21 //////////////////////////////////////////////////////
  22
  23 #include "TChain.h"
  24 #include "TList.h"
  25 #include "TH1D.h"
  26 #include "TH2D.h"
  27 #include "TH3D.h"
  28 #include "TProfile.h"
  29 #include "TVector3.h"
  30
  31 #include "AliAnalysisTaskSE.h"
  32 #include "AliAnalysisManager.h"
  33
  34 #include "AliVParticle.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliESDv0.h"
  38 #include "AliESDcascade.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDtrackCuts.h"
  41 #include "AliCentrality.h"
  42 #include "AliVVertex.h"
  43 #include "AliESDVZERO.h"
  44 #include "AliESDUtils.h"
  45
  46 #include "AliTPCPIDResponse.h"
  47 #include "AliTOFPIDResponse.h"
  48 #include "AliPIDResponse.h"
  49
  50 #include "AliAODEvent.h"
  51 #include "AliAODInputHandler.h"
  52 #include "AliAODTrack.h"
  53 #include "AliAODVZERO.h"
  54 #include "AliAODcascade.h"
  55
  56 #include "TMath.h"
  57 #include "TObjArray.h"
  58 #include "AliFlowCandidateTrack.h"
  59
  60 #include "AliFlowTrackCuts.h"
  61 #include "AliFlowEventCuts.h"
  62 #include "AliFlowEvent.h"
  63 #include "AliFlowCommonConstants.h"
  64
  65 #include "AliAnalysisTaskFlowCascade.h"
  66
  67 ClassImp(AliAnalysisTaskFlowCascade)
  68
  69 //_____________________________________________________________________________
  70   AliAnalysisTaskFlowCascade::AliAnalysisTaskFlowCascade() :
  71     AliAnalysisTaskSE(),
  72     //    fMinCent(0), fMaxCent(0),
  73     fSpecie(0),
  74     fMassBins(0),
  75     fMinMass(0.0),
  76     fMaxMass(0.0),
  77     fCutsEvent(NULL),
  78     fCutsRPTPC(NULL),
  79     fCutsRPVZE(NULL),
  80     fCutsPOI(NULL),
  81     fCutsDau(NULL),
  82     fPIDResponse(NULL),
  83     fFlowEventTPC(NULL),
  84     fFlowEventVZE(NULL),
  85     fCandidates(NULL),
  86     fQAList(NULL)
  87 {
  88   //ctor
  89   for (Int_t i=0; i!=8; ++i)
  90     fCascadeCuts[i] = 0;
  91
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliAnalysisTaskFlowCascade
  96 ::AliAnalysisTaskFlowCascade(const char *name,
  97                              AliFlowEventCuts *cutsEvent,
  98                              AliFlowTrackCuts *cutsRPTPC,
  99                              AliFlowTrackCuts *cutsRPVZE,
 100                              /* AliESDtrackCuts */ AliFlowTrackCuts *cutsDau ) :
 101   AliAnalysisTaskSE(name),
 102   //fMinCent(minCent), fMaxCent(maxCent),
 103   fSpecie(0),
 104   fMassBins(0),
 105   fMinMass(0.0),
 106   fMaxMass(0.0),
 107   fCutsEvent(cutsEvent),
 108   fCutsRPTPC(cutsRPTPC),
 109   fCutsRPVZE(cutsRPVZE),
 110   fCutsPOI(NULL),
 111   fCutsDau(cutsDau),
 112   fPIDResponse(NULL),
 113   fFlowEventTPC(NULL),
 114   fFlowEventVZE(NULL),
 115   fCandidates(NULL),
 116   fQAList(NULL)
 117 {
 118   //ctor
 119   for (Int_t i=0; i!=8; ++i)
 120     fCascadeCuts[i] = 0;
 121
 122   DefineInput( 0,TChain::Class());
 123   DefineOutput(1,AliFlowEventSimple::Class()); // TPC object
 124   DefineOutput(2,AliFlowEventSimple::Class()); // VZE object
 125   DefineOutput(3,TList::Class());
 126 }
 127
 128 //_____________________________________________________________________________
 129 AliAnalysisTaskFlowCascade::~AliAnalysisTaskFlowCascade()
 130 {
 131   if(fQAList) delete fQAList;
 132   if (fFlowEventTPC) delete fFlowEventTPC;
 133   if (fFlowEventVZE) delete fFlowEventVZE;
 134   if (fCandidates)   delete fCandidates;
 135   if (fCutsDau)      delete fCutsDau;
 136   if (fCutsPOI)      delete fCutsPOI;
 137   if (fCutsRPTPC)   delete fCutsRPTPC;
 138   if (fCutsRPVZE)   delete fCutsRPVZE;
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliAnalysisTaskFlowCascade::UserCreateOutputObjects()
 143 {
 144
 145   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 146   AliInputEventHandler* inputHandler
 147     = (AliInputEventHandler*) (man->GetInputEventHandler());
 148   fPIDResponse = inputHandler->GetPIDResponse();
 149
 150   fQAList = new TList();
 151   fQAList->SetOwner();
 152   AddQAEvents();
 153   AddQACandidates();
 154   //  PostData(3,fQAList);
 155
 156   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 157   cc->SetNbinsMult(1);
 158   cc->SetMultMin(0);
 159   cc->SetMultMax(1);
 160
 161   cc->SetNbinsPt(20);
 162   cc->SetPtMin(0.0);
 163   cc->SetPtMax(10.0);
 164
 165   cc->SetNbinsPhi(1);
 166   cc->SetPhiMin(0.0);
 167   cc->SetPhiMax(TMath::TwoPi());
 168
 169   cc->SetNbinsEta(1);
 170   cc->SetEtaMin(-2.0);
 171   cc->SetEtaMax(+2.0);
 172
 173   cc->SetNbinsQ(3);
 174   cc->SetQMin(0.0);
 175   cc->SetQMax(3.0);
 176
 177   cc->SetNbinsMass(fMassBins);
 178   cc->SetMassMin(fMinMass);
 179   cc->SetMassMax(fMaxMass);
 180
 181   fCutsPOI = new AliFlowTrackCuts("null_cuts");
 182   fCutsPOI->SetParamType(fCutsRPTPC->GetParamType());
 183   fCutsPOI->SetPtRange(+1,-1); // select nothing QUICK
 184   fCutsPOI->SetEtaRange(+1,-1); // select nothing VZERO
 185
 186   fFlowEventTPC = new AliFlowEvent(3000);
 187   fFlowEventVZE = new AliFlowEvent(1000);
 188   fCandidates = new TObjArray(100);
 189   fCandidates->SetOwner();
 190
 191   PostData(1,fFlowEventTPC);
 192   PostData(2,fFlowEventVZE);
 193   PostData(3,fQAList);
 194 }
 195
 196 //_____________________________________________________________________________
 197 void AliAnalysisTaskFlowCascade::AddQAEvents()
 198 {
 199   TList *tQAEvents = new TList();
 200   tQAEvents->SetName("Events");
 201   tQAEvents->SetOwner();
 202   TH1I* tEvent = new TH1I("Event","Number of Events",   3,0,3);
 203   tQAEvents->Add(tEvent);
 204
 205   TH1D *tTPCRFP = new TH1D("RFPTPC",
 206                            "TPC Reference Flow Particles;multiplicity",
 207                            100, 0, 3000);
 208   tQAEvents->Add(tTPCRFP);
 209   TH1D *tVZERFP = new TH1D("RFPVZE",
 210                            "VZERO Reference Flow Particles;multiplicity",
 211                            100, 0, 30000);
 212   tQAEvents->Add(tVZERFP);
 213
 214   TProfile *tCuts = new TProfile("Cuts","Analysis Cuts",10,0,10);
 215   tCuts->Fill(0.5,fCascadeCuts[0],1);
 216   tCuts->GetXaxis()->SetBinLabel(1,"dcaXiDau");
 217   tCuts->Fill(1.5,fCascadeCuts[1],1);
 218   tCuts->GetXaxis()->SetBinLabel(2,"XiCPA");
 219   tCuts->Fill(2.5,fCascadeCuts[2],1);
 220   tCuts->GetXaxis()->SetBinLabel(3,"dcaV0Vtx");
 221   tCuts->Fill(3.5,fCascadeCuts[3],1);
 222   tCuts->GetXaxis()->SetBinLabel(4,"dcaBachVtx");
 223   tCuts->Fill(4.5,fCascadeCuts[4],1);
 224   tCuts->GetXaxis()->SetBinLabel(5,"dcaV0Dau");
 225   tCuts->Fill(5.5,fCascadeCuts[5],1);
 226   tCuts->GetXaxis()->SetBinLabel(6,"V0CPA");
 227   tCuts->Fill(6.5,fCascadeCuts[6],1);
 228   tCuts->GetXaxis()->SetBinLabel(7,"dcaV0DauVtx");
 229   tCuts->Fill(7.5,fCascadeCuts[7],1);
 230   tCuts->GetXaxis()->SetBinLabel(8,"V0Mass");
 231   tQAEvents->Add(tCuts);
 232
 233   fQAList->Add(tQAEvents);
 234 }
 235
 236 //_____________________________________________________________________________
 237 void AliAnalysisTaskFlowCascade::AddQACandidates()
 238 {
 239   TList *tQACandidates;
 240
 241   tQACandidates = new TList();
 242   tQACandidates->SetOwner();
 243   tQACandidates->SetName("Candidates");
 244
 245   TH1F* tChi2Xi = new TH1F("Chi2Xi",
 246                        "Cascade #chi^{2}; #chi^{2}; Number of Cascades",
 247                        160, 0, 160);
 248   tQACandidates->Add(tChi2Xi);
 249
 250   TH1F* tDCAXiDaughters
 251     = new TH1F( "DcaXiDaughters",
 252                 "DCA between Xi Daughters; DCA (cm); Number of Cascades",
 253                 100, 0., 0.5);
 254   tQACandidates->Add(tDCAXiDaughters);
 255
 256   TH1F * tDCABachToPrimVertex
 257     = new TH1F("DcaBachToPrimVertex",
 258                "DCA of Bach. to Prim. Vertex; DCA (cm);Number of Cascades",
 259                250, 0., 2.5);
 260   tQACandidates->Add(tDCABachToPrimVertex);
 261
 262   TH1F * tXiCosOfPointingAngle
 263     = new TH1F("XiCosOfPointingAngle",
 264                "Cos of Xi Pointing Angle; Cos (Xi Point.Angl);Number of Xis",
 265                200, 0.99, 1.0);
 266   tQACandidates->Add(tXiCosOfPointingAngle);
 267
 268   TH1F * tXiRadius = new TH1F("XiRadius",
 269                               "Casc. decay transv. radius; r (cm); Counts" ,
 270                               1050, 0., 105.0 );
 271   tQACandidates->Add(tXiRadius);
 272
 273   TH1F *tMassLambda
 274     = new TH1F("MassLambdaAsCascDghter",
 275                "#Lambda assoc. to Casc. candidates; Eff. Mass (GeV/c^{2}); Counts",
 276                300,1.00,1.3);
 277   tQACandidates->Add(tMassLambda);
 278
 279   TH1F *tV0Chi2 = new TH1F("V0Chi2Xi",
 280                        "V0 #chi^{2}, in cascade; #chi^{2};Counts",
 281                        160, 0, 40);
 282   tQACandidates->Add(tV0Chi2);
 283
 284   TH1F * tV0CosOfPointingAngle
 285     = new TH1F("V0CosOfPointingAngleXi",
 286                "Cos of V0 Pointing Angle, in cascade;Cos(V0 Point. Angl); Counts",
 287                200, 0.98, 1.0);
 288   tQACandidates->Add(tV0CosOfPointingAngle);
 289
 290   TH1F *tV0Radius  = new TH1F("V0RadiusXi",
 291                           "V0 decay radius, in cascade; radius (cm); Counts",
 292                           1050, 0., 105.0);
 293   tQACandidates->Add(tV0Radius);
 294
 295   TH1F * tDcaV0DaughtersXi
 296     = new TH1F("DcaV0DaughtersXi",
 297                "DCA between V0 daughters, in cascade;DCA (cm);Number of V0s",
 298                120, 0., 0.6);
 299   tQACandidates->Add(tDcaV0DaughtersXi);
 300
 301   TH1F * tDcaV0ToPrimVertex
 302     = new TH1F("DcaV0ToPrimVertexXi",
 303                "DCA of V0 to Prim. Vertex, in cascade;DCA (cm);Number of Cascades", 200, 0., 1.);
 304   tQACandidates->Add(tDcaV0ToPrimVertex);
 305
 306   TH1F * tDCAPosToPrimVertex =
 307     new TH1F("DcaPosToPrimVertexXi",
 308              "DCA of V0 pos daughter to Prim. Vertex;DCA (cm);Counts",
 309              300, 0, 3);
 310   tQACandidates->Add(tDCAPosToPrimVertex);
 311
 312   TH1F * tDCANegToPrimVertex
 313     =  new TH1F("DcaNegToPrimVertexXi",
 314                 "DCA of V0 neg daughter to Prim. Vertex;DCA (cm);Counts",
 315                 300, 0, 3);
 316   tQACandidates->Add(tDCANegToPrimVertex);
 317
 318   TH1F *tV0toXiCosOfPointingAngle
 319     = new TH1F("V0toXiCosOfPointingAngle",
 320                "Cos. of V0 Ptng Angl Xi vtx; Cos(V0 Point. Angl / Xi vtx); Counts",
 321                100, 0.99, 1.0);
 322   tQACandidates->Add(tV0toXiCosOfPointingAngle);
 323
 324   TH2F *th2Armenteros
 325     = new TH2F("Armenteros",
 326                "#alpha_{Arm}(casc. cand.) Vs Pt_{Arm}(casc. cand.); #alpha_{Arm} ; Pt_{Arm} (GeV/c)",
 327                140, -1.2, 1.2, 300, 0., 0.3);
 328   tQACandidates->Add(th2Armenteros);
 329
 330   TH2F *th2TPCdEdxOfCascDghters
 331     = new TH2F( "TPCdEdxOfCascDghters",
 332                 "TPC dE/dx of the cascade daughters; charge x || #vec{p}_{TPC inner wall}(Casc. daughter) || (GeV/c); TPC signal (ADC) ",
 333                 200, -10.0, 10.0, 450, 0., 900.);
 334   tQACandidates->Add(th2TPCdEdxOfCascDghters);
 335
 336   TH2F *th2MassVsPtAll
 337     = new TH2F("MassVsPtAll",
 338                "M_{candidates} vs Pt; Pt (GeV/c); M (GeV/c^{2})",
 339                100, 0., 10., fMassBins, fMinMass, fMaxMass);
 340   tQACandidates->Add(th2MassVsPtAll);
 341
 342   TH1F *tDistToVtxZAfter
 343     = new TH1F("DistToVtxZAfter",
 344                "Distance to vtx z after propagation to vtx; z [cm]",
 345                100, -5., 5.);
 346   tQACandidates->Add(tDistToVtxZAfter);
 347
 348   TH1F * tDistToVtxXYAfter
 349     = new TH1F("DistToVtxXYAfter",
 350                "Distance to vtx xy after propagation to vtx",
 351                500, 0., 50.);
 352   tQACandidates->Add(tDistToVtxXYAfter);
 353
 354   TH2F *th2DistToVtxZBeforeVsAfter
 355     = new TH2F("DistToVtxZBeforeVsAfter",
 356                "Distance to vtx z before vs after propagation; Distance before [cm]; Distance after [cm]",
 357                500, -50., 50., 100, -5., 5.);
 358   tQACandidates->Add(th2DistToVtxZBeforeVsAfter);
 359
 360   TH2F *th2DistToVtxXYBeforeVsAfter
 361     = new TH2F("DistToVtxXYBeforeVsAfter",
 362                "Distance to vtx xy before vs after propagation; Distance before [cm]; Distance after [cm]",
 363                500, 0., 50, 500, 0., 50);
 364   tQACandidates->Add(th2DistToVtxXYBeforeVsAfter);
 365
 366   TH2F * th2PxBeforeVsAfter
 367     = new TH2F("PxBeforeVsAfter",
 368                "Px before vs after propagation; Px [GeV/c]; Px' [GeV/c]",
 369                200, -10., 10, 200, -10., 10.);
 370   tQACandidates->Add(th2PxBeforeVsAfter);
 371
 372   TH2F * th2PyBeforeVsAfter
 373     = new TH2F("PyBeforeVsAfter",
 374                "Py before vs after propagation; Py [GeV/c]; Py' [GeV/c]",
 375                200, -10., 10, 200, -10., 10.);
 376   tQACandidates->Add(th2PyBeforeVsAfter);
 377
 378   TH2F * th2PhiPosBeforeVsAfter
 379     = new TH2F("PhiPosBeforeVsAfter",
 380                "Phi for positively charged candidates before vs after propagation; #phi; #phi'",
 381                360, 0., 2.0*TMath::Pi(), 360, 0., 2.0*TMath::Pi());
 382   tQACandidates->Add(th2PhiPosBeforeVsAfter);
 383
 384   TH2F *th2PhiNegBeforeVsAfter
 385     = new TH2F("PhiNegBeforeVsAfter",
 386                "Phi for negatively charged candidates before vs after propagation; #phi; #phi'",
 387                360, 0., 2.0*TMath::Pi(), 360, 0., 2.0*TMath::Pi());
 388   tQACandidates->Add(th2PhiNegBeforeVsAfter);
 389
 390   fQAList->Add(tQACandidates);
 391 }
 392
 393 //_____________________________________________________________________________
 394 void AliAnalysisTaskFlowCascade::NotifyRun()
 395 {
 396 }
 397 //_____________________________________________________________________________
 398 void AliAnalysisTaskFlowCascade::UserExec(Option_t *)
 399 {
 400   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 401   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 402   Bool_t acceptEvent=kFALSE;
 403   fCandidates->SetLast(-1);
 404
 405   if(fESD) {
 406     // recorrecting VZERO (for pass 1 only)
 407     /*
 408     Float_t *vChCorr = new Float_t[64];
 409     Float_t dummy;
 410     AliESDUtils::GetCorrV0(fESD,dummy,NULL,vChCorr);
 411     AliESDVZERO *vzero = (AliESDVZERO*) fESD->GetVZEROData();
 412     vzero->SetMultiplicity( vChCorr );
 413     delete [] vChCorr;
 414     */
 415     //
 416
 417     ((TH1I*)((TList*)fQAList->FindObject("Events"))->FindObject("Event"))->Fill(0);
 418
 419     const AliVVertex *vtxGlb = fESD->GetPrimaryVertexTracks();
 420     const AliVVertex *vtxSPD = fESD->GetPrimaryVertexSPD();
 421     if(!vtxGlb || !vtxSPD) return;
 422     if( fCutsEvent->IsSelected(fESD, 0) && (TMath::Abs(vtxSPD->GetZ()-vtxGlb->GetZ()) <= 0.5) ) {
 423       acceptEvent = kTRUE;
 424     ((TH1I*)((TList*)fQAList->FindObject("Events"))->FindObject("Event"))->Fill(2);
 425       ReadFromESDv0(fESD);
 426     }
 427   } else if(fAOD) {
 428     const AliVVertex *vtxGlb = fAOD->GetPrimaryVertex();
 429     const AliVVertex *vtxSPD = fAOD->GetPrimaryVertexSPD();
 430     if(!vtxGlb || !vtxSPD) return;
 431
 432     ((TH1I*)((TList*)fQAList->FindObject("Events"))->FindObject("Event"))->Fill(0);
 433
 434     if(fCutsEvent->IsSelected(fAOD, 0) && (TMath::Abs(vtxSPD->GetZ()-vtxGlb->GetZ()) <= 0.5) ) {
 435       acceptEvent = kTRUE;
 436       ((TH1I*)((TList*)fQAList->FindObject("Events"))->FindObject("Event"))->Fill(2);
 437       ReadFromAODv0(fAOD);
 438     }
 439
 440
 441     /*
 442
 443     AliAODHeader *aodHeader = fAOD->GetHeader();
 444     if(!aodHeader) return;
 445     AliCentrality *centrality = aodHeader->GetCentralityP();
 446     if(!centrality) return;
 447     Double_t cent = centrality->GetCentralityPercentile("V0M" );
 448     Double_t cent1 = centrality->GetCentralityPercentile("TRK" );
 449     if(TMath::Abs(cent-cent1) >= 5.) return;
 450
 451     if(cent<fMinCent||cent>=fMaxCent) return; //centrality cut
 452
 453     Double_t zvtx = fAOD->GetPrimaryVertex()->GetZ();
 454     if(TMath::Abs(zvtx)>10.) return; //vertex cut
 455
 456     ((TH1I*)((TList*)fQAList->FindObject("Events"))->FindObject("Event"))->Fill(2);
 457     ReadFromAODv0(fAOD);
 458     */
 459   }
 460
 461   if(!acceptEvent) return;
 462
 463   ((TH1D*)((TList*)fQAList->FindObject("Events"))->FindObject("RFPTPC"))
 464     ->Fill(fFlowEventTPC->GetNumberOfRPs() );
 465   Double_t mult=0;
 466   for(Int_t i=0; i != fFlowEventVZE->GetNumberOfRPs(); ++i) {
 467     AliFlowTrackSimple *pTrack = fFlowEventVZE->GetTrack(i);
 468     mult += pTrack->Weight();
 469   }
 470   ((TH1D*)((TList*)fQAList->FindObject("Events"))->FindObject("RFPVZE"))
 471     ->Fill( mult );
 472
 473   //  if(fDebug) printf("TPCevent %d | VZEevent %d\n",
 474   //                  fFlowEventTPC->NumberOfTracks(),
 475   //                  fFlowEventVZE->NumberOfTracks() );
 476   AddCandidates();
 477
 478   PostData(1,fFlowEventTPC);
 479   PostData(2,fFlowEventVZE);
 480   PostData(3,fQAList);
 481
 482   return;
 483 }
 484
 485 //_____________________________________________________________________________
 486 void AliAnalysisTaskFlowCascade::AddCandidates(){
 487
 488   //  if(fDebug) printf("I received %d candidates\n",
 489   //                fCandidates->GetEntriesFast());
 490   for(int iCand=0; iCand!=fCandidates->GetEntriesFast(); ++iCand ) {
 491     AliFlowCandidateTrack *cand
 492       = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
 493     if(!cand) continue;
 494     // if(fDebug)
 495     //  printf(" >Checking at candidate %d with %d daughters: mass %f\n",
 496     //       iCand, cand->GetNDaughters(), cand->Mass());
 497
 498     // untagging ===>
 499     for(int iDau=0; iDau != cand->GetNDaughters(); ++iDau) {
 500       // if(fDebug)
 501       // printf(" >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
 502       for(int iRPs=0; iRPs != fFlowEventTPC->NumberOfTracks(); ++iRPs ) {
 503         AliFlowTrack *iRP
 504           = dynamic_cast<AliFlowTrack*>(fFlowEventTPC->GetTrack( iRPs ));
 505         if (!iRP) continue;
 506         if( !iRP->InRPSelection() ) continue;
 507         if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
 508           //if(fDebug) printf(" was in RP set");
 509           iRP->SetForRPSelection(kFALSE);
 510           fFlowEventTPC->SetNumberOfRPs( fFlowEventTPC->GetNumberOfRPs() -1 );
 511         }
 512       }
 513       //if(fDebug) printf("\n");
 514     }
 515     // <=== untagging
 516     cand->SetForPOISelection(kTRUE);
 517     fFlowEventTPC->InsertTrack( ((AliFlowTrack*) cand) );
 518     fFlowEventVZE->InsertTrack( ((AliFlowTrack*) cand) );
 519   }
 520
 521   //  if(fDebug) printf("TPCevent %d | VZEevent %d\n",
 522   //                fFlowEventTPC->NumberOfTracks(),
 523   //                fFlowEventVZE->NumberOfTracks() );
 524
 525 }
 526
 527 //______________________________________________________________________________
 528 void AliAnalysisTaskFlowCascade::ReadFromESDv0(AliESDEvent *fESD)
 529 {
 530   //AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("null_cuts");
 531   //cutsPOI->SetParamType( fCutsRP->GetParamType() );
 532   //cutsPOI->SetParamType( AliFlowTrackCuts::kGlobal );
 533   // cutsPOI->SetPtRange(+1,-1); // select nothing
 534   //cutsPOI->SetEtaRange(+1,-1); // select nothing VZERO
 535
 536   fCutsRPTPC->SetEvent(fESD,MCEvent());
 537   fCutsRPVZE->SetEvent(fESD,MCEvent());
 538
 539   fCutsPOI->SetEvent(fESD,MCEvent());
 540
 541   fFlowEventTPC->Fill(fCutsRPTPC,fCutsPOI);
 542   fFlowEventVZE->Fill(fCutsRPVZE,fCutsPOI);
 543
 544   Double_t trkPrimaryVtxPos[3] = {-100., -100., -100.};
 545   Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
 546   int nCascades=fESD->GetNumberOfCascades();
 547
 548   const AliESDVertex *primaryTrackingESDVtx = fESD->GetPrimaryVertexTracks();
 549   primaryTrackingESDVtx->GetXYZ(trkPrimaryVtxPos);
 550
 551   const AliESDVertex *primaryBestESDVtx = fESD->GetPrimaryVertex();
 552   primaryBestESDVtx->GetXYZ(bestPrimaryVtxPos);
 553
 554   Double_t b = fESD->GetMagneticField();
 555
 556   for(int i = 0; i != nCascades; ++i) {
 557
 558     // Double_t trkPrimaryVtxRadius3D = -500.;
 559     // Double_t bestPrimaryVtxRadius3D = -500.;
 560     Double_t effMassXi = 0.;
 561     Double_t chi2Xi = -1.;
 562     Double_t dcaXiDaughters = -1.;
 563     Double_t XiCosOfPointingAngle = -1.;
 564     Double_t posXi[3] = {-1000., -1000., -1000.};
 565     Double_t XiRadius = -1000.;
 566
 567     // Double_t innerWallMomCascDghters[3] = {-100., -100., -100.};
 568     //Double_t tpcSignalCascDghters[3] = {-100., -100., -100.};
 569
 570     Double_t invMassLambdaAsCascDghter = 0.;
 571     Double_t V0Chi2Xi = -1.;
 572     Double_t dcaV0DaughtersXi = -1.;
 573
 574     Double_t dcaBachToPrimaryVtxXi = -1.;
 575     Double_t dcaV0ToPrimaryVtxXi = -1.;
 576     Double_t dcaPosToPrimaryVtxXi = -1.;
 577     Double_t dcaNegToPrimaryVtxXi = -1.;
 578     Double_t V0CosOfPointingAngleXi = -1.;
 579     Double_t posV0Xi[3] = {-1000., -1000., -1000.};
 580     Double_t V0RadiusXi = -1000.;
 581     Double_t V0quality = 0.;
 582
 583     Double_t invMassXiMinus = 0.;
 584     Double_t invMassXiPlus = 0.;
 585     Double_t invMassOmegaMinus = 0.;
 586     Double_t invMassOmegaPlus = 0.;
 587
 588     /*
 589     Bool_t isPosInXiProton = kFALSE;
 590     Bool_t isPosInXiPion = kFALSE;
 591     Bool_t isPosInOmegaProton = kFALSE;
 592     Bool_t isPosInOmegaPion = kFALSE;
 593
 594     Bool_t isNegInXiProton = kFALSE;
 595     Bool_t isNegInXiPion = kFALSE;
 596     Bool_t isNegInOmegaProton = kFALSE;
 597     Bool_t isNegInOmegaPion = kFALSE;
 598
 599     Bool_t isBachelorKaon = kFALSE;
 600     Bool_t isBachelorPion = kFALSE;
 601     */
 602
 603     Bool_t isBachelorKaonForTPC = kFALSE;
 604     Bool_t isBachelorPionForTPC = kFALSE;
 605     Bool_t isNegPionForTPC = kFALSE;
 606     Bool_t isPosPionForTPC = kFALSE;
 607     Bool_t isNegProtonForTPC = kFALSE;
 608     Bool_t isPosProtonForTPC = kFALSE;
 609
 610     Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
 611     Double_t XiPt = 0.;
 612     Double_t XiPtot = 0.;
 613
 614     Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
 615     Double_t bachPt = 0.;
 616     Double_t bachPtot = 0.;
 617
 618     //Short_t chargeXi = -2;
 619     Double_t V0toXiCosOfPointingAngle = 0.;
 620
 621     Double_t rapXi = -20.;
 622     Double_t rapOmega = -20.;
 623     Double_t phi = 6.3;
 624     Double_t alphaXi = -200.;
 625     Double_t ptArmXi = -200.;
 626     //    TLorentzVector lv1, lv2, lv3, lv12, lvXi;
 627
 628     Double_t distToVtxZBefore = -999.;
 629     Double_t distToVtxZAfter = -999.;
 630     Double_t distToVtxXYBefore = -999.;
 631     Double_t distToVtxXYAfter = -999.;
 632     Double_t XiPAfter[3] = {-999., -999., -999.};
 633     Double_t phiAfter = -999.;
 634
 635     AliESDcascade *xi = fESD->GetCascade(i);
 636     if(!xi) continue;
 637
 638     if(xi->Charge()<0)
 639       xi->ChangeMassHypothesis(V0quality, 3312); // Xi- hypothesis
 640     else if(xi->Charge() > 0)
 641       xi->ChangeMassHypothesis(V0quality, -3312);
 642     else continue;
 643
 644     effMassXi = xi->GetEffMassXi();
 645     chi2Xi = xi->GetChi2Xi();
 646     dcaXiDaughters = xi->GetDcaXiDaughters();
 647     XiCosOfPointingAngle
 648       = xi->GetCascadeCosineOfPointingAngle(bestPrimaryVtxPos[0],
 649                                             bestPrimaryVtxPos[1],
 650                                             bestPrimaryVtxPos[2]);
 651     xi->GetXYZcascade(posXi[0], posXi[1], posXi[2]);
 652     XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
 653                            +posXi[1]*posXi[1]
 654                            +posXi[2]*posXi[2]);
 655
 656     UInt_t idxPosXi = (UInt_t)TMath::Abs(xi->GetPindex());
 657     UInt_t idxNegXi = (UInt_t)TMath::Abs(xi->GetNindex());
 658     UInt_t idxBach = (UInt_t)TMath::Abs(xi->GetBindex());
 659
 660     if(idxBach == idxPosXi || idxBach == idxNegXi) continue;
 661
 662     AliESDtrack *pTrkXi = fESD->GetTrack(idxPosXi);
 663     AliESDtrack *nTrkXi = fESD->GetTrack(idxNegXi);
 664     AliESDtrack *bTrkXi = fESD->GetTrack(idxBach);
 665
 666     if( !pTrkXi || !nTrkXi || !bTrkXi ) continue;
 667
 668     if( !fCutsDau->IsSelected(pTrkXi)
 669         || !fCutsDau->IsSelected(nTrkXi)
 670         || !fCutsDau->IsSelected(bTrkXi) ) continue;
 671
 672     const AliExternalTrackParam *pExtTrk = pTrkXi->GetInnerParam();
 673     const AliExternalTrackParam *nExtTrk = nTrkXi->GetInnerParam();
 674     const AliExternalTrackParam *bExtTrk = bTrkXi->GetInnerParam();
 675
 676     if(pExtTrk && pTrkXi->IsOn(AliESDtrack::kTPCin)){
 677       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(pExtTrk->GetP()*pExtTrk->Charge(), pTrkXi->GetTPCsignal());
 678     }
 679     if(nExtTrk && nTrkXi->IsOn(AliESDtrack::kTPCin)){
 680       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(nExtTrk->GetP()*nExtTrk->Charge(), nTrkXi->GetTPCsignal());
 681     }
 682     if(bExtTrk && bTrkXi->IsOn(AliESDtrack::kTPCin)){
 683       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(bExtTrk->GetP()*bExtTrk->Charge(), bTrkXi->GetTPCsignal());
 684     }
 685
 686     invMassLambdaAsCascDghter = xi->GetEffMass(); // from V0
 687     dcaV0DaughtersXi = xi->GetDcaV0Daughters();
 688     V0Chi2Xi = xi->GetChi2V0();
 689     V0CosOfPointingAngleXi
 690       = xi->GetV0CosineOfPointingAngle(bestPrimaryVtxPos[0],
 691                                        bestPrimaryVtxPos[1],
 692                                        bestPrimaryVtxPos[2]);
 693     dcaV0ToPrimaryVtxXi = xi->GetD(bestPrimaryVtxPos[0],
 694                                    bestPrimaryVtxPos[1],
 695                                    bestPrimaryVtxPos[2]);
 696     dcaBachToPrimaryVtxXi = TMath::Abs(bTrkXi->GetD(bestPrimaryVtxPos[0],
 697                                                     bestPrimaryVtxPos[1],
 698                                                     bestPrimaryVtxPos[2]));
 699
 700     //V0
 701     xi->GetXYZ(posV0Xi[0], posV0Xi[1], posV0Xi[2]);
 702     V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
 703                              +posV0Xi[1]*posV0Xi[1]
 704                              +posV0Xi[2]*posV0Xi[2]);
 705     dcaPosToPrimaryVtxXi = TMath::Abs(pTrkXi->GetD(bestPrimaryVtxPos[0],
 706                                                    bestPrimaryVtxPos[1],
 707                                                    bestPrimaryVtxPos[2]));
 708     dcaNegToPrimaryVtxXi = TMath::Abs(nTrkXi->GetD(bestPrimaryVtxPos[0],
 709                                                    bestPrimaryVtxPos[1],
 710                                                    bestPrimaryVtxPos[2]));
 711
 712     //apply cuts
 713     //if(XiRadius < 0.9 || XiRadius > 100.) continue;
 714     //if(dcaXiDaughters > 0.2) continue;
 715     //if(XiCosOfPointingAngle < 0.99) continue;
 716     //if(dcaV0ToPrimaryVtxXi < 0.03) continue;
 717     //if(dcaBachToPrimaryVtxXi < 0.01) continue;
 718
 719     if(dcaXiDaughters > fCascadeCuts[0]) continue;
 720     if(XiCosOfPointingAngle < fCascadeCuts[1]) continue;
 721     if(dcaV0ToPrimaryVtxXi < fCascadeCuts[2]) continue;
 722     if(dcaBachToPrimaryVtxXi < fCascadeCuts[3]) continue;
 723
 724     //V0 mass cut?
 725     //  if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.01) continue;
 726     if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > fCascadeCuts[7])
 727       continue;
 728
 729     //if(dcaV0DaughtersXi > 1.) continue;
 730     //if(V0CosOfPointingAngleXi > 0.9999) continue;
 731     //if(dcaPosToPrimaryVtxXi < 0.1) continue;
 732     //if(dcaNegToPrimaryVtxXi < 0.1) continue;
 733
 734     if(dcaV0DaughtersXi > fCascadeCuts[4]) continue;
 735     if(V0CosOfPointingAngleXi > fCascadeCuts[5]) continue;
 736     if(dcaPosToPrimaryVtxXi < fCascadeCuts[6]) continue;
 737     if(dcaNegToPrimaryVtxXi < fCascadeCuts[6]) continue;
 738
 739     //if(V0RadiusXi < 1.0 || V0RadiusXi > 100) continue;
 740
 741     //other cuts?
 742     // change mass hypothesis to cover all the possibilities
 743     if(bTrkXi->Charge()<0){
 744       V0quality = 0.;
 745       xi->ChangeMassHypothesis(V0quality, 3312); //Xi- hyp.
 746       invMassXiMinus = xi->GetEffMassXi();
 747
 748       V0quality = 0.;
 749       xi->ChangeMassHypothesis(V0quality, 3334); //Omega- hyp.
 750       invMassOmegaMinus = xi->GetEffMassXi();
 751
 752       V0quality = 0.;
 753       xi->ChangeMassHypothesis(V0quality, 3312); //back to default hyp.
 754     }
 755
 756     if(bTrkXi->Charge() > 0){
 757       V0quality = 0.;
 758       xi->ChangeMassHypothesis(V0quality, -3312); //anti-Xi- hyp.
 759       invMassXiPlus = xi->GetEffMassXi();
 760
 761       V0quality = 0.;
 762       xi->ChangeMassHypothesis(V0quality, -3334); //anti-Omega- hyp.
 763       invMassOmegaPlus = xi->GetEffMassXi();
 764
 765       V0quality = 0.;
 766       xi->ChangeMassHypothesis(V0quality, -3312); //back to default hyp.
 767     }
 768
 769     //PID on the daughter tracks
 770     /*
 771     //A - Combined PID
 772     //Resonable guess the priors for the cascade track sample
 773     //(e, mu, pi, K, p)
 774     Double_t priorsGuessXi[5] = {0, 0, 2, 0, 1};
 775     Double_t priorsGuessOmega[5] = {0, 0, 1, 1, 1};
 776
 777     //Combined bachelor-daughter PID
 778     AliPID pidXi;
 779     pidXi.SetPriors(priorsGuessXi);
 780     AliPID pidOmega;
 781     pidOmega.SetPriors(priorsGuessOmega);
 782
 783     if(pTrkXi->IsOn(AliESDtrack::kESDpid)){// combined PID exists
 784       Double_t r[10] = {0.};
 785       pTrkXi->GetESDpid(r);
 786       pidXi.SetProbabilities(r);
 787       pidOmega.SetProbabilities(r);
 788
 789       //Check if the V0 postive track is proton (case for Xi-)
 790       Double_t pProton = pidXi.GetProbability(AliPID::kProton);
 791       if(pProton > pidXi.GetProbability(AliPID::kElectron)
 792          && pProton > pidXi.GetProbability(AliPID::kMuon)
 793          && pProton > pidXi.GetProbability(AliPID::kPion)
 794          && pProton > pidXi.GetProbability(AliPID::kKaon))
 795         isPosInXiProton = kTRUE;
 796
 797       //Check if the V0 postive track is a pi+ (case for Xi+)
 798       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
 799       if(pPion > pidXi.GetProbability(AliPID::kElectron)
 800          && pPion > pidXi.GetProbability(AliPID::kMuon)
 801          && pPion > pidXi.GetProbability(AliPID::kKaon)
 802          && pPion > pidXi.GetProbability(AliPID::kProton))
 803         isPosInXiPion = kTRUE;
 804       // Check if the V0 positive track is a proton (case for Omega-)
 805       pProton = pidOmega.GetProbability(AliPID::kProton);
 806       if(pProton > pidOmega.GetProbability(AliPID::kElectron)
 807          && pProton > pidOmega.GetProbability(AliPID::kMuon)
 808          && pProton > pidOmega.GetProbability(AliPID::kPion)
 809          && pProton > pidOmega.GetProbability(AliPID::kKaon))
 810         isPosInOmegaProton = kTRUE;
 811
 812       // Check if the V0 positive track is a pi+ (case for Omega+)
 813       pPion =  pidOmega.GetProbability(AliPID::kPion);
 814       if(pPion > pidOmega.GetProbability(AliPID::kElectron)
 815          && pPion > pidOmega.GetProbability(AliPID::kMuon)
 816          && pPion > pidOmega.GetProbability(AliPID::kKaon)
 817          && pPion > pidOmega.GetProbability(AliPID::kProton))
 818         isPosInOmegaPion = kTRUE;
 819     }
 820
 821     //Combined V0-negative-daughter PID
 822     pidXi.SetPriors(priorsGuessXi);
 823     pidOmega.SetPriors(priorsGuessOmega);
 824     if(nTrkXi->IsOn(AliESDtrack::kESDpid)){
 825       Double_t r[10] = {0.};
 826       nTrkXi->GetESDpid(r);
 827       pidXi.SetProbabilities(r);
 828       pidOmega.SetProbabilities(r);
 829
 830       // Check if the V0 negative track is a pi- (case for Xi-)
 831       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
 832       if(pPion > pidXi.GetProbability(AliPID::kElectron)
 833          && pPion > pidXi.GetProbability(AliPID::kMuon)
 834          && pPion >  pidXi.GetProbability(AliPID::kKaon)
 835          && pPion > pidXi.GetProbability(AliPID::kProton))
 836         isNegInXiPion = kTRUE;
 837
 838       // Check if the V0 negative track is an anti-proton (case for Xi+)
 839       Double_t pProton = pidXi.GetProbability(AliPID::kProton);
 840       if(pProton > pidXi.GetProbability(AliPID::kElectron)
 841          && pProton > pidXi.GetProbability(AliPID::kMuon)
 842          && pProton > pidXi.GetProbability(AliPID::kPion)
 843          && pProton > pidXi.GetProbability(AliPID::kKaon))
 844         isNegInXiProton = kTRUE;
 845
 846       // Check if the V0 negative track is a pi- (case for Omega-)
 847       pPion = pidOmega.GetProbability(AliPID::kPion);
 848       if(pPion > pidOmega.GetProbability(AliPID::kElectron)
 849          && pPion > pidOmega.GetProbability(AliPID::kMuon)
 850          && pPion > pidOmega.GetProbability(AliPID::kKaon)
 851          && pPion > pidOmega.GetProbability(AliPID::kProton))
 852         isNegInOmegaPion = kTRUE;
 853
 854       // Check if the V0 negative track is an anti-proton (case for Omega+)
 855       pProton =  pidOmega.GetProbability(AliPID::kProton);
 856       if(pProton > pidOmega.GetProbability(AliPID::kElectron)
 857          && pProton > pidOmega.GetProbability(AliPID::kMuon)
 858          && pProton > pidOmega.GetProbability(AliPID::kPion)
 859          && pProton > pidOmega.GetProbability(AliPID::kKaon))
 860         isNegInOmegaProton = kTRUE;
 861
 862     }
 863
 864     // Combined bachelor PID
 865     pidXi.SetPriors(priorsGuessXi);
 866     pidOmega.SetPriors(priorsGuessOmega);
 867     if(bTrkXi->IsOn(AliESDtrack::kESDpid)){//Combined PID exists
 868       Double_t r[10] = {0.};
 869       bTrkXi->GetESDpid(r);
 870       pidXi.SetProbabilities(r);
 871       pidOmega.SetProbabilities(r);
 872
 873       //Check if the bachelor track is a pion
 874       Double_t pPion = pidXi.GetProbability(AliPID::kPion);
 875       if(pPion >  pidXi.GetProbability(AliPID::kElectron)
 876          && pPion >  pidXi.GetProbability(AliPID::kMuon)
 877          && pPion > pidXi.GetProbability(AliPID::kKaon)
 878          && pPion > pidXi.GetProbability(AliPID::kProton))
 879         isBachelorPion = kTRUE;
 880
 881       // Check if the bachelor track is a kaon
 882       Double_t pKaon = pidOmega.GetProbability(AliPID::kKaon);
 883       if(pKaon > pidOmega.GetProbability(AliPID::kElectron)
 884          && pKaon > pidOmega.GetProbability(AliPID::kMuon)
 885          && pKaon > pidOmega.GetProbability(AliPID::kPion)
 886          && pKaon > pidOmega.GetProbability(AliPID::kProton))
 887         isBachelorKaon = kTRUE;
 888     }
 889     */
 890
 891
 892     //B - TPC PID: 3-sigma bands on Bethe-Bloch curve
 893     //Bachelor
 894     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
 895       isBachelorKaonForTPC = kTRUE;
 896     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
 897       isBachelorPionForTPC = kTRUE;
 898
 899     //Negative V0 daughter
 900     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
 901       isNegPionForTPC = kTRUE;
 902     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
 903       isNegProtonForTPC = kTRUE;
 904
 905     //Positive V0 daughter
 906     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
 907       isPosPionForTPC = kTRUE;
 908     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
 909       isPosProtonForTPC = kTRUE;
 910
 911
 912     //Extra QA information
 913     xi->GetPxPyPz(XiPx, XiPy, XiPz);
 914     XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
 915     XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
 916
 917     XiPAfter[0] = XiPx;
 918     XiPAfter[1] = XiPy;
 919     XiPAfter[2] = XiPz;
 920
 921     xi->GetBPxPyPz(bachPx, bachPy, bachPz);
 922     bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
 923     bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
 924
 925     //chargeXi = xi->Charge();
 926
 927     V0toXiCosOfPointingAngle
 928       = xi->GetV0CosineOfPointingAngle(posXi[0], posXi[1], posXi[2]);
 929     rapXi = xi->RapXi();
 930     rapOmega = xi->RapOmega();
 931     phi = xi->Phi();
 932     alphaXi = xi->AlphaXi();
 933     ptArmXi = xi->PtArmXi();
 934
 935     distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
 936     distToVtxXYBefore
 937       = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
 938                     *(posXi[0] - bestPrimaryVtxPos[0])
 939                     +(posXi[1] - bestPrimaryVtxPos[1])
 940                     *(posXi[1] - bestPrimaryVtxPos[1]));
 941
 942     //propagation to the best primary vertex to determine the momentum
 943     Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->Charge());
 944     distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
 945     distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
 946                                    *(posXi[0] - bestPrimaryVtxPos[0])
 947                                    +(posXi[1] - bestPrimaryVtxPos[1])
 948                                    *(posXi[1] - bestPrimaryVtxPos[1]));
 949     phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
 950
 951     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxZAfter"))->Fill(distToVtxZAfter);
 952     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxXYAfter"))->Fill(distToVtxXYAfter);
 953     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxZBeforeVsAfter"))->Fill(distToVtxZBefore, distToVtxZAfter);
 954     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxXYBeforeVsAfter"))->Fill(distToVtxXYBefore, distToVtxXYAfter);
 955     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PxBeforeVsAfter"))->Fill(XiPx, XiPAfter[0]);
 956     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PyBeforeVsAfter"))->Fill(XiPy, XiPAfter[1]);
 957     if(xi->Charge()>0)
 958       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PhiPosBeforeVsAfter"))->Fill(phi, phiAfter);
 959     else if(xi->Charge()<0)
 960       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PhiNegBeforeVsAfter"))->Fill(phi, phiAfter);
 961
 962
 963     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("Chi2Xi"))->Fill(chi2Xi);
 964     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaXiDaughters"))->Fill(dcaXiDaughters);
 965     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaBachToPrimVertex"))->Fill(dcaBachToPrimaryVtxXi);
 966     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("XiCosOfPointingAngle"))->Fill(XiCosOfPointingAngle);
 967     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("XiRadius"))->Fill(XiRadius);
 968
 969     //V0
 970     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassLambdaAsCascDghter"))->Fill(invMassLambdaAsCascDghter);
 971     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0Chi2Xi"))->Fill(V0Chi2Xi);
 972     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaV0DaughtersXi"))->Fill(dcaV0DaughtersXi);
 973     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0CosOfPointingAngleXi"))->Fill(V0CosOfPointingAngleXi);
 974     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0RadiusXi"))->Fill(V0RadiusXi);
 975     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaV0ToPrimVertexXi"))->Fill(dcaV0ToPrimaryVtxXi);
 976     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaPosToPrimVertexXi"))->Fill(dcaPosToPrimaryVtxXi);
 977     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaNegToPrimVertexXi"))->Fill(dcaNegToPrimaryVtxXi);
 978
 979     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0toXiCosOfPointingAngle"))->Fill(V0toXiCosOfPointingAngle);
 980
 981     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("Armenteros"))->Fill(alphaXi, ptArmXi);
 982
 983     //PID cuts with TPC cuts
 984     if(xi->Charge() < 0){
 985       if(isPosProtonForTPC
 986          && isNegPionForTPC){
 987
 988         switch(fSpecie) {
 989         case 0:
 990           if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
 991             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassXiMinus);
 992
 993           //candidate inserting
 994             MakeTrack(invMassXiMinus, XiPt, /*xi->Phi()*/
 995                       phiAfter, xi->Eta(),  pTrkXi->GetID(),
 996                       nTrkXi->GetID(), bTrkXi->GetID());
 997           }
 998           break;
 999
1000         case 1:
1001           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1002             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassOmegaMinus);
1003             MakeTrack(invMassOmegaMinus, XiPt, /*xi->Phi()*/
1004                       phiAfter, xi->Eta(),
1005                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1006           }
1007           break;
1008         }
1009       }
1010     }
1011
1012     if(xi->Charge() > 0){
1013       if(isNegProtonForTPC
1014          && isPosPionForTPC){
1015
1016         switch (fSpecie){
1017         case 0:
1018           if(isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1019             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassXiPlus);
1020
1021           //candidate inserting
1022             MakeTrack(invMassXiPlus, XiPt, /*xi->Phi()*/
1023                       phiAfter, xi->Eta(),
1024                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1025           }
1026           break;
1027
1028         case 1:
1029           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8){
1030             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassOmegaPlus);
1031             MakeTrack(invMassOmegaPlus, XiPt, /*xi->Phi()*/
1032                       phiAfter, xi->Eta(),
1033                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1034           }
1035           break;
1036         }
1037       }
1038     }
1039
1040   }
1041
1042   return;
1043 }
1044
1045 //______________________________________________________________________________
1046 void AliAnalysisTaskFlowCascade::ReadFromAODv0(AliAODEvent *fAOD)
1047 {
1048
1049   fCutsRPTPC->SetEvent(fAOD, MCEvent());
1050   fCutsRPVZE->SetEvent(fAOD, MCEvent());
1051   fCutsPOI->SetEvent(fAOD, MCEvent());
1052   fFlowEventTPC->Fill(fCutsRPTPC, fCutsPOI);
1053   fFlowEventVZE->Fill(fCutsRPVZE, fCutsPOI);
1054
1055   //  Double_t trkPrimaryVtxPos[3] = {-100., -100., -100.};
1056   Double_t bestPrimaryVtxPos[3] = {-100., -100., -100.};
1057
1058   Double_t b = fAOD->GetMagneticField();
1059
1060   int nCascades=fAOD->GetNumberOfCascades();
1061   const AliAODVertex *primaryBestAODVtx = fAOD->GetPrimaryVertex();
1062   primaryBestAODVtx->GetXYZ(bestPrimaryVtxPos);
1063
1064   // calculation part dedicated to Xi vertices
1065   for(Int_t iXi = 0; iXi < nCascades; iXi++){
1066     Double_t effMassXi = 0.;
1067     Double_t chi2Xi = -1.;
1068     Double_t dcaXiDaughters = -1.;
1069     Double_t XiCosOfPointingAngle = -1.;
1070     Double_t posXi[3] = {-1000., -1000., -1000.};
1071     Double_t XiRadius = -1000.;
1072
1073     Double_t invMassLambdaAsCascDghter = 0.;
1074     Double_t V0Chi2Xi = -1.;
1075     Double_t dcaV0DaughtersXi = -1.;
1076
1077     Double_t dcaBachToPrimaryVtxXi = -1.;
1078     Double_t dcaV0ToPrimaryVtxXi = -1.;
1079     Double_t dcaPosToPrimaryVtxXi = -1.;
1080     Double_t dcaNegToPrimaryVtxXi = -1.;
1081     Double_t V0CosOfPointingAngleXi = -1.;
1082     Double_t posV0Xi[3] = {-1000., -1000., -1000.};
1083     Double_t V0RadiusXi = -1000.;
1084     //    Double_t V0quality = 0.;
1085
1086     Double_t invMassXiMinus = 0.;
1087     Double_t invMassXiPlus = 0.;
1088     Double_t invMassOmegaMinus = 0.;
1089     Double_t invMassOmegaPlus = 0.;
1090
1091     /*
1092     Bool_t isPosInXiProton = kFALSE;
1093     Bool_t isPosInXiPion = kFALSE;
1094     Bool_t isPosInOmegaProton = kFALSE;
1095     Bool_t isPosInOmegaPion = kFALSE;
1096
1097     Bool_t isNegInXiProton = kFALSE;
1098     Bool_t isNegInXiPion = kFALSE;
1099     Bool_t isNegInOmegaProton = kFALSE;
1100     Bool_t isNegInOmegaPion = kFALSE;
1101
1102     Bool_t isBachelorKaon = kFALSE;
1103     Bool_t isBachelorPion = kFALSE;
1104     */
1105
1106
1107     Bool_t isBachelorKaonForTPC = kFALSE;
1108     Bool_t isBachelorPionForTPC = kFALSE;
1109     Bool_t isNegPionForTPC = kFALSE;
1110     Bool_t isPosPionForTPC = kFALSE;
1111     Bool_t isNegProtonForTPC = kFALSE;
1112     Bool_t isPosProtonForTPC = kFALSE;
1113
1114     Double_t XiPx = 0., XiPy = 0., XiPz = 0.;
1115     Double_t XiPt = 0.;
1116     Double_t XiPtot = 0.;
1117
1118     Double_t bachPx = 0., bachPy = 0., bachPz = 0.;
1119     Double_t bachPt = 0.;
1120     Double_t bachPtot = 0.;
1121
1122     //Short_t chargeXi = -2;
1123     Double_t V0toXiCosOfPointingAngle = 0.;
1124
1125     Double_t rapXi = -20.;
1126     Double_t rapOmega = -20.;
1127     Double_t phi = 6.3;
1128     Double_t alphaXi = -200.;
1129     Double_t ptArmXi = -200.;
1130
1131     Double_t distToVtxZBefore = -999.;
1132     Double_t distToVtxZAfter = -999.;
1133     Double_t distToVtxXYBefore = -999.;
1134     Double_t distToVtxXYAfter = -999.;
1135     Double_t XiPAfter[3] = {-999., -999., -999.};
1136     Double_t phiAfter = -999.;
1137
1138     const AliAODcascade *xi = fAOD->GetCascade(iXi);
1139     if (!xi) continue;
1140
1141     effMassXi = xi->MassXi(); //default working hypothesis: Xi- decay
1142     chi2Xi = xi->Chi2Xi();
1143     dcaXiDaughters = xi->DcaXiDaughters();
1144     XiCosOfPointingAngle = xi->CosPointingAngleXi(bestPrimaryVtxPos[0],
1145                                                   bestPrimaryVtxPos[1],
1146                                                   bestPrimaryVtxPos[2]);
1147     posXi[0] = xi->DecayVertexXiX();
1148     posXi[1] = xi->DecayVertexXiY();
1149     posXi[2] = xi->DecayVertexXiZ();
1150     XiRadius = TMath::Sqrt(posXi[0]*posXi[0]
1151                            +posXi[1]*posXi[1]
1152                            +posXi[2]*posXi[2]);
1153
1154     AliAODTrack *pTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
1155     AliAODTrack *nTrkXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
1156     AliAODTrack *bTrkXi
1157       = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
1158
1159     if(!pTrkXi || !nTrkXi || !bTrkXi) continue;
1160
1161     UInt_t idxPosXi  = (UInt_t) TMath::Abs( pTrkXi->GetID() );
1162     UInt_t idxNegXi  = (UInt_t) TMath::Abs( nTrkXi->GetID() );
1163     UInt_t idxBach   = (UInt_t) TMath::Abs( bTrkXi->GetID() );
1164
1165     if(idxBach == idxNegXi || idxBach == idxPosXi) continue;
1166
1167     if( !fCutsDau->IsSelected(pTrkXi)
1168         || !fCutsDau->IsSelected(nTrkXi)
1169         || !fCutsDau->IsSelected(bTrkXi) ) continue;
1170
1171
1172     if(pTrkXi->IsOn(AliESDtrack::kTPCin)){
1173       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(pTrkXi->P()*pTrkXi->Charge(), pTrkXi->GetTPCsignal());
1174     }
1175     if( nTrkXi->IsOn(AliESDtrack::kTPCin) ){
1176       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(nTrkXi->P()*nTrkXi->Charge(), nTrkXi->GetTPCsignal());
1177     }
1178     if(bTrkXi->IsOn(AliESDtrack::kTPCin)){
1179       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("TPCdEdxOfCascDghters"))->Fill(bTrkXi->P()*bTrkXi->Charge(), bTrkXi->GetTPCsignal());
1180     }
1181
1182     if(xi->ChargeXi() < 0)
1183       invMassLambdaAsCascDghter = xi->MassLambda();
1184     else
1185       invMassLambdaAsCascDghter = xi->MassAntiLambda();
1186
1187     dcaV0DaughtersXi = xi->DcaV0Daughters();
1188     V0Chi2Xi = xi->Chi2V0();
1189     V0CosOfPointingAngleXi
1190       = xi->CosPointingAngle(bestPrimaryVtxPos);
1191     dcaV0ToPrimaryVtxXi = xi->DcaV0ToPrimVertex();
1192     dcaBachToPrimaryVtxXi = xi->DcaBachToPrimVertex();
1193
1194     //V0
1195     posV0Xi[0] = xi->DecayVertexV0X();
1196     posV0Xi[1] = xi->DecayVertexV0Y();
1197     posV0Xi[2] = xi->DecayVertexV0Z();
1198     V0RadiusXi = TMath::Sqrt(posV0Xi[0]*posV0Xi[0]
1199                              +posV0Xi[1]*posV0Xi[1]
1200                              +posV0Xi[2]*posV0Xi[2]);
1201     dcaPosToPrimaryVtxXi = xi->DcaPosToPrimVertex();
1202     dcaNegToPrimaryVtxXi = xi->DcaNegToPrimVertex();
1203
1204     //apply cuts ?
1205     // if(XiRadius < 1. || XiRadius > 100.) continue;
1206     //if(dcaXiDaughters > 0.1) continue;
1207     //if(XiCosOfPointingAngle < 0.999) continue;
1208     //if(dcaV0ToPrimaryVtxXi < 0.05) continue;
1209     //if(dcaBachToPrimaryVtxXi < 0.03) continue;
1210
1211     if(dcaXiDaughters > fCascadeCuts[0]) continue;
1212     if(XiCosOfPointingAngle < fCascadeCuts[1]) continue;
1213     if(dcaV0ToPrimaryVtxXi < fCascadeCuts[2]) continue;
1214     if(dcaBachToPrimaryVtxXi < fCascadeCuts[3]) continue;
1215
1216     //V0 mass cut?
1217     //if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > 0.006) continue;
1218     if(TMath::Abs(invMassLambdaAsCascDghter-1.11568) > fCascadeCuts[7])
1219       continue;
1220
1221     //if(dcaV0DaughtersXi > 1.) continue;
1222     //if(V0CosOfPointingAngleXi > 0.9999) continue;
1223     //if(dcaPosToPrimaryVtxXi < 0.1) continue;
1224     //if(dcaNegToPrimaryVtxXi < 0.1) continue;
1225     if(dcaV0DaughtersXi > fCascadeCuts[4]) continue;
1226     if(V0CosOfPointingAngleXi > fCascadeCuts[5]) continue;
1227     if(dcaPosToPrimaryVtxXi < fCascadeCuts[6]) continue;
1228     if(dcaNegToPrimaryVtxXi < fCascadeCuts[6]) continue;
1229
1230     // if(V0RadiusXi < 1.0 || V0RadiusXi > 100) continue;
1231
1232     //other cuts?
1233
1234
1235     //???
1236     if(xi->ChargeXi()<0){
1237       invMassXiMinus = xi->MassXi();
1238       invMassOmegaMinus = xi->MassOmega();
1239     }else{
1240       invMassXiPlus = xi->MassXi();
1241       invMassOmegaPlus = xi->MassOmega();
1242     }
1243
1244     /*
1245     if(pTrkXi->GetMostProbablePID() == AliAODTrack::kProton) {
1246       isPosInXiProton = kTRUE;
1247       isPosInOmegaProton = kTRUE;
1248     }
1249     if(pTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1250       isPosInXiPion = kTRUE;
1251       isPosInOmegaPion = kTRUE;
1252     }
1253
1254     if(nTrkXi->GetMostProbablePID() == AliAODTrack::kPion){
1255       isNegInXiPion = kTRUE;
1256       isNegInOmegaPion = kTRUE;
1257     }
1258     if(nTrkXi->GetMostProbablePID() == AliAODTrack::kProton){
1259       isNegInXiProton = kTRUE;
1260       isNegInOmegaProton = kTRUE;
1261     }
1262
1263     if(bTrkXi->GetMostProbablePID() == AliAODTrack::kPion)
1264       isBachelorPion = kTRUE;
1265     if(bTrkXi->GetMostProbablePID() == AliAODTrack::kKaon)
1266       isBachelorKaon = kTRUE;
1267     */
1268
1269     //PID with TPC only:
1270     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kKaon))<3.)
1271       isBachelorKaonForTPC = kTRUE;
1272     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bTrkXi, AliPID::kPion))<3.)
1273       isBachelorPionForTPC = kTRUE;
1274
1275     //Negative V0 daughter
1276     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kPion))<3.)
1277       isNegPionForTPC = kTRUE;
1278     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrkXi, AliPID::kProton))<3.)
1279       isNegProtonForTPC = kTRUE;
1280
1281     //Positive V0 daughter
1282     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kPion))<3.)
1283       isPosPionForTPC = kTRUE;
1284     if(TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrkXi, AliPID::kProton))<3.)
1285       isPosProtonForTPC = kTRUE;
1286
1287     //Extra QA information
1288     XiPx = xi->MomXiX();
1289     XiPy = xi->MomXiY();
1290     XiPz = xi->MomXiZ();
1291     XiPt = TMath::Sqrt(XiPx*XiPx + XiPy*XiPy);
1292     XiPtot= TMath::Sqrt(XiPx*XiPx + XiPy*XiPy + XiPz*XiPz);
1293
1294     bachPx = xi->MomBachX();
1295     bachPy = xi->MomBachY();
1296     bachPz = xi->MomBachZ();
1297
1298     bachPt = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy);
1299     bachPtot = TMath::Sqrt(bachPx*bachPx + bachPy*bachPy + bachPz*bachPz);
1300
1301     V0toXiCosOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
1302
1303     rapXi = xi->RapXi();
1304     rapOmega = xi->RapOmega();
1305     phi = xi->Phi();
1306     alphaXi = xi->AlphaXi();
1307     ptArmXi = xi->PtArmXi();
1308
1309     distToVtxZBefore = posXi[2]-bestPrimaryVtxPos[2];
1310     distToVtxXYBefore
1311       = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1312                     *(posXi[0] - bestPrimaryVtxPos[0])
1313                     +(posXi[1] - bestPrimaryVtxPos[1])
1314                     *(posXi[1] - bestPrimaryVtxPos[1]));
1315
1316
1317     XiPAfter[0] = XiPx;
1318     XiPAfter[1] = XiPy;
1319     XiPAfter[2] = XiPz;
1320     //propagation to the best primary vertex to determine the momentum
1321     Propagate(bestPrimaryVtxPos, posXi, XiPAfter, b, xi->ChargeXi());
1322     distToVtxZAfter = posXi[2] - bestPrimaryVtxPos[2];
1323     distToVtxXYAfter = TMath::Sqrt((posXi[0] - bestPrimaryVtxPos[0])
1324                                    *(posXi[0] - bestPrimaryVtxPos[0])
1325                                    +(posXi[1] - bestPrimaryVtxPos[1])
1326                                    *(posXi[1] - bestPrimaryVtxPos[1]));
1327     phiAfter = TMath::Pi() + TMath::ATan2(-XiPAfter[1],-XiPAfter[0]);
1328
1329     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxZAfter"))->Fill(distToVtxZAfter);
1330     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxXYAfter"))->Fill(distToVtxXYAfter);
1331     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxZBeforeVsAfter"))->Fill(distToVtxZBefore, distToVtxZAfter);
1332     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DistToVtxXYBeforeVsAfter"))->Fill(distToVtxXYBefore, distToVtxXYAfter);
1333     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PxBeforeVsAfter"))->Fill(XiPx, XiPAfter[0]);
1334     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PyBeforeVsAfter"))->Fill(XiPy, XiPAfter[1]);
1335     if(xi->ChargeXi()>0)
1336       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PhiPosBeforeVsAfter"))->Fill(phi, phiAfter);
1337     else if(xi->ChargeXi()<0)
1338       ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("PhiNegBeforeVsAfter"))->Fill(phi, phiAfter);
1339
1340     //for default hypothesis
1341     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("Chi2Xi"))->Fill(chi2Xi);
1342     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaXiDaughters"))->Fill(dcaXiDaughters);
1343     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaBachToPrimVertex"))->Fill(dcaBachToPrimaryVtxXi);
1344     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("XiCosOfPointingAngle"))->Fill(XiCosOfPointingAngle);
1345     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("XiRadius"))->Fill(XiRadius);
1346
1347     //V0
1348     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassLambdaAsCascDghter"))->Fill(invMassLambdaAsCascDghter);
1349     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0Chi2Xi"))->Fill(V0Chi2Xi);
1350     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaV0DaughtersXi"))->Fill(dcaV0DaughtersXi);
1351     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0CosOfPointingAngleXi"))->Fill(V0CosOfPointingAngleXi);
1352     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0RadiusXi"))->Fill(V0RadiusXi);
1353     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaV0ToPrimVertexXi"))->Fill(dcaV0ToPrimaryVtxXi);
1354     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaPosToPrimVertexXi"))->Fill(dcaPosToPrimaryVtxXi);
1355     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("DcaNegToPrimVertexXi"))->Fill(dcaNegToPrimaryVtxXi);
1356
1357     ((TH1F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("V0toXiCosOfPointingAngle"))->Fill(V0toXiCosOfPointingAngle);
1358
1359     ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("Armenteros"))->Fill(alphaXi, ptArmXi);
1360
1361     //with PID cuts
1362     if(xi->ChargeXi()<0){
1363       if(isPosProtonForTPC && isNegPionForTPC){
1364         switch (fSpecie){
1365         case 0:
1366           if( isBachelorPionForTPC && TMath::Abs(rapXi) < 0.8){
1367             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassXiMinus);
1368             MakeTrack(invMassXiMinus, XiPt, /*xi->Phi(),*/
1369                       phiAfter, xi->Eta(),
1370                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1371           }// endif
1372
1373           break;
1374
1375         case 1:
1376           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8
1377              && (invMassXiMinus > 1.32486 || invMassXiMinus < 1.30486)){
1378             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassOmegaMinus);
1379
1380             MakeTrack(invMassOmegaMinus, XiPt, /* xi->Phi(),*/
1381                       phiAfter,   xi->Eta(),
1382                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1383           }//endif
1384           break;
1385         }
1386       }
1387     }//end if ChargeXi()<0
1388
1389     if(xi->ChargeXi() > 0){
1390       if(isNegProtonForTPC && isPosPionForTPC){
1391         switch(fSpecie){
1392         case 0:
1393           if (isBachelorPionForTPC  && TMath::Abs(rapXi) < 0.8){
1394             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassXiPlus);
1395
1396             //candidate inserting
1397             MakeTrack(invMassXiPlus, XiPt, /* xi->Phi(),*/
1398                       phiAfter, xi->Eta(),
1399                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1400           }//endif particle id
1401           break;
1402
1403         case 1:
1404           if(isBachelorKaonForTPC && TMath::Abs(rapOmega) < 0.8
1405              && (invMassXiPlus > 1.32486 || invMassXiPlus < 1.30486)){
1406             ((TH2F*)((TList*)fQAList->FindObject("Candidates"))->FindObject("MassVsPtAll"))->Fill(XiPt, invMassOmegaPlus);
1407             MakeTrack(invMassOmegaPlus, XiPt, /* xi->Phi(),*/
1408                       phiAfter, xi->Eta(),
1409                       pTrkXi->GetID(), nTrkXi->GetID(), bTrkXi->GetID());
1410           }//endif particle id
1411         }
1412       }
1413     }//endif ChargeXi()>0
1414   }//for Xi candidate loop
1415
1416   return;
1417
1418 }
1419
1420
1421 void AliAnalysisTaskFlowCascade::MakeTrack( double mass,
1422                                             double pt,
1423                                             double phi,
1424                                             double eta,
1425                                             int iid,
1426                                             int jid,
1427                                             int kid) {
1428   // create track for flow tasks
1429   if(fCandidates->GetLast()+1>=fCandidates->GetSize()) {
1430     fCandidates->Expand( 2*fCandidates->GetSize() );
1431   }
1432   Bool_t overwrite = kTRUE;
1433
1434   AliFlowCandidateTrack *sTrack
1435     = (static_cast<AliFlowCandidateTrack*> (fCandidates->At( fCandidates->GetLast()+1 )));
1436   if( !sTrack ) { // creates new
1437     sTrack = new AliFlowCandidateTrack();
1438     overwrite = kFALSE;
1439   } else { // overwrites
1440     sTrack->ClearMe();
1441   }
1442
1443
1444   sTrack->SetMass(mass);
1445   sTrack->SetPt(pt);
1446   sTrack->SetPhi(phi);
1447   sTrack->SetEta(eta);
1448   sTrack->AddDaughter(iid);
1449   sTrack->AddDaughter(jid);
1450   sTrack->AddDaughter(kid);
1451   sTrack->SetForPOISelection(kTRUE);
1452   sTrack->SetForRPSelection(kFALSE);
1453
1454   if(overwrite) {
1455     fCandidates->SetLast( fCandidates->GetLast()+1 );
1456   } else {
1457     fCandidates->AddLast(sTrack);
1458   }
1459
1460   return;
1461 }
1462 //_____________________________________________________________________________
1463 void AliAnalysisTaskFlowCascade::Terminate(Option_t *)
1464 {
1465
1466 }
1467
1468 void AliAnalysisTaskFlowCascade::Propagate(Double_t vv[3],
1469                                            Double_t x[3],
1470                                            Double_t p[3],
1471                                            Double_t bz,
1472                                            Short_t sign){
1473   //Propagation to the primary vertex to determine the px and py
1474   //x, p are the position and momentum as input and output
1475   //bz is the magnetic field along z direction
1476   //sign is the charge of particle for propagation
1477
1478   Double_t pp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
1479   Double_t len = (vv[2]-x[2])*pp/p[2];
1480   Double_t a = -kB2C*bz*sign;
1481
1482   Double_t rho = a/pp;
1483   x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
1484   x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
1485   x[2] += p[2]*len/pp;
1486
1487   Double_t p0=p[0];
1488   p[0] = p0  *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
1489   p[1] = p[1]*TMath::Cos(rho*len) + p0  *TMath::Sin(rho*len);
1490 }
1491
1492
1493 //=====================================================================
1494 void AliAnalysisTaskFlowCascade::SetCommonConstants(Int_t massBins,
1495                                                     Double_t minMass,
1496                                                     Double_t maxMass)
1497 {
1498   // setter for mass bins
1499
1500   fMassBins = massBins;
1501   fMinMass = minMass;
1502   fMaxMass = maxMass;
1503 }
1504
1505
1506 //====================================================================
1507 void AliAnalysisTaskFlowCascade::SetCuts2010(int set) {
1508
1509   // fCascadeCuts[0]: DcaXiDaughter; fCascadeCuts[1]: XiCosOfPointingAngle
1510   // fCascadeCuts[2]: DcaV0ToPrimaryVtxXi; fCascadeCuts[3]: DcaBachToPrimaryVtxXi
1511   // fCascadeCuts[4]: DcaV0DaughersXi; fCascadeCuts[5]: V0CosOfPointingAngleXi
1512   // fCascadeCuts[6]: DcaV0DaughterToPrimaryVtxXi; fCascadeCuts[7]: V0MassWidth
1513
1514   switch(set){
1515
1516   case 0: //tighter
1517     fCascadeCuts[0] = 0.2; fCascadeCuts[1] = 0.999;
1518     fCascadeCuts[2] = 0.03; fCascadeCuts[3] = 0.05;
1519     fCascadeCuts[4] = .5; fCascadeCuts[5] = 0.9998;
1520     fCascadeCuts[6] = 0.1; fCascadeCuts[7] = 0.006;
1521     break;
1522
1523   case 1: //middle
1524     fCascadeCuts[0] = 0.3; fCascadeCuts[1] = 0.99;
1525     fCascadeCuts[2] = 0.01; fCascadeCuts[3] = 0.03;
1526     fCascadeCuts[4] = .6; fCascadeCuts[5] = 0.9999;
1527     fCascadeCuts[6] = 0.1; fCascadeCuts[7] = 0.008;
1528     break;
1529
1530   case 2: //looser
1531     fCascadeCuts[0] = 0.3; fCascadeCuts[1] = 0.99;
1532     fCascadeCuts[2] = 0.01; fCascadeCuts[3] = 0.03;
1533     fCascadeCuts[4] = 1.; fCascadeCuts[5] = 1.;
1534     fCascadeCuts[6] = 0.1; fCascadeCuts[7] = 0.01;
1535     break;
1536   }
1537
1538 }