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