PWGLF/STRANGENESS/LambdaK0/AliAnalysisTaskExtractPerformanceV0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  17 //
  18 // Modified version of AliAnalysisTaskCheckCascade.cxx.
  19 // This is a 'hybrid' output version, in that it uses a classic TTree
  20 // ROOT object to store the candidates, plus a couple of histograms filled on
  21 // a per-event basis for storing variables too numerous to put in a tree.
  22 //
  23 // --- Adapted to look for lambdas as well, using code from
  24 //        AliAnalysisTaskCheckPerformanceStrange.cxx
  25 //
  26 //  --- Algorithm Description
  27 //   1. Loop over primaries in stack to acquire generated charged Xi
  28 //   2. Loop over stack to find V0s, fill TH3Fs "PrimRawPt"s for Efficiency
  29 //   3. Perform Physics Selection
  30 //   4. Perform Primary Vertex |z|<10cm selection
  31 //   5. Perform Primary Vertex NoTPCOnly vertexing selection (>0 contrib.)
  32 //   6. Perform Pileup Rejection
  33 //   7. Analysis Loops:
  34 //    7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
  35 //    7b. Fill TTree object with V0 information, candidates
  36 //
  37 //  Please Report Any Bugs!
  38 //
  39 //   --- David Dobrigkeit Chinellato
  40 //        (david.chinellato@gmail.com)
  41 //
  42 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  43
  44 class TTree;
  45 class TParticle;
  46 class TVector3;
  47
  48 //class AliMCEventHandler;
  49 //class AliMCEvent;
  50 //class AliStack;
  51
  52 class AliESDVertex;
  53 class AliAODVertex;
  54 class AliESDv0;
  55 class AliAODv0;
  56
  57 #include <Riostream.h>
  58 #include "TList.h"
  59 #include "TH1.h"
  60 #include "TH2.h"
  61 #include "TH3.h"
  62 #include "TFile.h"
  63 #include "THnSparse.h"
  64 #include "TVector3.h"
  65 #include "TCanvas.h"
  66 #include "TMath.h"
  67 #include "TLegend.h"
  68 #include "AliLog.h"
  69
  70 #include "AliESDEvent.h"
  71 #include "AliAODEvent.h"
  72 #include "AliV0vertexer.h"
  73 #include "AliCascadeVertexer.h"
  74 #include "AliESDpid.h"
  75 #include "AliESDtrack.h"
  76 #include "AliESDtrackCuts.h"
  77 #include "AliInputEventHandler.h"
  78 #include "AliAnalysisManager.h"
  79 #include "AliMCEventHandler.h"
  80 #include "AliMCEvent.h"
  81 #include "AliStack.h"
  82
  83 #include "AliCFContainer.h"
  84 #include "AliMultiplicity.h"
  85 #include "AliAODMCParticle.h"
  86 #include "AliESDcascade.h"
  87 #include "AliAODcascade.h"
  88 #include "AliESDUtils.h"
  89 #include "AliGenEventHeader.h"
  90
  91 #include "AliAnalysisTaskExtractPerformanceV0.h"
  92
  93 using std::cout;
  94 using std::endl;
  95
  96 ClassImp(AliAnalysisTaskExtractPerformanceV0)
  97
  98 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0()
  99   : AliAnalysisTaskSE(), fListHistV0(0), fTree(0), fPIDResponse(0), fESDtrackCuts(0),
 100    fkIsNuclear   ( kFALSE ),
 101    fkLowEnergyPP ( kFALSE ),
 102    fkUseOnTheFly ( kFALSE ),
 103
 104 //------------------------------------------------
 105 // Tree Variables
 106
 107   fTreeVariablePrimaryStatus(0),
 108   fTreeVariablePrimaryStatusMother(0),
 109   fTreeVariableChi2V0(0),
 110         fTreeVariableDcaV0Daughters(0),
 111         fTreeVariableDcaV0ToPrimVertex(0),
 112         fTreeVariableDcaPosToPrimVertex(0),
 113         fTreeVariableDcaNegToPrimVertex(0),
 114         fTreeVariableV0CosineOfPointingAngle(0),
 115         fTreeVariableV0Radius(0),
 116         fTreeVariablePt(0),
 117         fTreeVariablePtMC(0),
 118         fTreeVariableRapK0Short(0),
 119         fTreeVariableRapLambda(0),
 120         fTreeVariableRapMC(0),
 121         fTreeVariableInvMassK0s(0),
 122         fTreeVariableInvMassLambda(0),
 123         fTreeVariableInvMassAntiLambda(0),
 124         fTreeVariableAlphaV0(0),
 125         fTreeVariablePtArmV0(0),
 126         fTreeVariableNegTotMomentum(0),
 127         fTreeVariablePosTotMomentum(0),
 128         fTreeVariableNegTransvMomentum(0),
 129         fTreeVariablePosTransvMomentum(0),
 130         fTreeVariableNegTransvMomentumMC(0),
 131         fTreeVariablePosTransvMomentumMC(0),
 132
 133         fTreeVariableNSigmasPosProton(0),
 134         fTreeVariableNSigmasPosPion(0),
 135         fTreeVariableNSigmasNegProton(0),
 136         fTreeVariableNSigmasNegPion(0),
 137
 138         fTreeVariablePtMother(0),
 139         fTreeVariableV0CreationRadius(0),
 140   fTreeVariablePID(0),
 141   fTreeVariablePIDPositive(0),
 142   fTreeVariablePIDNegative(0),
 143   fTreeVariablePIDMother(0),
 144   fTreeVariableIndexStatus(0),
 145   fTreeVariableIndexStatusMother(0),
 146
 147   fTreeVariableRunNumber(0),
 148   fTreeVariableEventNumber(0),
 149
 150         fTreeVariableDistOverTotMom(0),
 151
 152         fTreeVariablePosEta(0),
 153         fTreeVariableNegEta(0),
 154
 155         fTreeVariableVertexZ(0),
 156
 157   fTreeVariableLeastNbrCrossedRows(0),
 158   fTreeVariableLeastRatioCrossedRowsOverFindable(0),
 159   fTreeVariableMultiplicity(0),
 160   fTreeVariableMultiplicityMC(0),
 161
 162   fTreeVariableV0x(0),
 163   fTreeVariableV0y(0),
 164   fTreeVariableV0z(0),
 165
 166   fTreeVariableV0Px(0),
 167   fTreeVariableV0Py(0),
 168   fTreeVariableV0Pz(0),
 169
 170   fTreeVariableMCV0x(0),
 171   fTreeVariableMCV0y(0),
 172   fTreeVariableMCV0z(0),
 173
 174   fTreeVariableMCV0Px(0),
 175   fTreeVariableMCV0Py(0),
 176   fTreeVariableMCV0Pz(0),
 177
 178   fTreeVariablePVx(0),
 179   fTreeVariablePVy(0),
 180   fTreeVariablePVz(0),
 181
 182   fTreeVariableMCPVx(0),
 183   fTreeVariableMCPVy(0),
 184   fTreeVariableMCPVz(0),
 185
 186 //------------------------------------------------
 187 // HISTOGRAMS
 188 // --- Filled on an Event-by-event basis
 189 //------------------------------------------------
 190    fHistV0MultiplicityBeforeTrigSel(0),
 191    fHistV0MultiplicityForTrigEvt(0),
 192    fHistV0MultiplicityForSelEvt(0),
 193    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 194    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 195    fHistMultiplicityBeforeTrigSel(0),
 196    fHistMultiplicityForTrigEvt(0),
 197    fHistMultiplicity(0),
 198    fHistMultiplicityNoTPCOnly(0),
 199    fHistMultiplicityNoTPCOnlyNoPileup(0),
 200
 201         f2dHistMultiplicityVsTrueBeforeTrigSel(0),
 202         f2dHistMultiplicityVsTrueForTrigEvt(0),
 203         f2dHistMultiplicityVsTrue(0),
 204         f2dHistMultiplicityVsTrueNoTPCOnly(0),
 205         f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup(0),
 206
 207   //Raw Data for Vertex Z position estimator change
 208         f2dHistMultiplicityVsVertexZBeforeTrigSel(0),
 209         f2dHistMultiplicityVsVertexZForTrigEvt(0),
 210         f2dHistMultiplicityVsVertexZ(0),
 211         f2dHistMultiplicityVsVertexZNoTPCOnly(0),
 212         f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup(0),
 213
 214   fHistGenVertexZBeforeTrigSel(0),
 215   fHistGenVertexZForTrigEvt(0),
 216   fHistGenVertexZ(0),
 217   fHistGenVertexZNoTPCOnly(0),
 218   fHistGenVertexZNoTPCOnlyNoPileup(0),
 219
 220 //------------------------------------------------
 221 // PARTICLE HISTOGRAMS
 222 // --- Filled on a Particle-by-Particle basis
 223 //------------------------------------------------
 224    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
 225    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
 226    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
 227    f3dHistPrimRawPtVsYVsMultLambda(0),
 228    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
 229    f3dHistPrimRawPtVsYVsMultK0Short(0),
 230    f3dHistPrimRawPtVsYVsMultMCLambda(0),
 231    f3dHistPrimRawPtVsYVsMultMCAntiLambda(0),
 232    f3dHistPrimRawPtVsYVsMultMCK0Short(0),
 233    f3dHistPrimRawPtVsYVsVertexZLambda(0),
 234    f3dHistPrimRawPtVsYVsVertexZAntiLambda(0),
 235    f3dHistPrimRawPtVsYVsVertexZK0Short(0),
 236    f3dHistPrimCloseToPVPtVsYVsMultLambda(0),
 237    f3dHistPrimCloseToPVPtVsYVsMultAntiLambda(0),
 238    f3dHistPrimCloseToPVPtVsYVsMultK0Short(0),
 239    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
 240    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
 241    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
 242    f3dHistGenPtVsYVsMultXiMinus(0),
 243    f3dHistGenPtVsYVsMultXiPlus(0),
 244    f3dHistGenPtVsYVsMultOmegaMinus(0),
 245    f3dHistGenPtVsYVsMultOmegaPlus(0),
 246    f3dHistGenSelectedPtVsYVsMultXiMinus(0),
 247    f3dHistGenSelectedPtVsYVsMultXiPlus(0),
 248    f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
 249    f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
 250    fHistPVx(0),
 251    fHistPVy(0),
 252    fHistPVz(0),
 253    fHistPVxAnalysis(0),
 254    fHistPVyAnalysis(0),
 255    fHistPVzAnalysis(0),
 256    fHistPVxAnalysisHasHighPtLambda(0),
 257    fHistPVyAnalysisHasHighPtLambda(0),
 258    fHistPVzAnalysisHasHighPtLambda(0),
 259    fHistSwappedV0Counter(0)
 260 {
 261   // Dummy Constructor
 262 }
 263
 264 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0(const char *name)
 265   : AliAnalysisTaskSE(name), fListHistV0(0), fTree(0), fPIDResponse(0), fESDtrackCuts(0),
 266    fkIsNuclear   ( kFALSE ),
 267    fkLowEnergyPP ( kFALSE ),
 268    fkUseOnTheFly ( kFALSE ),
 269
 270 //------------------------------------------------
 271 // Tree Variables
 272
 273   fTreeVariablePrimaryStatus(0),
 274   fTreeVariablePrimaryStatusMother(0),
 275   fTreeVariableChi2V0(0),
 276         fTreeVariableDcaV0Daughters(0),
 277         fTreeVariableDcaV0ToPrimVertex(0),
 278         fTreeVariableDcaPosToPrimVertex(0),
 279         fTreeVariableDcaNegToPrimVertex(0),
 280         fTreeVariableV0CosineOfPointingAngle(0),
 281         fTreeVariableV0Radius(0),
 282         fTreeVariablePt(0),
 283         fTreeVariablePtMC(0),
 284         fTreeVariableRapK0Short(0),
 285         fTreeVariableRapLambda(0),
 286         fTreeVariableRapMC(0),
 287         fTreeVariableInvMassK0s(0),
 288         fTreeVariableInvMassLambda(0),
 289         fTreeVariableInvMassAntiLambda(0),
 290         fTreeVariableAlphaV0(0),
 291         fTreeVariablePtArmV0(0),
 292         fTreeVariableNegTotMomentum(0),
 293         fTreeVariablePosTotMomentum(0),
 294         fTreeVariableNegTransvMomentum(0),
 295         fTreeVariablePosTransvMomentum(0),
 296         fTreeVariableNegTransvMomentumMC(0),
 297         fTreeVariablePosTransvMomentumMC(0),
 298
 299         fTreeVariableNSigmasPosProton(0),
 300         fTreeVariableNSigmasPosPion(0),
 301         fTreeVariableNSigmasNegProton(0),
 302         fTreeVariableNSigmasNegPion(0),
 303
 304         fTreeVariablePtMother(0),
 305         fTreeVariableV0CreationRadius(0),
 306   fTreeVariablePID(0),
 307   fTreeVariablePIDPositive(0),
 308   fTreeVariablePIDNegative(0),
 309   fTreeVariablePIDMother(0),
 310   fTreeVariableIndexStatus(0),
 311   fTreeVariableIndexStatusMother(0),
 312
 313   fTreeVariableRunNumber(0),
 314   fTreeVariableEventNumber(0),
 315
 316         fTreeVariableDistOverTotMom(0),
 317
 318         fTreeVariablePosEta(0),
 319         fTreeVariableNegEta(0),
 320
 321         fTreeVariableVertexZ(0),
 322
 323   fTreeVariableLeastNbrCrossedRows(0),
 324   fTreeVariableLeastRatioCrossedRowsOverFindable(0),
 325   fTreeVariableMultiplicity(0),
 326   fTreeVariableMultiplicityMC(0),
 327
 328   fTreeVariableV0x(0),
 329   fTreeVariableV0y(0),
 330   fTreeVariableV0z(0),
 331
 332   fTreeVariableV0Px(0),
 333   fTreeVariableV0Py(0),
 334   fTreeVariableV0Pz(0),
 335
 336   fTreeVariableMCV0x(0),
 337   fTreeVariableMCV0y(0),
 338   fTreeVariableMCV0z(0),
 339
 340   fTreeVariableMCV0Px(0),
 341   fTreeVariableMCV0Py(0),
 342   fTreeVariableMCV0Pz(0),
 343
 344   fTreeVariablePVx(0),
 345   fTreeVariablePVy(0),
 346   fTreeVariablePVz(0),
 347
 348   fTreeVariableMCPVx(0),
 349   fTreeVariableMCPVy(0),
 350   fTreeVariableMCPVz(0),
 351
 352 //------------------------------------------------
 353 // HISTOGRAMS
 354 // --- Filled on an Event-by-event basis
 355 //------------------------------------------------
 356    fHistV0MultiplicityBeforeTrigSel(0),
 357    fHistV0MultiplicityForTrigEvt(0),
 358    fHistV0MultiplicityForSelEvt(0),
 359    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 360    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 361    fHistMultiplicityBeforeTrigSel(0),
 362    fHistMultiplicityForTrigEvt(0),
 363    fHistMultiplicity(0),
 364    fHistMultiplicityNoTPCOnly(0),
 365    fHistMultiplicityNoTPCOnlyNoPileup(0),
 366
 367         f2dHistMultiplicityVsTrueBeforeTrigSel(0),
 368         f2dHistMultiplicityVsTrueForTrigEvt(0),
 369         f2dHistMultiplicityVsTrue(0),
 370         f2dHistMultiplicityVsTrueNoTPCOnly(0),
 371         f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup(0),
 372
 373   //Raw Data for Vertex Z position estimator change
 374         f2dHistMultiplicityVsVertexZBeforeTrigSel(0),
 375         f2dHistMultiplicityVsVertexZForTrigEvt(0),
 376         f2dHistMultiplicityVsVertexZ(0),
 377         f2dHistMultiplicityVsVertexZNoTPCOnly(0),
 378         f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup(0),
 379
 380   fHistGenVertexZBeforeTrigSel(0),
 381   fHistGenVertexZForTrigEvt(0),
 382   fHistGenVertexZ(0),
 383   fHistGenVertexZNoTPCOnly(0),
 384   fHistGenVertexZNoTPCOnlyNoPileup(0),
 385
 386 //------------------------------------------------
 387 // PARTICLE HISTOGRAMS
 388 // --- Filled on a Particle-by-Particle basis
 389 //------------------------------------------------
 390    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
 391    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
 392    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
 393    f3dHistPrimRawPtVsYVsMultLambda(0),
 394    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
 395    f3dHistPrimRawPtVsYVsMultK0Short(0),
 396    f3dHistPrimRawPtVsYVsMultMCLambda(0),
 397    f3dHistPrimRawPtVsYVsMultMCAntiLambda(0),
 398    f3dHistPrimRawPtVsYVsMultMCK0Short(0),
 399    f3dHistPrimRawPtVsYVsVertexZLambda(0),
 400    f3dHistPrimRawPtVsYVsVertexZAntiLambda(0),
 401    f3dHistPrimRawPtVsYVsVertexZK0Short(0),
 402    f3dHistPrimCloseToPVPtVsYVsMultLambda(0),
 403    f3dHistPrimCloseToPVPtVsYVsMultAntiLambda(0),
 404    f3dHistPrimCloseToPVPtVsYVsMultK0Short(0),
 405    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
 406    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
 407    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
 408    f3dHistGenPtVsYVsMultXiMinus(0),
 409    f3dHistGenPtVsYVsMultXiPlus(0),
 410    f3dHistGenPtVsYVsMultOmegaMinus(0),
 411    f3dHistGenPtVsYVsMultOmegaPlus(0),
 412    f3dHistGenSelectedPtVsYVsMultXiMinus(0),
 413    f3dHistGenSelectedPtVsYVsMultXiPlus(0),
 414    f3dHistGenSelectedPtVsYVsMultOmegaMinus(0),
 415    f3dHistGenSelectedPtVsYVsMultOmegaPlus(0),
 416    fHistPVx(0),
 417    fHistPVy(0),
 418    fHistPVz(0),
 419    fHistPVxAnalysis(0),
 420    fHistPVyAnalysis(0),
 421    fHistPVzAnalysis(0),
 422    fHistPVxAnalysisHasHighPtLambda(0),
 423    fHistPVyAnalysisHasHighPtLambda(0),
 424    fHistPVzAnalysisHasHighPtLambda(0),
 425    fHistSwappedV0Counter(0)
 426 {
 427    // Constructor
 428    // Output slot #0 writes into a TList container (Cascade)
 429    DefineOutput(1, TList::Class());
 430    DefineOutput(2, TTree::Class());
 431 }
 432
 433
 434 AliAnalysisTaskExtractPerformanceV0::~AliAnalysisTaskExtractPerformanceV0()
 435 {
 436 //------------------------------------------------
 437 // DESTRUCTOR
 438 //------------------------------------------------
 439
 440    if (fListHistV0){
 441       delete fListHistV0;
 442       fListHistV0 = 0x0;
 443    }
 444    if (fTree){
 445       delete fTree;
 446       fTree = 0x0;
 447    }
 448     //cleanup esd track cuts object too...
 449    if (fESDtrackCuts){
 450     delete fESDtrackCuts;
 451     fESDtrackCuts = 0x0;
 452   }
 453 }
 454
 455 //________________________________________________________________________
 456 void AliAnalysisTaskExtractPerformanceV0::UserCreateOutputObjects()
 457 {
 458
 459    OpenFile(2);
 460    // Called once
 461
 462 //------------------------------------------------
 463
 464    fTree = new TTree("fTree","V0Candidates");
 465
 466 //------------------------------------------------
 467 // fTree Branch definitions - V0 Tree
 468 //------------------------------------------------
 469
 470 //-----------BASIC-INFO---------------------------
 471 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatus",&fTreeVariablePrimaryStatus,"fTreeVariablePrimaryStatus/I");
 472 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatusMother",&fTreeVariablePrimaryStatusMother,"fTreeVariablePrimaryStatusMother/I");
 473 /* 2*/   fTree->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"Chi2V0/F");
 474 /* 3*/   fTree->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
 475 /* 4*/   fTree->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
 476 /* 5*/   fTree->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
 477 /* 6*/   fTree->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
 478 /* 7*/   fTree->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
 479 /* 7*/   fTree->Branch("fTreeVariablePtMC",&fTreeVariablePtMC,"fTreeVariablePtMC/F");
 480 /* 8*/   fTree->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
 481 /* 9*/   fTree->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
 482 /*10*/   fTree->Branch("fTreeVariableRapMC",&fTreeVariableRapMC,"fTreeVariableRapMC/F");
 483 /*11*/   fTree->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
 484 /*12*/   fTree->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
 485 /*13*/   fTree->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
 486 /*14*/   fTree->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
 487 /*15*/   fTree->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
 488 /*16*/   fTree->Branch("fTreeVariableNegTransvMomentum",&fTreeVariableNegTransvMomentum,"fTreeVariableNegTransvMomentum/F");
 489 /*17*/   fTree->Branch("fTreeVariablePosTransvMomentum",&fTreeVariablePosTransvMomentum,"fTreeVariablePosTransvMomentum/F");
 490 /*18*/   fTree->Branch("fTreeVariableNegTransvMomentumMC",&fTreeVariableNegTransvMomentumMC,"fTreeVariableNegTransvMomentumMC/F");
 491 /*19*/   fTree->Branch("fTreeVariablePosTransvMomentumMC",&fTreeVariablePosTransvMomentumMC,"fTreeVariablePosTransvMomentumMC/F");
 492 /*20*/   fTree->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
 493 /*21*/   fTree->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
 494 /*22*/   fTree->Branch("fTreeVariablePID",&fTreeVariablePID,"fTreeVariablePID/I");
 495 /*23*/   fTree->Branch("fTreeVariablePIDPositive",&fTreeVariablePIDPositive,"fTreeVariablePIDPositive/I");
 496 /*24*/   fTree->Branch("fTreeVariablePIDNegative",&fTreeVariablePIDNegative,"fTreeVariablePIDNegative/I");
 497 /*25*/   fTree->Branch("fTreeVariablePIDMother",&fTreeVariablePIDMother,"fTreeVariablePIDMother/I");
 498 /*26*/   fTree->Branch("fTreeVariablePtXiMother",&fTreeVariablePtMother,"fTreeVariablePtMother/F");
 499 /*27*/   fTree->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
 500 //-----------MULTIPLICITY-INFO--------------------
 501 /*28*/   fTree->Branch("fTreeVariableMultiplicity",&fTreeVariableMultiplicity,"fTreeVariableMultiplicity/I");
 502 /*28*/   fTree->Branch("fTreeVariableMultiplicityMC",&fTreeVariableMultiplicityMC,"fTreeVariableMultiplicityMC/I");
 503 //------------------------------------------------
 504 /*29*/   fTree->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
 505 /*30*/   fTree->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
 506 /*31*/   fTree->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
 507 /*32*/   fTree->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
 508 /*33*/   fTree->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
 509 //------------------------------------------------
 510 /*34*/   fTree->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
 511 /*35*/   fTree->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
 512 /*36*/   fTree->Branch("fTreeVariableV0CreationRadius",&fTreeVariableV0CreationRadius,"fTreeVariableV0CreationRadius/F");
 513 /*37*/   fTree->Branch("fTreeVariableIndexStatus",&fTreeVariableIndexStatus,"fTreeVariableIndexStatus/I");
 514 /*38*/   fTree->Branch("fTreeVariableIndexStatusMother",&fTreeVariableIndexStatusMother,"fTreeVariableIndexStatusMother/I");
 515
 516 /*39*/   fTree->Branch("fTreeVariableRunNumber",&fTreeVariableRunNumber,"fTreeVariableRunNumber/I");
 517 /*40*/   fTree->Branch("fTreeVariableEventNumber",&fTreeVariableEventNumber,"fTreeVariableEventNumber/l");
 518
 519 /*34*/   fTree->Branch("fTreeVariableVertexZ",&fTreeVariableVertexZ,"fTreeVariableVertexZ/F");
 520
 521 //-----------FOR CTAU DEBUGGING: Full Phase Space + Decay Position Info
 522         fTree->Branch("fTreeVariableV0x",&fTreeVariableV0x,"fTreeVariableV0x/F");
 523         fTree->Branch("fTreeVariableV0y",&fTreeVariableV0y,"fTreeVariableV0y/F");
 524         fTree->Branch("fTreeVariableV0z",&fTreeVariableV0z,"fTreeVariableV0z/F");
 525
 526         fTree->Branch("fTreeVariableV0Px",&fTreeVariableV0Px,"fTreeVariableV0Px/F");
 527         fTree->Branch("fTreeVariableV0Py",&fTreeVariableV0Py,"fTreeVariableV0Py/F");
 528         fTree->Branch("fTreeVariableV0Pz",&fTreeVariableV0Pz,"fTreeVariableV0Pz/F");
 529
 530 //-----------FOR CTAU DEBUGGING: Full Phase Space + Decay Position Info, perfect info from MC
 531         fTree->Branch("fTreeVariableMCV0x",&fTreeVariableMCV0x,"fTreeVariableMCV0x/F");
 532         fTree->Branch("fTreeVariableMCV0y",&fTreeVariableMCV0y,"fTreeVariableMCV0y/F");
 533         fTree->Branch("fTreeVariableMCV0z",&fTreeVariableMCV0z,"fTreeVariableMCV0z/F");
 534
 535         fTree->Branch("fTreeVariableMCV0Px",&fTreeVariableMCV0Px,"fTreeVariableMCV0Px/F");
 536         fTree->Branch("fTreeVariableMCV0Py",&fTreeVariableMCV0Py,"fTreeVariableMCV0Py/F");
 537         fTree->Branch("fTreeVariableMCV0Pz",&fTreeVariableMCV0Pz,"fTreeVariableMCV0Pz/F");
 538
 539 //-----------FOR CTAU DEBUGGING: Primary vertex info
 540         fTree->Branch("fTreeVariablePVx",&fTreeVariablePVx,"fTreeVariablePVx/F");
 541         fTree->Branch("fTreeVariablePVy",&fTreeVariablePVy,"fTreeVariablePVy/F");
 542         fTree->Branch("fTreeVariablePVz",&fTreeVariablePVz,"fTreeVariablePVz/F");
 543
 544         fTree->Branch("fTreeVariableMCPVx",&fTreeVariableMCPVx,"fTreeVariableMCPVx/F");
 545         fTree->Branch("fTreeVariableMCPVy",&fTreeVariableMCPVy,"fTreeVariableMCPVy/F");
 546         fTree->Branch("fTreeVariableMCPVz",&fTreeVariableMCPVz,"fTreeVariableMCPVz/F");
 547
 548 //------------------------------------------------
 549 // Particle Identification Setup
 550 //------------------------------------------------
 551
 552    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 553    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 554    fPIDResponse = inputHandler->GetPIDResponse();
 555
 556   // Multiplicity
 557
 558     if(! fESDtrackCuts ){
 559           fESDtrackCuts = new AliESDtrackCuts();
 560     }
 561
 562 //------------------------------------------------
 563 // V0 Multiplicity Histograms
 564 //------------------------------------------------
 565
 566    // Create histograms
 567    OpenFile(1);
 568    fListHistV0 = new TList();
 569    fListHistV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
 570
 571
 572    if(! fHistV0MultiplicityBeforeTrigSel) {
 573       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
 574          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
 575          25, 0, 25);
 576       fListHistV0->Add(fHistV0MultiplicityBeforeTrigSel);
 577    }
 578
 579    if(! fHistV0MultiplicityForTrigEvt) {
 580       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
 581          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
 582          25, 0, 25);
 583       fListHistV0->Add(fHistV0MultiplicityForTrigEvt);
 584    }
 585
 586    if(! fHistV0MultiplicityForSelEvt) {
 587       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
 588          "V0s per event;Nbr of V0s/Evt;Events",
 589          25, 0, 25);
 590       fListHistV0->Add(fHistV0MultiplicityForSelEvt);
 591    }
 592
 593    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
 594       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
 595          "V0s per event;Nbr of V0s/Evt;Events",
 596          25, 0, 25);
 597       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
 598    }
 599    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
 600       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
 601          "V0s per event;Nbr of V0s/Evt;Events",
 602          25, 0, 25);
 603       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
 604    }
 605
 606 //------------------------------------------------
 607 // Track Multiplicity Histograms
 608 //------------------------------------------------
 609
 610    if(! fHistMultiplicityBeforeTrigSel) {
 611       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
 612          "Tracks per event;Nbr of Tracks;Events",
 613          200, 0, 200);
 614       fListHistV0->Add(fHistMultiplicityBeforeTrigSel);
 615    }
 616    if(! fHistMultiplicityForTrigEvt) {
 617       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
 618          "Tracks per event;Nbr of Tracks;Events",
 619          200, 0, 200);
 620       fListHistV0->Add(fHistMultiplicityForTrigEvt);
 621    }
 622    if(! fHistMultiplicity) {
 623       fHistMultiplicity = new TH1F("fHistMultiplicity",
 624          "Tracks per event;Nbr of Tracks;Events",
 625          200, 0, 200);
 626       fListHistV0->Add(fHistMultiplicity);
 627    }
 628    if(! fHistMultiplicityNoTPCOnly) {
 629       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
 630          "Tracks per event;Nbr of Tracks;Events",
 631          200, 0, 200);
 632       fListHistV0->Add(fHistMultiplicityNoTPCOnly);
 633    }
 634    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
 635       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
 636          "Tracks per event;Nbr of Tracks;Events",
 637          200, 0, 200);
 638       fListHistV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
 639    }
 640
 641   //Raw Data for J/Psi paper Technique
 642         //TH2F    *f2dHistMultiplicityVsTrueBeforeTrigSel;              //! multiplicity distribution
 643         //TH2F    *f2dHistMultiplicityVsTrueForTrigEvt;                         //! multiplicity distribution
 644         //TH2F    *f2dHistMultiplicityVsTrue;                                                   //! multiplicity distribution
 645         //TH2F    *f2dHistMultiplicityVsTrueNoTPCOnly;                          //! multiplicity distribution
 646         //TH2F    *f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup;                  //! multiplicity distribution
 647
 648    if(! f2dHistMultiplicityVsTrueBeforeTrigSel) {
 649       f2dHistMultiplicityVsTrueBeforeTrigSel = new TH2F("f2dHistMultiplicityVsTrueBeforeTrigSel",
 650          "Tracks per event", 200, 0, 200, 200, 0, 200);
 651       fListHistV0->Add(f2dHistMultiplicityVsTrueBeforeTrigSel);
 652    }
 653    if(! f2dHistMultiplicityVsTrueForTrigEvt) {
 654       f2dHistMultiplicityVsTrueForTrigEvt = new TH2F("f2dHistMultiplicityVsTrueForTrigEvt",
 655          "Tracks per event", 200, 0, 200, 200, 0, 200);
 656       fListHistV0->Add(f2dHistMultiplicityVsTrueForTrigEvt);
 657    }
 658    if(! f2dHistMultiplicityVsTrue) {
 659       f2dHistMultiplicityVsTrue = new TH2F("f2dHistMultiplicityVsTrue",
 660          "Tracks per event", 200, 0, 200, 200, 0, 200);
 661       fListHistV0->Add(f2dHistMultiplicityVsTrue);
 662    }
 663    if(! f2dHistMultiplicityVsTrueNoTPCOnly) {
 664       f2dHistMultiplicityVsTrueNoTPCOnly = new TH2F("f2dHistMultiplicityVsTrueNoTPCOnly",
 665          "Tracks per event", 200, 0, 200, 200, 0, 200);
 666       fListHistV0->Add(f2dHistMultiplicityVsTrueNoTPCOnly);
 667    }
 668    if(! f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup) {
 669       f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup = new TH2F("f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup",
 670          "Tracks per event", 200, 0, 200, 200, 0, 200);
 671       fListHistV0->Add(f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup);
 672    }
 673
 674
 675   //Raw Data for Vertex Z position estimator change
 676         //TH2F    *f2dHistMultiplicityVsVertexZBeforeTrigSel;           //! multiplicity distribution
 677         //TH2F    *f2dHistMultiplicityVsVertexZForTrigEvt;                      //! multiplicity distribution
 678         //TH2F    *f2dHistMultiplicityVsVertexZ;                                                //! multiplicity distribution
 679         //TH2F    *f2dHistMultiplicityVsVertexZNoTPCOnly;                               //! multiplicity distribution
 680         //TH2F    *f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup;                       //! multiplicity distribution
 681
 682    if(! f2dHistMultiplicityVsVertexZBeforeTrigSel) {
 683       f2dHistMultiplicityVsVertexZBeforeTrigSel = new TH2F("f2dHistMultiplicityVsVertexZBeforeTrigSel",
 684          "Tracks per event", 200, 0, 200,400, -20, 20);
 685       fListHistV0->Add(f2dHistMultiplicityVsVertexZBeforeTrigSel);
 686    }
 687    if(! f2dHistMultiplicityVsVertexZForTrigEvt) {
 688       f2dHistMultiplicityVsVertexZForTrigEvt = new TH2F("f2dHistMultiplicityVsVertexZForTrigEvt",
 689          "Tracks per event", 200, 0, 200, 400, -20, 20);
 690       fListHistV0->Add(f2dHistMultiplicityVsVertexZForTrigEvt);
 691    }
 692    if(! f2dHistMultiplicityVsVertexZ) {
 693       f2dHistMultiplicityVsVertexZ = new TH2F("f2dHistMultiplicityVsVertexZ",
 694          "Tracks per event", 200, 0, 200, 400, -20, 20);
 695       fListHistV0->Add(f2dHistMultiplicityVsVertexZ);
 696    }
 697    if(! f2dHistMultiplicityVsVertexZNoTPCOnly) {
 698       f2dHistMultiplicityVsVertexZNoTPCOnly = new TH2F("f2dHistMultiplicityVsVertexZNoTPCOnly",
 699          "Tracks per event", 200, 0, 200, 400, -20, 20);
 700       fListHistV0->Add(f2dHistMultiplicityVsVertexZNoTPCOnly);
 701    }
 702    if(! f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup) {
 703       f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup = new TH2F("f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup",
 704          "Tracks per event", 200, 0, 200, 400, -20, 20);
 705       fListHistV0->Add(f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup);
 706    }
 707
 708 //Generated PVz test
 709
 710 //   TH1F      *fHistGenVertexZBeforeTrigSel;     //! multiplicity distribution
 711 //   TH1F      *fHistGenVertexZForTrigEvt;        //! multiplicity distribution
 712 //   TH1F      *fHistGenVertexZ;                  //! multiplicity distribution
 713 //   TH1F      *fHistGenVertexZNoTPCOnly;         //! multiplicity distribution
 714 //   TH1F      *fHistGenVertexZNoTPCOnlyNoPileup; //! multiplicity distribution
 715
 716    if(! fHistGenVertexZBeforeTrigSel) {
 717          fHistGenVertexZBeforeTrigSel = new TH1F("fHistGenVertexZBeforeTrigSel",
 718             "PV z position;Nbr of Evts;z",
 719             400, -20, 20);
 720       fListHistV0->Add(fHistGenVertexZBeforeTrigSel);
 721    }
 722    if(! fHistGenVertexZForTrigEvt) {
 723          fHistGenVertexZForTrigEvt = new TH1F("fHistGenVertexZForTrigEvt",
 724             "PV z position;Nbr of Evts;z",
 725             400, -20, 20);
 726       fListHistV0->Add(fHistGenVertexZForTrigEvt);
 727    }
 728    if(! fHistGenVertexZ) {
 729          fHistGenVertexZ = new TH1F("fHistGenVertexZ",
 730             "PV z position;Nbr of Evts;z",
 731             400, -20, 20);
 732       fListHistV0->Add(fHistGenVertexZ);
 733    }
 734    if(! fHistGenVertexZNoTPCOnly) {
 735          fHistGenVertexZNoTPCOnly = new TH1F("fHistGenVertexZNoTPCOnly",
 736             "PV z position;Nbr of Evts;z",
 737             400, -20, 20);
 738       fListHistV0->Add(fHistGenVertexZNoTPCOnly);
 739    }
 740    if(! fHistGenVertexZNoTPCOnlyNoPileup) {
 741          fHistGenVertexZNoTPCOnlyNoPileup = new TH1F("fHistGenVertexZNoTPCOnlyNoPileup",
 742             "PV z position;Nbr of Evts;z",
 743             400, -20, 20);
 744       fListHistV0->Add(fHistGenVertexZNoTPCOnlyNoPileup);
 745    }
 746
 747
 748 //------------------------------------------------
 749 // Generated Particle Histograms
 750 //------------------------------------------------
 751
 752    Int_t lCustomNBins = 200;
 753    Double_t lCustomPtUpperLimit = 20;
 754    Int_t lCustomNBinsMultiplicity = 100;
 755
 756 //----------------------------------
 757 // Raw Generated (Pre-physics-selection)
 758 //----------------------------------
 759
 760 //--- 3D Histo (Pt, Y, Multiplicity)
 761
 762    if(! f3dHistPrimRawPtVsYVsMultLambda) {
 763       f3dHistPrimRawPtVsYVsMultLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 764       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultLambda);
 765    }
 766    if(! f3dHistPrimRawPtVsYVsMultAntiLambda) {
 767       f3dHistPrimRawPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 768       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultAntiLambda);
 769    }
 770    if(! f3dHistPrimRawPtVsYVsMultK0Short) {
 771       f3dHistPrimRawPtVsYVsMultK0Short = new TH3F( "f3dHistPrimRawPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 772       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultK0Short);
 773    }
 774
 775 //--- 3D Histo (Pt, Y, MultiplicityMC)
 776
 777    if(! f3dHistPrimRawPtVsYVsMultMCLambda) {
 778       f3dHistPrimRawPtVsYVsMultMCLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultMCLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 779       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCLambda);
 780    }
 781    if(! f3dHistPrimRawPtVsYVsMultMCAntiLambda) {
 782       f3dHistPrimRawPtVsYVsMultMCAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsMultMCAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 783       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCAntiLambda);
 784    }
 785    if(! f3dHistPrimRawPtVsYVsMultMCK0Short) {
 786       f3dHistPrimRawPtVsYVsMultMCK0Short = new TH3F( "f3dHistPrimRawPtVsYVsMultMCK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; MultMC", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 787       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultMCK0Short);
 788    }
 789
 790 //--- 3D Histo (Pt, Y, VertexZ)
 791
 792    if(! f3dHistPrimRawPtVsYVsVertexZLambda) {
 793       f3dHistPrimRawPtVsYVsVertexZLambda = new TH3F( "f3dHistPrimRawPtVsYVsVertexZLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs VertexZiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
 794       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZLambda);
 795    }
 796    if(! f3dHistPrimRawPtVsYVsVertexZAntiLambda) {
 797       f3dHistPrimRawPtVsYVsVertexZAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsVertexZAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs VertexZiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
 798       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZAntiLambda);
 799    }
 800    if(! f3dHistPrimRawPtVsYVsVertexZK0Short) {
 801       f3dHistPrimRawPtVsYVsVertexZK0Short = new TH3F( "f3dHistPrimRawPtVsYVsVertexZK0Short", "Pt_{K0S} Vs Y_{K0S} Vs VertexZiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; VertexZ", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,40,-10,10);
 802       fListHistV0->Add(f3dHistPrimRawPtVsYVsVertexZK0Short);
 803    }
 804
 805 //--- 3D Histo (Pt, Y, Multiplicity), close to PV criterion
 806
 807    if(! f3dHistPrimCloseToPVPtVsYVsMultLambda) {
 808       f3dHistPrimCloseToPVPtVsYVsMultLambda = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 809       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultLambda);
 810    }
 811    if(! f3dHistPrimCloseToPVPtVsYVsMultAntiLambda) {
 812       f3dHistPrimCloseToPVPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 813       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultAntiLambda);
 814    }
 815    if(! f3dHistPrimCloseToPVPtVsYVsMultK0Short) {
 816       f3dHistPrimCloseToPVPtVsYVsMultK0Short = new TH3F( "f3dHistPrimCloseToPVPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 817       fListHistV0->Add(f3dHistPrimCloseToPVPtVsYVsMultK0Short);
 818    }
 819
 820
 821 //--- 3D Histo (Pt, Y, Proper Decay Length)
 822
 823    if(! f3dHistPrimRawPtVsYVsDecayLengthLambda) {
 824       f3dHistPrimRawPtVsYVsDecayLengthLambda = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs DecayLength; Pt_{lambda} (GeV/c); Y_{#Lambda} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
 825       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthLambda);
 826    }
 827    if(! f3dHistPrimRawPtVsYVsDecayLengthAntiLambda) {
 828       f3dHistPrimRawPtVsYVsDecayLengthAntiLambda = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs DecayLength; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
 829       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthAntiLambda);
 830    }
 831    if(! f3dHistPrimRawPtVsYVsDecayLengthK0Short) {
 832       f3dHistPrimRawPtVsYVsDecayLengthK0Short = new TH3F( "f3dHistPrimRawPtVsYVsDecayLengthK0Short", "Pt_{K0S} Vs Y_{K0S} Vs DecayLength; Pt_{K0S} (GeV/c); Y_{K0S} ; DecayLength", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,200,0,50);
 833       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthK0Short);
 834    }
 835
 836 //--------------------------------------------------------------------------------------
 837 //--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, all generated
 838
 839    if(! f3dHistGenPtVsYVsMultXiMinus) {
 840       f3dHistGenPtVsYVsMultXiMinus = new TH3F( "f3dHistGenPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 841       fListHistV0->Add(f3dHistGenPtVsYVsMultXiMinus);
 842    }
 843    if(! f3dHistGenPtVsYVsMultXiPlus) {
 844       f3dHistGenPtVsYVsMultXiPlus = new TH3F( "f3dHistGenPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 845       fListHistV0->Add(f3dHistGenPtVsYVsMultXiPlus);
 846    }
 847 //--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
 848
 849    if(! f3dHistGenPtVsYVsMultOmegaMinus) {
 850       f3dHistGenPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 851       fListHistV0->Add(f3dHistGenPtVsYVsMultOmegaMinus);
 852    }
 853    if(! f3dHistGenPtVsYVsMultOmegaPlus) {
 854       f3dHistGenPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 855       fListHistV0->Add(f3dHistGenPtVsYVsMultOmegaPlus);
 856    }
 857
 858 //--------------------------------------------------------------------------------------
 859 //--- 3D Histo (Pt, Y, Multiplicity) for generated XiMinus/Plus, at selected analysis evts
 860
 861    if(! f3dHistGenSelectedPtVsYVsMultXiMinus) {
 862       f3dHistGenSelectedPtVsYVsMultXiMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiMinus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 863       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultXiMinus);
 864    }
 865    if(! f3dHistGenSelectedPtVsYVsMultXiPlus) {
 866       f3dHistGenSelectedPtVsYVsMultXiPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultXiPlus", "Pt_{#Xi} Vs Y_{#Xi} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Xi} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 867       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultXiPlus);
 868    }
 869 //--- 3D Histo (Pt, Y, Multiplicity) for generated OmegaMinus/Plus
 870
 871    if(! f3dHistGenSelectedPtVsYVsMultOmegaMinus) {
 872       f3dHistGenSelectedPtVsYVsMultOmegaMinus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaMinus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 873       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultOmegaMinus);
 874    }
 875    if(! f3dHistGenSelectedPtVsYVsMultOmegaPlus) {
 876       f3dHistGenSelectedPtVsYVsMultOmegaPlus = new TH3F( "f3dHistGenSelectedPtVsYVsMultOmegaPlus", "Pt_{#Omega} Vs Y_{#Omega} Vs Multiplicity; Pt_{cascade} (GeV/c); Y_{#Omega} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 877       fListHistV0->Add(f3dHistGenSelectedPtVsYVsMultOmegaPlus);
 878    }
 879
 880
 881 //----------------------------------
 882 // Histos at analysis level
 883 //----------------------------------
 884
 885    if(! f3dHistPrimAnalysisPtVsYVsMultLambda) {
 886       f3dHistPrimAnalysisPtVsYVsMultLambda = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultLambda", "Pt_{lambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{lambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 887       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultLambda);
 888    }
 889    if(! f3dHistPrimAnalysisPtVsYVsMultAntiLambda) {
 890       f3dHistPrimAnalysisPtVsYVsMultAntiLambda = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultAntiLambda", "Pt_{antilambda} Vs Y_{#Lambda} Vs Multiplicity; Pt_{antilambda} (GeV/c); Y_{#Lambda} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 891       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultAntiLambda);
 892    }
 893    if(! f3dHistPrimAnalysisPtVsYVsMultK0Short) {
 894       f3dHistPrimAnalysisPtVsYVsMultK0Short = new TH3F( "f3dHistPrimAnalysisPtVsYVsMultK0Short", "Pt_{K0S} Vs Y_{K0S} Vs Multiplicity; Pt_{K0S} (GeV/c); Y_{K0S} ; Mult", lCustomNBins, 0., lCustomPtUpperLimit, 48, -1.2,1.2,lCustomNBinsMultiplicity,0,lCustomNBinsMultiplicity);
 895       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultK0Short);
 896    }
 897
 898 //----------------------------------
 899 // Primary Vertex Position Histos
 900 //----------------------------------
 901
 902    if(! fHistPVx) {
 903          fHistPVx = new TH1F("fHistPVx",
 904             "PV x position;Nbr of Evts;x",
 905             2000, -0.5, 0.5);
 906       fListHistV0->Add(fHistPVx);
 907    }
 908    if(! fHistPVy) {
 909          fHistPVy = new TH1F("fHistPVy",
 910             "PV y position;Nbr of Evts;y",
 911             2000, -0.5, 0.5);
 912       fListHistV0->Add(fHistPVy);
 913    }
 914    if(! fHistPVz) {
 915          fHistPVz = new TH1F("fHistPVz",
 916             "PV z position;Nbr of Evts;z",
 917             400, -20, 20);
 918       fListHistV0->Add(fHistPVz);
 919    }
 920
 921    if(! fHistPVxAnalysis) {
 922          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
 923             "PV x position;Nbr of Evts;x",
 924             2000, -0.5, 0.5);
 925       fListHistV0->Add(fHistPVxAnalysis);
 926    }
 927    if(! fHistPVyAnalysis) {
 928          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
 929             "PV y position;Nbr of Evts;y",
 930             2000, -0.5, 0.5);
 931       fListHistV0->Add(fHistPVyAnalysis);
 932    }
 933    if(! fHistPVzAnalysis) {
 934          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
 935             "PV z position;Nbr of Evts;z",
 936             400, -20, 20);
 937       fListHistV0->Add(fHistPVzAnalysis);
 938    }
 939
 940    if(! fHistPVxAnalysisHasHighPtLambda) {
 941          fHistPVxAnalysisHasHighPtLambda = new TH1F("fHistPVxAnalysisHasHighPtLambda",
 942             "PV x position;Nbr of Evts;x",
 943             2000, -0.5, 0.5);
 944       fListHistV0->Add(fHistPVxAnalysisHasHighPtLambda);
 945    }
 946    if(! fHistPVyAnalysisHasHighPtLambda) {
 947          fHistPVyAnalysisHasHighPtLambda = new TH1F("fHistPVyAnalysisHasHighPtLambda",
 948             "PV y position;Nbr of Evts;y",
 949             2000, -0.5, 0.5);
 950       fListHistV0->Add(fHistPVyAnalysisHasHighPtLambda);
 951    }
 952    if(! fHistPVzAnalysisHasHighPtLambda) {
 953          fHistPVzAnalysisHasHighPtLambda = new TH1F("fHistPVzAnalysisHasHighPtLambda",
 954             "PV z position;Nbr of Evts;z",
 955             400, -20, 20);
 956       fListHistV0->Add(fHistPVzAnalysisHasHighPtLambda);
 957    }
 958    if(! fHistSwappedV0Counter) {
 959       fHistSwappedV0Counter = new TH1F("fHistSwappedV0Counter",
 960          "Swap or not histo;Swapped (1) or not (0); count",
 961          2, 0, 2);
 962       fListHistV0->Add(fHistSwappedV0Counter);
 963    }
 964
 965    //List of Histograms: Normal
 966    PostData(1, fListHistV0);
 967
 968    //TTree Object: Saved to base directory. Should cache to disk while saving.
 969    //(Important to avoid excessive memory usage, particularly when merging)
 970    PostData(2, fTree);
 971
 972 }// end UserCreateOutputObjects
 973
 974
 975 //________________________________________________________________________
 976 void AliAnalysisTaskExtractPerformanceV0::UserExec(Option_t *)
 977 {
 978   // Main loop
 979   // Called for each event
 980
 981    AliESDEvent *lESDevent = 0x0;
 982    AliMCEvent  *lMCevent  = 0x0;
 983    AliStack    *lMCstack  = 0x0;
 984
 985    Int_t    lNumberOfV0s                      = -1;
 986    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 987    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 988    Double_t lMagneticField                 = -10.;
 989
 990   // Connect to the InputEvent
 991   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
 992
 993    // Appropriate for ESD analysis!
 994
 995    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 996    if (!lESDevent) {
 997       AliWarning("ERROR: lESDevent not available \n");
 998       return;
 999    }
1000
1001    fTreeVariableRunNumber = lESDevent->GetRunNumber();
1002    fTreeVariableEventNumber =
1003     ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
1004       ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
1005         ((ULong64_t)lESDevent->GetBunchCrossNumber() )  );
1006
1007    lMCevent = MCEvent();
1008    if (!lMCevent) {
1009       Printf("ERROR: Could not retrieve MC event \n");
1010       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1011       return;
1012    }
1013
1014    lMCstack = lMCevent->Stack();
1015    if (!lMCstack) {
1016       Printf("ERROR: Could not retrieve MC stack \n");
1017       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
1018       return;
1019    }
1020    TArrayF mcPrimaryVtx;
1021    AliGenEventHeader* mcHeader=lMCevent->GenEventHeader();
1022    if(!mcHeader) return;
1023    mcHeader->PrimaryVertex(mcPrimaryVtx);
1024
1025 //------------------------------------------------
1026 // Multiplicity Information Acquistion
1027 //------------------------------------------------
1028
1029    //REVISED multiplicity estimator after 'multiplicity day' (2011)
1030    Int_t lMultiplicity = -100;
1031
1032    //testing purposes
1033    if(fkIsNuclear == kFALSE) lMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,0.5);
1034
1035    //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
1036    //---> Warning: Experimental
1037    if(fkIsNuclear == kTRUE){
1038       AliCentrality* centrality;
1039       centrality = lESDevent->GetCentrality();
1040       lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( "V0M" ) ) );
1041       if (centrality->GetQuality()>1) {
1042         PostData(1, fListHistV0);
1043         PostData(2, fTree);
1044         return;
1045       }
1046    }
1047
1048    //Set variable for filling tree afterwards!
1049    //---> pp case......: GetReferenceMultiplicity
1050    //---> Pb-Pb case...: Centrality by V0M
1051    fTreeVariableMultiplicity = lMultiplicity;
1052
1053    fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
1054    fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
1055
1056 //------------------------------------------------
1057 // MC Information Acquistion
1058 //------------------------------------------------
1059
1060    Int_t iNumberOfPrimaries = -1;
1061    iNumberOfPrimaries = lMCstack->GetNprimary();
1062    if(iNumberOfPrimaries < 1) return;
1063    Bool_t lHasHighPtLambda = kFALSE;
1064
1065 //------------------------------------------------
1066 // Variable Definition
1067 //------------------------------------------------
1068
1069    Int_t lNbMCPrimary        = 0;
1070
1071    Int_t lPdgcodeCurrentPart = 0;
1072    Double_t lRapCurrentPart  = 0;
1073    Double_t lPtCurrentPart   = 0;
1074
1075    //Int_t lComeFromSigma      = 0;
1076
1077    // current mc particle 's mother
1078    //Int_t iCurrentMother  = 0;
1079    lNbMCPrimary = lMCstack->GetNprimary();
1080
1081 //------------------------------------------------
1082 // Pre-Physics Selection
1083 //------------------------------------------------
1084
1085 //----- Loop on primary Xi, Omega --------------------------------------------------------------
1086    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++)
1087    {// This is the begining of the loop on primaries
1088
1089       TParticle* lCurrentParticlePrimary = 0x0;
1090       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
1091       if(!lCurrentParticlePrimary){
1092          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1093          continue;
1094       }
1095       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) {
1096          Double_t lRapXiMCPrimary = -100;
1097          if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) {
1098            if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
1099              lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
1100            }
1101          }
1102
1103          //=================================================================================
1104          // Xi Histograms
1105          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){
1106             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1107             f3dHistGenPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1108          }
1109          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){
1110             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1111             f3dHistGenPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1112          }
1113          // Omega Histograms
1114          if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){
1115             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1116             f3dHistGenPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1117          }
1118          if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){
1119             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1120             f3dHistGenPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1121          }
1122       }
1123    }
1124 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
1125
1126 //--------- GENERATED NUMBER OF CHARGED PARTICLES
1127 // ---> Set Variables to Zero again
1128 // ---> Variable Definition
1129
1130   Long_t lNumberOfCharged = 0;
1131
1132 //----- Loop on Stack ----------------------------------------------------------------
1133    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++)
1134    {// This is the begining of the loop on tracks
1135       TParticle* particleOne = lMCstack->Particle(iCurrentLabelStack);
1136       if(!particleOne) continue;
1137       if(!particleOne->GetPDG()) continue;
1138       Double_t lThisCharge = particleOne->GetPDG()->Charge()/3.;
1139       if(TMath::Abs(lThisCharge)<0.001) continue;
1140       if(! (lMCstack->IsPhysicalPrimary(iCurrentLabelStack)) ) continue;
1141
1142       Double_t gpt = particleOne -> Pt();
1143       Double_t geta = particleOne -> Eta();
1144
1145       if( TMath::Abs(geta) < 0.5) lNumberOfCharged++;
1146    }//End of loop on tracks
1147 //----- End Loop on Stack ------------------------------------------------------------
1148
1149
1150 //----- Loop on Lambda, K0Short ----------------------------------------------------------------
1151    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++)
1152    {// This is the begining of the loop on tracks
1153
1154       TParticle* lCurrentParticleForLambdaCheck = 0x0;
1155       lCurrentParticleForLambdaCheck = lMCstack->Particle( iCurrentLabelStack );
1156       if(!lCurrentParticleForLambdaCheck){
1157          Printf("V0s loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1158          continue;
1159       }
1160
1161       //=================================================================================
1162       //Single-Strange checks
1163       // Keep only K0s, Lambda and AntiLambda:
1164       lPdgcodeCurrentPart = lCurrentParticleForLambdaCheck->GetPdgCode();
1165
1166       if ( (lCurrentParticleForLambdaCheck->GetPdgCode() == 310   ) ||
1167            (lCurrentParticleForLambdaCheck->GetPdgCode() == 3122  ) ||
1168            (lCurrentParticleForLambdaCheck->GetPdgCode() == -3122 ) )
1169            {
1170          lRapCurrentPart   = MyRapidity(lCurrentParticleForLambdaCheck->Energy(),lCurrentParticleForLambdaCheck->Pz());
1171          lPtCurrentPart    = lCurrentParticleForLambdaCheck->Pt();
1172
1173           //Use Close to PV for filling CloseToPV histograms!
1174          Double_t dx, dy, dz;
1175
1176          dx = ( (mcPrimaryVtx.At(0)) - (lCurrentParticleForLambdaCheck->Vx()) );
1177          dy = ( (mcPrimaryVtx.At(1)) - (lCurrentParticleForLambdaCheck->Vy()) );
1178          dz = ( (mcPrimaryVtx.At(2)) - (lCurrentParticleForLambdaCheck->Vz()) );
1179          Double_t lDistToPV = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
1180          if( lDistToPV <= 0.001){
1181            if( lPdgcodeCurrentPart == 3122 ){
1182               f3dHistPrimCloseToPVPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1183            }
1184            if( lPdgcodeCurrentPart == -3122 ){
1185               f3dHistPrimCloseToPVPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1186            }
1187            if( lPdgcodeCurrentPart == 310 ){
1188               f3dHistPrimCloseToPVPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1189            }
1190          }
1191
1192          //Use Physical Primaries only for filling PrimRaw Histograms!
1193          if ( lMCstack->IsPhysicalPrimary(iCurrentLabelStack)!=kTRUE ) continue;
1194
1195          if( lPdgcodeCurrentPart == 3122 ){
1196             f3dHistPrimRawPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1197             f3dHistPrimRawPtVsYVsMultMCLambda->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1198             f3dHistPrimRawPtVsYVsVertexZLambda->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1199             if( TMath::Abs( lCurrentParticleForLambdaCheck->Eta() )<1.2 && lPtCurrentPart>2 ){
1200                lHasHighPtLambda = kTRUE; //Keep track of events with Lambda within |eta|<1.2 and pt>2
1201             }
1202          }
1203          if( lPdgcodeCurrentPart == -3122 ){
1204             f3dHistPrimRawPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1205             f3dHistPrimRawPtVsYVsMultMCAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1206             f3dHistPrimRawPtVsYVsVertexZAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1207          }
1208          if( lPdgcodeCurrentPart == 310 ){
1209             f3dHistPrimRawPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1210             f3dHistPrimRawPtVsYVsMultMCK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lNumberOfCharged);
1211             f3dHistPrimRawPtVsYVsVertexZK0Short->Fill(lPtCurrentPart, lRapCurrentPart, mcPrimaryVtx.At(2));
1212          }
1213          //Decay Length Acquisition=====================================================
1214          Double_t decaylength = -1;
1215          Double_t lV0Mass = -1;
1216
1217          if( !(lCurrentParticleForLambdaCheck->GetDaughter(0) < 0) ) {
1218             TParticle* lDght0ofV0 = lMCstack->Particle(  lCurrentParticleForLambdaCheck->GetDaughter(0) ); //get first daughter
1219             if(lDght0ofV0){ // skip if not defined.
1220                decaylength = TMath::Sqrt(
1221                                         TMath::Power( lCurrentParticleForLambdaCheck->Vx() - lDght0ofV0->Vx() , 2) +
1222                                         TMath::Power( lCurrentParticleForLambdaCheck->Vy() - lDght0ofV0->Vy() , 2) +
1223                                         TMath::Power( lCurrentParticleForLambdaCheck->Vz() - lDght0ofV0->Vz() , 2)
1224                );
1225                //Need to correct for relativitity! Involves multiplying by mass and dividing by momentum.
1226                if(TMath::Abs( lPdgcodeCurrentPart ) == 3122 ) { lV0Mass = 1.115683; }
1227                if(TMath::Abs( lPdgcodeCurrentPart ) == 310 ) { lV0Mass = 0.497614; }
1228                if( lCurrentParticleForLambdaCheck->P() + 1e-10 != 0 ) decaylength = ( lV0Mass * decaylength ) / ( lCurrentParticleForLambdaCheck->P() + 1e-10 );
1229                if( lCurrentParticleForLambdaCheck->P() + 1e-10 == 0 ) decaylength = 1e+5;
1230             }
1231          }
1232          if( lPdgcodeCurrentPart == 3122) f3dHistPrimRawPtVsYVsDecayLengthLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
1233          if( lPdgcodeCurrentPart == -3122) f3dHistPrimRawPtVsYVsDecayLengthAntiLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
1234          if( lPdgcodeCurrentPart == 310) f3dHistPrimRawPtVsYVsDecayLengthK0Short ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
1235       }
1236    }//End of loop on tracks
1237 //----- End Loop on Lambda, K0Short ------------------------------------------------------------
1238
1239
1240    f2dHistMultiplicityVsTrueBeforeTrigSel->Fill ( lMultiplicity , lNumberOfCharged );
1241
1242     fTreeVariableMultiplicityMC = lNumberOfCharged;
1243
1244    fHistGenVertexZBeforeTrigSel->Fill( (mcPrimaryVtx.At(2)) );
1245
1246    lPdgcodeCurrentPart = 0;
1247    lRapCurrentPart  = 0;
1248    lPtCurrentPart   = 0;
1249
1250 //------------------------------------------------
1251 // Physics Selection
1252 //------------------------------------------------
1253
1254    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1255    Bool_t isSelected = 0;
1256    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
1257
1258    //pp at 2.76TeV: special case, ignore FastOnly
1259    if ( (fkLowEnergyPP == kTRUE) && (maskIsSelected& AliVEvent::kFastOnly) ){
1260       PostData(1, fListHistV0);
1261       PostData(2, fTree);
1262       return;
1263    }
1264    //Standard Min-Bias Selection
1265    if ( ! isSelected ) {
1266       PostData(1, fListHistV0);
1267       PostData(2, fTree);
1268       return;
1269    }
1270
1271    f2dHistMultiplicityVsTrueForTrigEvt->Fill ( lMultiplicity , lNumberOfCharged );
1272    fHistGenVertexZForTrigEvt->Fill( mcPrimaryVtx.At(2) );
1273 //------------------------------------------------
1274 // After Trigger Selection
1275 //------------------------------------------------
1276
1277    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
1278
1279    //Set variable for filling tree afterwards!
1280    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
1281    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
1282
1283 //------------------------------------------------
1284 // Getting: Primary Vertex + MagField Info
1285 //------------------------------------------------
1286
1287    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
1288    // get the vtx stored in ESD found with tracks
1289    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
1290
1291    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
1292    // get the best primary vertex available for the event
1293    // As done in AliCascadeVertexer, we keep the one which is the best one available.
1294    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
1295    // This one will be used for next calculations (DCA essentially)
1296    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
1297
1298    Double_t lPrimaryVtxPosition[3];
1299    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
1300    lPrimaryVtxPosition[0] = primaryVtx->GetX();
1301    lPrimaryVtxPosition[1] = primaryVtx->GetY();
1302    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
1303    fHistPVx->Fill( lPrimaryVtxPosition[0] );
1304    fHistPVy->Fill( lPrimaryVtxPosition[1] );
1305    fHistPVz->Fill( lPrimaryVtxPosition[2] );
1306
1307    f2dHistMultiplicityVsVertexZForTrigEvt->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1308
1309 //------------------------------------------------
1310 // Primary Vertex Z position: SKIP
1311 //------------------------------------------------
1312
1313    if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
1314       AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
1315       PostData(1, fListHistV0);
1316       PostData(2, fTree);
1317       return;
1318    }
1319
1320    f2dHistMultiplicityVsVertexZ->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1321
1322    lMagneticField = lESDevent->GetMagneticField( );
1323    fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
1324    fHistMultiplicity->Fill(lMultiplicity);
1325    f2dHistMultiplicityVsTrue->Fill ( lMultiplicity , lNumberOfCharged );
1326    fHistGenVertexZ->Fill( (mcPrimaryVtx.At(2)) );
1327 //------------------------------------------------
1328 // SKIP: Events with well-established PVtx
1329 //------------------------------------------------
1330
1331    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
1332    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
1333    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
1334       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
1335       PostData(1, fListHistV0);
1336       PostData(2, fTree);
1337       return;
1338    }
1339
1340    f2dHistMultiplicityVsVertexZNoTPCOnly->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1341    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
1342    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
1343    f2dHistMultiplicityVsTrueNoTPCOnly->Fill ( lMultiplicity , lNumberOfCharged );
1344    fHistGenVertexZNoTPCOnly->Fill( (mcPrimaryVtx.At(2)) );
1345 //------------------------------------------------
1346 // Pileup Rejection Studies
1347 //------------------------------------------------
1348
1349    // FIXME : quality selection regarding pile-up rejection
1350    if(lESDevent->IsPileupFromSPD() && !fkIsNuclear ){// minContributors=3, minZdist=0.8, nSigmaZdist=3., nSigmaDiamXY=2., nSigmaDiamZ=5.  -> see http://alisoft.cern.ch/viewvc/trunk/STEER/AliESDEvent.h?root=AliRoot&r1=41914&r2=42199&pathrev=42199
1351       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
1352       PostData(1, fListHistV0);
1353       PostData(2, fTree);
1354       return;
1355    }
1356    f2dHistMultiplicityVsVertexZNoTPCOnlyNoPileup->Fill( lMultiplicity, lPrimaryVtxPosition[2] );
1357    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
1358    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
1359    f2dHistMultiplicityVsTrueNoTPCOnlyNoPileup->Fill ( lMultiplicity , lNumberOfCharged );
1360    fHistGenVertexZNoTPCOnlyNoPileup->Fill( (mcPrimaryVtx.At(2)) );
1361    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
1362    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){
1363       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
1364       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
1365       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
1366       if ( lHasHighPtLambda == kTRUE ){
1367          fHistPVxAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[0] );
1368          fHistPVyAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[1] );
1369          fHistPVzAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[2] );
1370       }
1371    }
1372
1373   fTreeVariableVertexZ = lPrimaryVtxPosition[2];
1374
1375   fTreeVariablePVx = lPrimaryVtxPosition[0];
1376   fTreeVariablePVy = lPrimaryVtxPosition[1];
1377   fTreeVariablePVz = lPrimaryVtxPosition[2];
1378
1379   fTreeVariableMCPVx = (mcPrimaryVtx.At(0));
1380   fTreeVariableMCPVy = (mcPrimaryVtx.At(1));
1381   fTreeVariableMCPVz = (mcPrimaryVtx.At(2));
1382
1383 //------------------------------------------------
1384 // stack loop starts here
1385 //------------------------------------------------
1386
1387 //---> Loop over ALL PARTICLES
1388
1389    for (Int_t iMc = 0; iMc < (lMCstack->GetNtrack()); iMc++) {
1390       TParticle *p0 = lMCstack->Particle(iMc);
1391       if (!p0) {
1392          //Printf("ERROR: particle with label %d not found in lMCstack (mc loop)", iMc);
1393          continue;
1394       }
1395       lPdgcodeCurrentPart = p0->GetPdgCode();
1396
1397       // Keep only K0s, Lambda and AntiLambda:
1398       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
1399
1400       lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
1401       lPtCurrentPart    = p0->Pt();
1402
1403         //Use Physical Primaries only for filling PrimRaw Histograms!
1404       if ( lMCstack->IsPhysicalPrimary(iMc)!=kTRUE ) continue;
1405
1406       if( lPdgcodeCurrentPart == 3122 ){
1407          f3dHistPrimAnalysisPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1408       }
1409       if( lPdgcodeCurrentPart == -3122 ){
1410          f3dHistPrimAnalysisPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1411       }
1412       if( lPdgcodeCurrentPart == 310 ){
1413          f3dHistPrimAnalysisPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
1414       }
1415    }
1416
1417 //----- Loop on primary Xi, Omega --------------------------------------------------------------
1418    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++)
1419    {// This is the begining of the loop on primaries
1420
1421       TParticle* lCurrentParticlePrimary = 0x0;
1422       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
1423       if(!lCurrentParticlePrimary){
1424          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
1425          continue;
1426       }
1427       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 || TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3334 ) {
1428          Double_t lRapXiMCPrimary = -100;
1429          if( (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) != 0 ) {
1430            if ( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) !=0 ){
1431              lRapXiMCPrimary = 0.5*TMath::Log( (lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13) );
1432            }
1433          }
1434
1435          //=================================================================================
1436          // Xi Histograms
1437          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){
1438             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1439             f3dHistGenSelectedPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1440          }
1441          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){
1442             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1443             f3dHistGenSelectedPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1444          }
1445          // Omega Histograms
1446          if( lCurrentParticlePrimary->GetPdgCode() == 3334 ){
1447             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1448             f3dHistGenSelectedPtVsYVsMultOmegaMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1449          }
1450          if( lCurrentParticlePrimary->GetPdgCode() == -3334 ){
1451             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
1452             f3dHistGenSelectedPtVsYVsMultOmegaPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
1453          }
1454       }
1455    }
1456 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
1457
1458 //------------------------------------------------
1459 // MAIN LAMBDA LOOP STARTS HERE
1460 //------------------------------------------------
1461
1462    //Variable definition
1463    Int_t    lOnFlyStatus = 0;
1464    Double_t lChi2V0 = 0;
1465    Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
1466    Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
1467    Double_t lV0CosineOfPointingAngle = 0;
1468    Double_t lV0Radius = 0, lPt = 0;
1469    Double_t lRapK0Short = 0, lRapLambda = 0;
1470    Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
1471    Double_t lAlphaV0 = 0, lPtArmV0 = 0;
1472    Double_t fMinV0Pt = 0;
1473    Double_t fMaxV0Pt = 100;
1474
1475    Int_t nv0s = 0;
1476    nv0s = lESDevent->GetNumberOfV0s();
1477
1478    for (Int_t iV0 = 0; iV0 < nv0s; iV0++)
1479         {// This is the begining of the V0 loop
1480       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
1481       if (!v0) continue;
1482
1483       //---> Fix On-the-Fly candidates
1484       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
1485         fHistSwappedV0Counter -> Fill( 1 );
1486       }else{
1487         fHistSwappedV0Counter -> Fill( 0 );
1488       }
1489       if ( fkUseOnTheFly ) CheckChargeV0(v0);
1490
1491
1492       Double_t tV0mom[3];
1493       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] );
1494       Double_t lV0TotalMomentum = TMath::Sqrt(
1495          tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
1496
1497       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]);
1498       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
1499       lPt = v0->Pt();
1500       lRapK0Short = v0->RapK0Short();
1501       lRapLambda  = v0->RapLambda();
1502
1503       //Set Variables for later filling
1504       fTreeVariableV0x = tDecayVertexV0[0];
1505       fTreeVariableV0y = tDecayVertexV0[1];
1506       fTreeVariableV0z = tDecayVertexV0[2];
1507
1508       //Set Variables for later filling
1509       fTreeVariableV0Px = tV0mom[0];
1510       fTreeVariableV0Py = tV0mom[1];
1511       fTreeVariableV0Pz = tV0mom[2];
1512
1513       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
1514
1515       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
1516       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
1517
1518       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
1519       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
1520
1521       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
1522       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
1523       if (!pTrack || !nTrack) {
1524          Printf("ERROR: Could not retreive one of the daughter track");
1525          continue;
1526       }
1527
1528       fTreeVariableNegEta = nTrack->Eta();
1529       fTreeVariablePosEta = pTrack->Eta();
1530
1531       // Filter like-sign V0 (next: add counter and distribution)
1532       if ( pTrack->GetSign() == nTrack->GetSign()){
1533          continue;
1534       }
1535
1536       //________________________________________________________________________
1537       // Track quality cuts
1538       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
1539       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
1540       fTreeVariableLeastNbrCrossedRows = (Int_t) lPosTrackCrossedRows;
1541       if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
1542          fTreeVariableLeastNbrCrossedRows = (Int_t) lNegTrackCrossedRows;
1543
1544       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
1545       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
1546       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
1547
1548       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
1549
1550       //GetKinkIndex condition
1551       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
1552
1553       //Findable clusters > 0 condition
1554       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
1555
1556       //Compute ratio Crossed Rows / Findable clusters
1557       //Note: above test avoids division by zero!
1558       Float_t lPosTrackCrossedRowsOverFindable = -1;
1559       Float_t lNegTrackCrossedRowsOverFindable = -1;
1560       if ( ((double)(pTrack->GetTPCNclsF()) ) != 0 ) lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF()));
1561       if ( ((double)(nTrack->GetTPCNclsF()) ) != 0 ) lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF()));
1562
1563       fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
1564       if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
1565          fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
1566
1567       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
1568       if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
1569
1570       //End track Quality Cuts
1571       //________________________________________________________________________
1572
1573       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(lPrimaryVtxPosition[0],
1574                                                         lPrimaryVtxPosition[1],
1575                                                         lMagneticField) );
1576
1577       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(lPrimaryVtxPosition[0],
1578                                                         lPrimaryVtxPosition[1],
1579                                                         lMagneticField) );
1580
1581       lOnFlyStatus = v0->GetOnFlyStatus();
1582       lChi2V0 = v0->GetChi2V0();
1583       lDcaV0Daughters = v0->GetDcaV0Daughters();
1584       lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1585       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
1586       fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
1587
1588       // Getting invariant mass infos directly from ESD
1589       v0->ChangeMassHypothesis(310);
1590       lInvMassK0s = v0->GetEffMass();
1591       v0->ChangeMassHypothesis(3122);
1592       lInvMassLambda = v0->GetEffMass();
1593       v0->ChangeMassHypothesis(-3122);
1594       lInvMassAntiLambda = v0->GetEffMass();
1595       lAlphaV0 = v0->AlphaV0();
1596       lPtArmV0 = v0->PtArmV0();
1597
1598       //fTreeVariableOnFlyStatus = lOnFlyStatus;
1599       //fHistV0OnFlyStatus->Fill(lOnFlyStatus);
1600
1601 //===============================================
1602 // Monte Carlo Association starts here
1603 //===============================================
1604
1605       //---> Set Everything to "I don't know" before starting
1606
1607       fTreeVariablePIDPositive = 0;
1608       fTreeVariablePIDNegative = 0;
1609
1610       fTreeVariableIndexStatus = 0;
1611       fTreeVariableIndexStatusMother = 0;
1612
1613       fTreeVariablePtMother = -1;
1614       fTreeVariablePtMC = -1;
1615       fTreeVariableRapMC = -100;
1616
1617       fTreeVariablePID = -1;
1618       fTreeVariablePIDMother = -1;
1619
1620       fTreeVariablePrimaryStatus = 0;
1621       fTreeVariablePrimaryStatusMother = 0;
1622       fTreeVariableV0CreationRadius = -1;
1623
1624       Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrack->GetLabel() );
1625       Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrack->GetLabel() );
1626
1627       TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
1628       TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
1629
1630       fTreeVariablePosTransvMomentumMC = mcPosV0Dghter->Pt();
1631       fTreeVariableNegTransvMomentumMC = mcNegV0Dghter->Pt();
1632
1633       Int_t lPIDPositive = mcPosV0Dghter -> GetPdgCode();
1634       Int_t lPIDNegative = mcNegV0Dghter -> GetPdgCode();
1635
1636       fTreeVariablePIDPositive = lPIDPositive;
1637       fTreeVariablePIDNegative = lPIDNegative;
1638
1639       Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
1640       Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
1641
1642       if( lblMotherPosV0Dghter == lblMotherNegV0Dghter && lblMotherPosV0Dghter > -1 ){
1643          //either label is fine, they're equal at this stage
1644          TParticle* pThisV0 = lMCstack->Particle( lblMotherPosV0Dghter );
1645          //Set tree variables
1646          fTreeVariablePID   = pThisV0->GetPdgCode(); //PDG Code
1647          fTreeVariablePtMC  = pThisV0->Pt(); //Perfect Pt
1648
1649          //Set Variables for later filling
1650          //Be careful: Vx, Vy, Vz: Creation vertex. So decay position is the
1651          //Creation vertex of any one of the daughters!
1652          fTreeVariableMCV0x = mcPosV0Dghter->Vx();
1653          fTreeVariableMCV0y = mcPosV0Dghter->Vy();
1654          fTreeVariableMCV0z = mcPosV0Dghter->Vz();
1655
1656          //Set Variables for later filling
1657          fTreeVariableMCV0Px = pThisV0->Px();
1658          fTreeVariableMCV0Py = pThisV0->Py();
1659          fTreeVariableMCV0Pz = pThisV0->Pz();
1660
1661          //Only Interested if it's a Lambda, AntiLambda or K0s
1662          //Avoid the Junction Bug! PYTHIA has particles with Px=Py=Pz=E=0 occasionally,
1663          //having particle code 88 (unrecognized by PDG), for documentation purposes.
1664          //Even ROOT's TParticle::Y() is not prepared to deal with that exception!
1665          //Note that TParticle::Pt() is immune (that would just return 0)...
1666          //Though granted that that should be extremely rare in this precise condition...
1667          if( TMath::Abs(fTreeVariablePID) == 3122 || fTreeVariablePID==310 ){
1668             fTreeVariableRapMC = pThisV0->Y(); //Perfect Y
1669          }
1670          fTreeVariableV0CreationRadius = TMath::Sqrt(
1671           TMath::Power(  ( (mcPrimaryVtx.At(0)) - (pThisV0->Vx()) ) , 2) +
1672           TMath::Power(  ( (mcPrimaryVtx.At(1)) - (pThisV0->Vy()) ) , 2) +
1673           TMath::Power(  ( (mcPrimaryVtx.At(2)) - (pThisV0->Vz()) ) , 2)
1674          );
1675          if( lblMotherPosV0Dghter  < lNbMCPrimary ) fTreeVariableIndexStatus = 1; //looks primary
1676          if( lblMotherPosV0Dghter >= lNbMCPrimary ) fTreeVariableIndexStatus = 2; //looks secondary
1677          if( lMCstack->IsPhysicalPrimary       (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 1; //Is Primary!
1678          if( lMCstack->IsSecondaryFromWeakDecay(lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 2; //Weak Decay!
1679          if( lMCstack->IsSecondaryFromMaterial (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 3; //Material Int!
1680
1681          //Now we try to acquire the V0 parent particle, if possible
1682          Int_t lblThisV0Parent = pThisV0->GetFirstMother();
1683          if ( lblThisV0Parent > -1 ){ //if it has a parent, get it and store specs
1684             TParticle* pThisV0Parent = lMCstack->Particle( lblThisV0Parent );
1685             fTreeVariablePIDMother   = pThisV0Parent->GetPdgCode(); //V0 Mother PDG
1686             fTreeVariablePtMother    = pThisV0Parent->Pt();         //V0 Mother Pt
1687             //Primary Status for the V0 Mother particle
1688             if( lblThisV0Parent  < lNbMCPrimary ) fTreeVariableIndexStatusMother = 1; //looks primary
1689             if( lblThisV0Parent >= lNbMCPrimary ) fTreeVariableIndexStatusMother = 2; //looks secondary
1690             if( lMCstack->IsPhysicalPrimary       (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 1; //Is Primary!
1691             if( lMCstack->IsSecondaryFromWeakDecay(lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 2; //Weak Decay!
1692             if( lMCstack->IsSecondaryFromMaterial (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 3; //Material Int!
1693          }
1694       }
1695
1696       fTreeVariablePt = v0->Pt();
1697       fTreeVariableChi2V0 = lChi2V0;
1698       fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
1699       fTreeVariableDcaV0Daughters = lDcaV0Daughters;
1700       fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle;
1701       fTreeVariableV0Radius = lV0Radius;
1702       fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
1703       fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
1704       fTreeVariableInvMassK0s = lInvMassK0s;
1705       fTreeVariableInvMassLambda = lInvMassLambda;
1706       fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
1707       fTreeVariableRapK0Short = lRapK0Short;
1708
1709       fTreeVariableRapLambda = lRapLambda;
1710       fTreeVariableAlphaV0 = lAlphaV0;
1711       fTreeVariablePtArmV0 = lPtArmV0;
1712
1713       //Official means of acquiring N-sigmas
1714       fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
1715       fTreeVariableNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
1716       fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
1717       fTreeVariableNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
1718
1719 //tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]
1720       Double_t lDistanceTravelled = TMath::Sqrt(
1721                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
1722                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
1723                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
1724                                         );
1725       fTreeVariableDistOverTotMom = 1e+5;
1726       if( lV0TotalMomentum + 1e-10 != 0 ) fTreeVariableDistOverTotMom = lDistanceTravelled / (lV0TotalMomentum + 1e-10); //avoid division by zero, to be sure
1727
1728       Double_t lMomentumPosTemp[3];
1729       pTrack->GetPxPyPz(lMomentumPosTemp);
1730       Double_t lPtPosTemporary = sqrt(pow(lMomentumPosTemp[0],2) + pow(lMomentumPosTemp[1],2));
1731
1732       Double_t lMomentumNegTemp[3];
1733       nTrack->GetPxPyPz(lMomentumNegTemp);
1734       Double_t lPtNegTemporary = sqrt(pow(lMomentumNegTemp[0],2) + pow(lMomentumNegTemp[1],2));
1735
1736       fTreeVariablePosTransvMomentum = lPtPosTemporary;
1737       fTreeVariableNegTransvMomentum = lPtNegTemporary;
1738
1739
1740 //------------------------------------------------
1741 // Fill Tree!
1742 //------------------------------------------------
1743
1744       // The conditionals are meant to decrease excessive
1745       // memory usage!
1746
1747       //Modified version: Keep only OnFlyStatus == 0
1748       //Keep only if included in a parametric InvMass Region 20 sigmas away from peak
1749
1750       //First Selection: Reject OnFly
1751       if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){
1752          //Second Selection: rough 20-sigma band, parametric.
1753          //K0Short: Enough to parametrize peak broadening with linear function.
1754          Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt;
1755          Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
1756          //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
1757          //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
1758          Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt);
1759          Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
1760          //Do Selection
1761          if( (fTreeVariableInvMassLambda     < lUpperLimitLambda  && fTreeVariableInvMassLambda     > lLowerLimitLambda     ) ||
1762              (fTreeVariableInvMassAntiLambda < lUpperLimitLambda  && fTreeVariableInvMassAntiLambda > lLowerLimitLambda     ) ||
1763              (fTreeVariableInvMassK0s        < lUpperLimitK0Short && fTreeVariableInvMassK0s        > lLowerLimitK0Short    ) ){
1764              //Pre-selection in case this is AA...
1765              if( fkIsNuclear == kFALSE ) fTree->Fill();
1766              if( fkIsNuclear == kTRUE){
1767              //If this is a nuclear collision___________________
1768              // ... pre-filter with daughter eta selection only (not TPC)
1769                if ( TMath::Abs(fTreeVariableNegEta)<0.8 && TMath::Abs(fTreeVariablePosEta)<0.8 ) fTree->Fill();
1770              }//end nuclear_____________________________________
1771          }
1772       }
1773
1774 //------------------------------------------------
1775 // Fill tree over.
1776 //------------------------------------------------
1777
1778
1779    }// This is the end of the V0 loop
1780
1781 //------------------------------------------------
1782
1783    // Post output data.
1784    PostData(1, fListHistV0);
1785    PostData(2, fTree);
1786 }
1787
1788 //________________________________________________________________________
1789 void AliAnalysisTaskExtractPerformanceV0::Terminate(Option_t *)
1790 {
1791    // Draw result to the screen
1792    // Called once at the end of the query
1793
1794    TList *cRetrievedList = 0x0;
1795    cRetrievedList = (TList*)GetOutputData(1);
1796    if(!cRetrievedList){
1797       Printf("ERROR - AliAnalysisTaskExtractV0 : ouput data container list not available\n");
1798       return;
1799    }
1800
1801    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1802    if (!fHistV0MultiplicityForTrigEvt) {
1803       Printf("ERROR - AliAnalysisTaskExtractV0 : fHistV0MultiplicityForTrigEvt not available");
1804       return;
1805    }
1806
1807    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractV0","V0 Multiplicity",10,10,510,510);
1808    canCheck->cd(1)->SetLogy();
1809
1810    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1811    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1812 }
1813
1814 //----------------------------------------------------------------------------
1815
1816 Double_t AliAnalysisTaskExtractPerformanceV0::MyRapidity(Double_t rE, Double_t rPz) const
1817 {
1818    // Local calculation for rapidity
1819    Double_t ReturnValue = -100;
1820    if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){
1821       ReturnValue =  0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1822    }
1823    return ReturnValue;
1824 }
1825
1826 //________________________________________________________________________
1827 void AliAnalysisTaskExtractPerformanceV0::CheckChargeV0(AliESDv0 *v0)
1828 {
1829    // This function checks charge of negative and positive daughter tracks.
1830    // If incorrectly defined (onfly vertexer), swaps out.
1831    if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){
1832       //V0 daughter track swapping is required! Note: everything is swapped here... P->N, N->P
1833       Long_t lCorrectNidx = v0->GetPindex();
1834       Long_t lCorrectPidx = v0->GetNindex();
1835       Double32_t        lCorrectNmom[3];
1836       Double32_t        lCorrectPmom[3];
1837       v0->GetPPxPyPz( lCorrectNmom[0], lCorrectNmom[1], lCorrectNmom[2] );
1838       v0->GetNPxPyPz( lCorrectPmom[0], lCorrectPmom[1], lCorrectPmom[2] );
1839
1840       AliExternalTrackParam     lCorrectParamN(
1841         v0->GetParamP()->GetX() ,
1842         v0->GetParamP()->GetAlpha() ,
1843         v0->GetParamP()->GetParameter() ,
1844         v0->GetParamP()->GetCovariance()
1845       );
1846       AliExternalTrackParam     lCorrectParamP(
1847         v0->GetParamN()->GetX() ,
1848         v0->GetParamN()->GetAlpha() ,
1849         v0->GetParamN()->GetParameter() ,
1850         v0->GetParamN()->GetCovariance()
1851       );
1852       lCorrectParamN.SetMostProbablePt( v0->GetParamP()->GetMostProbablePt() );
1853       lCorrectParamP.SetMostProbablePt( v0->GetParamN()->GetMostProbablePt() );
1854
1855       //Get Variables___________________________________________________
1856       Double_t lDcaV0Daughters = v0 -> GetDcaV0Daughters();
1857       Double_t lCosPALocal     = v0 -> GetV0CosineOfPointingAngle();
1858       Bool_t lOnFlyStatusLocal = v0 -> GetOnFlyStatus();
1859
1860       //Create Replacement Object_______________________________________
1861       AliESDv0 *v0correct = new AliESDv0(lCorrectParamN,lCorrectNidx,lCorrectParamP,lCorrectPidx);
1862       v0correct->SetDcaV0Daughters          ( lDcaV0Daughters   );
1863       v0correct->SetV0CosineOfPointingAngle ( lCosPALocal       );
1864       v0correct->ChangeMassHypothesis       ( kK0Short          );
1865       v0correct->SetOnFlyStatus             ( lOnFlyStatusLocal );
1866
1867       //Reverse Cluster info..._________________________________________
1868       v0correct->SetClusters( v0->GetClusters( 1 ), v0->GetClusters ( 0 ) );
1869
1870       *v0 = *v0correct;
1871       //Proper cleanup..._______________________________________________
1872       v0correct->Delete();
1873       v0correct = 0x0;
1874
1875       //Just another cross-check and output_____________________________
1876       if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ) {
1877         AliWarning("Found Swapped Charges, tried to correct but something FAILED!");
1878       }else{
1879         //AliWarning("Found Swapped Charges and fixed.");
1880       }
1881       //________________________________________________________________
1882    }else{
1883       //Don't touch it! ---
1884       //Printf("Ah, nice. Charges are already ordered...");
1885    }
1886    return;
1887 }
1888