PWG/FLOW/Tasks/AliAnalysisTaskFlowCascade.cxx

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