PWGLF/STRANGENESS/Cascades/AliAnalysisTaskExtractCascade.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 //
  24 //  --- Algorithm Description
  25 //   1. Loop over primaries in stack to acquire generated charged Xi
  26 //   2. Loop over stack to find Cascades, fill TH3Fs "PrimRawPt"s for Efficiency
  27 //   3. Perform Physics Selection
  28 //   4. Perform Primary Vertex |z|<10cm selection
  29 //   5. Perform Primary Vertex NoTPCOnly vertexing selection
  30 //   6. Perform Pileup Rejection
  31 //   7. Analysis Loops:
  32 //    7a. Fill TH3Fs "PrimAnalysisPt" for control purposes only
  33 //
  34 //  Please Report Any Bugs!
  35 //
  36 //   --- David Dobrigkeit Chinellato
  37 //        (david.chinellato@gmail.com)
  38 //
  39 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
  40
  41 class TTree;
  42 class TParticle;
  43 class TVector3;
  44
  45 //class AliMCEventHandler;
  46 //class AliMCEvent;
  47 //class AliStack;
  48
  49 class AliESDVertex;
  50 class AliAODVertex;
  51 class AliESDv0;
  52 class AliAODv0;
  53
  54 #include <Riostream.h>
  55 #include "TList.h"
  56 #include "TH1.h"
  57 #include "TH2.h"
  58 #include "TH3.h"
  59 #include "TFile.h"
  60 #include "THnSparse.h"
  61 #include "TVector3.h"
  62 #include "TCanvas.h"
  63 #include "TMath.h"
  64 #include "TLegend.h"
  65 #include "AliLog.h"
  66
  67 #include "AliESDEvent.h"
  68 #include "AliAODEvent.h"
  69 #include "AliV0vertexer.h"
  70 #include "AliCascadeVertexer.h"
  71 #include "AliESDpid.h"
  72 #include "AliESDtrack.h"
  73 #include "AliESDtrackCuts.h"
  74 #include "AliInputEventHandler.h"
  75 #include "AliAnalysisManager.h"
  76 #include "AliMCEventHandler.h"
  77 #include "AliMCEvent.h"
  78 #include "AliStack.h"
  79
  80 #include "AliV0vertexer.h"
  81 #include "AliCascadeVertexer.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 #include "AliAnalysisUtils.h"
  91
  92 #include "AliAnalysisTaskExtractCascade.h"
  93
  94 using std::cout;
  95 using std::endl;
  96
  97 ClassImp(AliAnalysisTaskExtractCascade)
  98
  99 AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade()
 100   : AliAnalysisTaskSE(), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), fUtils(0),
 101    fkIsNuclear   ( kFALSE ),
 102    fkSwitchINT7  ( kFALSE ),
 103    fCentralityEstimator("V0M"),
 104    fkpAVertexSelection( kFALSE ),
 105    fEtaRefMult ( 0.5 ),
 106 //------------------------------------------------
 107 // Tree Variables
 108 //------------------------------------------------
 109
 110    fTreeCascVarCharge(0),
 111    fTreeCascVarMassAsXi(0),
 112    fTreeCascVarMassAsOmega(0),
 113    fTreeCascVarPt(0),
 114    fTreeCascVarPtMC(0),
 115    fTreeCascVarRapMC(0),
 116    fTreeCascVarRapXi(0),
 117    fTreeCascVarRapOmega(0),
 118    fTreeCascVarNegEta(0),
 119    fTreeCascVarPosEta(0),
 120    fTreeCascVarBachEta(0),
 121    fTreeCascVarDCACascDaughters(0),
 122    fTreeCascVarDCABachToPrimVtx(0),
 123    fTreeCascVarDCAV0Daughters(0),
 124    fTreeCascVarDCAV0ToPrimVtx(0),
 125    fTreeCascVarDCAPosToPrimVtx(0),
 126    fTreeCascVarDCANegToPrimVtx(0),
 127    fTreeCascVarCascCosPointingAngle(0),
 128    fTreeCascVarCascRadius(0),
 129    fTreeCascVarV0Mass(0),
 130    fTreeCascVarV0CosPointingAngle(0),
 131    fTreeCascVarV0Radius(0),
 132    fTreeCascVarLeastNbrClusters(0),
 133    fTreeCascVarMultiplicity(0),
 134    fTreeCascVarMultiplicityV0A(0),
 135    fTreeCascVarMultiplicityZNA(0),
 136    fTreeCascVarMultiplicityTRK(0),
 137    fTreeCascVarMultiplicitySPD(0),
 138    fTreeCascVarDistOverTotMom(0),
 139    fTreeCascVarPID(0),
 140    fTreeCascVarPIDBachelor(0),
 141    fTreeCascVarPIDNegative(0),
 142    fTreeCascVarPIDPositive(0),
 143    fTreeCascVarPosTransMom(0),
 144    fTreeCascVarNegTransMom(0),
 145    fTreeCascVarPosTransMomMC(0),
 146    fTreeCascVarNegTransMomMC(0),
 147    fTreeCascVarNegNSigmaPion(0),
 148    fTreeCascVarNegNSigmaProton(0),
 149    fTreeCascVarPosNSigmaPion(0),
 150    fTreeCascVarPosNSigmaProton(0),
 151    fTreeCascVarBachNSigmaPion(0),
 152    fTreeCascVarBachNSigmaKaon(0),
 153
 154 //------------------------------------------------
 155 // HISTOGRAMS
 156 // --- Filled on an Event-by-event basis
 157 //------------------------------------------------
 158    fHistV0MultiplicityBeforeTrigSel(0),
 159    fHistV0MultiplicityForTrigEvt(0),
 160    fHistV0MultiplicityForSelEvt(0),
 161    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 162    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 163    fHistMultiplicityBeforeTrigSel(0),
 164    fHistMultiplicityForTrigEvt(0),
 165    fHistMultiplicity(0),
 166    fHistMultiplicityNoTPCOnly(0),
 167    fHistMultiplicityNoTPCOnlyNoPileup(0),
 168
 169 //V0A Centrality
 170 fHistMultiplicityV0ABeforeTrigSel(0),
 171 fHistMultiplicityV0AForTrigEvt(0),
 172 fHistMultiplicityV0A(0),
 173 fHistMultiplicityV0ANoTPCOnly(0),
 174 fHistMultiplicityV0ANoTPCOnlyNoPileup(0),
 175
 176 //ZNA Centrality
 177 fHistMultiplicityZNABeforeTrigSel(0),
 178 fHistMultiplicityZNAForTrigEvt(0),
 179 fHistMultiplicityZNA(0),
 180 fHistMultiplicityZNANoTPCOnly(0),
 181 fHistMultiplicityZNANoTPCOnlyNoPileup(0),
 182
 183 //TRK Centrality
 184 fHistMultiplicityTRKBeforeTrigSel(0),
 185 fHistMultiplicityTRKForTrigEvt(0),
 186 fHistMultiplicityTRK(0),
 187 fHistMultiplicityTRKNoTPCOnly(0),
 188 fHistMultiplicityTRKNoTPCOnlyNoPileup(0),
 189
 190 //SPD Centrality
 191 fHistMultiplicitySPDBeforeTrigSel(0),
 192 fHistMultiplicitySPDForTrigEvt(0),
 193 fHistMultiplicitySPD(0),
 194 fHistMultiplicitySPDNoTPCOnly(0),
 195 fHistMultiplicitySPDNoTPCOnlyNoPileup(0),
 196
 197    fHistPVx(0),
 198    fHistPVy(0),
 199    fHistPVz(0),
 200    fHistPVxAnalysis(0),
 201    fHistPVyAnalysis(0),
 202    fHistPVzAnalysis(0)
 203 {
 204   // Dummy Constructor
 205 }
 206
 207 AliAnalysisTaskExtractCascade::AliAnalysisTaskExtractCascade(const char *name)
 208   : AliAnalysisTaskSE(name), fListHist(0), fTreeCascade(0), fPIDResponse(0), fESDtrackCuts(0), fUtils(0),
 209    fkIsNuclear   ( kFALSE ),
 210    fkSwitchINT7  ( kFALSE ),
 211    fCentralityEstimator("V0M"),
 212    fkpAVertexSelection( kFALSE ),
 213    fEtaRefMult ( 0.5 ),
 214 //------------------------------------------------
 215 // Tree Variables
 216 //------------------------------------------------
 217
 218    fTreeCascVarCharge(0),
 219    fTreeCascVarMassAsXi(0),
 220    fTreeCascVarMassAsOmega(0),
 221    fTreeCascVarPt(0),
 222    fTreeCascVarPtMC(0),
 223    fTreeCascVarRapMC(0),
 224    fTreeCascVarRapXi(0),
 225    fTreeCascVarRapOmega(0),
 226    fTreeCascVarNegEta(0),
 227    fTreeCascVarPosEta(0),
 228    fTreeCascVarBachEta(0),
 229    fTreeCascVarDCACascDaughters(0),
 230    fTreeCascVarDCABachToPrimVtx(0),
 231    fTreeCascVarDCAV0Daughters(0),
 232    fTreeCascVarDCAV0ToPrimVtx(0),
 233    fTreeCascVarDCAPosToPrimVtx(0),
 234    fTreeCascVarDCANegToPrimVtx(0),
 235    fTreeCascVarCascCosPointingAngle(0),
 236    fTreeCascVarCascRadius(0),
 237    fTreeCascVarV0Mass(0),
 238    fTreeCascVarV0CosPointingAngle(0),
 239    fTreeCascVarV0Radius(0),
 240    fTreeCascVarLeastNbrClusters(0),
 241    fTreeCascVarMultiplicity(0),
 242    fTreeCascVarMultiplicityV0A(0),
 243    fTreeCascVarMultiplicityZNA(0),
 244    fTreeCascVarMultiplicityTRK(0),
 245    fTreeCascVarMultiplicitySPD(0),
 246    fTreeCascVarDistOverTotMom(0),
 247    fTreeCascVarPID(0),
 248    fTreeCascVarPIDBachelor(0),
 249    fTreeCascVarPIDNegative(0),
 250    fTreeCascVarPIDPositive(0),
 251    fTreeCascVarPosTransMom(0),
 252    fTreeCascVarNegTransMom(0),
 253    fTreeCascVarPosTransMomMC(0),
 254    fTreeCascVarNegTransMomMC(0),
 255    fTreeCascVarNegNSigmaPion(0),
 256    fTreeCascVarNegNSigmaProton(0),
 257    fTreeCascVarPosNSigmaPion(0),
 258    fTreeCascVarPosNSigmaProton(0),
 259    fTreeCascVarBachNSigmaPion(0),
 260    fTreeCascVarBachNSigmaKaon(0),
 261
 262 //------------------------------------------------
 263 // HISTOGRAMS
 264 // --- Filled on an Event-by-event basis
 265 //------------------------------------------------
 266    fHistV0MultiplicityBeforeTrigSel(0),
 267    fHistV0MultiplicityForTrigEvt(0),
 268    fHistV0MultiplicityForSelEvt(0),
 269    fHistV0MultiplicityForSelEvtNoTPCOnly(0),
 270    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),
 271    fHistMultiplicityBeforeTrigSel(0),
 272    fHistMultiplicityForTrigEvt(0),
 273    fHistMultiplicity(0),
 274    fHistMultiplicityNoTPCOnly(0),
 275    fHistMultiplicityNoTPCOnlyNoPileup(0),
 276
 277 //V0A Centrality
 278 fHistMultiplicityV0ABeforeTrigSel(0),
 279 fHistMultiplicityV0AForTrigEvt(0),
 280 fHistMultiplicityV0A(0),
 281 fHistMultiplicityV0ANoTPCOnly(0),
 282 fHistMultiplicityV0ANoTPCOnlyNoPileup(0),
 283
 284 //ZNA Centrality
 285 fHistMultiplicityZNABeforeTrigSel(0),
 286 fHistMultiplicityZNAForTrigEvt(0),
 287 fHistMultiplicityZNA(0),
 288 fHistMultiplicityZNANoTPCOnly(0),
 289 fHistMultiplicityZNANoTPCOnlyNoPileup(0),
 290
 291 //TRK Centrality
 292 fHistMultiplicityTRKBeforeTrigSel(0),
 293 fHistMultiplicityTRKForTrigEvt(0),
 294 fHistMultiplicityTRK(0),
 295 fHistMultiplicityTRKNoTPCOnly(0),
 296 fHistMultiplicityTRKNoTPCOnlyNoPileup(0),
 297
 298 //SPD Centrality
 299 fHistMultiplicitySPDBeforeTrigSel(0),
 300 fHistMultiplicitySPDForTrigEvt(0),
 301 fHistMultiplicitySPD(0),
 302 fHistMultiplicitySPDNoTPCOnly(0),
 303 fHistMultiplicitySPDNoTPCOnlyNoPileup(0),
 304
 305    fHistPVx(0),
 306    fHistPVy(0),
 307    fHistPVz(0),
 308    fHistPVxAnalysis(0),
 309    fHistPVyAnalysis(0),
 310    fHistPVzAnalysis(0)
 311 {
 312    // Constructor
 313
 314    //Set Variables for re-running the cascade vertexers (as done for MS paper)
 315
 316         // New Loose : 1st step for the 7 TeV pp analysis
 317         /*
 318         fV0Sels[0] =  33.  ;  // max allowed chi2
 319         fV0Sels[1] =   0.02;  // min allowed impact parameter for the 1st daughter (LHC09a4 : 0.05)
 320         fV0Sels[2] =   0.02;  // min allowed impact parameter for the 2nd daughter (LHC09a4 : 0.05)
 321         fV0Sels[3] =   2.0 ;  // max allowed DCA between the daughter tracks       (LHC09a4 : 0.5)
 322         fV0Sels[4] =   0.95;  // min allowed cosine of V0's pointing angle         (LHC09a4 : 0.99)
 323         fV0Sels[5] =   1.0 ;  // min radius of the fiducial volume                 (LHC09a4 : 0.2)
 324         fV0Sels[6] = 100.  ;  // max radius of the fiducial volume                 (LHC09a4 : 100.0)
 325
 326         fCascSels[0] =  33.   ;  // max allowed chi2 (same as PDC07)
 327         fCascSels[1] =   0.05 ;  // min allowed V0 impact parameter                    (PDC07 : 0.05   / LHC09a4 : 0.025 )
 328         fCascSels[2] =   0.010;  // "window" around the Lambda mass                    (PDC07 : 0.008  / LHC09a4 : 0.010 )
 329         fCascSels[3] =   0.03 ;  // min allowed bachelor's impact parameter            (PDC07 : 0.035  / LHC09a4 : 0.025 )
 330         fCascSels[4] =   2.0  ;  // max allowed DCA between the V0 and the bachelor    (PDC07 : 0.1    / LHC09a4 : 0.2   )
 331         fCascSels[5] =   0.95 ;  // min allowed cosine of the cascade pointing angle   (PDC07 : 0.9985 / LHC09a4 : 0.998 )
 332         fCascSels[6] =   0.4  ;  // min radius of the fiducial volume                  (PDC07 : 0.9    / LHC09a4 : 0.2   )
 333         fCascSels[7] = 100.   ;  // max radius of the fiducial volume                  (PDC07 : 100    / LHC09a4 : 100   )
 334         */
 335    // Output slot #0 writes into a TList container (Cascade)
 336    DefineOutput(1, TList::Class());
 337    DefineOutput(2, TTree::Class());
 338 }
 339
 340
 341 AliAnalysisTaskExtractCascade::~AliAnalysisTaskExtractCascade()
 342 {
 343 //------------------------------------------------
 344 // DESTRUCTOR
 345 //------------------------------------------------
 346
 347    if (fListHist){
 348       delete fListHist;
 349       fListHist = 0x0;
 350    }
 351    if (fTreeCascade){
 352       delete fTreeCascade;
 353       fTreeCascade = 0x0;
 354    }
 355     //cleanup esd track cuts object too...
 356    if (fESDtrackCuts){
 357     delete fESDtrackCuts;
 358     fESDtrackCuts = 0x0;
 359   }
 360   if (fUtils){
 361     delete fUtils;
 362     fUtils = 0x0;
 363   }
 364
 365 }
 366
 367 //________________________________________________________________________
 368 void AliAnalysisTaskExtractCascade::UserCreateOutputObjects()
 369 {
 370    OpenFile(2);
 371    // Called once
 372
 373 //------------------------------------------------
 374
 375    fTreeCascade = new TTree("fTreeCascade","CascadeCandidates");
 376
 377 //------------------------------------------------
 378 // fTreeCascade Branch definitions - Cascade Tree
 379 //------------------------------------------------
 380
 381 //------------------------------------------------
 382 // fTreeCascade Branch definitions
 383 //------------------------------------------------
 384
 385 //-----------BASIC-INFO---------------------------
 386 /* 1*/          fTreeCascade->Branch("fTreeCascVarCharge",&fTreeCascVarCharge,"fTreeCascVarCharge/I");
 387 /* 2*/          fTreeCascade->Branch("fTreeCascVarMassAsXi",&fTreeCascVarMassAsXi,"fTreeCascVarMassAsXi/F");
 388 /* 3*/          fTreeCascade->Branch("fTreeCascVarMassAsOmega",&fTreeCascVarMassAsOmega,"fTreeCascVarMassAsOmega/F");
 389 /* 4*/          fTreeCascade->Branch("fTreeCascVarPt",&fTreeCascVarPt,"fTreeCascVarPt/F");
 390 /* 5*/          fTreeCascade->Branch("fTreeCascVarRapXi",&fTreeCascVarRapXi,"fTreeCascVarRapXi/F");
 391 /* 6*/          fTreeCascade->Branch("fTreeCascVarRapOmega",&fTreeCascVarRapOmega,"fTreeCascVarRapOmega/F");
 392 /* 7*/          fTreeCascade->Branch("fTreeCascVarNegEta",&fTreeCascVarNegEta,"fTreeCascVarNegEta/F");
 393 /* 8*/          fTreeCascade->Branch("fTreeCascVarPosEta",&fTreeCascVarPosEta,"fTreeCascVarPosEta/F");
 394 /* 9*/          fTreeCascade->Branch("fTreeCascVarBachEta",&fTreeCascVarBachEta,"fTreeCascVarBachEta/F");
 395 //-----------INFO-FOR-CUTS------------------------
 396 /*10*/          fTreeCascade->Branch("fTreeCascVarDCACascDaughters",&fTreeCascVarDCACascDaughters,"fTreeCascVarDCACascDaughters/F");
 397 /*11*/          fTreeCascade->Branch("fTreeCascVarDCABachToPrimVtx",&fTreeCascVarDCABachToPrimVtx,"fTreeCascVarDCABachToPrimVtx/F");
 398 /*12*/          fTreeCascade->Branch("fTreeCascVarDCAV0Daughters",&fTreeCascVarDCAV0Daughters,"fTreeCascVarDCAV0Daughters/F");
 399 /*13*/          fTreeCascade->Branch("fTreeCascVarDCAV0ToPrimVtx",&fTreeCascVarDCAV0ToPrimVtx,"fTreeCascVarDCAV0ToPrimVtx/F");
 400 /*14*/          fTreeCascade->Branch("fTreeCascVarDCAPosToPrimVtx",&fTreeCascVarDCAPosToPrimVtx,"fTreeCascVarDCAPosToPrimVtx/F");
 401 /*15*/          fTreeCascade->Branch("fTreeCascVarDCANegToPrimVtx",&fTreeCascVarDCANegToPrimVtx,"fTreeCascVarDCANegToPrimVtx/F");
 402 /*16*/          fTreeCascade->Branch("fTreeCascVarCascCosPointingAngle",&fTreeCascVarCascCosPointingAngle,"fTreeCascVarCascCosPointingAngle/F");
 403 /*17*/          fTreeCascade->Branch("fTreeCascVarCascRadius",&fTreeCascVarCascRadius,"fTreeCascVarCascRadius/F");
 404 /*18*/          fTreeCascade->Branch("fTreeCascVarV0Mass",&fTreeCascVarV0Mass,"fTreeCascVarV0Mass/F");
 405 /*19*/          fTreeCascade->Branch("fTreeCascVarV0CosPointingAngle",&fTreeCascVarV0CosPointingAngle,"fTreeCascVarV0CosPointingAngle/F");
 406 /*20*/          fTreeCascade->Branch("fTreeCascVarV0Radius",&fTreeCascVarV0Radius,"fTreeCascVarV0Radius/F");
 407 /*21*/          fTreeCascade->Branch("fTreeCascVarLeastNbrClusters",&fTreeCascVarLeastNbrClusters,"fTreeCascVarLeastNbrClusters/I");
 408 //-----------MULTIPLICITY-INFO--------------------
 409 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicity",&fTreeCascVarMultiplicity,"fTreeCascVarMultiplicity/I");
 410 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityV0A",&fTreeCascVarMultiplicityV0A,"fTreeCascVarMultiplicityV0A/I");
 411 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityZNA",&fTreeCascVarMultiplicityZNA,"fTreeCascVarMultiplicityZNA/I");
 412 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicityTRK",&fTreeCascVarMultiplicityTRK,"fTreeCascVarMultiplicityTRK/I");
 413 /*22*/          fTreeCascade->Branch("fTreeCascVarMultiplicitySPD",&fTreeCascVarMultiplicitySPD,"fTreeCascVarMultiplicitySPD/I");
 414 //-----------DECAY-LENGTH-INFO--------------------
 415 /*23*/          fTreeCascade->Branch("fTreeCascVarDistOverTotMom",&fTreeCascVarDistOverTotMom,"fTreeCascVarDistOverTotMom/F");
 416 //------------------------------------------------
 417 /*24*/          fTreeCascade->Branch("fTreeCascVarNegNSigmaPion",&fTreeCascVarNegNSigmaPion,"fTreeCascVarNegNSigmaPion/F");
 418 /*25*/          fTreeCascade->Branch("fTreeCascVarNegNSigmaProton",&fTreeCascVarNegNSigmaProton,"fTreeCascVarNegNSigmaProton/F");
 419 /*26*/          fTreeCascade->Branch("fTreeCascVarPosNSigmaPion",&fTreeCascVarPosNSigmaPion,"fTreeCascVarPosNSigmaPion/F");
 420 /*27*/          fTreeCascade->Branch("fTreeCascVarPosNSigmaProton",&fTreeCascVarPosNSigmaProton,"fTreeCascVarPosNSigmaProton/F");
 421 /*28*/          fTreeCascade->Branch("fTreeCascVarBachNSigmaPion",&fTreeCascVarBachNSigmaPion,"fTreeCascVarBachNSigmaPion/F");
 422 /*29*/          fTreeCascade->Branch("fTreeCascVarBachNSigmaKaon",&fTreeCascVarBachNSigmaKaon,"fTreeCascVarBachNSigmaKaon/F");
 423
 424 //------------------------------------------------
 425 // Particle Identification Setup
 426 //------------------------------------------------
 427
 428    AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 429    AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 430    fPIDResponse = inputHandler->GetPIDResponse();
 431
 432 // Multiplicity
 433
 434   if(! fESDtrackCuts ){
 435     fESDtrackCuts = new AliESDtrackCuts();
 436   }
 437   if(! fUtils ){
 438     fUtils = new AliAnalysisUtils();
 439   }
 440
 441 //------------------------------------------------
 442 // V0 Multiplicity Histograms
 443 //------------------------------------------------
 444
 445    // Create histograms
 446    OpenFile(1);
 447    fListHist = new TList();
 448    fListHist->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner
 449
 450
 451    if(! fHistV0MultiplicityBeforeTrigSel) {
 452       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel",
 453          "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events",
 454          25, 0, 25);
 455       fListHist->Add(fHistV0MultiplicityBeforeTrigSel);
 456    }
 457
 458    if(! fHistV0MultiplicityForTrigEvt) {
 459       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt",
 460          "V0s per event (for triggered evt);Nbr of V0s/Evt;Events",
 461          25, 0, 25);
 462       fListHist->Add(fHistV0MultiplicityForTrigEvt);
 463    }
 464
 465    if(! fHistV0MultiplicityForSelEvt) {
 466       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt",
 467          "V0s per event;Nbr of V0s/Evt;Events",
 468          25, 0, 25);
 469       fListHist->Add(fHistV0MultiplicityForSelEvt);
 470    }
 471
 472    if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {
 473       fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly",
 474          "V0s per event;Nbr of V0s/Evt;Events",
 475          25, 0, 25);
 476       fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);
 477    }
 478    if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {
 479       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup",
 480          "V0s per event;Nbr of V0s/Evt;Events",
 481          25, 0, 25);
 482       fListHist->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);
 483    }
 484
 485 //------------------------------------------------
 486 // Track Multiplicity Histograms
 487 //------------------------------------------------
 488
 489    if(! fHistMultiplicityBeforeTrigSel) {
 490       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel",
 491          "Tracks per event;Nbr of Tracks;Events",
 492          200, 0, 200);
 493       fListHist->Add(fHistMultiplicityBeforeTrigSel);
 494    }
 495    if(! fHistMultiplicityForTrigEvt) {
 496       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt",
 497          "Tracks per event;Nbr of Tracks;Events",
 498          200, 0, 200);
 499       fListHist->Add(fHistMultiplicityForTrigEvt);
 500    }
 501    if(! fHistMultiplicity) {
 502       fHistMultiplicity = new TH1F("fHistMultiplicity",
 503          "Tracks per event;Nbr of Tracks;Events",
 504          200, 0, 200);
 505       fListHist->Add(fHistMultiplicity);
 506    }
 507    if(! fHistMultiplicityNoTPCOnly) {
 508       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly",
 509          "Tracks per event;Nbr of Tracks;Events",
 510          200, 0, 200);
 511       fListHist->Add(fHistMultiplicityNoTPCOnly);
 512    }
 513    if(! fHistMultiplicityNoTPCOnlyNoPileup) {
 514       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup",
 515          "Tracks per event;Nbr of Tracks;Events",
 516          200, 0, 200);
 517       fListHist->Add(fHistMultiplicityNoTPCOnlyNoPileup);
 518    }
 519
 520
 521   //V0A Centrality (if PbPb / pPb)
 522   if(! fHistMultiplicityV0ABeforeTrigSel) {
 523     fHistMultiplicityV0ABeforeTrigSel = new TH1F("fHistMultiplicityV0ABeforeTrigSel",
 524                                                  "Centrality Distribution: V0A;V0A Centrality;Events",
 525                                                  200, 0, 200);
 526     fListHist->Add(fHistMultiplicityV0ABeforeTrigSel);
 527   }
 528   if(! fHistMultiplicityV0AForTrigEvt) {
 529     fHistMultiplicityV0AForTrigEvt = new TH1F("fHistMultiplicityV0AForTrigEvt",
 530                                               "Centrality Distribution: V0A;V0A Centrality;Events",
 531                                               200, 0, 200);
 532     fListHist->Add(fHistMultiplicityV0AForTrigEvt);
 533   }
 534   if(! fHistMultiplicityV0A) {
 535     fHistMultiplicityV0A = new TH1F("fHistMultiplicityV0A",
 536                                     "Centrality Distribution: V0A;V0A Centrality;Events",
 537                                     200, 0, 200);
 538     fListHist->Add(fHistMultiplicityV0A);
 539   }
 540   if(! fHistMultiplicityV0ANoTPCOnly) {
 541     fHistMultiplicityV0ANoTPCOnly = new TH1F("fHistMultiplicityV0ANoTPCOnly",
 542                                              "Centrality Distribution: V0A;V0A Centrality;Events",
 543                                              200, 0, 200);
 544     fListHist->Add(fHistMultiplicityV0ANoTPCOnly);
 545   }
 546   if(! fHistMultiplicityV0ANoTPCOnlyNoPileup) {
 547     fHistMultiplicityV0ANoTPCOnlyNoPileup = new TH1F("fHistMultiplicityV0ANoTPCOnlyNoPileup",
 548                                                      "Centrality Distribution: V0A;V0A Centrality;Events",
 549                                                      200, 0, 200);
 550     fListHist->Add(fHistMultiplicityV0ANoTPCOnlyNoPileup);
 551   }
 552
 553   //ZNA Centrality (if PbPb / pPb)
 554   if(! fHistMultiplicityZNABeforeTrigSel) {
 555     fHistMultiplicityZNABeforeTrigSel = new TH1F("fHistMultiplicityZNABeforeTrigSel",
 556                                                  "Centrality Distribution: ZNA;ZNA Centrality;Events",
 557                                                  200, 0, 200);
 558     fListHist->Add(fHistMultiplicityZNABeforeTrigSel);
 559   }
 560   if(! fHistMultiplicityZNAForTrigEvt) {
 561     fHistMultiplicityZNAForTrigEvt = new TH1F("fHistMultiplicityZNAForTrigEvt",
 562                                               "Centrality Distribution: ZNA;ZNA Centrality;Events",
 563                                               200, 0, 200);
 564     fListHist->Add(fHistMultiplicityZNAForTrigEvt);
 565   }
 566   if(! fHistMultiplicityZNA) {
 567     fHistMultiplicityZNA = new TH1F("fHistMultiplicityZNA",
 568                                     "Centrality Distribution: ZNA;ZNA Centrality;Events",
 569                                     200, 0, 200);
 570     fListHist->Add(fHistMultiplicityZNA);
 571   }
 572   if(! fHistMultiplicityZNANoTPCOnly) {
 573     fHistMultiplicityZNANoTPCOnly = new TH1F("fHistMultiplicityZNANoTPCOnly",
 574                                              "Centrality Distribution: ZNA;ZNA Centrality;Events",
 575                                              200, 0, 200);
 576     fListHist->Add(fHistMultiplicityZNANoTPCOnly);
 577   }
 578   if(! fHistMultiplicityZNANoTPCOnlyNoPileup) {
 579     fHistMultiplicityZNANoTPCOnlyNoPileup = new TH1F("fHistMultiplicityZNANoTPCOnlyNoPileup",
 580                                                      "Centrality Distribution: ZNA;ZNA Centrality;Events",
 581                                                      200, 0, 200);
 582     fListHist->Add(fHistMultiplicityZNANoTPCOnlyNoPileup);
 583   }
 584
 585   //TRK Centrality (if PbPb / pPb)
 586   if(! fHistMultiplicityTRKBeforeTrigSel) {
 587     fHistMultiplicityTRKBeforeTrigSel = new TH1F("fHistMultiplicityTRKBeforeTrigSel",
 588                                                  "Centrality Distribution: TRK;TRK Centrality;Events",
 589                                                  200, 0, 200);
 590     fListHist->Add(fHistMultiplicityTRKBeforeTrigSel);
 591   }
 592   if(! fHistMultiplicityTRKForTrigEvt) {
 593     fHistMultiplicityTRKForTrigEvt = new TH1F("fHistMultiplicityTRKForTrigEvt",
 594                                               "Centrality Distribution: TRK;TRK Centrality;Events",
 595                                               200, 0, 200);
 596     fListHist->Add(fHistMultiplicityTRKForTrigEvt);
 597   }
 598   if(! fHistMultiplicityTRK) {
 599     fHistMultiplicityTRK = new TH1F("fHistMultiplicityTRK",
 600                                     "Centrality Distribution: TRK;TRK Centrality;Events",
 601                                     200, 0, 200);
 602     fListHist->Add(fHistMultiplicityTRK);
 603   }
 604   if(! fHistMultiplicityTRKNoTPCOnly) {
 605     fHistMultiplicityTRKNoTPCOnly = new TH1F("fHistMultiplicityTRKNoTPCOnly",
 606                                              "Centrality Distribution: TRK;TRK Centrality;Events",
 607                                              200, 0, 200);
 608     fListHist->Add(fHistMultiplicityTRKNoTPCOnly);
 609   }
 610   if(! fHistMultiplicityTRKNoTPCOnlyNoPileup) {
 611     fHistMultiplicityTRKNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityTRKNoTPCOnlyNoPileup",
 612                                                      "Centrality Distribution: TRK;TRK Centrality;Events",
 613                                                      200, 0, 200);
 614     fListHist->Add(fHistMultiplicityTRKNoTPCOnlyNoPileup);
 615   }
 616
 617   //SPD Centrality (if PbPb / pPb)
 618   if(! fHistMultiplicitySPDBeforeTrigSel) {
 619     fHistMultiplicitySPDBeforeTrigSel = new TH1F("fHistMultiplicitySPDBeforeTrigSel",
 620                                                  "Centrality Distribution: SPD;SPD Centrality;Events",
 621                                                  200, 0, 200);
 622     fListHist->Add(fHistMultiplicitySPDBeforeTrigSel);
 623   }
 624   if(! fHistMultiplicitySPDForTrigEvt) {
 625     fHistMultiplicitySPDForTrigEvt = new TH1F("fHistMultiplicitySPDForTrigEvt",
 626                                               "Centrality Distribution: SPD;SPD Centrality;Events",
 627                                               200, 0, 200);
 628     fListHist->Add(fHistMultiplicitySPDForTrigEvt);
 629   }
 630   if(! fHistMultiplicitySPD) {
 631     fHistMultiplicitySPD = new TH1F("fHistMultiplicitySPD",
 632                                     "Centrality Distribution: SPD;SPD Centrality;Events",
 633                                     200, 0, 200);
 634     fListHist->Add(fHistMultiplicitySPD);
 635   }
 636   if(! fHistMultiplicitySPDNoTPCOnly) {
 637     fHistMultiplicitySPDNoTPCOnly = new TH1F("fHistMultiplicitySPDNoTPCOnly",
 638                                              "Centrality Distribution: SPD;SPD Centrality;Events",
 639                                              200, 0, 200);
 640     fListHist->Add(fHistMultiplicitySPDNoTPCOnly);
 641   }
 642   if(! fHistMultiplicitySPDNoTPCOnlyNoPileup) {
 643     fHistMultiplicitySPDNoTPCOnlyNoPileup = new TH1F("fHistMultiplicitySPDNoTPCOnlyNoPileup",
 644                                                      "Centrality Distribution: SPD;SPD Centrality;Events",
 645                                                      200, 0, 200);
 646     fListHist->Add(fHistMultiplicitySPDNoTPCOnlyNoPileup);
 647   }
 648
 649
 650 //----------------------------------
 651 // Primary Vertex Position Histos
 652 //----------------------------------
 653
 654    if(! fHistPVx) {
 655          fHistPVx = new TH1F("fHistPVx",
 656             "PV x position;Nbr of Evts;x",
 657             2000, -0.5, 0.5);
 658       fListHist->Add(fHistPVx);
 659    }
 660    if(! fHistPVy) {
 661          fHistPVy = new TH1F("fHistPVy",
 662             "PV y position;Nbr of Evts;y",
 663             2000, -0.5, 0.5);
 664       fListHist->Add(fHistPVy);
 665    }
 666    if(! fHistPVz) {
 667          fHistPVz = new TH1F("fHistPVz",
 668             "PV z position;Nbr of Evts;z",
 669             400, -20, 20);
 670       fListHist->Add(fHistPVz);
 671    }
 672
 673    if(! fHistPVxAnalysis) {
 674          fHistPVxAnalysis = new TH1F("fHistPVxAnalysis",
 675             "PV x position;Nbr of Evts;x",
 676             2000, -0.5, 0.5);
 677       fListHist->Add(fHistPVxAnalysis);
 678    }
 679    if(! fHistPVyAnalysis) {
 680          fHistPVyAnalysis = new TH1F("fHistPVyAnalysis",
 681             "PV y position;Nbr of Evts;y",
 682             2000, -0.5, 0.5);
 683       fListHist->Add(fHistPVyAnalysis);
 684    }
 685    if(! fHistPVzAnalysis) {
 686          fHistPVzAnalysis = new TH1F("fHistPVzAnalysis",
 687             "PV z position;Nbr of Evts;z",
 688             400, -20, 20);
 689       fListHist->Add(fHistPVzAnalysis);
 690    }
 691
 692    //List of Histograms: Normal
 693    PostData(1, fListHist);
 694
 695    //TTree Object: Saved to base directory. Should cache to disk while saving.
 696    //(Important to avoid excessive memory usage, particularly when merging)
 697    PostData(2, fTreeCascade);
 698
 699 }// end UserCreateOutputObjects
 700
 701
 702 //________________________________________________________________________
 703 void AliAnalysisTaskExtractCascade::UserExec(Option_t *)
 704 {
 705   // Main loop
 706   // Called for each event
 707
 708    AliESDEvent *lESDevent = 0x0;
 709
 710    Int_t    lNumberOfV0s                      = -1;
 711    Double_t lTrkgPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 712    Double_t lBestPrimaryVtxPos[3]          = {-100.0, -100.0, -100.0};
 713    Double_t lMagneticField                 = -10.;
 714
 715   // Connect to the InputEvent
 716   // After these lines, we should have an ESD/AOD event + the number of V0s in it.
 717
 718    // Appropriate for ESD analysis!
 719
 720    lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 721    if (!lESDevent) {
 722       AliWarning("ERROR: lESDevent not available \n");
 723       return;
 724    }
 725
 726 /* --- Acquisition of exact event ID
 727    fTreeVariableRunNumber = lESDevent->GetRunNumber();
 728    fTreeVariableEventNumber =
 729     ( ( ((ULong64_t)lESDevent->GetPeriodNumber() ) << 36 ) |
 730       ( ((ULong64_t)lESDevent->GetOrbitNumber () ) << 12 ) |
 731         ((ULong64_t)lESDevent->GetBunchCrossNumber() )  );
 732 */
 733
 734 //------------------------------------------------
 735 // Multiplicity Information Acquistion
 736 //------------------------------------------------
 737
 738    //REVISED multiplicity estimator after 'multiplicity day' (2011)
 739    Int_t lMultiplicity = -100;
 740   Int_t lMultiplicityV0A = -100;
 741   Int_t lMultiplicityZNA = -100;
 742   Int_t lMultiplicityTRK = -100;
 743   Int_t lMultiplicitySPD = -100;
 744
 745    //testing purposes
 746    if(fkIsNuclear == kFALSE) lMultiplicity =  fESDtrackCuts->GetReferenceMultiplicity(lESDevent, AliESDtrackCuts::kTrackletsITSTPC,  fEtaRefMult );
 747
 748    //---> If this is a nuclear collision, then go nuclear on "multiplicity" variable...
 749    //---> Warning: Experimental
 750    if(fkIsNuclear == kTRUE){
 751       AliCentrality* centrality;
 752       centrality = lESDevent->GetCentrality();
 753       lMultiplicity = ( ( Int_t ) ( centrality->GetCentralityPercentile( fCentralityEstimator.Data() ) ) );
 754       lMultiplicityV0A = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "V0A" ) ) );
 755       lMultiplicityZNA = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "ZNA" ) ) );
 756       lMultiplicityTRK = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "TRK" ) ) );
 757       lMultiplicitySPD = ( ( Int_t ) ( centrality->GetCentralityPercentile(   "CL1" ) ) );
 758       if (centrality->GetQuality()>1) {
 759         PostData(1, fListHist);
 760         PostData(2, fTreeCascade);
 761         return;
 762       }
 763    }
 764
 765    //Set variable for filling tree afterwards!
 766    //---> pp case......: GetReferenceMultiplicity
 767    //---> Pb-Pb case...: Centrality by V0M
 768
 769   fTreeCascVarMultiplicity = lMultiplicity;
 770   fTreeCascVarMultiplicityV0A = lMultiplicityV0A;
 771   fTreeCascVarMultiplicityZNA = lMultiplicityZNA;
 772   fTreeCascVarMultiplicityTRK = lMultiplicityTRK;
 773   fTreeCascVarMultiplicitySPD = lMultiplicitySPD;
 774
 775   fHistV0MultiplicityBeforeTrigSel->Fill ( lESDevent->GetNumberOfV0s() );
 776   fHistMultiplicityBeforeTrigSel->Fill ( lMultiplicity );
 777   fHistMultiplicityV0ABeforeTrigSel->Fill ( lMultiplicityV0A );
 778   fHistMultiplicityZNABeforeTrigSel->Fill ( lMultiplicityZNA );
 779   fHistMultiplicityTRKBeforeTrigSel->Fill ( lMultiplicityTRK );
 780   fHistMultiplicitySPDBeforeTrigSel->Fill ( lMultiplicitySPD );
 781
 782 //------------------------------------------------
 783 // Physics Selection
 784 //------------------------------------------------
 785
 786    UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 787    Bool_t isSelected = 0;
 788    isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
 789
 790    //pA triggering: CINT7
 791    if( fkSwitchINT7 ) isSelected = (maskIsSelected & AliVEvent::kINT7) == AliVEvent::kINT7;
 792
 793    //Standard Min-Bias Selection
 794    if ( ! isSelected ) {
 795         PostData(1, fListHist);
 796         PostData(2, fTreeCascade);
 797       return;
 798    }
 799
 800 //------------------------------------------------
 801 // Rerun cascade vertexer!
 802 //------------------------------------------------
 803 /*
 804   lESDevent->ResetCascades();
 805   lESDevent->ResetV0s();
 806
 807   AliV0vertexer lV0vtxer;
 808   AliCascadeVertexer lCascVtxer;
 809
 810   lV0vtxer.SetDefaultCuts(fV0Sels);
 811   lCascVtxer.SetDefaultCuts(fCascSels);
 812
 813   lV0vtxer.Tracks2V0vertices(lESDevent);
 814   lCascVtxer.V0sTracks2CascadeVertices(lESDevent);
 815 */
 816 //------------------------------------------------
 817 // After Trigger Selection
 818 //------------------------------------------------
 819
 820    lNumberOfV0s          = lESDevent->GetNumberOfV0s();
 821
 822    //Set variable for filling tree afterwards!
 823    fHistV0MultiplicityForTrigEvt->Fill(lNumberOfV0s);
 824    fHistMultiplicityForTrigEvt->Fill ( lMultiplicity );
 825    fHistMultiplicityV0AForTrigEvt       ->Fill( lMultiplicityV0A  );
 826    fHistMultiplicityZNAForTrigEvt       ->Fill( lMultiplicityZNA  );
 827    fHistMultiplicityTRKForTrigEvt       ->Fill( lMultiplicityTRK  );
 828    fHistMultiplicitySPDForTrigEvt       ->Fill( lMultiplicitySPD  );
 829
 830 //------------------------------------------------
 831 // Getting: Primary Vertex + MagField Info
 832 //------------------------------------------------
 833
 834    const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
 835    // get the vtx stored in ESD found with tracks
 836    lPrimaryTrackingESDVtx->GetXYZ( lTrkgPrimaryVtxPos );
 837
 838    const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
 839    // get the best primary vertex available for the event
 840    // As done in AliCascadeVertexer, we keep the one which is the best one available.
 841    // between : Tracking vertex > SPD vertex > TPC vertex > default SPD vertex
 842    // This one will be used for next calculations (DCA essentially)
 843    lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
 844
 845    Double_t lPrimaryVtxPosition[3];
 846    const AliVVertex *primaryVtx = lESDevent->GetPrimaryVertex();
 847    lPrimaryVtxPosition[0] = primaryVtx->GetX();
 848    lPrimaryVtxPosition[1] = primaryVtx->GetY();
 849    lPrimaryVtxPosition[2] = primaryVtx->GetZ();
 850    fHistPVx->Fill( lPrimaryVtxPosition[0] );
 851    fHistPVy->Fill( lPrimaryVtxPosition[1] );
 852    fHistPVz->Fill( lPrimaryVtxPosition[2] );
 853
 854   //------------------------------------------------
 855   // Primary Vertex Requirements Section:
 856   //  ---> pp and PbPb: Only requires |z|<10cm
 857   //  ---> pPb: all requirements checked at this stage
 858   //------------------------------------------------
 859
 860   //Roberto's PV selection criteria, implemented 17th April 2013
 861
 862   /* vertex selection */
 863   Bool_t fHasVertex = kFALSE;
 864
 865   const AliESDVertex *vertex = lESDevent->GetPrimaryVertexTracks();
 866   if (vertex->GetNContributors() < 1) {
 867     vertex = lESDevent->GetPrimaryVertexSPD();
 868     if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
 869     else fHasVertex = kTRUE;
 870     TString vtxTyp = vertex->GetTitle();
 871     Double_t cov[6]={0};
 872     vertex->GetCovarianceMatrix(cov);
 873     Double_t zRes = TMath::Sqrt(cov[5]);
 874     if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fHasVertex = kFALSE;
 875   }
 876   else fHasVertex = kTRUE;
 877
 878   //Is First event in chunk rejection: Still present!
 879   if(fkpAVertexSelection==kTRUE && fHasVertex == kFALSE) {
 880     AliWarning("Pb / | PV does not satisfy selection criteria!");
 881     PostData(1, fListHist);
 882     PostData(2, fTreeCascade);
 883     return;
 884   }
 885
 886   //Is First event in chunk rejection: Still present!
 887   if(fkpAVertexSelection==kTRUE && fUtils->IsFirstEventInChunk(lESDevent)) {
 888     AliWarning("Pb / | This is the first event in the chunk!");
 889     PostData(1, fListHist);
 890     PostData(2, fTreeCascade);
 891     return;
 892   }
 893
 894   //17 April Fix: Always do primary vertex Z selection, after pA vertex selection from Roberto
 895   if(TMath::Abs(lBestPrimaryVtxPos[2]) > 10.0) {
 896     AliWarning("Pb / | Z position of Best Prim Vtx | > 10.0 cm ... return !");
 897     PostData(1, fListHist);
 898     PostData(2, fTreeCascade);
 899     return;
 900   }
 901
 902
 903   lMagneticField = lESDevent->GetMagneticField( );
 904   fHistV0MultiplicityForSelEvt ->Fill( lNumberOfV0s );
 905   fHistMultiplicity->Fill(lMultiplicity);
 906   fHistMultiplicityV0A->Fill(lMultiplicityV0A);
 907   fHistMultiplicityZNA->Fill(lMultiplicityZNA);
 908   fHistMultiplicityTRK->Fill(lMultiplicityTRK);
 909   fHistMultiplicitySPD->Fill(lMultiplicitySPD);
 910
 911 //------------------------------------------------
 912 // SKIP: Events with well-established PVtx
 913 //------------------------------------------------
 914
 915    const AliESDVertex *lPrimaryTrackingESDVtxCheck = lESDevent->GetPrimaryVertexTracks();
 916    const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
 917    if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtxCheck->GetStatus() && fkpAVertexSelection==kFALSE ){
 918       AliWarning("Pb / No SPD prim. vertex nor prim. Tracking vertex ... return !");
 919         PostData(1, fListHist);
 920         PostData(2, fTreeCascade);
 921       return;
 922    }
 923    fHistV0MultiplicityForSelEvtNoTPCOnly ->Fill( lNumberOfV0s );
 924    fHistMultiplicityNoTPCOnly->Fill(lMultiplicity);
 925    fHistMultiplicityV0ANoTPCOnly->Fill(lMultiplicityV0A);
 926    fHistMultiplicityZNANoTPCOnly->Fill(lMultiplicityZNA);
 927    fHistMultiplicityTRKNoTPCOnly->Fill(lMultiplicityTRK);
 928    fHistMultiplicitySPDNoTPCOnly->Fill(lMultiplicitySPD);
 929
 930 //------------------------------------------------
 931 // Pileup Rejection Studies
 932 //------------------------------------------------
 933
 934    // FIXME : quality selection regarding pile-up rejection
 935    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
 936       AliWarning("Pb / This is tagged as Pileup from SPD... return !");
 937         PostData(1, fListHist);
 938         PostData(2, fTreeCascade);
 939       return;
 940    }
 941    fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( lNumberOfV0s );
 942    fHistMultiplicityNoTPCOnlyNoPileup->Fill(lMultiplicity);
 943    fHistMultiplicityV0ANoTPCOnlyNoPileup->Fill(lMultiplicityV0A);
 944    fHistMultiplicityZNANoTPCOnlyNoPileup->Fill(lMultiplicityZNA);
 945    fHistMultiplicityTRKNoTPCOnlyNoPileup->Fill(lMultiplicityTRK);
 946    fHistMultiplicitySPDNoTPCOnlyNoPileup->Fill(lMultiplicitySPD);
 947
 948    //Do control histograms without the IsFromVertexerZ events, but consider them in analysis...
 949    if( ! (lESDevent->GetPrimaryVertex()->IsFromVertexerZ() )     ){
 950       fHistPVxAnalysis->Fill( lPrimaryVtxPosition[0] );
 951       fHistPVyAnalysis->Fill( lPrimaryVtxPosition[1] );
 952       fHistPVzAnalysis->Fill( lPrimaryVtxPosition[2] );
 953    }
 954
 955 //------------------------------------------------
 956 // MAIN CASCADE LOOP STARTS HERE
 957 //------------------------------------------------
 958 // Code Credit: Antonin Maire (thanks^100)
 959 // ---> This is an adaptation
 960
 961   Long_t ncascades = 0;
 962         ncascades = lESDevent->GetNumberOfCascades();
 963
 964   for (Int_t iXi = 0; iXi < ncascades; iXi++){
 965     //------------------------------------------------
 966     // Initializations
 967     //------------------------------------------------
 968           //Double_t lTrkgPrimaryVtxRadius3D = -500.0;
 969           //Double_t lBestPrimaryVtxRadius3D = -500.0;
 970
 971           // - 1st part of initialisation : variables needed to store AliESDCascade data members
 972           Double_t lEffMassXi      = 0. ;
 973           //Double_t lChi2Xi         = -1. ;
 974           Double_t lDcaXiDaughters = -1. ;
 975           Double_t lXiCosineOfPointingAngle = -1. ;
 976           Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
 977           Double_t lXiRadius = -1000. ;
 978
 979           // - 2nd part of initialisation : Nbr of clusters within TPC for the 3 daughter cascade tracks
 980           Int_t    lPosTPCClusters    = -1; // For ESD only ...//FIXME : wait for availability in AOD
 981           Int_t    lNegTPCClusters    = -1; // For ESD only ...
 982           Int_t    lBachTPCClusters   = -1; // For ESD only ...
 983
 984           // - 3rd part of initialisation : about V0 part in cascades
 985           Double_t lInvMassLambdaAsCascDghter = 0.;
 986           //Double_t lV0Chi2Xi         = -1. ;
 987           Double_t lDcaV0DaughtersXi = -1.;
 988
 989           Double_t lDcaBachToPrimVertexXi = -1., lDcaV0ToPrimVertexXi = -1.;
 990           Double_t lDcaPosToPrimVertexXi  = -1.;
 991           Double_t lDcaNegToPrimVertexXi  = -1.;
 992           Double_t lV0CosineOfPointingAngleXi = -1. ;
 993           Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; // Position of VO coming from cascade
 994           Double_t lV0RadiusXi = -1000.0;
 995           Double_t lV0quality  = 0.;
 996
 997           // - 4th part of initialisation : Effective masses
 998           Double_t lInvMassXiMinus    = 0.;
 999           Double_t lInvMassXiPlus     = 0.;
1000           Double_t lInvMassOmegaMinus = 0.;
1001           Double_t lInvMassOmegaPlus  = 0.;
1002
1003           // - 6th part of initialisation : extra info for QA
1004           Double_t lXiMomX       = 0. , lXiMomY = 0., lXiMomZ = 0.;
1005           Double_t lXiTransvMom  = 0. ;
1006           Double_t lXiTransvMomMC= 0. ;
1007           Double_t lXiTotMom     = 0. ;
1008
1009           Double_t lBachMomX       = 0., lBachMomY  = 0., lBachMomZ   = 0.;
1010           //Double_t lBachTransvMom  = 0.;
1011           //Double_t lBachTotMom     = 0.;
1012
1013     fTreeCascVarNegNSigmaPion   = -100;
1014     fTreeCascVarNegNSigmaProton = -100;
1015     fTreeCascVarPosNSigmaPion   = -100;
1016     fTreeCascVarPosNSigmaProton = -100;
1017     fTreeCascVarBachNSigmaPion  = -100;
1018     fTreeCascVarBachNSigmaKaon  = -100;
1019
1020           Short_t  lChargeXi = -2;
1021           //Double_t lV0toXiCosineOfPointingAngle = 0. ;
1022
1023           Double_t lRapXi   = -20.0, lRapOmega = -20.0, lRapMC = -20.0; //  lEta = -20.0, lTheta = 360., lPhi = 720. ;
1024           //Double_t lAlphaXi = -200., lPtArmXi  = -200.0;
1025
1026     // -------------------------------------
1027     // II.ESD - Calculation Part dedicated to Xi vertices (ESD)
1028
1029           AliESDcascade *xi = lESDevent->GetCascade(iXi);
1030           if (!xi) continue;
1031
1032
1033                   // - II.Step 1 : around primary vertex
1034                   //-------------
1035           //lTrkgPrimaryVtxRadius3D = TMath::Sqrt(  lTrkgPrimaryVtxPos[0] * lTrkgPrimaryVtxPos[0] +
1036           //                                        lTrkgPrimaryVtxPos[1] * lTrkgPrimaryVtxPos[1] +
1037           //                                        lTrkgPrimaryVtxPos[2] * lTrkgPrimaryVtxPos[2] );
1038
1039           //lBestPrimaryVtxRadius3D = TMath::Sqrt(  lBestPrimaryVtxPos[0] * lBestPrimaryVtxPos[0] +
1040           //                                        lBestPrimaryVtxPos[1] * lBestPrimaryVtxPos[1] +
1041           //                                        lBestPrimaryVtxPos[2] * lBestPrimaryVtxPos[2] );
1042
1043                 // - II.Step 2 : Assigning the necessary variables for specific AliESDcascade data members (ESD)
1044                 //-------------
1045           lV0quality = 0.;
1046           xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
1047
1048           lEffMassXi                    = xi->GetEffMassXi();
1049           //lChi2Xi                         = xi->GetChi2Xi();
1050           lDcaXiDaughters       = xi->GetDcaXiDaughters();
1051           lXiCosineOfPointingAngle                  = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0],
1052                                                                                  lBestPrimaryVtxPos[1],
1053                                                                                  lBestPrimaryVtxPos[2] );
1054                   // Take care : the best available vertex should be used (like in AliCascadeVertexer)
1055
1056           xi->GetXYZcascade( lPosXi[0],  lPosXi[1], lPosXi[2] );
1057           lXiRadius                     = TMath::Sqrt( lPosXi[0]*lPosXi[0]  +  lPosXi[1]*lPosXi[1] );
1058
1059                 // - II.Step 3 : around the tracks : Bach + V0 (ESD)
1060                 // ~ Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
1061                 //-------------
1062
1063         UInt_t lIdxPosXi        = (UInt_t) TMath::Abs( xi->GetPindex() );
1064         UInt_t lIdxNegXi        = (UInt_t) TMath::Abs( xi->GetNindex() );
1065         UInt_t lBachIdx         = (UInt_t) TMath::Abs( xi->GetBindex() );
1066                 // Care track label can be negative in MC production (linked with the track quality)
1067                 // However = normally, not the case for track index ...
1068
1069           // FIXME : rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
1070           if(lBachIdx == lIdxNegXi) {
1071                   AliWarning("Pb / Idx(Bach. track) = Idx(Neg. track) ... continue!"); continue;
1072           }
1073     if(lBachIdx == lIdxPosXi) {
1074         AliWarning("Pb / Idx(Bach. track) = Idx(Pos. track) ... continue!"); continue;
1075           }
1076
1077           AliESDtrack *pTrackXi         = lESDevent->GetTrack( lIdxPosXi );
1078           AliESDtrack *nTrackXi         = lESDevent->GetTrack( lIdxNegXi );
1079           AliESDtrack *bachTrackXi      = lESDevent->GetTrack( lBachIdx );
1080
1081           if (!pTrackXi || !nTrackXi || !bachTrackXi ) {
1082                   AliWarning("ERROR: Could not retrieve one of the 3 ESD daughter tracks of the cascade ...");
1083                   continue;
1084           }
1085
1086    fTreeCascVarPosEta = pTrackXi->Eta();
1087    fTreeCascVarNegEta = nTrackXi->Eta();
1088    fTreeCascVarBachEta = bachTrackXi->Eta();
1089
1090     //------------------------------------------------
1091     // TPC dEdx information
1092     //------------------------------------------------
1093     fTreeCascVarNegNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kPion   );
1094     fTreeCascVarNegNSigmaProton = fPIDResponse->NumberOfSigmasTPC( nTrackXi, AliPID::kProton );
1095     fTreeCascVarPosNSigmaPion   = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kPion );
1096     fTreeCascVarPosNSigmaProton = fPIDResponse->NumberOfSigmasTPC( pTrackXi, AliPID::kProton );
1097     fTreeCascVarBachNSigmaPion  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kPion );
1098     fTreeCascVarBachNSigmaKaon  = fPIDResponse->NumberOfSigmasTPC( bachTrackXi, AliPID::kKaon );
1099
1100     //------------------------------------------------
1101     // TPC Number of clusters info
1102     // --- modified to save the smallest number
1103     // --- of TPC clusters for the 3 tracks
1104     //------------------------------------------------
1105
1106           lPosTPCClusters   = pTrackXi->GetTPCNcls();
1107           lNegTPCClusters   = nTrackXi->GetTPCNcls();
1108           lBachTPCClusters  = bachTrackXi->GetTPCNcls();
1109
1110     // 1 - Poor quality related to TPCrefit
1111           ULong_t pStatus    = pTrackXi->GetStatus();
1112           ULong_t nStatus    = nTrackXi->GetStatus();
1113           ULong_t bachStatus = bachTrackXi->GetStatus();
1114     if ((pStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Pos. track has no TPCrefit ... continue!"); continue; }
1115     if ((nStatus&AliESDtrack::kTPCrefit)    == 0) { AliWarning("Pb / V0 Neg. track has no TPCrefit ... continue!"); continue; }
1116     if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning("Pb / Bach.   track has no TPCrefit ... continue!"); continue; }
1117           // 2 - Poor quality related to TPC clusters: lowest cut of 70 clusters
1118     if(lPosTPCClusters  < 70) { AliWarning("Pb / V0 Pos. track has less than 70 TPC clusters ... continue!"); continue; }
1119           if(lNegTPCClusters  < 70) { AliWarning("Pb / V0 Neg. track has less than 70 TPC clusters ... continue!"); continue; }
1120           if(lBachTPCClusters < 70) { AliWarning("Pb / Bach.   track has less than 70 TPC clusters ... continue!"); continue; }
1121           Int_t leastnumberofclusters = 1000;
1122           if( lPosTPCClusters < leastnumberofclusters ) leastnumberofclusters = lPosTPCClusters;
1123           if( lNegTPCClusters < leastnumberofclusters ) leastnumberofclusters = lNegTPCClusters;
1124           if( lBachTPCClusters < leastnumberofclusters ) leastnumberofclusters = lBachTPCClusters;
1125
1126           lInvMassLambdaAsCascDghter    = xi->GetEffMass();
1127           // This value shouldn't change, whatever the working hyp. is : Xi-, Xi+, Omega-, Omega+
1128           lDcaV0DaughtersXi             = xi->GetDcaV0Daughters();
1129           //lV0Chi2Xi                   = xi->GetChi2V0();
1130
1131           lV0CosineOfPointingAngleXi    = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0],
1132                                                                             lBestPrimaryVtxPos[1],
1133                                                                             lBestPrimaryVtxPos[2] );
1134
1135           lDcaV0ToPrimVertexXi          = xi->GetD( lBestPrimaryVtxPos[0],
1136                                                       lBestPrimaryVtxPos[1],
1137                                                       lBestPrimaryVtxPos[2] );
1138
1139           lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD(       lBestPrimaryVtxPos[0],
1140                                                                 lBestPrimaryVtxPos[1],
1141                                                                 lMagneticField  ) );
1142                                           // Note : AliExternalTrackParam::GetD returns an algebraic value ...
1143
1144           xi->GetXYZ( lPosV0Xi[0],  lPosV0Xi[1], lPosV0Xi[2] );
1145           lV0RadiusXi           = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0]  +  lPosV0Xi[1]*lPosV0Xi[1] );
1146
1147           lDcaPosToPrimVertexXi         = TMath::Abs( pTrackXi  ->GetD( lBestPrimaryVtxPos[0],
1148                                                                 lBestPrimaryVtxPos[1],
1149                                                                 lMagneticField  )     );
1150
1151           lDcaNegToPrimVertexXi         = TMath::Abs( nTrackXi  ->GetD( lBestPrimaryVtxPos[0],
1152                                                                 lBestPrimaryVtxPos[1],
1153                                                                 lMagneticField  )     );
1154
1155           // - II.Step 4 : around effective masses (ESD)
1156           // ~ change mass hypotheses to cover all the possibilities :  Xi-/+, Omega -/+
1157
1158           if( bachTrackXi->Charge() < 0 )       {
1159                   lV0quality = 0.;
1160                   xi->ChangeMassHypothesis(lV0quality , 3312);
1161                           // Calculate the effective mass of the Xi- candidate.
1162                           // pdg code 3312 = Xi-
1163                   lInvMassXiMinus = xi->GetEffMassXi();
1164
1165                   lV0quality = 0.;
1166                   xi->ChangeMassHypothesis(lV0quality , 3334);
1167                           // Calculate the effective mass of the Xi- candidate.
1168                           // pdg code 3334 = Omega-
1169                   lInvMassOmegaMinus = xi->GetEffMassXi();
1170
1171                   lV0quality = 0.;
1172                   xi->ChangeMassHypothesis(lV0quality , 3312);  // Back to default hyp.
1173           }// end if negative bachelor
1174
1175
1176           if( bachTrackXi->Charge() >  0 ){
1177                   lV0quality = 0.;
1178                   xi->ChangeMassHypothesis(lV0quality , -3312);
1179                           // Calculate the effective mass of the Xi+ candidate.
1180                           // pdg code -3312 = Xi+
1181                   lInvMassXiPlus = xi->GetEffMassXi();
1182
1183                   lV0quality = 0.;
1184                   xi->ChangeMassHypothesis(lV0quality , -3334);
1185                           // Calculate the effective mass of the Xi+ candidate.
1186                           // pdg code -3334  = Omega+
1187                   lInvMassOmegaPlus = xi->GetEffMassXi();
1188
1189                   lV0quality = 0.;
1190                   xi->ChangeMassHypothesis(lV0quality , -3312);         // Back to "default" hyp.
1191           }// end if positive bachelor
1192                   // - II.Step 6 : extra info for QA (ESD)
1193                   // miscellaneous pieces of info that may help regarding data quality assessment.
1194                   //-------------
1195
1196           xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
1197                   lXiTransvMom          = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY );
1198                   lXiTotMom     = TMath::Sqrt( lXiMomX*lXiMomX   + lXiMomY*lXiMomY   + lXiMomZ*lXiMomZ );
1199
1200           xi->GetBPxPyPz(  lBachMomX,  lBachMomY,  lBachMomZ );
1201                   //lBachTransvMom  = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY );
1202                   //lBachTotMom         = TMath::Sqrt( lBachMomX*lBachMomX   + lBachMomY*lBachMomY  +  lBachMomZ*lBachMomZ  );
1203
1204           lChargeXi = xi->Charge();
1205
1206           //lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
1207
1208           lRapXi    = xi->RapXi();
1209           lRapOmega = xi->RapOmega();
1210           //lEta      = xi->Eta();
1211           //lTheta    = xi->Theta() *180.0/TMath::Pi();
1212           //lPhi      = xi->Phi()   *180.0/TMath::Pi();
1213           //lAlphaXi  = xi->AlphaXi();
1214           //lPtArmXi  = xi->PtArmXi();
1215
1216   //------------------------------------------------
1217   // Set Variables for adding to tree
1218   //------------------------------------------------
1219
1220 /* 1*/          fTreeCascVarCharge      = lChargeXi;
1221 /* 2*/          if(lInvMassXiMinus!=0)    fTreeCascVarMassAsXi = lInvMassXiMinus;
1222 /* 2*/          if(lInvMassXiPlus!=0)     fTreeCascVarMassAsXi = lInvMassXiPlus;
1223 /* 3*/          if(lInvMassOmegaMinus!=0) fTreeCascVarMassAsOmega = lInvMassOmegaMinus;
1224 /* 3*/          if(lInvMassOmegaPlus!=0)  fTreeCascVarMassAsOmega = lInvMassOmegaPlus;
1225 /* 4*/          fTreeCascVarPt = lXiTransvMom;
1226 /* 4*/          fTreeCascVarPtMC = lXiTransvMomMC;
1227 /* 5*/          fTreeCascVarRapXi = lRapXi ;
1228 /* 5*/          fTreeCascVarRapMC = lRapMC ;
1229 /* 6*/          fTreeCascVarRapOmega = lRapOmega ;
1230 /* 7*/          fTreeCascVarDCACascDaughters = lDcaXiDaughters;
1231 /* 8*/          fTreeCascVarDCABachToPrimVtx = lDcaBachToPrimVertexXi;
1232 /* 9*/          fTreeCascVarDCAV0Daughters = lDcaV0DaughtersXi;
1233 /*10*/          fTreeCascVarDCAV0ToPrimVtx = lDcaV0ToPrimVertexXi;
1234 /*11*/          fTreeCascVarDCAPosToPrimVtx = lDcaPosToPrimVertexXi;
1235 /*12*/          fTreeCascVarDCANegToPrimVtx = lDcaNegToPrimVertexXi;
1236 /*13*/          fTreeCascVarCascCosPointingAngle = lXiCosineOfPointingAngle;
1237 /*14*/          fTreeCascVarCascRadius = lXiRadius;
1238 /*15*/          fTreeCascVarV0Mass = lInvMassLambdaAsCascDghter;
1239 /*16*/          fTreeCascVarV0CosPointingAngle = lV0CosineOfPointingAngleXi;
1240 /*17*/          fTreeCascVarV0Radius = lV0RadiusXi;
1241 /*20*/          fTreeCascVarLeastNbrClusters = leastnumberofclusters;
1242 /*21*/          fTreeCascVarMultiplicity = lMultiplicity; //multiplicity, whatever that may be
1243
1244 /*23*/          fTreeCascVarDistOverTotMom = TMath::Sqrt(
1245                                                 TMath::Power( lPosXi[0] - lBestPrimaryVtxPos[0] , 2) +
1246                                                 TMath::Power( lPosXi[1] - lBestPrimaryVtxPos[1] , 2) +
1247                                                 TMath::Power( lPosXi[2] - lBestPrimaryVtxPos[2] , 2)
1248                                         );
1249 /*23*/          fTreeCascVarDistOverTotMom /= (lXiTotMom+1e-13);
1250
1251 //All vars not specified here: specified elsewhere!
1252
1253 //------------------------------------------------
1254 // Fill Tree!
1255 //------------------------------------------------
1256
1257 // The conditional is meant to decrease excessive
1258 // memory usage! Be careful when loosening the
1259 // cut!
1260
1261   //Xi    Mass window: 150MeV wide
1262   //Omega mass window: 150MeV wide
1263
1264   if( (fTreeCascVarMassAsXi<1.32+0.075&&fTreeCascVarMassAsXi>1.32-0.075) ||
1265       (fTreeCascVarMassAsOmega<1.68+0.075&&fTreeCascVarMassAsOmega>1.68-0.075) ){
1266       fTreeCascade->Fill();
1267   }
1268
1269 //------------------------------------------------
1270 // Fill tree over.
1271 //------------------------------------------------
1272
1273         }// end of the Cascade loop (ESD or AOD)
1274
1275    // Post output data.
1276    PostData(1, fListHist);
1277    PostData(2, fTreeCascade);
1278 }
1279
1280 //________________________________________________________________________
1281 void AliAnalysisTaskExtractCascade::Terminate(Option_t *)
1282 {
1283    // Draw result to the screen
1284    // Called once at the end of the query
1285
1286    TList *cRetrievedList = 0x0;
1287    cRetrievedList = (TList*)GetOutputData(1);
1288    if(!cRetrievedList){
1289       Printf("ERROR - AliAnalysisTaskExtractCascade : ouput data container list not available\n");
1290       return;
1291    }
1292
1293    fHistV0MultiplicityForTrigEvt = dynamic_cast<TH1F*> (  cRetrievedList->FindObject("fHistV0MultiplicityForTrigEvt")  );
1294    if (!fHistV0MultiplicityForTrigEvt) {
1295       Printf("ERROR - AliAnalysisTaskExtractCascade : fHistV0MultiplicityForTrigEvt not available");
1296       return;
1297    }
1298
1299    TCanvas *canCheck = new TCanvas("AliAnalysisTaskExtractCascade","V0 Multiplicity",10,10,510,510);
1300    canCheck->cd(1)->SetLogy();
1301
1302    fHistV0MultiplicityForTrigEvt->SetMarkerStyle(22);
1303    fHistV0MultiplicityForTrigEvt->DrawCopy("E");
1304 }
1305
1306 //----------------------------------------------------------------------------
1307
1308 Double_t AliAnalysisTaskExtractCascade::MyRapidity(Double_t rE, Double_t rPz) const
1309 {
1310    // Local calculation for rapidity
1311    Double_t ReturnValue = -100;
1312    if( (rE-rPz+1.e-13) != 0 && (rE+rPz) != 0 ){
1313       ReturnValue =  0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
1314    }
1315    return ReturnValue;
1316 }