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