PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskPerformanceStrange.cxx

   1 /* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   2  *                                                                        *
   3  * Author: The ALICE Off-line Project.                                    *
   4  * Contributors are mentioned in the code where appropriate.              *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //--------------------------------------------------------------------------
  16 //              AliAnalysisTaskPerformanceStrange class
  17 //    This task is for a performance study of V0 identification.
  18 //                It works with MC info and ESD tree.
  19 //                 Author: Peter Kalinak  pkalinak@cern.ch kalinak@saske.sk
  20 //--------------------------------------------------------------------------
  21
  22 #include <Riostream.h>
  23
  24 #include <stdio.h>
  25 #include <iostream>
  26 #include "TChain.h"
  27 #include "TTree.h"
  28 #include "TH1F.h"
  29 #include "TH2F.h"
  30 #include "TH3F.h"
  31 #include "TF1.h"
  32 #include "TList.h"
  33 #include "TMath.h"
  34 #include "TCanvas.h"
  35
  36 #include "AliAnalysisManager.h"
  37
  38 #include "AliPhysicsSelection.h"
  39 #include "AliBackgroundSelection.h"
  40
  41 #include "AliESDVertex.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDv0.h"
  46 #include "AliESDcascade.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliESDpid.h"
  49 #include "AliMultiplicity.h"
  50
  51 #include "AliAODEvent.h"
  52 #include "AliAODVertex.h"
  53 #include "AliAODTrack.h"
  54 #include "AliAODv0.h"
  55 #include "AliAODMCHeader.h"
  56 #include "AliAODInputHandler.h"
  57
  58 //#include "AliV0vertexer.h"
  59
  60 #include "AliAODMCParticle.h"
  61
  62 #include "AliMCEventHandler.h"
  63 #include "AliMCEvent.h"
  64 #include "AliStack.h"
  65 #include "AliGenEventHeader.h"
  66
  67 #include "AliLog.h"
  68
  69 #include "AliKFVertex.h"
  70 #include "AliVertexerTracks.h"
  71
  72 #include "AliAnalysisTaskPerformanceStrange.h"
  73 #include "AliAnalysisCentralitySelector.h"
  74 #include "AliPIDResponse.h"
  75 #include "AliCentrality.h"
  76
  77
  78
  79 ClassImp(AliAnalysisTaskPerformanceStrange)
  80
  81
  82 //________________________________________________________________________
  83 AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange()
  84 : AliAnalysisTaskSE(), fAnalysisMC(0), fAnalysisType("infoType"),  fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infoCut"),fDown(0),fUp(0), fESD(0), fListHist(0),fCentrSelector(0),fTracksCuts(0),fPIDResponse(0),fQASelector(0),
  85
  86
  87     /////// primary vertex
  88     fHistMCPrimaryVertexX(0),   fHistMCPrimaryVertexY(0),   fHistMCPrimaryVertexZ(0),
  89     ////// tracks & multiplicity
  90     fHistPtTracks(0),
  91     fHistMCMultiplicityPrimary(0),  fHistMCMultiplicityTracks(0),  fHistTPCTracks(0),
  92     ///// Transverse Momentum
  93     fHistMCPtAllK0s(0),
  94     fHistMCPtAllLambda(0),   fHistMCPtAllAntiLambda(0),
  95     fHistMCPtAllXi(0),  fHistMCPtAllAntiXi(0),
  96     fHistMCPtAllOmega(0), fHistMCPtAllAntiOmega(0),
  97     /// Rapidity
  98     fHistMCRapK0s(0),
  99     fHistMCRapLambda(0), fHistMCRapAntiLambda(0),
 100     fHistMCRapXi(0),
 101
 102     ///// Transverse Momentum primary
 103     fHistMCPtK0s(0),
 104     fHistMCPtLambda(0),  fHistMCPtAntiLambda(0),
 105
 106     ///////////////////////////////////////////
 107     ////   ESD
 108
 109     fHistNumberEvents(0),
 110     fHistTrackPerEvent(0),
 111     fHistTPCMult(0),
 112     fHistTrackletPerEvent(0),
 113     fHistSPDPrimaryVertexZ(0),
 114     fHistPrimaryVertexX(0),
 115     fHistPrimaryVertexY(0),
 116     fHistPrimaryVertexZ(0),
 117
 118     fHistV0Multiplicity(0),
 119     ///// inv. mass
 120     fHistMassK0(0),
 121     fHistMassLambda(0),
 122     fHistMassAntiLambda(0),
 123     fHistMassXi(0),
 124     fHistMassAntiXi(0),
 125     fHistMassOmega(0),
 126     fHistMassAntiOmega(0),
 127
 128     //inv mass vs PID
 129     fHistMassXiVsPID(0),
 130
 131     ///////////////////////////////////////
 132     fHistPtVsMassK0(0),
 133     fHistPtVsMassLambda(0),
 134     fHistPtVsMassAntiLambda(0),
 135
 136     ////////////////////////////////////////
 137
 138     fHistArmenterosPodolanski(0),
 139     fHistK0sMassVsLambdaMass(0),
 140     fHistTPCsigPLambda(0),
 141     fHistTPCsigPAntiLambda(0),
 142     fHistNSigmaProton(0),
 143
 144     /// Associated histos
 145     ///rapidity
 146     fHistAsMcRapK0(0),
 147     fHistAsMcRapLambda(0),   fHistAsMcRapAntiLambda(0),
 148
 149     // pt distribution    /////////////////////
 150     fHistAsMcPtK0(0),
 151     fHistAsMcPtLambda(0),   fHistAsMcPtAntiLambda(0),
 152
 153     /////////////////////////////////////
 154     fHistPidMcMassK0(0),
 155     fHistPidMcMassLambda(0),
 156     fHistPidMcMassAntiLambda(0),
 157
 158     ///inv. mass
 159     fHistAsMcMassK0(0),
 160     fHistAsMcMassLambda(0), fHistAsMcMassAntiLambda(0),
 161     ///transv. momentum
 162     fHistAsMcPtVsMassK0(0),
 163     fHistAsMcPtVsMassLambda(0), fHistAsMcPtVsMassAntiLambda(0),
 164
 165     /////background composition
 166     fHistCompositionXi(0),
 167     fHistCompositionAntiXi(0),
 168     fHistCompositionOmega(0),
 169   fHistCompositionAntiOmega(0),
 170
 171   fHistMCIndexes(0)
 172 {
 173   // Constructor
 174 }
 175
 176 //________________________________________________________________________
 177 AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange(const char *name)
 178   : AliAnalysisTaskSE(name), fAnalysisMC(0), fAnalysisType("infoType"), fCollidingSystems(0), fUsePID("infoPID"), fUseCut("infocut"),fDown(0),fUp(0), fESD(0), fListHist(),fCentrSelector(0), fTracksCuts(0),fPIDResponse(0),fQASelector(0),
 179
 180     /////// primary vertex
 181     fHistMCPrimaryVertexX(0),   fHistMCPrimaryVertexY(0),   fHistMCPrimaryVertexZ(0),
 182     ////// tracks & multiplicity
 183     fHistPtTracks(0),
 184     fHistMCMultiplicityPrimary(0),  fHistMCMultiplicityTracks(0),  fHistTPCTracks(0),
 185     ///// Transverse Momentum
 186     fHistMCPtAllK0s(0),
 187     fHistMCPtAllLambda(0),   fHistMCPtAllAntiLambda(0),
 188     fHistMCPtAllXi(0),  fHistMCPtAllAntiXi(0),
 189     fHistMCPtAllOmega(0), fHistMCPtAllAntiOmega(0),
 190     /// Rapidity
 191     fHistMCRapK0s(0),
 192     fHistMCRapLambda(0), fHistMCRapAntiLambda(0),
 193     fHistMCRapXi(0),
 194
 195     ///// Transverse Momentum primary
 196     fHistMCPtK0s(0),
 197     fHistMCPtLambda(0),  fHistMCPtAntiLambda(0),
 198
 199     ///////////////////////////////////////////
 200     ////   ESD
 201
 202     fHistNumberEvents(0),
 203     fHistTrackPerEvent(0),
 204     fHistTPCMult(0),
 205     fHistTrackletPerEvent(0),
 206     fHistSPDPrimaryVertexZ(0),
 207     fHistPrimaryVertexX(0),
 208     fHistPrimaryVertexY(0),
 209     fHistPrimaryVertexZ(0),
 210
 211     fHistV0Multiplicity(0),
 212     ///// inv. mass
 213     fHistMassK0(0),
 214     fHistMassLambda(0),
 215     fHistMassAntiLambda(0),
 216     fHistMassXi(0),
 217     fHistMassAntiXi(0),
 218     fHistMassOmega(0),
 219     fHistMassAntiOmega(0),
 220
 221     //inv mass vs PID
 222     fHistMassXiVsPID(0),
 223
 224     ///////////////////////////////////////
 225     fHistPtVsMassK0(0),
 226     fHistPtVsMassLambda(0),
 227     fHistPtVsMassAntiLambda(0),
 228
 229     ////////////////////////////////////////
 230
 231     fHistArmenterosPodolanski(0),
 232     fHistK0sMassVsLambdaMass(0),
 233     fHistTPCsigPLambda(0),
 234     fHistTPCsigPAntiLambda(0),
 235     fHistNSigmaProton(0),
 236
 237     /// Associated histos
 238     ///rapidity
 239     fHistAsMcRapK0(0),
 240     fHistAsMcRapLambda(0),   fHistAsMcRapAntiLambda(0),
 241
 242     // pt distribution    /////////////////////
 243     fHistAsMcPtK0(0),
 244     fHistAsMcPtLambda(0),   fHistAsMcPtAntiLambda(0),
 245
 246     /////////////////////////////////////
 247     fHistPidMcMassK0(0),
 248     fHistPidMcMassLambda(0),
 249     fHistPidMcMassAntiLambda(0),
 250
 251     ///inv. mass
 252     fHistAsMcMassK0(0),
 253     fHistAsMcMassLambda(0), fHistAsMcMassAntiLambda(0),
 254     ///transv. momentum
 255     fHistAsMcPtVsMassK0(0),
 256     fHistAsMcPtVsMassLambda(0), fHistAsMcPtVsMassAntiLambda(0),
 257
 258     /////background composition
 259     fHistCompositionXi(0),
 260     fHistCompositionAntiXi(0),
 261     fHistCompositionOmega(0),
 262     fHistCompositionAntiOmega(0),
 263
 264     fHistMCIndexes(0)
 265 {
 266   // Constructor
 267
 268   // Define output slots only here
 269   // Output slot #1 writes into a TList container
 270   DefineOutput(1, TList::Class());
 271   DefineOutput(2, AliAnalysisCentralitySelector::Class());
 272   DefineOutput(3, AliESDtrackCuts::Class());
 273 }
 274 AliAnalysisTaskPerformanceStrange::~AliAnalysisTaskPerformanceStrange() {
 275   //
 276   // Destructor
 277   //
 278   if (fListHist && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())  { delete fListHist;     fListHist = 0x0;    }
 279   if (fCentrSelector && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())  { delete fCentrSelector;    fCentrSelector = 0x0;    }
 280   if (fTracksCuts && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())  { delete fTracksCuts;     fTracksCuts = 0x0;    }
 281
 282
 283 }
 284 //________________________________________________________________________
 285 void AliAnalysisTaskPerformanceStrange::UserCreateOutputObjects()
 286 {
 287
 288   //******************
 289   // Create histograms
 290   //*******************
 291   fListHist = new TList();
 292   fListHist->SetOwner();
 293   //fListHistCuts = new TList();
 294   //fListHistCuts->SetOwner();
 295
 296   // Bo: tbd: condition before allocation (i.e. if (!fHistMCMultiplicityPrimary){...} for each histo...
 297
 298   //***************
 299   // MC histograms
 300   //***************
 301
 302   // Primary Vertex:
 303   fHistMCPrimaryVertexX          = new TH1F("h1MCPrimaryVertexX", "MC Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 304   fListHist->Add(fHistMCPrimaryVertexX);
 305
 306   fHistMCPrimaryVertexY          = new TH1F("h1MCPrimaryVertexY", "MC Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 307   fListHist->Add(fHistMCPrimaryVertexY);
 308
 309   fHistMCPrimaryVertexZ          = new TH1F("h1MCPrimaryVertexZ", "MC Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 310   fListHist->Add(fHistMCPrimaryVertexZ);
 311
 312   // Pt track distribution
 313
 314   fHistPtTracks            = new TH1F("h1PtTracks", "Pt tracks;p_{t} (GeV/c);Counts", 240,0,12);
 315   fListHist->Add(fHistPtTracks);
 316
 317   // Multiplicity
 318   fHistMCMultiplicityPrimary           = new TH1F("h1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
 319   fListHist->Add(fHistMCMultiplicityPrimary);
 320
 321   fHistMCMultiplicityTracks            = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
 322   fListHist->Add(fHistMCMultiplicityTracks);
 323
 324   // Pt Distribution of non-primary particles:
 325   fHistMCPtAllK0s                      = new TH1F("h1MCPtAllK0s", "Non-primary MC K^{0};p_{t} (GeV/c);Counts",240,0,12);
 326   fListHist->Add(fHistMCPtAllK0s);
 327
 328   fHistMCPtAllLambda                   = new TH1F("h1MCPtAllLambda", "Non-primary MC #Lambda^{0};p_{t} (GeV/c);Counts",240,0,12);
 329   fListHist->Add(fHistMCPtAllLambda);
 330
 331   fHistMCPtAllAntiLambda               = new TH1F("h1MCPtAllAntiLambda", "Non-primary MC #bar{#Lambda}^{0};p_{t} (GeV/c);Counts",240,0,12);
 332   fListHist->Add(fHistMCPtAllAntiLambda);
 333
 334   fHistMCPtAllXi               = new TH1F("h1MCPtAllXi", "Non-primary MC #Xi;p_{t} (GeV/c);Counts",240,0,12);
 335   fListHist->Add(fHistMCPtAllXi);
 336
 337   fHistMCPtAllAntiXi               = new TH1F("h1MCPtAllAntiXi", "Non-primary MC #bar{#Xi};p_{t} (GeV/c);Counts",240,0,12);
 338   fListHist->Add(fHistMCPtAllAntiXi);
 339
 340   fHistMCPtAllOmega               = new TH1F("h1MCPtAllOmega", "Non-primary MC #Omega;p_{t} (GeV/c);Counts",240,0,12);
 341   fListHist->Add(fHistMCPtAllOmega);
 342
 343   fHistMCPtAllAntiOmega               = new TH1F("h1MCPtAllAntiOmega", "Non-primary MC #bar{#Omega};p_{t} (GeV/c);Counts",240,0,12);
 344   fListHist->Add(fHistMCPtAllAntiOmega);
 345
 346
 347   // Rapidity distribution:
 348   fHistMCRapK0s                 = new TH1F("h1MCRapK0s", "K^{0};y",160,-4,4);
 349   fListHist->Add(fHistMCRapK0s);
 350
 351   fHistMCRapLambda              = new TH1F("h1MCRapLambda", "#Lambda;y",160,-4,4);
 352   fListHist->Add(fHistMCRapLambda);
 353
 354
 355   fHistMCRapAntiLambda          = new TH1F("h1MCRapAntiLambda", "#bar{#Lambda};y",160,-4,4);
 356   fListHist->Add(fHistMCRapAntiLambda);
 357
 358
 359   // Pt distribution:
 360   fHistMCPtK0s               = new TH1F("h1MCPtK0s", "K^{0};p_{t} (GeV/c)",240,0,12);
 361   fListHist->Add(fHistMCPtK0s);
 362
 363   fHistMCPtLambda            = new TH1F("h1MCPtLambda", "#Lambda^{0};p_{t} (GeV/c)",240,0,12);
 364   fListHist->Add(fHistMCPtLambda);
 365
 366   fHistMCPtAntiLambda            = new TH1F("h1MCPtAntiLambda", "#AntiLambda^{0};p_{t} (GeV/c)",240,0,12);
 367   fListHist->Add(fHistMCPtAntiLambda);
 368
 369   //***********************************
 370   // Reconstructed particles histograms
 371   //***********************************
 372
 373   // Number of events;
 374   fHistNumberEvents           = new TH1F("h1NumberEvents", "Number of events; index;Number of Events",10,0,10);
 375   fListHist->Add(fHistNumberEvents);
 376
 377   // multiplicity
 378   fHistTrackPerEvent           = new TH1F("h1TrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",10000,0,10000);
 379   fListHist->Add(fHistTrackPerEvent);
 380
 381   fHistTPCMult           = new TH1F("h1HistTPCMult", "TPC tracks per event;Number of Tracks;Number of Events",10000,0,10000);
 382   fListHist->Add(fHistTPCMult);
 383
 384
 385   fHistTrackletPerEvent       = new TH1F("h1TrackletPerEvent", "Number of tracklets;Number of tracklets per events;Number of events",1000,0,1000);
 386   fListHist->Add(fHistTrackletPerEvent);
 387
 388   fHistTPCTracks               = new TH1F("h1TPCTracks","Distribution of TPC tracks;Number of TPC tracks:Number of events",1000,0,10000);
 389   fListHist->Add(fHistTPCTracks);
 390
 391   // Primary Vertex:
 392   fHistSPDPrimaryVertexZ          = new TH1F("h1SPDPrimaryVertexZ", "SPD Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 393   fListHist->Add(fHistSPDPrimaryVertexZ);
 394
 395   fHistPrimaryVertexX          = new TH1F("h1PrimaryVertexX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 396   fListHist->Add(fHistPrimaryVertexX);
 397
 398   fHistPrimaryVertexY          = new TH1F("h1PrimaryVertexY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 399   fListHist->Add(fHistPrimaryVertexY);
 400
 401   fHistPrimaryVertexZ          = new TH1F("h1PrimaryVertexZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-20,20);
 402   fListHist->Add(fHistPrimaryVertexZ);
 403
 404   //////K0s///////////////// 2D histos: cut vs on fly status////
 405
 406   // V0 Multiplicity
 407   if (!fHistV0Multiplicity) {
 408     if (fCollidingSystems)
 409       fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 200, 0, 40000);
 410     else
 411       fHistV0Multiplicity = new TH1F("fHistV0Multiplicity", "Multiplicity distribution;Number of Offline V0s;Events", 10, 0, 10);
 412     fListHist->Add(fHistV0Multiplicity);
 413   }
 414
 415
 416   // Mass:
 417   fHistMassK0                   = new TH1F("h1MassK0", "K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 200, 0.4, 0.6);
 418   fListHist->Add(fHistMassK0);
 419
 420   fHistMassLambda               = new TH1F("h1MassLambda", "#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});Counts", 150, 1.05, 1.2);
 421   fListHist->Add(fHistMassLambda);
 422
 423   fHistMassAntiLambda           = new TH1F("h1MassAntiLambda", "#bar{#Lambda}^{0} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 150, 1.05, 1.2);
 424   fListHist->Add(fHistMassAntiLambda);
 425
 426   fHistMassXi           = new TH1F("h1MassXi", "#Xi candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 150, 1.25, 1.4);  // 1321.71
 427   fListHist->Add(fHistMassXi);
 428
 429   fHistMassAntiXi           = new TH1F("h1MassAntiXi", "#bar{#Xi} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 150, 1.25, 1.4);
 430   fListHist->Add(fHistMassAntiXi);
 431
 432   fHistMassOmega           = new TH1F("h1MassOmega", "#Omega candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 150, 1.6, 1.75); //1672.45
 433   fListHist->Add(fHistMassOmega);
 434
 435   fHistMassAntiOmega           = new TH1F("h1MassAntiOmega", "#bar{#Omega} candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 150, 1.6, 1.75);
 436   fListHist->Add(fHistMassAntiOmega);
 437
 438   //PID vs Mass
 439   fHistMassXiVsPID           = new TH2F("h1MassXiVsPID", "#Xi candidates;M(#bar{p}#pi^{+}) (GeV/c^{2});PID", 150, 1.25, 1.4,5,0,5);  // 1321.71
 440   fListHist->Add(fHistMassXiVsPID);
 441
 442
 443
 444   // Pt Vs Mass
 445   fHistPtVsMassK0               = new TH2F("h2PtVsMassK0","K^{0} candidates;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",400, 0.4, 0.6,240,0,12);
 446   fListHist->Add(fHistPtVsMassK0);
 447
 448   fHistPtVsMassLambda           = new TH2F("h2PtVsMassLambda","#Lambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",280, 1.06, 1.2,240,0,12);
 449   fListHist->Add(fHistPtVsMassLambda);
 450
 451   fHistPtVsMassAntiLambda           = new TH2F("h2PtVsMassAntiLambda","#AntiLambda^{0} candidates;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",280, 1.06, 1.2,240,0,12);
 452   fListHist->Add(fHistPtVsMassAntiLambda);
 453
 454
 455   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 456
 457   ///Armenteros Podolansky
 458   fHistArmenterosPodolanski     = new TH2F("h2ArmenterosPodolanski","Armenteros-Podolanski phase space;#alpha;p_{t} arm",100,-1.0,1.0,50,0,0.5);
 459   fListHist->Add(fHistArmenterosPodolanski);
 460
 461   ///Inv. Mass K0s vs Inv. Mass. Lambda
 462   fHistK0sMassVsLambdaMass      = new TH2F("h2HistK0sMassVsLambdaMass","K^{0} vs #Lambda^{0} candidates; M(#pi^{+}#pi^{-}) (GeV/c^{2}); M(p#pi^{-}) (GeV/c^{2})",200, 0.4, 0.6,140, 1.06, 1.2);
 463   fListHist->Add(fHistK0sMassVsLambdaMass);
 464
 465   //dE/dx vs P daughters
 466   fHistTPCsigPLambda                            = new TH2F("h2TPCsignalVsPLambda","TPC signal Vs p_{t} daughters;  p (GeV/c);TPC signal",1000,0,2,1000,0,1000);
 467   fListHist->Add(fHistTPCsigPLambda);
 468
 469
 470   fHistTPCsigPAntiLambda                            = new TH2F("h2TPCsignalVsPAntiLambda","TPC signal Vs p_{t} daughters;  p (GeV/c);TPC signal",1000,0,2,1000,0,1000);
 471   fListHist->Add(fHistTPCsigPAntiLambda);
 472
 473
 474   fHistNSigmaProton                          =new TH1F("h1NSigmaProton","Number of sigmas for proton;;Count",600,0,6);
 475   fListHist->Add(fHistNSigmaProton);
 476
 477
 478   //********************************
 479   // Associated particles histograms
 480   //********************************
 481
 482   // Rap distribution
 483   fHistAsMcRapK0                = new TH1F("h1AsMcRapK0", "K^{0} associated;eta;Counts", 60, -1.5, 1.5);
 484   fListHist->Add(fHistAsMcRapK0);
 485
 486   fHistAsMcRapLambda            = new TH1F("h1AsMcRapLambda", "#Lambda^{0} associated;eta;Counts", 60, -1.5, 1.5);
 487   fListHist->Add(fHistAsMcRapLambda);
 488
 489   fHistAsMcRapAntiLambda        = new TH1F("h1AsMcRapAntiLambda", "#bar{#Lambda}^{0} associated;eta;Counts", 60, -1.5, 1.5);
 490   fListHist->Add(fHistAsMcRapAntiLambda);
 491
 492   //Pt distribution
 493   fHistAsMcPtK0                = new TH1F("h1AsMcPtK0", "K^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
 494   fListHist->Add(fHistAsMcPtK0);
 495
 496   fHistAsMcPtLambda            = new TH1F("h1AsMcPtLambda", "#Lambda^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
 497   fListHist->Add(fHistAsMcPtLambda);
 498
 499   fHistAsMcPtAntiLambda            = new TH1F("h1AsMcPtAntiLambda", "#AntiLambda^{0} associated;p_{t} (GeV/c);Counts", 240,0,12);
 500   fListHist->Add(fHistAsMcPtAntiLambda);
 501
 502   // Mass
 503   fHistPidMcMassK0             = new TH1F("h1PidMcMassK0", "K^{0} MC PId checked;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 504   fListHist->Add(fHistPidMcMassK0);
 505
 506   fHistPidMcMassLambda         = new TH1F("h1PidMcMassLambda", "#Lambda^{0} MC PId checked;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 507   fListHist->Add(fHistPidMcMassLambda);
 508
 509   fHistPidMcMassAntiLambda     = new TH1F("h1PidMcMassAntiLambda", "#bar{#Lambda}^{0} MC PId checked;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 510   fListHist->Add(fHistPidMcMassAntiLambda);
 511
 512   fHistAsMcMassK0              = new TH1F("h1AsMcMassK0", "K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});Counts", 100, 0.4, 0.6);
 513   fListHist->Add(fHistAsMcMassK0);
 514
 515   fHistAsMcMassLambda          = new TH1F("h1AsMcMassLambda", "#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 516   fListHist->Add(fHistAsMcMassLambda);
 517
 518   fHistAsMcMassAntiLambda      = new TH1F("h1AsMcMassAntiLambda", "#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});Counts", 75, 1.05, 1.2);
 519   fListHist->Add(fHistAsMcMassAntiLambda);
 520
 521   //Pt versus Mass
 522   fHistAsMcPtVsMassK0               = new TH2F("h2AsMcPtVsMassK0","K^{0} associated;M(#pi^{+}#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",200, 0.4, 0.6,240,0,12);
 523   fListHist->Add(fHistAsMcPtVsMassK0);
 524
 525   fHistAsMcPtVsMassLambda           = new TH2F("h2AsMcPtVsMassLambda","#Lambda^{0} associated;M(p#pi^{-}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
 526   fListHist->Add(fHistAsMcPtVsMassLambda);
 527
 528   fHistAsMcPtVsMassAntiLambda       = new TH2F("h2AsMcPtVsMassAntiLambda","#bar{#Lambda}^{0} associated;M(#bar{p}#pi^{+}) (GeV/c^{2});p_{t} (GeV/c)",140, 1.06, 1.2,240,0,12);
 529   fListHist->Add(fHistAsMcPtVsMassAntiLambda);
 530
 531   fHistCompositionXi = new TH1F("h1CompisitionXi","background composition of Xi;Part ID;PID",8000,-4000,4000);
 532   fListHist->Add(fHistCompositionXi);
 533
 534   fHistCompositionAntiXi = new TH1F("h1CompisitionAntiXi","background composition of AntiXi;Part ID;Counts",8000,-4000,4000);
 535   fListHist->Add(fHistCompositionAntiXi);
 536
 537   fHistCompositionOmega = new TH1F("h1CompisitionOmega","background composition of Omega;Part ID;Counts",8000,-4000,4000);
 538   fListHist->Add(fHistCompositionOmega);
 539
 540   fHistCompositionAntiOmega = new TH1F("h1CompisitionAntiOmega","background composition of AntiOmega;Part ID;Counts",8000,-4000,4000);
 541   fListHist->Add(fHistCompositionAntiOmega);
 542
 543   fHistMCIndexes = new TH1I("h1MCIndexes","MC Indexes;Index;Counts",210,-10,200);
 544   fListHist->Add(fHistMCIndexes);
 545
 546   PostData(1, fListHist);
 547   PostData(2, fCentrSelector);
 548   PostData(3, fTracksCuts);
 549 }
 550
 551 //________________________________________________________________________
 552 void AliAnalysisTaskPerformanceStrange::UserExec(Option_t *)
 553 {
 554   // Main loop
 555   // Called for each event
 556
 557   AliStack* stack = NULL;
 558   //  TClonesArray *mcArray = NULL;
 559   TArrayF mcPrimaryVtx;
 560
 561   fESD=(AliESDEvent *)InputEvent();
 562
 563   if (!fESD) {
 564     Printf("ERROR: fESD not available");
 565     return;
 566   }
 567
 568   AliVEvent* lEvent = InputEvent();
 569
 570   if (!lEvent) {
 571     Printf("ERROR: Event not available");
 572     return;
 573   }
 574
 575   ///////
 576   // PID
 577   ///////
 578
 579   if (fUsePID.Contains("withPID")) {
 580     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 581     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 582     fPIDResponse = inputHandler->GetPIDResponse();
 583   }
 584
 585   fHistNumberEvents->Fill(0.5);
 586   fHistNumberEvents->Fill(1.5);
 587   // Centrality selection
 588
 589   static AliESDtrackCuts * trackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
 590   Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,trackCuts);
 591   if(!isCentralitySelected) return;
 592
 593   // Remove Events with no tracks
 594   if (!(fESD->GetNumberOfTracks()))  return;
 595
 596   //*************************************
 597   // Cut used:
 598   //*************************************
 599
 600   // Cut Rapidity:
 601   Double_t lCutRap  = 0.5;
 602
 603   // cut Pseudorapidity:
 604
 605   Double_t lCutPseudorap = 0.8;
 606
 607   // If PID is used:
 608   Double_t lLimitPPID    = 0.7;
 609   Float_t cutNSigmaLowP  = 1E3;
 610   Float_t cutNSigmaHighP = 1E3;
 611   if (fUsePID.Contains("withPID")) {
 612     cutNSigmaLowP  = 4.0;
 613     cutNSigmaHighP = 4.0;
 614   }
 615   // Cut Daughters pt (GeV/c):
 616   //Double_t cutMinPtDaughter = 0.160;
 617
 618   // Cut primary vertex:
 619   Double_t cutPrimVertex = 10.0;
 620
 621   // cut ctau
 622   Double_t cutcTauL = 3*7.89;
 623   Double_t cutcTauK0s = 3*2.68;
 624
 625   // Min number of TPC clusters:
 626   // Int_t nbMinTPCclusters = 80;
 627   // PID flags:
 628   Int_t LambdaPID = 0;
 629   Int_t AntiLambdaPID = 0;
 630
 631   //////******************************************
 632   ////// Access MC: ******************************
 633   //////******************************************
 634
 635   if (fAnalysisMC) {
 636     if(fAnalysisType == "ESD") {
 637       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 638       if (!eventHandler) {
 639         Printf("ERROR: Could not retrieve MC event handler");
 640         return;
 641       }
 642       AliMCEvent* mcEvent = eventHandler->MCEvent();
 643       if (!mcEvent) {
 644         Printf("ERROR: Could not retrieve MC event");
 645         return;
 646       }
 647       stack = mcEvent->Stack();
 648       if (!stack) {
 649         Printf("ERROR: Could not retrieve stack");
 650         return;
 651       }
 652
 653       AliGenEventHeader* mcHeader=mcEvent->GenEventHeader();
 654       if(!mcHeader) return;
 655       mcHeader->PrimaryVertex(mcPrimaryVtx);
 656
 657       if (TMath::Abs(mcPrimaryVtx.At(2)) > cutPrimVertex) return;  /// cut on z of prim. vertex !!!!!!
 658     }
 659
 660     /*
 661     // PID parameters for MC simulations:
 662     fAlephParameters[0] = 2.15898e+00/50.;
 663     fAlephParameters[1] = 1.75295e+01;
 664     fAlephParameters[2] = 3.40030e-09;
 665     fAlephParameters[3] = 1.96178e+00;
 666     fAlephParameters[4] = 3.91720e+00;
 667     */
 668   }
 669
 670   //**********************************************
 671   // MC loop
 672   //**********************************************
 673
 674   //  Double_t lmcPrimVtxR      = 0;
 675
 676   Int_t lNbMCPrimary        = 0;
 677   Int_t lNbMCPart           = 0;
 678
 679   Int_t lPdgcodeCurrentPart = 0;
 680   Double_t lRapCurrentPart  = 0;
 681   Double_t lPtCurrentPart   = 0;
 682
 683   Int_t lComeFromSigma      = 0;
 684
 685   // Production Radius
 686   Double_t mcPosX     = 0.0,  mcPosY      = 0.0,  mcPosZ      = 0.0;
 687   Double_t mcPosR     = 0.0;
 688
 689   // Decay Radius
 690   Double_t mcDecayPosX = 0, mcDecayPosY = 0, mcDecayPosR = 0;
 691
 692   // current mc particle 's mother
 693   Int_t iCurrentMother  = 0, lPdgCurrentMother    = 0;
 694   //  Bool_t lCurrentMotherIsPrimary;
 695
 696   // variables for multiple reconstruction studies:
 697   Int_t id0           = 0, id1          = 0;
 698   Int_t lNtimesReconstructedK0s   = 0, lNtimesReconstructedLambda   = 0, lNtimesReconstructedAntiLambda   = 0;
 699   // Int_t lNtimesReconstructedK0sMI = 0, lNtimesReconstructedLambdaMI = 0, lNtimesReconstructedAntiLambdaMI = 0;
 700
 701   ////********************************
 702   ////Start loop over MC particles****
 703   ////********************************
 704
 705   if (fAnalysisMC) {
 706
 707     // Primary vertex
 708     fHistMCPrimaryVertexX->Fill(mcPrimaryVtx.At(0));
 709     fHistMCPrimaryVertexY->Fill(mcPrimaryVtx.At(1));
 710     fHistMCPrimaryVertexZ->Fill(mcPrimaryVtx.At(2));
 711
 712     //    lmcPrimVtxR = TMath::Sqrt(mcPrimaryVtx.At(0)*mcPrimaryVtx.At(0)+mcPrimaryVtx.At(1)*mcPrimaryVtx.At(1));
 713
 714     if(fAnalysisType == "ESD") {
 715
 716       lNbMCPrimary = stack->GetNprimary();
 717      lNbMCPart    = stack->GetNtrack();
 718
 719       fHistMCMultiplicityPrimary->Fill(lNbMCPrimary);
 720       fHistMCMultiplicityTracks->Fill(lNbMCPart);
 721
 722       for (Int_t iMc = 0; iMc < (stack->GetNtrack()); iMc++) {
 723         TParticle *p0 = stack->Particle(iMc);
 724         if (!p0) {
 725           //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
 726           continue;
 727         }
 728         lPdgcodeCurrentPart = p0->GetPdgCode();
 729
 730         // Keep only K0s, Lambda and AntiLambda, Xi-, antiXi+ and Omega-, antiOmega+:
 731         if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) && (lPdgcodeCurrentPart != 3312 ) && (lPdgcodeCurrentPart != -3312) && (lPdgcodeCurrentPart != 3334)  && (lPdgcodeCurrentPart != -3334) ) continue;
 732
 733         lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
 734         //lEtaCurrentPart   = p0->Eta();
 735         lPtCurrentPart    = p0->Pt();
 736         iCurrentMother    = p0->GetFirstMother();
 737
 738         //      lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();
 739         if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}
 740
 741         mcPosX = p0->Vx();
 742         mcPosY = p0->Vy();
 743         mcPosZ = p0->Vz();
 744         mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
 745
 746         id0  = p0->GetDaughter(0);
 747         id1  = p0->GetDaughter(1);
 748
 749         // Decay Radius and Production Radius
 750         if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {
 751           TParticle *pDaughter0 = stack->Particle(id0);
 752           //      TParticle *pDaughter1 = stack->Particle(id1);
 753           //      lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();
 754           //      lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();
 755
 756           mcDecayPosX = pDaughter0->Vx();
 757           mcDecayPosY = pDaughter0->Vy();
 758           mcDecayPosR = TMath::Sqrt(mcDecayPosX*mcDecayPosX+mcDecayPosY*mcDecayPosY);
 759         }
 760         else  {
 761           //Printf("ERROR: particle with label %d and/or %d not found in stack (mc loop)", id0,id1);
 762           mcDecayPosR = -1.0;
 763         }
 764
 765         if (lPdgcodeCurrentPart==310)   {
 766           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllK0s->Fill(lPtCurrentPart);
 767         }
 768         else if (lPdgcodeCurrentPart==3122)  {
 769           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllLambda->Fill(lPtCurrentPart);
 770         }
 771         else if (lPdgcodeCurrentPart==-3122) {
 772           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiLambda->Fill(lPtCurrentPart);
 773         }
 774         else if (lPdgcodeCurrentPart==3312) { // Xi
 775           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllXi->Fill(lPtCurrentPart);
 776         }
 777         else if (lPdgcodeCurrentPart==-3312) { //Anti Xi
 778           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiXi->Fill(lPtCurrentPart);
 779         }
 780         else if (lPdgcodeCurrentPart==3334) {  //Omega
 781           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllXi->Fill(lPtCurrentPart);
 782         }
 783         else if (lPdgcodeCurrentPart==-3334) { //Anti Omega
 784           if (TMath::Abs(lRapCurrentPart) < lCutRap) fHistMCPtAllAntiXi->Fill(lPtCurrentPart);
 785         }
 786
 787         if ( ( ( TMath::Abs(lPdgCurrentMother) == 3212)  ||
 788                ( TMath::Abs(lPdgCurrentMother) == 3224)  ||
 789                ( TMath::Abs(lPdgCurrentMother) == 3214)  ||
 790                ( TMath::Abs(lPdgCurrentMother) == 3114) )
 791              && ( iCurrentMother <= lNbMCPrimary )
 792              ) lComeFromSigma = 1;
 793         else lComeFromSigma = 0;
 794
 795         //*********************************************
 796         // Now keep only primary particles
 797         // if ( ( iMc > lNbMCPrimary ) && (!lComeFromSigma) ) continue;
 798
 799         //*************************************//
 800         // new definition of primary particles //
 801         //*************************************//
 802
 803         Double_t dx = 0;
 804         Double_t dy = 0;
 805         Double_t dz = 0;
 806         Double_t ProdDistance = 0;
 807
 808         dx = ( (mcPrimaryVtx.At(0)) - (p0->Vx()) );
 809         dy = ( (mcPrimaryVtx.At(1)) - (p0->Vy()) );
 810         dz = ( (mcPrimaryVtx.At(2)) - (p0->Vz()) );
 811
 812         ProdDistance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 813
 814         if (ProdDistance > 0.001) continue; // secondary V0
 815
 816         lNtimesReconstructedK0s   = 0; lNtimesReconstructedLambda   = 0; lNtimesReconstructedAntiLambda   = 0;
 817
 818         // Rap distribution
 819         if (lPdgcodeCurrentPart==310) {
 820           fHistMCRapK0s->Fill(lRapCurrentPart);
 821         }
 822
 823         if (lPdgcodeCurrentPart==3122) {
 824           fHistMCRapLambda->Fill(lRapCurrentPart);
 825         }
 826
 827         if (lPdgcodeCurrentPart==-3122) {
 828           fHistMCRapAntiLambda->Fill(lRapCurrentPart);
 829         }
 830         /*
 831         if (lPdgcodeCurrentPart==3312) {
 832           fHistMCRapXi->Fill(lRapCurrentPart);
 833         }
 834         */
 835         // Rapidity Cut
 836         if (TMath::Abs(lRapCurrentPart) > lCutRap) continue;
 837
 838         if (lPdgcodeCurrentPart==310) {
 839           fHistMCPtK0s->Fill(lPtCurrentPart);
 840         }
 841         else
 842           if (lPdgcodeCurrentPart==3122) {
 843             fHistMCPtLambda->Fill(lPtCurrentPart);
 844           }
 845           else
 846             if (lPdgcodeCurrentPart==-3122) {
 847               fHistMCPtAntiLambda->Fill(lPtCurrentPart);
 848             }
 849
 850       } // end loop ESD MC
 851     } // end ESD condition
 852   } // End Loop over MC condition
 853
 854   //************************************
 855   // ESD loop
 856   //************************************
 857
 858   //  Double_t  lLambdaMass = 1.115683;  //PDG
 859   Double_t lPLambda = 0;
 860   // Double_t lPAntiLambda = 0;
 861   //Double_t lPK0s = 0;
 862   Double_t lMagneticField = 999;
 863
 864   //Multiplcity:
 865   Int_t    nv0sTot= 0, nv0s = 0;
 866   Int_t    nCasTot= 0;
 867
 868   // Qualities
 869   Double_t v0q = 0, kinCasQual = 0;
 870   //  Int_t nv0sMI =0;
 871   // Variables:
 872   Double_t  lV0Position[3];
 873
 874   Double_t lDcaPosToPrimVertex = 0;
 875   Double_t lDcaNegToPrimVertex = 0;
 876   Double_t lDcaV0Daughters     = 0;
 877   Double_t lV0cosPointAngle    = 0;
 878   Double_t lChi2V0             = 0;
 879   Double_t lV0DecayLength      = 0;
 880   Double_t lV0tDecayLength     = 0; //transverse decay length in xy plain
 881   Double_t lV0Radius           = 0;
 882   //  Double_t lDcaV0ToPrimVertex  = 0;
 883   Double_t lcTauLambda         = 0;
 884   Double_t lcTauAntiLambda     = 0;
 885   Double_t lcTauK0s            = 0;
 886   Int_t    lOnFlyStatus        = 0;
 887
 888
 889   Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
 890   Double_t lInvMassXi = 0,lInvMassAntiXi = 0,lInvMassOmega = 0,lInvMassAntiOmega = 0;
 891   Double_t lPtK0s      = 0, lPtLambda      = 0, lPtAntiLambda      = 0;
 892   Double_t lRapK0s     = 0, lRapLambda     = 0, lRapAntiLambda     = 0;
 893   //  Double_t lEtaK0s     = 0, lEtaLambda     = 0, lEtaAntiLambda     = 0;
 894   Double_t lAlphaV0      = 0, lPtArmV0       = 0;
 895
 896   // to study Associated V0s:
 897   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg         = 0;
 898   UInt_t   lLabelTrackPos       = 0, lLabelTrackNeg         = 0;
 899   Int_t    lCheckPIdK0Short     = 0, lCheckMcK0Short        = 0;
 900   Int_t    lCheckPIdLambda      = 0, lCheckMcLambda         = 0;
 901   Int_t    lCheckPIdAntiLambda  = 0, lCheckMcAntiLambda     = 0;
 902   Int_t    lCheckSecondaryK0s   = 0, lCheckSecondaryLambda  = 0, lCheckSecondaryAntiLambda  = 0;
 903   //  Int_t    lCheckGamma          = 0;
 904   Double_t mcPosMotherX         = 0, mcPosMotherY           = 0, mcPosMotherZ  = 0;
 905   Double_t mcPosMotherR         = 0;
 906   //  Double_t mcMotherPt           = 0;
 907
 908   Int_t lIndexPosMother        = 0;
 909   Int_t lIndexNegMother        = 0;
 910   Int_t lIndexMotherOfMother   = 0;
 911   Int_t lPDGCodePosDaughter    = 0;
 912   Int_t lPDGCodeNegDaughter    = 0;
 913   Int_t lPdgcodeMother         = 0;
 914   Int_t lPdgcodeMotherOfMother = 0;
 915
 916   // Reconstructed position
 917   Double_t rcPosXK0s        = 0,  rcPosYK0s        = 0, rcPosZK0s        = 0;
 918   Double_t rcPosRK0s        = 0;
 919   Double_t rcPosXLambda     = 0,  rcPosYLambda     = 0, rcPosZLambda     = 0;
 920   Double_t rcPosRLambda     = 0;
 921   Double_t rcPosXAntiLambda = 0,  rcPosYAntiLambda = 0, rcPosZAntiLambda = 0;
 922   Double_t rcPosRAntiLambda = 0;
 923
 924   // Pt resolution
 925   Double_t deltaPtK0s  = 0, deltaPtLambda  = 0, deltaPtAntiLambda  = 0;
 926
 927   // Daughters
 928   AliESDtrack  *myTrackPos  = NULL;
 929   AliESDtrack  *myTrackNeg  = NULL;
 930
 931   // Daughters' momentum:
 932   Double_t  lMomPos[3] = {999,999,999};
 933   Double_t  lMomNeg[3] = {999,999,999};
 934
 935   // Inner Wall parameters:
 936   Double_t  lMomInnerWallPos =999, lMomInnerWallNeg = 999;
 937
 938   // PID
 939   Float_t nSigmaPosPion   = 0;
 940   Float_t nSigmaNegPion   = 0;
 941
 942   Float_t nSigmaPosProton = 0;
 943   Float_t nSigmaNegProton = 0;
 944
 945 // Bachelor
 946 Bool_t lIsBachelorKaonForTPC = kFALSE;
 947 Bool_t lIsBachelorPionForTPC = kFALSE;
 948
 949 // Negative V0 daughter
 950 Bool_t lIsNegPionForTPC   = kFALSE;
 951 Bool_t lIsNegProtonForTPC = kFALSE;
 952
 953 // Positive V0 daughter
 954 Bool_t lIsPosPionForTPC   = kFALSE;
 955 Bool_t lIsPosProtonForTPC = kFALSE;
 956
 957  Double_t lEffMassXi      = 0. ;
 958  Double_t lChi2Xi         = -1. ;
 959  Double_t lDcaXiDaughters = -1. ;
 960  Double_t lXiCosineOfPointingAngle = -1. ;
 961  Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
 962  Double_t lXiRadius2D = -1000. ;
 963  Double_t lXiRadius3D = -1000. ;
 964
 965  Int_t lIDs = 0;
 966  Int_t lIndexBachelorMother = 0;
 967  Int_t lIndexBachelor = 0;
 968   //  Int_t lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
 969   Int_t lCheckPIDLambdaPosDaughter     = 0;// lCheckPIDLambdaNegDaughter     = 0;
 970   //  Int_t lCheckPIDAntiLambdaPosDaughter = 0;
 971   Int_t lCheckPIDAntiLambdaNegDaughter = 0;
 972
 973   //****************************************************
 974   // Primary Vertex cuts &
 975   // Magnetic field and Quality tracks cuts
 976   //****************************************************
 977   Double_t  lPrimaryVtxPosition[3];
 978   Double_t  lPrimaryVtxCov[6];
 979   Double_t  lPrimaryVtxChi2 = 999;
 980   Double_t  lResPrimaryVtxX = 999;
 981   Double_t  lResPrimaryVtxY = 999;
 982   Double_t  lResPrimaryVtxZ = 999;
 983
 984   AliAODVertex *myPrimaryVertex = NULL;
 985   //const AliVVertex *mySPDPrimaryVertex = NULL;
 986
 987   const AliMultiplicity *myMultiplicty = ((AliESDEvent*)fESD)->GetMultiplicity();
 988
 989   if(fAnalysisType == "ESD") {
 990     ////////////////////////////////////////////////////////////////////////////////////
 991     //////   Best Primary Vertex:
 992     if(fCollidingSystems ==0 || fCollidingSystems == 1){ //pp, PbPb Analysis
 993     const AliESDVertex *myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
 994     myBestPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertex();
 995     if (!myBestPrimaryVertex) return;
 996     if (!myBestPrimaryVertex->GetStatus()) return;
 997
 998     fHistNumberEvents->Fill(3.5);
 999
1000     myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1001     myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1002     if ( ( TMath::Abs(lPrimaryVtxPosition[2]) ) > cutPrimVertex) return ; //// cut on z of prim. vertex!!!!!
1003     fHistNumberEvents->Fill(4.5);
1004     lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1005     lResPrimaryVtxX = myBestPrimaryVertex->GetXRes();
1006     lResPrimaryVtxY = myBestPrimaryVertex->GetYRes();
1007     lResPrimaryVtxZ = myBestPrimaryVertex->GetZRes();
1008     // remove TPC-only primary vertex : retain only events with tracking + SPD vertex
1009     const AliESDVertex *mySPDPrimaryVertex = ((AliESDEvent*)fESD)->GetPrimaryVertexSPD();
1010     if (!mySPDPrimaryVertex) return;
1011     fHistSPDPrimaryVertexZ->Fill(mySPDPrimaryVertex->GetZ());
1012     const AliESDVertex *myPrimaryVertexTracking = ((AliESDEvent*)fESD)->GetPrimaryVertexTracks();
1013     if (!myPrimaryVertexTracking) return;
1014     if (!mySPDPrimaryVertex->GetStatus() && !myPrimaryVertexTracking->GetStatus() ) return;
1015     fHistNumberEvents->Fill(5.5);
1016     fHistTrackPerEvent->Fill(fESD->GetNumberOfTracks());
1017     myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1018     if (!myPrimaryVertex) return;
1019   }
1020     ///  pPb analysis
1021     //    if(fCollidingSystems == 2){  //twiky https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PAVertexSelectionStudies
1022     if(1){  //twiky https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PAVertexSelectionStudies
1023       const AliESDVertex* trkVtx = ((AliESDEvent*)fESD)->GetPrimaryVertex();
1024       if (!trkVtx || trkVtx->GetNContributors()<=0) return;
1025       TString vtxTtl = trkVtx->GetTitle();
1026       if (!vtxTtl.Contains("VertexerTracks")) return;
1027       Float_t zvtx = trkVtx->GetZ();
1028       const AliESDVertex* spdVtx = ((AliESDEvent*)fESD)->GetPrimaryVertexSPD();
1029       if (spdVtx->GetNContributors()<=0) return;
1030       TString vtxTyp = spdVtx->GetTitle();
1031       Double_t cov[6]={0};
1032       spdVtx->GetCovarianceMatrix(cov);
1033       Double_t zRes = TMath::Sqrt(cov[5]);
1034       if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1035       if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1036       if (TMath::Abs(zvtx) > cutPrimVertex) return;
1037       trkVtx->GetXYZ(lPrimaryVtxPosition);
1038     }
1039
1040     // Number of Tracklets:
1041     fHistTrackletPerEvent->Fill(myMultiplicty->GetNumberOfTracklets());
1042     lMagneticField = ((AliESDEvent*)fESD)->GetMagneticField();
1043     fHistTPCTracks->Fill(AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent*)fESD, kTRUE));
1044     //////simple chack for multiplicity////////////////////////////////////////////////////////
1045     Int_t i =0;
1046
1047     for (Int_t jTracks=0;jTracks<fESD->GetNumberOfTracks();jTracks++){
1048
1049       AliESDtrack* tPCtrack=fESD->GetTrack(jTracks);
1050       Float_t xy=0;
1051       Float_t z=0;
1052       tPCtrack->GetImpactParameters(xy,z);
1053       if ((fTracksCuts->IsSelected(tPCtrack))&&(xy<1.0)&&(z<1.0)) {i=i+1;}
1054     }
1055     fHistTPCMult->Fill(i);
1056     /////////////////////////////////////////////////////////////////////////////////////////
1057       }
1058
1059   fHistPrimaryVertexX->Fill(lPrimaryVtxPosition[0]);
1060   fHistPrimaryVertexY->Fill(lPrimaryVtxPosition[1]);
1061   fHistPrimaryVertexZ->Fill(lPrimaryVtxPosition[2]);
1062   //Double_t lrcPrimVtxR = TMath::Sqrt(lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]+lPrimaryVtxPosition[0]*lPrimaryVtxPosition[0]);
1063
1064   ////*************************
1065   //// Cascade loop ***********
1066   ////*************************
1067   nCasTot = fESD->GetNumberOfCascades();
1068
1069   for (Int_t iCas = 0; iCas < nCasTot; iCas++) {
1070
1071     AliESDcascade *Cas = ((AliESDEvent*)fESD)->GetCascade(iCas);
1072       if (!Cas) continue;
1073
1074       lIndexBachelor = Cas->GetBindex();
1075       fHistMCIndexes->Fill(lIndexBachelor);
1076
1077
1078       UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( Cas->GetPindex() );
1079       UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( Cas->GetNindex() );
1080       UInt_t lBachIdx         = (UInt_t) TMath::Abs( lIndexBachelor );
1081       // Care track label can be negative in MC production (linked with the track quality)
1082       // However = normally, not the case for track index ...
1083
1084       // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
1085       if(lBachIdx == lIdxNegXi) {
1086         AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
1087       }
1088       if(lBachIdx == lIdxPosXi) {
1089         AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
1090       }
1091
1092       lEffMassXi                      = Cas->GetEffMassXi();
1093       lChi2Xi                         = Cas->GetChi2Xi();
1094       lDcaXiDaughters                 = Cas->GetDcaXiDaughters();
1095       lXiCosineOfPointingAngle        = Cas->GetCascadeCosineOfPointingAngle( lPrimaryVtxPosition[0],
1096                                                                              lPrimaryVtxPosition[1],
1097                                                                              lPrimaryVtxPosition[2] );
1098       // Take care : the best available vertex should be used (like in AliCascadeVertexer)
1099
1100       Cas->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
1101       lXiRadius2D    = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
1102       lXiRadius3D    = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] +  lPosXi[2]*lPosXi[2]);
1103
1104
1105       AliESDtrack *pTrackXi           = fESD->GetTrack( lIdxPosXi );
1106       AliESDtrack *nTrackXi           = fESD->GetTrack( lIdxNegXi );
1107       AliESDtrack *bachTrackXi        = fESD->GetTrack( lBachIdx );
1108       if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1109         AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
1110         continue;
1111       }
1112
1113       // Bachelor
1114       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
1115       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
1116
1117       // Negative V0 daughter
1118       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion   )) < 4) lIsNegPionForTPC   = kTRUE;
1119       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
1120
1121       // Positive V0 daughter
1122       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion   )) < 4) lIsPosPionForTPC   = kTRUE;
1123       if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
1124
1125
1126       if( bachTrackXi->Charge() < 0 ) {
1127         v0q = 0;
1128       kinCasQual =  Cas->ChangeMassHypothesis(v0q,3312);
1129       lInvMassXi = Cas->GetEffMassXi();
1130         v0q = 0;
1131       kinCasQual =  Cas->ChangeMassHypothesis(v0q,3334);
1132       lInvMassOmega = Cas->GetEffMassXi();
1133       }
1134
1135       if( bachTrackXi->Charge() > 0 ) {
1136         v0q = 0;
1137       kinCasQual =  Cas->ChangeMassHypothesis(v0q,-3312);
1138       lInvMassAntiXi = Cas->GetEffMassXi();
1139         v0q = 0;
1140       kinCasQual =  Cas->ChangeMassHypothesis(v0q,-3334);
1141       lInvMassAntiOmega = Cas->GetEffMassXi();
1142       }
1143
1144
1145       if( lIsBachelorPionForTPC == kTRUE && lIsNegPionForTPC == kTRUE && lIsPosProtonForTPC == kTRUE){
1146       fHistMassXi->Fill(lInvMassXi);
1147       if(fAnalysisMC && (lInvMassXi > 1.25) && (lInvMassXi < 1.4)){
1148         TParticle *Bachelor = stack->Particle(lIndexBachelor);
1149                 //        lIndexBachelorMother = Bachelor->GetMother(1);
1150         lIndexBachelorMother = Bachelor->GetFirstMother();
1151         if(lIndexBachelorMother == -1) continue;
1152         if(lIndexBachelorMother == 1){  fHistMassXiVsPID->Fill(lInvMassXi,0);fHistCompositionXi->Fill(0); continue;} //combinatorial
1153         lIDs = stack->Particle(lIndexBachelorMother)->GetPdgCode();
1154         if(lIDs == 111) { fHistMassXiVsPID->Fill(lInvMassXi,1); } //pi0
1155         if(lIDs == 213) { fHistMassXiVsPID->Fill(lInvMassXi,2); } //rho+
1156         if(lIDs == 2212) { fHistMassXiVsPID->Fill(lInvMassXi,3); } //Proton
1157         if(lIDs == 2112) { fHistMassXiVsPID->Fill(lInvMassXi,4); } //Neutron
1158         fHistCompositionXi->Fill(lIDs);
1159        }
1160
1161       }
1162       if( lIsBachelorKaonForTPC == kTRUE && lIsNegPionForTPC == kTRUE && lIsPosProtonForTPC == kTRUE){
1163       fHistMassOmega->Fill(lInvMassOmega);
1164       if(fAnalysisMC && lInvMassOmega > 1.6 && lInvMassOmega < 1.75){
1165         TParticle *Bachelor = stack->Particle(lIndexBachelor);
1166                 //        lIndexBachelorMother = Bachelor->GetMother(1);
1167         lIndexBachelorMother = Bachelor->GetFirstMother();
1168         if(lIndexBachelorMother == -1) continue;
1169         if(lIndexBachelorMother == 1) { fHistCompositionOmega->Fill(0); continue;}
1170         lIDs = stack->Particle(lIndexBachelorMother)->GetPdgCode();
1171       fHistCompositionOmega->Fill(lIDs);
1172        }
1173
1174       }
1175       if( lIsBachelorPionForTPC == kTRUE && lIsPosPionForTPC == kTRUE && lIsNegProtonForTPC == kTRUE){
1176       fHistMassAntiXi->Fill(lInvMassAntiXi);
1177       if(fAnalysisMC && lInvMassAntiXi > 1.25 && lInvMassAntiXi < 1.4){
1178         TParticle *Bachelor = stack->Particle(lIndexBachelor);
1179                 //        lIndexBachelorMother = Bachelor->GetMother(1);
1180         lIndexBachelorMother = Bachelor->GetFirstMother();
1181         if(lIndexBachelorMother == -1) continue;
1182         if(lIndexBachelorMother == 1) { fHistCompositionAntiXi->Fill(0); continue;}
1183         lIDs = stack->Particle(lIndexBachelorMother)->GetPdgCode();
1184       fHistCompositionAntiXi->Fill(lIDs);
1185        }
1186
1187       }
1188       if( lIsBachelorKaonForTPC == kTRUE && lIsPosPionForTPC == kTRUE && lIsNegProtonForTPC == kTRUE){
1189       fHistMassAntiOmega->Fill(lInvMassAntiOmega);
1190       if(fAnalysisMC && lInvMassAntiOmega > 1.6 && lInvMassAntiOmega < 1.75){
1191         TParticle *Bachelor = stack->Particle(lIndexBachelor);
1192                 //        lIndexBachelorMother = Bachelor->GetMother(1);
1193         lIndexBachelorMother = Bachelor->GetFirstMother();
1194         if(lIndexBachelorMother == -1) continue;
1195         if(lIndexBachelorMother == 1) { fHistCompositionAntiOmega->Fill(0); continue;}
1196         lIDs = stack->Particle(lIndexBachelorMother)->GetPdgCode();
1197       fHistCompositionAntiOmega->Fill(lIDs);
1198        }
1199
1200       }
1201
1202
1203
1204
1205   } /// end of Cascade loop
1206
1207
1208   ////*************************
1209   //// V0 loop ****************
1210   ////*************************
1211
1212   nv0sTot = fESD->GetNumberOfV0s();
1213   if (!nv0sTot) fHistNumberEvents->Fill(6.5);
1214
1215   for (Int_t iV0 = 0; iV0 < nv0sTot; iV0++) {
1216
1217     lIndexPosMother     = 0; lIndexNegMother     = 0; lIndexMotherOfMother       = 0;
1218     lCheckPIdK0Short    = 0; lCheckMcK0Short     = 0; lCheckSecondaryK0s         = 0;
1219     lCheckPIdLambda     = 0; lCheckMcLambda      = 0; lCheckSecondaryLambda      = 0;
1220     lCheckPIdAntiLambda = 0; lCheckMcAntiLambda  = 0; lCheckSecondaryAntiLambda  = 0;
1221     lComeFromSigma      = 0; //lCheckGamma = 0;
1222
1223     if(fAnalysisType == "ESD") {
1224
1225       AliESDv0 *v0 = ((AliESDEvent*)fESD)->GetV0(iV0);
1226       if (!v0) continue;
1227       //      if ((v0->Pt())<0.6) continue;
1228
1229       // V0's Daughters
1230       lIndexTrackPos = TMath::Abs(v0->GetPindex());
1231       lIndexTrackNeg = TMath::Abs(v0->GetNindex());
1232       AliESDtrack *myTrackPosTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1233       AliESDtrack *myTrackNegTest = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1234       if (!myTrackPosTest || !myTrackNegTest) {
1235         Printf("strange analysis::UserExec:: Error:Could not retreive one of the daughter track\n");
1236         continue;
1237       }
1238       // Remove like-sign
1239       if ((Int_t)myTrackPosTest->GetSign() == (Int_t)myTrackNegTest->GetSign()){
1240         continue;
1241       }
1242
1243       // VO's main characteristics to check the reconstruction cuts
1244       lOnFlyStatus       = v0->GetOnFlyStatus();
1245       lChi2V0            = v0->GetChi2V0();
1246       lDcaV0Daughters    = v0->GetDcaV0Daughters();
1247       //      lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1248       lV0cosPointAngle   = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1], lPrimaryVtxPosition[2]);
1249
1250       v0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1251
1252       lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1253       lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1254                                    TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1255                                    TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1256       lV0tDecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1257                                     TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2));
1258
1259       if( myTrackPosTest->GetSign() ==1){
1260         myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1261         myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1262         // Daughters' momentum;
1263         v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1264         v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1265       }
1266
1267       if( myTrackPosTest->GetSign() ==-1){
1268         myTrackPos = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackNeg);
1269         myTrackNeg = ((AliESDEvent*)fESD)->GetTrack(lIndexTrackPos);
1270         // Daughters' momentum;
1271         v0->GetPPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1272         v0->GetNPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1273       }
1274
1275       lLabelTrackPos = (UInt_t)TMath::Abs(myTrackPos->GetLabel());
1276       lLabelTrackNeg = (UInt_t)TMath::Abs(myTrackNeg->GetLabel());
1277
1278       // Inner Wall parameter:
1279       const AliExternalTrackParam *myInnerWallTrackPos = myTrackPos->GetInnerParam();
1280
1281       if (myInnerWallTrackPos) {lMomInnerWallPos = myInnerWallTrackPos->GetP();}
1282       else {continue;}
1283
1284       const AliExternalTrackParam *myInnerWallTrackNeg = myTrackNeg->GetInnerParam();
1285
1286       if (myInnerWallTrackNeg) {lMomInnerWallNeg = myInnerWallTrackNeg->GetP();}
1287       else {continue;}
1288
1289       // DCA between daughter and Primary Vertex:
1290       if (myTrackPos) lDcaPosToPrimVertex = TMath::Abs(myTrackPos->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
1291
1292       if (myTrackNeg) lDcaNegToPrimVertex = TMath::Abs(myTrackNeg->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lMagneticField) );
1293
1294       // Quality tracks cuts:
1295       if ( !(fTracksCuts->IsSelected(myTrackPos)) || !(fTracksCuts->IsSelected(myTrackNeg)) )
1296         {         continue;}
1297
1298       if ( ( TMath::Abs(myTrackPos->Eta()) > lCutPseudorap ) || ( TMath::Abs(myTrackNeg->Eta()) > lCutPseudorap ) ) {continue;}
1299
1300       // Armenteros variables:
1301       lAlphaV0      =  v0->AlphaV0();
1302       lPtArmV0      =  v0->PtArmV0();
1303
1304       // Pseudorapidity:
1305       //      lV0Eta = v0->Eta();
1306       //////////////////////////////////////////////////////////////////////////
1307       // Invariant mass
1308       v0->ChangeMassHypothesis(310);
1309       lInvMassK0s = v0->GetEffMass();
1310       lPtK0s = v0->Pt();
1311       // lPzK0s = v0->Pz();
1312
1313       v0->ChangeMassHypothesis(3122);
1314       lInvMassLambda = v0->GetEffMass();
1315       lPtLambda = v0->Pt();
1316       //lPzLambda = v0->Pz();
1317
1318       v0->ChangeMassHypothesis(-3122);
1319       lInvMassAntiLambda = v0->GetEffMass();
1320       lPtAntiLambda = v0->Pt();
1321       //lPzAntiLambda = v0->Pz();
1322
1323       // Rapidity:
1324       lRapK0s    = v0->Y(310);
1325       lRapLambda = v0->Y(3122);
1326       lRapAntiLambda = v0->Y(-3122);
1327
1328       if (lPtK0s==0)    {continue;}
1329       if (lPtLambda==0)         {continue;}
1330
1331       if (lPtAntiLambda==0)     {continue;}
1332
1333       ///////////////////////////////////////////////////////////////////////
1334
1335       // PID  new method July 2011
1336       if (fUsePID.Contains("withPID")) {
1337         //      nSigmaPosPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kPion));
1338         nSigmaPosPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kPion));
1339         //      nSigmaNegPion   = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kPion));
1340         nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(myTrackNeg, AliPID::kPion));
1341         //      nSigmaPosProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackPos,AliPID::kProton));
1342         nSigmaPosProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(myTrackPos, AliPID::kProton));
1343         //      nSigmaNegProton = TMath::Abs(fESDpid->NumberOfSigmasTPC(myTrackNeg,AliPID::kProton));
1344         nSigmaNegProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(myTrackNeg, AliPID::kProton));
1345
1346       }
1347       else {
1348         nSigmaPosPion = 0; nSigmaNegPion =0; nSigmaPosProton = 0; nSigmaNegProton= 0;
1349       }
1350
1351       // Monte-Carlo particle associated to reconstructed particles:
1352       if (fAnalysisMC) {
1353         //if (lLabelTrackPos < 0 || lLabelTrackNeg < 0) continue;
1354         TParticle  *lMCESDPartPos  = stack->Particle(lLabelTrackPos);
1355         if(!lMCESDPartPos) {
1356           //  Printf("no MC particle for positive and/or negative daughter\n");
1357           continue;
1358         }
1359         TParticle  *lMCESDPartNeg  = stack->Particle(lLabelTrackNeg);
1360         if (!lMCESDPartNeg)     {  continue;}
1361         lPDGCodePosDaughter = lMCESDPartPos->GetPdgCode();
1362         lPDGCodeNegDaughter = lMCESDPartNeg->GetPdgCode();
1363         lIndexPosMother = lMCESDPartPos->GetFirstMother();
1364         lIndexNegMother = lMCESDPartNeg->GetFirstMother();
1365
1366         if (lIndexPosMother == -1) {
1367           lPdgcodeMother = 0;
1368           lIndexMotherOfMother = 0;
1369           mcPosX = 0;
1370           mcPosY = 0;
1371           mcPosZ = 0;
1372           mcPosR = 0;
1373           mcPosMotherX = 0;
1374           mcPosMotherY = 0;
1375           mcPosMotherZ = 0;
1376           mcPosMotherR = 0;
1377           // mcMotherPt = 1;
1378         }
1379
1380         else {
1381
1382           TParticle  *lMCESDMother    = stack->Particle(lIndexPosMother);
1383           if (!lMCESDMother)    { continue;}
1384           lPdgcodeMother         = lMCESDMother->GetPdgCode();
1385           lIndexMotherOfMother   = lMCESDMother->GetFirstMother();
1386           if (lIndexMotherOfMother ==-1) lPdgcodeMotherOfMother = 0;
1387           else {
1388             TParticle  *lMCESDMotherOfMother    = stack->Particle(lIndexMotherOfMother);
1389             if (!lMCESDMotherOfMother)  { continue;}
1390             lPdgcodeMotherOfMother = lMCESDMotherOfMother->GetPdgCode();
1391           }
1392
1393           mcPosX = lMCESDPartPos->Vx();
1394           mcPosY = lMCESDPartPos->Vy();
1395           mcPosZ = lMCESDPartPos->Vz();
1396           mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);
1397           mcPosMotherX = lMCESDMother->Vx();
1398           mcPosMotherY = lMCESDMother->Vy();
1399           mcPosMotherZ = lMCESDMother->Vz();
1400           mcPosMotherR = TMath::Sqrt(mcPosMotherX*mcPosMotherX+mcPosMotherY*mcPosMotherY);
1401
1402           //  mcMotherPt   = lMCESDMother->Pt();
1403         }
1404       }
1405     } // end ESD condition
1406
1407     // Multiplicity:
1408     if(!lOnFlyStatus) nv0s++;
1409     //    else  if(lOnFlyStatus) nv0sMI++;
1410
1411     // Daughter momentum cut: ! FIX it in case of AOD !
1412     //if ( (lPtPos  < cutMinPtDaughter ) ||
1413     //     (lPtNeg  < cutMinPtDaughter )
1414     //  )       { continue;}
1415
1416     // Look for associated particles:
1417     if (fAnalysisMC) {
1418       if( (lIndexPosMother==-1) || (lIndexNegMother==-1) ) {
1419         //      fHistMCDaughterTrack->Fill(1);
1420       }
1421
1422       else if( ( (lPDGCodePosDaughter==+211) && (lPDGCodeNegDaughter==-211) ) )
1423         {
1424           lCheckPIdK0Short    = 1;
1425           //      fHistMCDaughterTrack->Fill(3);
1426           if ( (lIndexPosMother==lIndexNegMother) &&
1427                (lPdgcodeMother==310) ) {
1428
1429
1430             //if (lIndexPosMother <= lNbMCPrimary) lCheckMcK0Short  = 1;
1431             //else lCheckSecondaryK0s = 1;
1432
1433             Double_t dx = 0;
1434             Double_t dy = 0;
1435             Double_t dz = 0;
1436             Double_t ProdDistance = 0;
1437
1438             dx = ( (mcPrimaryVtx.At(0)) - (mcPosMotherX) );
1439             dy = ( (mcPrimaryVtx.At(1)) - (mcPosMotherY) );
1440             dz = ( (mcPrimaryVtx.At(2)) - (mcPosMotherZ) );
1441
1442             ProdDistance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
1443
1444             if (ProdDistance < 0.001) lCheckMcK0Short  = 1;
1445             else lCheckSecondaryK0s = 1;
1446
1447           }
1448         }
1449       else if( ( (lPDGCodePosDaughter==+2212) && (lPDGCodeNegDaughter==-211)  )  )
1450         {
1451           lCheckPIdLambda     = 1;
1452           //      fHistMCDaughterTrack->Fill(5);
1453           if ( (lIndexPosMother==lIndexNegMother) &&
1454                (lPdgcodeMother==3122)  ){
1455             if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
1456                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
1457                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
1458                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
1459                  ) lComeFromSigma = 1;
1460             else lComeFromSigma = 0;
1461
1462
1463             // if ( (lIndexPosMother <= lNbMCPrimary) ||
1464             //     ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
1465             //      ) lCheckMcLambda  = 1;
1466             // else lCheckSecondaryLambda    = 1;
1467
1468             Double_t dx = 0;
1469             Double_t dy = 0;
1470             Double_t dz = 0;
1471             Double_t ProdDistance = 0;
1472
1473             dx = ( (mcPrimaryVtx.At(0)) - (mcPosMotherX) );
1474             dy = ( (mcPrimaryVtx.At(1)) - (mcPosMotherY) );
1475             dz = ( (mcPrimaryVtx.At(2)) - (mcPosMotherZ) );
1476
1477             ProdDistance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
1478
1479             if (ProdDistance < 0.001) lCheckMcLambda  = 1;
1480             else lCheckSecondaryLambda = 1;
1481
1482
1483           }
1484         }
1485       else if( ( (lPDGCodePosDaughter==211)   && (lPDGCodeNegDaughter==-2212) ) )
1486         {
1487           lCheckPIdAntiLambda = 1;
1488           //      fHistMCDaughterTrack->Fill(7);
1489           if ( (lIndexPosMother==lIndexNegMother) &&
1490                (lPdgcodeMother==-3122) ) {
1491             if ( ( TMath::Abs(lPdgcodeMotherOfMother) == 3212) ||
1492                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3224) ||
1493                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3214) ||
1494                  ( TMath::Abs(lPdgcodeMotherOfMother) == 3114)
1495                  ) lComeFromSigma = 1;
1496             else lComeFromSigma = 0;
1497
1498             //  if ( (lIndexPosMother <= lNbMCPrimary) ||
1499             //       ( ( lIndexPosMother > lNbMCPrimary) && (lComeFromSigma) )
1500             //       ) lCheckMcAntiLambda  = 1;
1501             //  else lCheckSecondaryAntiLambda = 1;
1502
1503             Double_t dx = 0;
1504             Double_t dy = 0;
1505             Double_t dz = 0;
1506             Double_t ProdDistance = 0;
1507
1508             dx = ( (mcPrimaryVtx.At(0)) - (mcPosMotherX) );
1509             dy = ( (mcPrimaryVtx.At(1)) - (mcPosMotherY) );
1510             dz = ( (mcPrimaryVtx.At(2)) - (mcPosMotherZ) );
1511
1512             ProdDistance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
1513
1514             if (ProdDistance < 0.001) lCheckMcAntiLambda = 1;
1515             else lCheckSecondaryAntiLambda = 1;
1516
1517           }
1518         }
1519
1520       // Gamma conversion
1521       //   else if ( (lPDGCodePosDaughter==-11) &&
1522       //        (lPDGCodeNegDaughter==11) &&
1523       //        (lPdgcodeMother==22 ) )
1524       //        lCheckGamma = 1;
1525
1526     } // end "look for associated particles
1527
1528     /////////////////////////////////////
1529     // PID condition for daughters tracks
1530     //////////////////////////////////////
1531
1532     //    lCheckPIDK0sPosDaughter        = 0, lCheckPIDK0sNegDaughter        = 0;
1533     lCheckPIDLambdaPosDaughter     = 0;//, lCheckPIDLambdaNegDaughter     = 0;
1534     //lCheckPIDAntiLambdaPosDaughter = 0;,
1535     lCheckPIDAntiLambdaNegDaughter = 0;
1536
1537     if (lMomInnerWallPos < lLimitPPID) {
1538       if (nSigmaPosPion < cutNSigmaLowP)   {
1539         //      lCheckPIDK0sPosDaughter        = 1;
1540         //lCheckPIDAntiLambdaPosDaughter = 1;
1541       }
1542       if (nSigmaPosProton < cutNSigmaLowP) lCheckPIDLambdaPosDaughter    = 1;
1543     }
1544
1545     else if (lMomInnerWallPos > lLimitPPID) {
1546       if (nSigmaPosPion < cutNSigmaHighP)   {
1547         //      lCheckPIDK0sPosDaughter        = 1;
1548         //      lCheckPIDAntiLambdaPosDaughter = 1;
1549       }
1550       if (nSigmaPosProton < cutNSigmaHighP) lCheckPIDLambdaPosDaughter    = 1;
1551     }
1552
1553     if (lMomInnerWallNeg < lLimitPPID) {
1554       if (nSigmaNegPion < cutNSigmaLowP)    {
1555         //      lCheckPIDK0sNegDaughter       = 1;
1556         //      lCheckPIDLambdaNegDaughter    = 1;
1557       }
1558       if (nSigmaNegProton < cutNSigmaLowP)  lCheckPIDAntiLambdaNegDaughter = 1;
1559
1560     }
1561     else if (lMomInnerWallNeg > lLimitPPID) {
1562       if (nSigmaNegPion < cutNSigmaHighP)   {
1563         //      lCheckPIDK0sNegDaughter       = 1;
1564         //      lCheckPIDLambdaNegDaughter    = 1;
1565       }
1566       if (nSigmaNegProton < cutNSigmaHighP) lCheckPIDAntiLambdaNegDaughter = 1;
1567     }
1568
1569
1570
1571     ///////////////values for cuts/////////////////////////////////////////////////////////////////////////////////////////
1572     if ((lDcaPosToPrimVertex < 0.1) || (lDcaNegToPrimVertex < 0.1) || (lDcaV0Daughters > 1.00) ||
1573         (lV0cosPointAngle < 0.998) || (lV0Radius < 0.9) || (lV0Radius > 100) )
1574
1575       {continue;}
1576
1577
1578     /////////////////////////////////
1579     //PID for Lambda and AntiLambda
1580     /////////////////////////////////
1581     if(fUsePID.Contains("withPID") && (lCheckPIDAntiLambdaNegDaughter==0) && (lCheckPIDLambdaPosDaughter==1)) LambdaPID = 1;
1582     else LambdaPID =0;
1583     if(fUsePID.Contains("withPID") && (lCheckPIDLambdaPosDaughter==0) && (lCheckPIDAntiLambdaNegDaughter==1)) AntiLambdaPID = 1;
1584     else AntiLambdaPID =0;
1585
1586     //ctau for lambda
1587     lcTauLambda     = (lV0tDecayLength*lInvMassLambda)/lPtLambda;
1588     //ctau for antilambda
1589     lcTauAntiLambda = (lV0tDecayLength*lInvMassAntiLambda)/lPtAntiLambda;
1590     //ctau for K0s
1591     lcTauK0s        = (lV0tDecayLength*lInvMassK0s)/lPtK0s;
1592
1593     //*****************************
1594     // filling histograms
1595     //*****************************
1596     if (lPLambda <1 && lOnFlyStatus==0 ){
1597       fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0);
1598     }
1599
1600     ////////////////
1601     //K0s particle
1602     ////////////////
1603     if (TMath::Abs(lRapK0s) < lCutRap ) {
1604       if (lcTauK0s< cutcTauK0s) {
1605
1606         if (lOnFlyStatus==0){
1607           fHistMassK0->Fill(lInvMassK0s);
1608           fHistPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
1609             } //end !lOnFlystatus
1610           } //end ctau cut
1611     } //end Rap condition
1612
1613     ///////////////////
1614     //Lambda particle
1615     ///////////////////
1616     if ((LambdaPID==1 && lMomInnerWallPos <=1 ) || (lMomInnerWallPos > 1) ||  !(fUsePID.Contains("withPID")  )){
1617       if (TMath::Abs(lRapLambda) < lCutRap) {
1618         if (lcTauLambda < cutcTauL){
1619           if (lOnFlyStatus==0){
1620             fHistMassLambda->Fill(lInvMassLambda);
1621             fHistPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
1622             if(lPtLambda <=1) fHistNSigmaProton->Fill(nSigmaPosProton);
1623           }
1624         }// end ctau condition
1625       } //end of Rap condition
1626     }// end of PID condition
1627
1628     //////////////////////////////
1629     // Anti Lambda ///////////////
1630     //////////////////////////////
1631     if ((AntiLambdaPID==1 && lMomInnerWallNeg <=1) || (lMomInnerWallNeg>1) ||  !(fUsePID.Contains("withPID"))){
1632       if (TMath::Abs(lRapAntiLambda) < lCutRap) {
1633         if (lcTauAntiLambda < cutcTauL){
1634           if (lOnFlyStatus==0){
1635             fHistMassAntiLambda->Fill(lInvMassAntiLambda);
1636             fHistPtVsMassAntiLambda->Fill(lInvMassAntiLambda,lPtAntiLambda);
1637           }
1638         } //end of Rap condition
1639       } // end of PID condition
1640     } //end ctau cut
1641
1642
1643
1644
1645     // Histo versus Rap and armenteros plot
1646     if (!lOnFlyStatus){
1647       if (lCheckMcK0Short) fHistAsMcRapK0->Fill(lRapK0s);
1648       if (lCheckMcLambda) fHistAsMcRapLambda->Fill(lRapLambda);
1649       if (lCheckMcAntiLambda) fHistAsMcRapLambda->Fill(lRapAntiLambda);
1650       //      fHistArmenterosPodolanski->Fill(lAlphaV0,lPtArmV0);
1651       if ((TMath::Abs(lRapK0s) < lCutRap)&&(TMath::Abs(lRapLambda) < lCutRap)) fHistK0sMassVsLambdaMass->Fill(lInvMassK0s,lInvMassLambda);
1652     }
1653
1654
1655     ///////////////////////////////////////////////////
1656     // K0s associated histograms in |rap| < lCutRap:
1657     ///////////////////////////////////////////////////
1658     if (TMath::Abs(lRapK0s) < lCutRap) {
1659       switch (lOnFlyStatus){
1660       case 0 :
1661         if (lcTauK0s< cutcTauK0s) {
1662           if(lCheckPIdK0Short) fHistPidMcMassK0->Fill(lInvMassK0s);
1663           if(lCheckMcK0Short) {
1664             fHistAsMcMassK0->Fill(lInvMassK0s);
1665             fHistAsMcPtK0->Fill(lPtK0s);
1666             fHistAsMcPtVsMassK0->Fill(lInvMassK0s,lPtK0s);
1667           }
1668         }
1669       } // end rapidity condition
1670     } // end ctau cut
1671
1672     ///////////////////////////////////////////////////
1673     // Associated Lambda histograms in |rap| < lCutRap
1674     ////////////////////////////////////////////////////
1675      if ((LambdaPID==1 && lMomInnerWallPos <=1) || (lMomInnerWallPos>1) ||  !(fUsePID.Contains("withPID"))){
1676     if (TMath::Abs(lRapLambda) < lCutRap) {
1677       switch (lOnFlyStatus){
1678       case 0 :
1679         if (lcTauLambda < cutcTauL){
1680           if(lCheckPIdLambda) fHistPidMcMassLambda->Fill(lInvMassLambda);
1681           if(lCheckMcLambda) {
1682             fHistAsMcMassLambda->Fill(lInvMassLambda);
1683             fHistAsMcPtLambda->Fill(lPtLambda);
1684             //    fHistCosPointAngleLvsMassVsPtsigL->Fill(lPtLambda,lV0cosPointAngle,lInvMassLambda);
1685             fHistAsMcPtVsMassLambda->Fill(lInvMassLambda,lPtLambda);
1686           }
1687         }
1688       } // end rapidity condition
1689      }// end PID condition
1690     }// end ctau condition
1691
1692     ////////////////////////////////////////////////////////
1693     // Associated AntiLambda histograms in |rap| < lCutRap
1694     ////////////////////////////////////////////////////////
1695       if ((AntiLambdaPID==1 && lMomInnerWallNeg <=1) || (lMomInnerWallNeg>1) ||  !(fUsePID.Contains("withPID"))){
1696     if (TMath::Abs(lRapAntiLambda) < lCutRap) {
1697       switch (lOnFlyStatus){
1698       case 0 :
1699         if (lcTauAntiLambda < cutcTauL){
1700           if(lCheckPIdAntiLambda) fHistPidMcMassAntiLambda->Fill(lInvMassAntiLambda);
1701           if(lCheckMcAntiLambda) {
1702             fHistAsMcMassAntiLambda->Fill(lInvMassAntiLambda);
1703             fHistAsMcPtAntiLambda->Fill(lPtAntiLambda);
1704             fHistAsMcPtVsMassAntiLambda->Fill(lInvMassAntiLambda,lPtAntiLambda);
1705           }
1706         }
1707       } // end rapidity condition
1708      }// end PID condition
1709     }// end ctau cut
1710
1711   } // end V0 loop
1712
1713   fHistV0Multiplicity->Fill(nv0s);
1714   //  fHistV0MultiplicityMI->Fill(nv0sMI);
1715
1716   if (fAnalysisType == "ESD") { if(myPrimaryVertex) delete myPrimaryVertex; }
1717   //  if (fAnalysisType == "ESD") { if(TestTrackCuts) delete TestTrackCuts; }
1718
1719
1720   // Post output data
1721 }
1722
1723 //________________________________________________________________________
1724 void AliAnalysisTaskPerformanceStrange::Terminate(Option_t *)
1725 {/*
1726  // Draw result to the screen
1727  // Called once at the end of the query
1728
1729  TList *cRetrievedList = 0x0;
1730  cRetrievedList = (TList*)GetOutputData(1);
1731
1732  if(!cRetrievedList){
1733  AliWarning("ERROR - AliAnalysisTaskPerformanceStrange: output data container list not available\n"); return;
1734  }
1735
1736
1737  fHistV0Multiplicity = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0Multiplicity"));
1738  if (!fHistV0Multiplicity) {
1739  Printf("ERROR: fHistV0Multiplicity not available");
1740  return;
1741  }
1742
1743  fHistV0MultiplicityMI = dynamic_cast<TH1F*> ( cRetrievedList->FindObject("fHistV0MultiplicityMI"));
1744  if (!fHistV0MultiplicityMI) {
1745  Printf("ERROR: fHistV0MultiplicityMI not available");
1746  return;
1747  }
1748
1749  TCanvas *canPerformanceStrange = new TCanvas("AliAnalysisTaskCheckV0","Multiplicity",10,10,510,510);
1750  canPerformanceStrange->Divide(2,1);
1751  if (fHistV0Multiplicity->GetMaximum() > 0.) canPerformanceStrange->cd(1)->SetLogy();
1752  fHistV0Multiplicity->SetMarkerStyle(25);
1753  fHistV0Multiplicity->DrawCopy("E");
1754  if (fHistV0MultiplicityMI->GetMaximum() > 0.) canPerformanceStrange->cd(2)->SetLogy();
1755  fHistV0MultiplicityMI->SetMarkerStyle(24);
1756  fHistV0MultiplicityMI->DrawCopy("E");
1757
1758
1759  */
1760 }
1761
1762 //----------------------------------------------------------------------------
1763
1764 Double_t AliAnalysisTaskPerformanceStrange::MyRapidity(Double_t rE, Double_t rPz) const
1765 {
1766   // Local calculation for rapidity
1767   return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1768 }
1769 //----------------------------------------------------------------------------
1770