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
  90 #include "AliAnalysisTaskExtractPerformanceV0.h"
  91
  92 ClassImp(AliAnalysisTaskExtractPerformanceV0)
  93
  94 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0()
  95   : AliAnalysisTaskSE(), fListHistV0(0), fTree(0),
  96
  97 //------------------------------------------------
  98 // HISTOGRAMS
  99 // --- Filled on an Event-by-event basis
 100 //------------------------------------------------
 101    fHistV0MultiplicityBeforeTrigSel(0),
 102    fHistV0MultiplicityForTrigEvt(0),
 103    fHistV0MultiplicityForSelEvt(0),
 104    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 105    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 106    fHistMultiplicityBeforeTrigSel(0),
 107    fHistMultiplicityForTrigEvt(0),
 108    fHistMultiplicity(0),
 109    fHistMultiplicityNoTPCOnly(0),
 110    fHistMultiplicityNoTPCOnlyNoPileup(0),
 111
 112 //------------------------------------------------
 113 // PARTICLE HISTOGRAMS
 114 // --- Filled on a Particle-by-Particle basis
 115 //------------------------------------------------
 116    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
 117    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
 118    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
 119    f3dHistPrimRawPtVsYVsMultLambda(0),
 120    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
 121    f3dHistPrimRawPtVsYVsMultK0Short(0),
 122    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
 123    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
 124    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
 125    f3dHistGenPtVsYVsMultXiMinus(0),
 126    f3dHistGenPtVsYVsMultXiPlus(0),
 127    fHistPVx(0),
 128    fHistPVy(0),
 129    fHistPVz(0),
 130    fHistPVxAnalysis(0),
 131    fHistPVyAnalysis(0),
 132    fHistPVzAnalysis(0),
 133    fHistPVxAnalysisHasHighPtLambda(0),
 134    fHistPVyAnalysisHasHighPtLambda(0),
 135    fHistPVzAnalysisHasHighPtLambda(0)
 136 {
 137   // Dummy Constructor
 138 }
 139
 140 AliAnalysisTaskExtractPerformanceV0::AliAnalysisTaskExtractPerformanceV0(const char *name)
 141   : AliAnalysisTaskSE(name), fListHistV0(0), fTree(0),
 142
 143 //------------------------------------------------
 144 // HISTOGRAMS
 145 // --- Filled on an Event-by-event basis
 146 //------------------------------------------------
 147    fHistV0MultiplicityBeforeTrigSel(0),
 148    fHistV0MultiplicityForTrigEvt(0),
 149    fHistV0MultiplicityForSelEvt(0),
 150    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 151    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 152    fHistMultiplicityBeforeTrigSel(0),
 153    fHistMultiplicityForTrigEvt(0),
 154    fHistMultiplicity(0),
 155    fHistMultiplicityNoTPCOnly(0),
 156    fHistMultiplicityNoTPCOnlyNoPileup(0),
 157
 158
 159 //------------------------------------------------
 160 // PARTICLE HISTOGRAMS
 161 // --- Filled on a Particle-by-Particle basis
 162 //------------------------------------------------
 163    f3dHistPrimAnalysisPtVsYVsMultLambda(0),
 164    f3dHistPrimAnalysisPtVsYVsMultAntiLambda(0),
 165    f3dHistPrimAnalysisPtVsYVsMultK0Short(0),
 166    f3dHistPrimRawPtVsYVsMultLambda(0),
 167    f3dHistPrimRawPtVsYVsMultAntiLambda(0),
 168    f3dHistPrimRawPtVsYVsMultK0Short(0),
 169    f3dHistPrimRawPtVsYVsDecayLengthLambda(0),
 170    f3dHistPrimRawPtVsYVsDecayLengthAntiLambda(0),
 171    f3dHistPrimRawPtVsYVsDecayLengthK0Short(0),
 172    f3dHistGenPtVsYVsMultXiMinus(0),
 173    f3dHistGenPtVsYVsMultXiPlus(0),
 174    fHistPVx(0),
 175    fHistPVy(0),
 176    fHistPVz(0),
 177    fHistPVxAnalysis(0),
 178    fHistPVyAnalysis(0),
 179    fHistPVzAnalysis(0),
 180    fHistPVxAnalysisHasHighPtLambda(0),
 181    fHistPVyAnalysisHasHighPtLambda(0),
 182    fHistPVzAnalysisHasHighPtLambda(0)
 183 {
 184    // Constructor
 185    // Output slot #0 writes into a TList container (Cascade)
 186    DefineOutput(1, TList::Class());
 187    DefineOutput(2, TTree::Class());
 188 }
 189
 190
 191 AliAnalysisTaskExtractPerformanceV0::~AliAnalysisTaskExtractPerformanceV0()
 192 {
 193 //------------------------------------------------
 194 // DESTRUCTOR
 195 //------------------------------------------------
 196
 197    if (fListHistV0){
 198       delete fListHistV0;
 199       fListHistV0 = 0x0;
 200    }
 201    if (fTree){
 202       delete fTree;
 203       fTree = 0x0;
 204    }
 205 }
 206
 207 //________________________________________________________________________
 208 void AliAnalysisTaskExtractPerformanceV0::UserCreateOutputObjects()
 209 {
 210
 211    OpenFile(2);
 212    // Called once
 213    fTree = new TTree("fTree","V0Candidates");
 214
 215 //------------------------------------------------
 216 // fTree Branch definitions
 217 //------------------------------------------------
 218
 219 //-----------BASIC-INFO---------------------------
 220 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatus",&fTreeVariablePrimaryStatus,"fTreeVariablePrimaryStatus/I");
 221 /* 1*/   fTree->Branch("fTreeVariablePrimaryStatusMother",&fTreeVariablePrimaryStatusMother,"fTreeVariablePrimaryStatusMother/I");
 222 /* 2*/   fTree->Branch("fTreeVariableChi2V0",&fTreeVariableChi2V0,"Chi2V0/F");
 223 /* 3*/   fTree->Branch("fTreeVariableDcaV0Daughters",&fTreeVariableDcaV0Daughters,"fTreeVariableDcaV0Daughters/F");
 224 /* 4*/   fTree->Branch("fTreeVariableDcaPosToPrimVertex",&fTreeVariableDcaPosToPrimVertex,"fTreeVariableDcaPosToPrimVertex/F");
 225 /* 5*/   fTree->Branch("fTreeVariableDcaNegToPrimVertex",&fTreeVariableDcaNegToPrimVertex,"fTreeVariableDcaNegToPrimVertex/F");
 226 /* 6*/   fTree->Branch("fTreeVariableV0Radius",&fTreeVariableV0Radius,"fTreeVariableV0Radius/F");
 227 /* 7*/   fTree->Branch("fTreeVariablePt",&fTreeVariablePt,"fTreeVariablePt/F");
 228 /* 7*/   fTree->Branch("fTreeVariablePtMC",&fTreeVariablePtMC,"fTreeVariablePtMC/F");
 229 /* 8*/   fTree->Branch("fTreeVariableRapK0Short",&fTreeVariableRapK0Short,"fTreeVariableRapK0Short/F");
 230 /* 9*/   fTree->Branch("fTreeVariableRapLambda",&fTreeVariableRapLambda,"fTreeVariableRapLambda/F");
 231 /*10*/   fTree->Branch("fTreeVariableRapMC",&fTreeVariableRapMC,"fTreeVariableRapMC/F");
 232 /*11*/   fTree->Branch("fTreeVariableInvMassK0s",&fTreeVariableInvMassK0s,"fTreeVariableInvMassK0s/F");
 233 /*12*/   fTree->Branch("fTreeVariableInvMassLambda",&fTreeVariableInvMassLambda,"fTreeVariableInvMassLambda/F");
 234 /*13*/   fTree->Branch("fTreeVariableInvMassAntiLambda",&fTreeVariableInvMassAntiLambda,"fTreeVariableInvMassAntiLambda/F");
 235 /*14*/   fTree->Branch("fTreeVariableAlphaV0",&fTreeVariableAlphaV0,"fTreeVariableAlphaV0/F");
 236 /*15*/   fTree->Branch("fTreeVariablePtArmV0",&fTreeVariablePtArmV0,"fTreeVariablePtArmV0/F");
 237 /*16*/   fTree->Branch("fTreeVariableNegTransvMomentum",&fTreeVariableNegTransvMomentum,"fTreeVariableNegTransvMomentum/F");
 238 /*17*/   fTree->Branch("fTreeVariablePosTransvMomentum",&fTreeVariablePosTransvMomentum,"fTreeVariablePosTransvMomentum/F");
 239 /*18*/   fTree->Branch("fTreeVariableNegTransvMomentumMC",&fTreeVariableNegTransvMomentumMC,"fTreeVariableNegTransvMomentumMC/F");
 240 /*19*/   fTree->Branch("fTreeVariablePosTransvMomentumMC",&fTreeVariablePosTransvMomentumMC,"fTreeVariablePosTransvMomentumMC/F");
 241 /*20*/   fTree->Branch("fTreeVariableLeastNbrCrossedRows",&fTreeVariableLeastNbrCrossedRows,"fTreeVariableLeastNbrCrossedRows/I");
 242 /*21*/   fTree->Branch("fTreeVariableLeastRatioCrossedRowsOverFindable",&fTreeVariableLeastRatioCrossedRowsOverFindable,"fTreeVariableLeastRatioCrossedRowsOverFindable/F");
 243 /*22*/   fTree->Branch("fTreeVariablePID",&fTreeVariablePID,"fTreeVariablePID/I");
 244 /*23*/   fTree->Branch("fTreeVariablePIDPositive",&fTreeVariablePIDPositive,"fTreeVariablePIDPositive/I");
 245 /*24*/   fTree->Branch("fTreeVariablePIDNegative",&fTreeVariablePIDNegative,"fTreeVariablePIDNegative/I");
 246 /*25*/   fTree->Branch("fTreeVariablePIDMother",&fTreeVariablePIDMother,"fTreeVariablePIDMother/I");
 247 /*26*/   fTree->Branch("fTreeVariablePtXiMother",&fTreeVariablePtMother,"fTreeVariablePtMother/F");
 248 /*27*/   fTree->Branch("fTreeVariableV0CosineOfPointingAngle",&fTreeVariableV0CosineOfPointingAngle,"fTreeVariableV0CosineOfPointingAngle/F");
 249 //-----------MULTIPLICITY-INFO--------------------
 250 /*28*/   fTree->Branch("fTreeVariableMultiplicity",&fTreeVariableMultiplicity,"fTreeVariableMultiplicity/I");
 251 //------------------------------------------------
 252 /*29*/   fTree->Branch("fTreeVariableDistOverTotMom",&fTreeVariableDistOverTotMom,"fTreeVariableDistOverTotMom/F");
 253 /*30*/   fTree->Branch("fTreeVariableNSigmasPosProton",&fTreeVariableNSigmasPosProton,"fTreeVariableNSigmasPosProton/F");
 254 /*31*/   fTree->Branch("fTreeVariableNSigmasPosPion",&fTreeVariableNSigmasPosPion,"fTreeVariableNSigmasPosPion/F");
 255 /*32*/   fTree->Branch("fTreeVariableNSigmasNegProton",&fTreeVariableNSigmasNegProton,"fTreeVariableNSigmasNegProton/F");
 256 /*33*/   fTree->Branch("fTreeVariableNSigmasNegPion",&fTreeVariableNSigmasNegPion,"fTreeVariableNSigmasNegPion/F");
 257 //------------------------------------------------
 258 /*34*/   fTree->Branch("fTreeVariableNegEta",&fTreeVariableNegEta,"fTreeVariableNegEta/F");
 259 /*35*/   fTree->Branch("fTreeVariablePosEta",&fTreeVariablePosEta,"fTreeVariablePosEta/F");
 260 /*36*/   fTree->Branch("fTreeVariableV0CreationRadius",&fTreeVariableV0CreationRadius,"fTreeVariableV0CreationRadius/F");
 261 /*37*/   fTree->Branch("fTreeVariableIndexStatus",&fTreeVariableIndexStatus,"fTreeVariableIndexStatus/I");
 262 /*38*/   fTree->Branch("fTreeVariableIndexStatusMother",&fTreeVariableIndexStatusMother,"fTreeVariableIndexStatusMother/I");
 263
 264 //------------------------------------------------
 265 // Particle Identification Setup
 266 //------------------------------------------------
 267
 268    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 269    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 270    fPIDResponse = inputHandler->GetPIDResponse();
 271
 272 //------------------------------------------------
 273 // V0 Multiplicity Histograms
 274 //------------------------------------------------
 275
 276    // Create histograms
 277    OpenFile(1);
 278    fListHistV0 = new TList();
 279    fListHistV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
 280
 281
 282    if(! fHistV0MultiplicityBeforeTrigSel) {
 283       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
 284          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
 285          25, 0, 25);
 286       fListHistV0->Add(fHistV0MultiplicityBeforeTrigSel);
 287    }
 288
 289    if(! fHistV0MultiplicityForTrigEvt) {
 290       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
 291          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
 292          25, 0, 25);
 293       fListHistV0->Add(fHistV0MultiplicityForTrigEvt);
 294    }
 295
 296    if(! fHistV0MultiplicityForSelEvt) {
 297       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
 298          "V0s per event;Nbr of V0s/Evt;Events",
 299          25, 0, 25);
 300       fListHistV0->Add(fHistV0MultiplicityForSelEvt);
 301    }
 302
 303    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
 304       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
 305          "V0s per event;Nbr of V0s/Evt;Events",
 306          25, 0, 25);
 307       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
 308    }
 309    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
 310       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
 311          "V0s per event;Nbr of V0s/Evt;Events",
 312          25, 0, 25);
 313       fListHistV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
 314    }
 315
 316 //------------------------------------------------
 317 // Track Multiplicity Histograms
 318 //------------------------------------------------
 319
 320    if(! fHistMultiplicityBeforeTrigSel) {
 321       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
 322          "Tracks per event;Nbr of Tracks;Events",
 323          200, 0, 200);
 324       fListHistV0->Add(fHistMultiplicityBeforeTrigSel);
 325    }
 326    if(! fHistMultiplicityForTrigEvt) {
 327       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
 328          "Tracks per event;Nbr of Tracks;Events",
 329          200, 0, 200);
 330       fListHistV0->Add(fHistMultiplicityForTrigEvt);
 331    }
 332    if(! fHistMultiplicity) {
 333       fHistMultiplicity = new TH1F("fHistMultiplicity",
 334          "Tracks per event;Nbr of Tracks;Events",
 335          200, 0, 200);
 336       fListHistV0->Add(fHistMultiplicity);
 337    }
 338    if(! fHistMultiplicityNoTPCOnly) {
 339       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
 340          "Tracks per event;Nbr of Tracks;Events",
 341          200, 0, 200);
 342       fListHistV0->Add(fHistMultiplicityNoTPCOnly);
 343    }
 344    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
 345       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
 346          "Tracks per event;Nbr of Tracks;Events",
 347          200, 0, 200);
 348       fListHistV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);
 349    }
 350
 351 //------------------------------------------------
 352 // Generated Particle Histograms
 353 //------------------------------------------------
 354
 355    Int_t lCustomNBins = 200;
 356    Double_t lCustomPtUpperLimit = 20;
 357    Int_t lCustomNBinsMultiplicity = 100;
 358
 359 //----------------------------------
 360 // Raw Generated (Pre-physics-selection)
 361 //----------------------------------
 362
 363 //--- 3D Histo (Pt, Y, Multiplicity)
 364
 365    if(! f3dHistPrimRawPtVsYVsMultLambda) {
 366       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);
 367       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultLambda);
 368    }
 369    if(! f3dHistPrimRawPtVsYVsMultAntiLambda) {
 370       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);
 371       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultAntiLambda);
 372    }
 373    if(! f3dHistPrimRawPtVsYVsMultK0Short) {
 374       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);
 375       fListHistV0->Add(f3dHistPrimRawPtVsYVsMultK0Short);
 376    }
 377
 378 //--- 3D Histo (Pt, Y, Proper Decay Length)
 379
 380    if(! f3dHistPrimRawPtVsYVsDecayLengthLambda) {
 381       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);
 382       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthLambda);
 383    }
 384    if(! f3dHistPrimRawPtVsYVsDecayLengthAntiLambda) {
 385       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);
 386       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthAntiLambda);
 387    }
 388    if(! f3dHistPrimRawPtVsYVsDecayLengthK0Short) {
 389       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);
 390       fListHistV0->Add(f3dHistPrimRawPtVsYVsDecayLengthK0Short);
 391    }
 392
 393 //--- 3D Histo (Pt, Y, Multiplicity) for generated charged Xi (feeddown)
 394
 395    if(! f3dHistGenPtVsYVsMultXiMinus) {
 396       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);
 397       fListHistV0->Add(f3dHistGenPtVsYVsMultXiMinus);
 398    }
 399    if(! f3dHistGenPtVsYVsMultXiPlus) {
 400       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);
 401       fListHistV0->Add(f3dHistGenPtVsYVsMultXiPlus);
 402    }
 403
 404 //----------------------------------
 405 // Histos at analysis level
 406 //----------------------------------
 407
 408    if(! f3dHistPrimAnalysisPtVsYVsMultLambda) {
 409       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);
 410       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultLambda);
 411    }
 412    if(! f3dHistPrimAnalysisPtVsYVsMultAntiLambda) {
 413       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);
 414       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultAntiLambda);
 415    }
 416    if(! f3dHistPrimAnalysisPtVsYVsMultK0Short) {
 417       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);
 418       fListHistV0->Add(f3dHistPrimAnalysisPtVsYVsMultK0Short);
 419    }
 420
 421 //----------------------------------
 422 // Primary Vertex Position Histos
 423 //----------------------------------
 424
 425    if(! fHistPVx) {
 426          fHistPVx = new TH1F("fHistPVx",
 427             "PV x position;Nbr of Evts;x",
 428             2000, -0.5, 0.5);
 429       fListHistV0->Add(fHistPVx);
 430    }
 431    if(! fHistPVy) {
 432          fHistPVy = new TH1F("fHistPVy",
 433             "PV y position;Nbr of Evts;y",
 434             2000, -0.5, 0.5);
 435       fListHistV0->Add(fHistPVy);
 436    }
 437    if(! fHistPVz) {
 438          fHistPVz = new TH1F("fHistPVz",
 439             "PV z position;Nbr of Evts;z",
 440             400, -20, 20);
 441       fListHistV0->Add(fHistPVz);
 442    }
 443
 444    if(! fHistPVxAnalysis) {
 445          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
 446             "PV x position;Nbr of Evts;x",
 447             2000, -0.5, 0.5);
 448       fListHistV0->Add(fHistPVxAnalysis);
 449    }
 450    if(! fHistPVyAnalysis) {
 451          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
 452             "PV y position;Nbr of Evts;y",
 453             2000, -0.5, 0.5);
 454       fListHistV0->Add(fHistPVyAnalysis);
 455    }
 456    if(! fHistPVzAnalysis) {
 457          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
 458             "PV z position;Nbr of Evts;z",
 459             400, -20, 20);
 460       fListHistV0->Add(fHistPVzAnalysis);
 461    }
 462
 463    if(! fHistPVxAnalysisHasHighPtLambda) {
 464          fHistPVxAnalysisHasHighPtLambda = new TH1F("fHistPVxAnalysisHasHighPtLambda",
 465             "PV x position;Nbr of Evts;x",
 466             2000, -0.5, 0.5);
 467       fListHistV0->Add(fHistPVxAnalysisHasHighPtLambda);
 468    }
 469    if(! fHistPVyAnalysisHasHighPtLambda) {
 470          fHistPVyAnalysisHasHighPtLambda = new TH1F("fHistPVyAnalysisHasHighPtLambda",
 471             "PV y position;Nbr of Evts;y",
 472             2000, -0.5, 0.5);
 473       fListHistV0->Add(fHistPVyAnalysisHasHighPtLambda);
 474    }
 475    if(! fHistPVzAnalysisHasHighPtLambda) {
 476          fHistPVzAnalysisHasHighPtLambda = new TH1F("fHistPVzAnalysisHasHighPtLambda",
 477             "PV z position;Nbr of Evts;z",
 478             400, -20, 20);
 479       fListHistV0->Add(fHistPVzAnalysisHasHighPtLambda);
 480    }
 481    //List of Histograms: Normal
 482    PostData(1, fListHistV0);
 483
 484    //TTree Object: Saved to base directory. Should cache to disk while saving.
 485    //(Important to avoid excessive memory usage, particularly when merging)
 486    PostData(2, fTree);
 487
 488 }// end UserCreateOutputObjects
 489
 490
 491 //________________________________________________________________________
 492 void AliAnalysisTaskExtractPerformanceV0::UserExec(Option_t *)
 493 {
 494   // Main loop
 495   // Called for each event
 496
 497    AliESDEvent *lESDevent = 0x0;
 498    AliMCEvent  *lMCevent  = 0x0;
 499    AliStack    *lMCstack  = 0x0;
 500
 501    Int_t    lNumberOfV0s                      = -1;
 502    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 503    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 504    Double_t lMagneticField                 = -10.;
 505
 506   // Connect to the InputEvent
 507   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
 508
 509    // Appropriate for ESD analysis!
 510
 511    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 512    if (!lESDevent) {
 513       AliWarning("ERROR: lESDevent not available \n");
 514       return;
 515    }
 516
 517    lMCevent = MCEvent();
 518    if (!lMCevent) {
 519       Printf("ERROR: Could not retrieve MC event \n");
 520       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
 521       return;
 522    }
 523
 524    lMCstack = lMCevent->Stack();
 525    if (!lMCstack) {
 526       Printf("ERROR: Could not retrieve MC stack \n");
 527       cout << "Name of the file with pb :" <<  fInputHandler->GetTree()->GetCurrentFile()->GetName() << endl;
 528       return;
 529    }
 530
 531 //------------------------------------------------
 532 // Multiplicity Information Acquistion
 533 //------------------------------------------------
 534
 535   //REVISED multiplicity estimator after 'multiplicity day' (2011)
 536    Int_t lMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity( lESDevent );
 537
 538    fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
 539    fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
 540
 541 //------------------------------------------------
 542 // MC Information Acquistion
 543 //------------------------------------------------
 544
 545    Int_t iNumberOfPrimaries = -1;
 546    iNumberOfPrimaries = lMCstack->GetNprimary();
 547    if(iNumberOfPrimaries < 1) return;
 548    Bool_t lHasHighPtLambda = kFALSE;
 549
 550 //------------------------------------------------
 551 // Variable Definition
 552 //------------------------------------------------
 553
 554    Int_t lNbMCPrimary        = 0;
 555
 556    Int_t lPdgcodeCurrentPart = 0;
 557    Double_t lRapCurrentPart  = 0;
 558    Double_t lPtCurrentPart   = 0;
 559
 560    //Int_t lComeFromSigma      = 0;
 561
 562    // current mc particle 's mother
 563    //Int_t iCurrentMother  = 0;
 564    lNbMCPrimary = lMCstack->GetNprimary();
 565
 566 //------------------------------------------------
 567 // Pre-Physics Selection
 568 //------------------------------------------------
 569
 570 //----- Loop on primary Xi, Omega --------------------------------------------------------------
 571    for (Int_t iCurrentLabelStack = 0; iCurrentLabelStack < lNbMCPrimary; iCurrentLabelStack++)
 572    {// This is the begining of the loop on primaries
 573
 574       TParticle* lCurrentParticlePrimary = 0x0;
 575       lCurrentParticlePrimary = lMCstack->Particle( iCurrentLabelStack );
 576       if(!lCurrentParticlePrimary){
 577          Printf("Cascade loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
 578          continue;
 579       }
 580       if ( TMath::Abs(lCurrentParticlePrimary->GetPdgCode()) == 3312 ){
 581          Double_t lRapXiMCPrimary = 0.5*TMath::Log((lCurrentParticlePrimary->Energy() + lCurrentParticlePrimary->Pz()) / (lCurrentParticlePrimary->Energy() - lCurrentParticlePrimary->Pz() +1.e-13));
 582
 583          //=================================================================================
 584          // Xi Histograms for Feeddown - Primary Charged Xis
 585          if( lCurrentParticlePrimary->GetPdgCode() == 3312 ){
 586             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
 587             f3dHistGenPtVsYVsMultXiMinus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
 588          }
 589          if( lCurrentParticlePrimary->GetPdgCode() == -3312 ){
 590             lPtCurrentPart    = lCurrentParticlePrimary->Pt();
 591             f3dHistGenPtVsYVsMultXiPlus->Fill(lPtCurrentPart, lRapXiMCPrimary, lMultiplicity);
 592          }
 593       }
 594    }
 595 //----- End Loop on primary Xi, Omega ----------------------------------------------------------
 596
 597 //----- Loop on Lambda, K0Short ----------------------------------------------------------------
 598    for (Int_t iCurrentLabelStack = 0;  iCurrentLabelStack < (lMCstack->GetNtrack()); iCurrentLabelStack++)
 599    {// This is the begining of the loop on tracks
 600
 601       TParticle* lCurrentParticleForLambdaCheck = 0x0;
 602       lCurrentParticleForLambdaCheck = lMCstack->Particle( iCurrentLabelStack );
 603       if(!lCurrentParticleForLambdaCheck){
 604          Printf("V0s loop %d - MC TParticle pointer to current stack particle = 0x0 ! Skip ...\n", iCurrentLabelStack );
 605          continue;
 606       }
 607
 608       //=================================================================================
 609       //Single-Strange checks
 610       // Keep only K0s, Lambda and AntiLambda:
 611       lPdgcodeCurrentPart = lCurrentParticleForLambdaCheck->GetPdgCode();
 612
 613       if ( (lCurrentParticleForLambdaCheck->GetPdgCode() == 310   ) ||
 614            (lCurrentParticleForLambdaCheck->GetPdgCode() == 3122  ) ||
 615            (lCurrentParticleForLambdaCheck->GetPdgCode() == -3122 ) )
 616            {
 617          lRapCurrentPart   = MyRapidity(lCurrentParticleForLambdaCheck->Energy(),lCurrentParticleForLambdaCheck->Pz());
 618          lPtCurrentPart    = lCurrentParticleForLambdaCheck->Pt();
 619
 620          //Use Physical Primaries only for filling PrimRaw Histograms!
 621          if ( lMCstack->IsPhysicalPrimary(iCurrentLabelStack)!=kTRUE ) continue;
 622
 623          if( lPdgcodeCurrentPart == 3122 ){
 624             f3dHistPrimRawPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 625             if( TMath::Abs( lCurrentParticleForLambdaCheck->Eta() )<1.2 && lPtCurrentPart>2 ){
 626                lHasHighPtLambda = kTRUE; //Keep track of events with Lambda within |eta|<1.2 and pt>2
 627             }
 628          }
 629          if( lPdgcodeCurrentPart == -3122 ){
 630             f3dHistPrimRawPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 631          }
 632          if( lPdgcodeCurrentPart == 310 ){
 633             f3dHistPrimRawPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 634          }
 635          //Decay Length Acquisition=====================================================
 636          Double_t decaylength = -1;
 637          Double_t lV0Mass = -1;
 638
 639          if( !(lCurrentParticleForLambdaCheck->GetDaughter(0) < 0) ) {
 640             TParticle* lDght0ofV0 = lMCstack->Particle(  lCurrentParticleForLambdaCheck->GetDaughter(0) ); //get first daughter
 641             if(lDght0ofV0){ // skip if not defined.
 642                decaylength = TMath::Sqrt(
 643                                         TMath::Power( lCurrentParticleForLambdaCheck->Vx() - lDght0ofV0->Vx() , 2) +
 644                                         TMath::Power( lCurrentParticleForLambdaCheck->Vy() - lDght0ofV0->Vy() , 2) +
 645                                         TMath::Power( lCurrentParticleForLambdaCheck->Vz() - lDght0ofV0->Vz() , 2)
 646                );
 647                //Need to correct for relativitity! Involves multiplying by mass and dividing by momentum.
 648                if(TMath::Abs( lPdgcodeCurrentPart ) == 3122 ) { lV0Mass = 1.115683; }
 649                if(TMath::Abs( lPdgcodeCurrentPart ) == 310 ) { lV0Mass = 0.497614; }
 650                decaylength = ( lV0Mass * decaylength ) / ( lCurrentParticleForLambdaCheck->P() + 1e-10 );
 651             }
 652          }
 653          if( lPdgcodeCurrentPart == 3122) f3dHistPrimRawPtVsYVsDecayLengthLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
 654          if( lPdgcodeCurrentPart == -3122) f3dHistPrimRawPtVsYVsDecayLengthAntiLambda ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
 655          if( lPdgcodeCurrentPart == 310) f3dHistPrimRawPtVsYVsDecayLengthK0Short ->Fill( lPtCurrentPart, lRapCurrentPart , decaylength );
 656       }
 657    }//End of loop on tracks
 658 //----- End Loop on Lambda, K0Short ------------------------------------------------------------
 659
 660 // ---> Set Variables to Zero again
 661 // ---> Variable Definition
 662
 663    lPdgcodeCurrentPart = 0;
 664    lRapCurrentPart  = 0;
 665    lPtCurrentPart   = 0;
 666
 667 //------------------------------------------------
 668 // Physics Selection
 669 //------------------------------------------------
 670
 671    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 672    Bool_t isSelected = 0;
 673    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
 674    if ( ! isSelected ) {
 675       PostData(1, fListHistV0);
 676       PostData(2, fTree);
 677       return;
 678    }
 679
 680 //------------------------------------------------
 681 // After Trigger Selection
 682 //------------------------------------------------
 683
 684    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
 685
 686    //Set variable for filling tree afterwards!
 687    fTreeVariableMultiplicity = lMultiplicity;
 688    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
 689    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
 690
 691 //------------------------------------------------
 692 // Getting: Primary Vertex + MagField Info
 693 //------------------------------------------------
 694
 695    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
 696    // get the vtx stored in ESD found with tracks
 697    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
 698
 699    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
 700    // get the best primary vertex available for the event
 701    // As done in AliCascadeVertexer, we keep the one which is the best one available.
 702    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
 703    // This one will be used for next calculations (DCA essentially)
 704    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
 705
 706    Double_t lPrimaryVtxPosition[3];
 707    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
 708    lPrimaryVtxPosition[0] = primaryVtx->GetX();
 709    lPrimaryVtxPosition[1] = primaryVtx->GetY();
 710    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
 711    fHistPVx->Fill( lPrimaryVtxPosition[0] );
 712    fHistPVy->Fill( lPrimaryVtxPosition[1] );
 713    fHistPVz->Fill( lPrimaryVtxPosition[2] );
 714
 715 //------------------------------------------------
 716 // Primary Vertex Z position: SKIP
 717 //------------------------------------------------
 718
 719    if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0 ) {
 720       AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
 721       PostData(1, fListHistV0);
 722       PostData(2, fTree);
 723       return;
 724    }
 725
 726    lMagneticField = lESDevent->GetMagneticField( );
 727    fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
 728    fHistMultiplicity->Fill(lMultiplicity);
 729
 730 //------------------------------------------------
 731 // SKIP: Events with well-established PVtx
 732 //------------------------------------------------
 733
 734    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
 735    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
 736    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() ){
 737       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
 738       PostData(1, fListHistV0);
 739       PostData(2, fTree);
 740       return;
 741    }
 742    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
 743    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
 744
 745 //------------------------------------------------
 746 // Pileup Rejection Studies
 747 //------------------------------------------------
 748
 749    // FIXME : quality selection regarding pile-up rejection
 750    if(lESDevent->IsPileupFromSPD() ){// 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
 751       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
 752       PostData(1, fListHistV0);
 753       PostData(2, fTree);
 754       return;
 755    }
 756    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
 757    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
 758
 759    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
 760    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){
 761       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
 762       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
 763       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
 764       if ( lHasHighPtLambda == kTRUE ){
 765          fHistPVxAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[0] );
 766          fHistPVyAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[1] );
 767          fHistPVzAnalysisHasHighPtLambda->Fill( lPrimaryVtxPosition[2] );
 768       }
 769    }
 770
 771 //------------------------------------------------
 772 // stack loop starts here
 773 //------------------------------------------------
 774
 775 //---> Loop over ALL PARTICLES
 776
 777    for (Int_t iMc = 0; iMc < (lMCstack->GetNtrack()); iMc++) {
 778       TParticle *p0 = lMCstack->Particle(iMc);
 779       if (!p0) {
 780          //Printf("ERROR: particle with label %d not found in lMCstack (mc loop)", iMc);
 781          continue;
 782       }
 783       lPdgcodeCurrentPart = p0->GetPdgCode();
 784
 785       // Keep only K0s, Lambda and AntiLambda:
 786       if ( (lPdgcodeCurrentPart != 310 ) && (lPdgcodeCurrentPart != 3122 ) && (lPdgcodeCurrentPart != -3122 ) ) continue;
 787
 788       lRapCurrentPart   = MyRapidity(p0->Energy(),p0->Pz());
 789       lPtCurrentPart    = p0->Pt();
 790
 791         //Use Physical Primaries only for filling PrimRaw Histograms!
 792       if ( lMCstack->IsPhysicalPrimary(iMc)!=kTRUE ) continue;
 793
 794       if( lPdgcodeCurrentPart == 3122 ){
 795          f3dHistPrimAnalysisPtVsYVsMultLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 796       }
 797       if( lPdgcodeCurrentPart == -3122 ){
 798          f3dHistPrimAnalysisPtVsYVsMultAntiLambda->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 799       }
 800       if( lPdgcodeCurrentPart == 310 ){
 801          f3dHistPrimAnalysisPtVsYVsMultK0Short->Fill(lPtCurrentPart, lRapCurrentPart, lMultiplicity);
 802       }
 803    }
 804
 805 //------------------------------------------------
 806 // MAIN LAMBDA LOOP STARTS HERE
 807 //------------------------------------------------
 808
 809    //Variable definition
 810    Int_t    lOnFlyStatus = 0;
 811    Double_t lChi2V0 = 0;
 812    Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;
 813    Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;
 814    Double_t lV0CosineOfPointingAngle = 0;
 815    Double_t lV0Radius = 0, lPt = 0;
 816    Double_t lRapK0Short = 0, lRapLambda = 0;
 817    Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;
 818    Double_t lAlphaV0 = 0, lPtArmV0 = 0;
 819    Double_t fMinV0Pt = 0;
 820    Double_t fMaxV0Pt = 100;
 821
 822    Int_t nv0s = 0;
 823    nv0s = lESDevent->GetNumberOfV0s();
 824
 825    for (Int_t iV0 = 0; iV0 < nv0s; iV0++)
 826         {// This is the begining of the V0 loop
 827       AliESDv0 *v0 = ((AliESDEvent*)lESDevent)->GetV0(iV0);
 828       if (!v0) continue;
 829
 830       Double_t tV0mom[3];
 831       v0->GetPxPyPz( tV0mom[0],tV0mom[1],tV0mom[2] );
 832       Double_t lV0TotalMomentum = TMath::Sqrt(
 833          tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );
 834
 835       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]);
 836       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);
 837       lPt = v0->Pt();
 838       lRapK0Short = v0->RapK0Short();
 839       lRapLambda  = v0->RapLambda();
 840       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;
 841
 842       UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPindex());
 843       UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetNindex());
 844
 845       Double_t lMomPos[3]; v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);
 846       Double_t lMomNeg[3]; v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);
 847
 848       AliESDtrack *pTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyPos);
 849       AliESDtrack *nTrack=((AliESDEvent*)lESDevent)->GetTrack(lKeyNeg);
 850       if (!pTrack || !nTrack) {
 851          Printf("ERROR: Could not retreive one of the daughter track");
 852          continue;
 853       }
 854
 855       fTreeVariableNegEta = nTrack->Eta();
 856       fTreeVariablePosEta = pTrack->Eta();
 857
 858       // Filter like-sign V0 (next: add counter and distribution)
 859       if ( pTrack->GetSign() == nTrack->GetSign()){
 860          continue;
 861       }
 862
 863       //________________________________________________________________________
 864       // Track quality cuts
 865       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);
 866       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);
 867       fTreeVariableLeastNbrCrossedRows = lPosTrackCrossedRows;
 868       if( lNegTrackCrossedRows < fTreeVariableLeastNbrCrossedRows )
 869          fTreeVariableLeastNbrCrossedRows = lNegTrackCrossedRows;
 870
 871       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)
 872       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 873       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;
 874
 875       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;
 876
 877       //GetKinkIndex condition
 878       if( pTrack->GetKinkIndex(0)>0 || nTrack->GetKinkIndex(0)>0 ) continue;
 879
 880       //Findable clusters > 0 condition
 881       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;
 882
 883       //Compute ratio Crossed Rows / Findable clusters
 884       //Note: above test avoids division by zero!
 885       Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF()));
 886       Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF()));
 887
 888       fTreeVariableLeastRatioCrossedRowsOverFindable = lPosTrackCrossedRowsOverFindable;
 889       if( lNegTrackCrossedRowsOverFindable < fTreeVariableLeastRatioCrossedRowsOverFindable )
 890          fTreeVariableLeastRatioCrossedRowsOverFindable = lNegTrackCrossedRowsOverFindable;
 891
 892       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here
 893       if ( fTreeVariableLeastRatioCrossedRowsOverFindable < 0.8 ) continue;
 894
 895       //End track Quality Cuts
 896       //________________________________________________________________________
 897
 898       lDcaPosToPrimVertex = TMath::Abs(pTrack->GetD(lPrimaryVtxPosition[0],
 899                                                         lPrimaryVtxPosition[1],
 900                                                         lMagneticField) );
 901
 902       lDcaNegToPrimVertex = TMath::Abs(nTrack->GetD(lPrimaryVtxPosition[0],
 903                                                         lPrimaryVtxPosition[1],
 904                                                         lMagneticField) );
 905
 906       lOnFlyStatus = v0->GetOnFlyStatus();
 907       lChi2V0 = v0->GetChi2V0();
 908       lDcaV0Daughters = v0->GetDcaV0Daughters();
 909       lDcaV0ToPrimVertex = v0->GetD(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
 910       lV0CosineOfPointingAngle = v0->GetV0CosineOfPointingAngle(lPrimaryVtxPosition[0],lPrimaryVtxPosition[1],lPrimaryVtxPosition[2]);
 911       fTreeVariableV0CosineOfPointingAngle=lV0CosineOfPointingAngle;
 912
 913       // Getting invariant mass infos directly from ESD
 914       v0->ChangeMassHypothesis(310);
 915       lInvMassK0s = v0->GetEffMass();
 916       v0->ChangeMassHypothesis(3122);
 917       lInvMassLambda = v0->GetEffMass();
 918       v0->ChangeMassHypothesis(-3122);
 919       lInvMassAntiLambda = v0->GetEffMass();
 920       lAlphaV0 = v0->AlphaV0();
 921       lPtArmV0 = v0->PtArmV0();
 922
 923       //fTreeVariableOnFlyStatus = lOnFlyStatus;
 924       //fHistV0OnFlyStatus->Fill(lOnFlyStatus);
 925
 926 //===============================================
 927 // Monte Carlo Association starts here
 928 //===============================================
 929
 930       //---> Set Everything to "I don't know" before starting
 931
 932       fTreeVariablePIDPositive = 0;
 933       fTreeVariablePIDNegative = 0;
 934
 935       fTreeVariableIndexStatus = 0;
 936       fTreeVariableIndexStatusMother = 0;
 937
 938       fTreeVariablePtMother = -1;
 939       fTreeVariablePtMC = -1;
 940       fTreeVariableRapMC = -100;
 941
 942       fTreeVariablePID = -1;
 943       fTreeVariablePIDMother = -1;
 944
 945       fTreeVariablePrimaryStatus = 0;
 946       fTreeVariablePrimaryStatusMother = 0;
 947       fTreeVariableV0CreationRadius = -1;
 948
 949       Int_t lblPosV0Dghter = (Int_t) TMath::Abs( pTrack->GetLabel() );
 950       Int_t lblNegV0Dghter = (Int_t) TMath::Abs( nTrack->GetLabel() );
 951
 952       TParticle* mcPosV0Dghter = lMCstack->Particle( lblPosV0Dghter );
 953       TParticle* mcNegV0Dghter = lMCstack->Particle( lblNegV0Dghter );
 954
 955       fTreeVariablePosTransvMomentumMC = mcPosV0Dghter->Pt();
 956       fTreeVariableNegTransvMomentumMC = mcNegV0Dghter->Pt();
 957
 958       Int_t lPIDPositive = mcPosV0Dghter -> GetPdgCode();
 959       Int_t lPIDNegative = mcNegV0Dghter -> GetPdgCode();
 960
 961       fTreeVariablePIDPositive = lPIDPositive;
 962       fTreeVariablePIDNegative = lPIDNegative;
 963
 964       Int_t lblMotherPosV0Dghter = mcPosV0Dghter->GetFirstMother() ;
 965       Int_t lblMotherNegV0Dghter = mcNegV0Dghter->GetFirstMother();
 966
 967       if( lblMotherPosV0Dghter == lblMotherNegV0Dghter && lblMotherPosV0Dghter > -1 ){
 968          //either label is fine, they're equal at this stage
 969          TParticle* pThisV0 = lMCstack->Particle( lblMotherPosV0Dghter );
 970          //Set tree variables
 971          fTreeVariablePID   = pThisV0->GetPdgCode(); //PDG Code
 972          fTreeVariablePtMC  = pThisV0->Pt(); //Perfect Pt
 973          //Only Interested if it's a Lambda, AntiLambda or K0s
 974          //Avoid the Junction Bug! PYTHIA has particles with Px=Py=Pz=E=0 occasionally,
 975          //having particle code 88 (unrecognized by PDG), for documentation purposes.
 976          //Even ROOT's TParticle::Y() is not prepared to deal with that exception!
 977          //Note that TParticle::Pt() is immune (that would just return 0)...
 978          //Though granted that that should be extremely rare in this precise condition...
 979          if( TMath::Abs(fTreeVariablePID) == 3122 || fTreeVariablePID==310 ){
 980             fTreeVariableRapMC = pThisV0->Y(); //Perfect Y
 981          }
 982          fTreeVariableV0CreationRadius = pThisV0->R(); // Creation Radius
 983          if( lblMotherPosV0Dghter  < lNbMCPrimary ) fTreeVariableIndexStatus = 1; //looks primary
 984          if( lblMotherPosV0Dghter >= lNbMCPrimary ) fTreeVariableIndexStatus = 2; //looks secondary
 985          if( lMCstack->IsPhysicalPrimary       (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 1; //Is Primary!
 986          if( lMCstack->IsSecondaryFromWeakDecay(lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 2; //Weak Decay!
 987          if( lMCstack->IsSecondaryFromMaterial (lblMotherPosV0Dghter) ) fTreeVariablePrimaryStatus = 3; //Material Int!
 988
 989          //Now we try to acquire the V0 parent particle, if possible
 990          Int_t lblThisV0Parent = pThisV0->GetFirstMother();
 991          if ( lblThisV0Parent > -1 ){ //if it has a parent, get it and store specs
 992             TParticle* pThisV0Parent = lMCstack->Particle( lblThisV0Parent );
 993             fTreeVariablePIDMother   = pThisV0Parent->GetPdgCode(); //V0 Mother PDG
 994             fTreeVariablePtMother    = pThisV0Parent->Pt();         //V0 Mother Pt
 995             //Primary Status for the V0 Mother particle
 996             if( lblThisV0Parent  < lNbMCPrimary ) fTreeVariableIndexStatusMother = 1; //looks primary
 997             if( lblThisV0Parent >= lNbMCPrimary ) fTreeVariableIndexStatusMother = 2; //looks secondary
 998             if( lMCstack->IsPhysicalPrimary       (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 1; //Is Primary!
 999             if( lMCstack->IsSecondaryFromWeakDecay(lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 2; //Weak Decay!
1000             if( lMCstack->IsSecondaryFromMaterial (lblThisV0Parent) ) fTreeVariablePrimaryStatusMother = 3; //Material Int!
1001          }
1002       }
1003
1004       fTreeVariablePt = v0->Pt();
1005       fTreeVariableChi2V0 = lChi2V0;
1006       fTreeVariableDcaV0ToPrimVertex = lDcaV0ToPrimVertex;
1007       fTreeVariableDcaV0Daughters = lDcaV0Daughters;
1008       fTreeVariableV0CosineOfPointingAngle = lV0CosineOfPointingAngle;
1009       fTreeVariableV0Radius = lV0Radius;
1010       fTreeVariableDcaPosToPrimVertex = lDcaPosToPrimVertex;
1011       fTreeVariableDcaNegToPrimVertex = lDcaNegToPrimVertex;
1012       fTreeVariableInvMassK0s = lInvMassK0s;
1013       fTreeVariableInvMassLambda = lInvMassLambda;
1014       fTreeVariableInvMassAntiLambda = lInvMassAntiLambda;
1015       fTreeVariableRapK0Short = lRapK0Short;
1016
1017       fTreeVariableRapLambda = lRapLambda;
1018       fTreeVariableAlphaV0 = lAlphaV0;
1019       fTreeVariablePtArmV0 = lPtArmV0;
1020
1021       //Official means of acquiring N-sigmas
1022       fTreeVariableNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );
1023       fTreeVariableNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );
1024       fTreeVariableNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );
1025       fTreeVariableNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );
1026
1027 //tDecayVertexV0[0],tDecayVertexV0[1],tDecayVertexV0[2]
1028       fTreeVariableDistOverTotMom = TMath::Sqrt(
1029                                                 TMath::Power( tDecayVertexV0[0] - lBestPrimaryVtxPos[0] , 2) +
1030                                                 TMath::Power( tDecayVertexV0[1] - lBestPrimaryVtxPos[1] , 2) +
1031                                                 TMath::Power( tDecayVertexV0[2] - lBestPrimaryVtxPos[2] , 2)
1032                                         );
1033       fTreeVariableDistOverTotMom /= (lV0TotalMomentum + 1e-10); //avoid division by zero, to be sure
1034
1035       Double_t lMomentumPosTemp[3];
1036       pTrack->GetPxPyPz(lMomentumPosTemp);
1037       Double_t lPtPosTemporary = sqrt(pow(lMomentumPosTemp[0],2) + pow(lMomentumPosTemp[1],2));
1038
1039       Double_t lMomentumNegTemp[3];
1040       nTrack->GetPxPyPz(lMomentumNegTemp);
1041       Double_t lPtNegTemporary = sqrt(pow(lMomentumNegTemp[0],2) + pow(lMomentumNegTemp[1],2));
1042
1043       fTreeVariablePosTransvMomentum = lPtPosTemporary;
1044       fTreeVariableNegTransvMomentum = lPtNegTemporary;
1045
1046
1047 //------------------------------------------------
1048 // Fill Tree!
1049 //------------------------------------------------
1050
1051       // The conditionals are meant to decrease excessive
1052       // memory usage!
1053
1054       //Modified version: Keep only OnFlyStatus == 0
1055       //Keep only if included in a parametric InvMass Region 20 sigmas away from peak
1056
1057       //First Selection: Reject OnFly
1058       if( lOnFlyStatus == 0 ){
1059          //Second Selection: rough 20-sigma band, parametric.
1060          //K0Short: Enough to parametrize peak broadening with linear function.
1061          Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*fTreeVariablePt;
1062          Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*fTreeVariablePt;
1063          //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
1064          //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
1065          Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*fTreeVariablePt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*fTreeVariablePt);
1066          Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*fTreeVariablePt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*fTreeVariablePt);
1067          //Do Selection
1068          if( (fTreeVariableInvMassLambda     < lUpperLimitLambda  && fTreeVariableInvMassLambda     > lLowerLimitLambda     ) ||
1069              (fTreeVariableInvMassAntiLambda < lUpperLimitLambda  && fTreeVariableInvMassAntiLambda > lLowerLimitLambda     ) ||
1070              (fTreeVariableInvMassK0s        < lUpperLimitK0Short && fTreeVariableInvMassK0s        > lLowerLimitK0Short    ) ){
1071             fTree->Fill();
1072          }
1073       }
1074
1075 //------------------------------------------------
1076 // Fill tree over.
1077 //------------------------------------------------
1078
1079
1080    }// This is the end of the V0 loop
1081
1082    // Post output data.
1083    PostData(1, fListHistV0);
1084    PostData(2, fTree);
1085 }
1086
1087 //________________________________________________________________________
1088 void AliAnalysisTaskExtractPerformanceV0::Terminate(Option_t *)
1089 {
1090    // Draw result to the screen
1091    // Called once at the end of the query
1092
1093    TList *cRetrievedList = 0x0;
1094    cRetrievedList = (TList*)GetOutputData(1);
1095    if(!cRetrievedList){
1096       Printf("ERROR - AliAnalysisTaskExtractV0 : ouput data container list not available\n");
1097       return;
1098    }
1099
1100    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1101    if (!fHistV0MultiplicityForTrigEvt) {
1102       Printf("ERROR - AliAnalysisTaskExtractV0 : fHistV0MultiplicityForTrigEvt not available");
1103       return;
1104    }
1105
1106    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractV0","V0 Multiplicity",10,10,510,510);
1107    canCheck->cd(1)->SetLogy();
1108
1109    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1110    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1111 }
1112
1113 //----------------------------------------------------------------------------
1114
1115 Double_t AliAnalysisTaskExtractPerformanceV0::MyRapidity(Double_t rE, Double_t rPz) const
1116 {
1117    // Local calculation for rapidity
1118    return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1119 }
1120 //----------------------------------------------------------------------------
1121